Comprehensive assessment of computational methods for cancer immunoediting

Abstract

Cancer immunoediting reflects the role of the immune system in eliminating tumor cells and shaping tumor immunogenicity, which leaves marks in the genome. In this study, we systematically evaluate four methods for quantifying immunoediting. In colorectal cancer samples from The Cancer Genome Atlas, we found that these methods identified 78.41%, 46.17%, 36.61%, and 4.92% of immunoedited samples, respectively, covering 92.90% of all colorectal cancer samples. Comparison of 10 patient-derived xenografts (PDXs) with their original tumors showed that different methods identified reduced immune selection in PDXs ranging from 44.44% to 60.0%. The proportion of such PDX-tumor pairs increases to 77.78% when considering the union of results from multiple methods, indicating the complementarity of these methods. We find that observed-to-expected ratios highly rely on neoantigen selection criteria and reference datasets. In contrast, HLA-binding mutation ratio, immune dN/dS, and enrichment score of cancer cell fraction were less affected by these factors. Our findings suggest integration of multiple methods may benefit future immunoediting analyses.

Keywords: CP: cancer biology; CP: genetics; computational methods; immune selection; immunoediting.

Graphical abstract

Figure 1

Schematic of the study and assessment based on TCGA multiple epithelial cancers (A) Study overview. First, genome data from TCGA, publicly available PDX models, and an additional cohort were collected. Second, immunoediting scores calculated by four types of methods (including five specific methods) are compared to assess the consistency of these methods in measuring immunoediting scores and the ability to detect immunoediting events in tumors. Finally, the potential impact of neoantigen selection strategies, mutation counts, HLA alleles, clinical features, and reference datasets on the performance of these methods is assessed. (B) Left: Spearman correlation between immunoediting scores derived from different immunoediting quantification methods across various cancer types. Each column represents a cancer type, and each row represents the comparison between two methods. Correlation coefficients with p < 0.05 are shown in each cell. Right: The distribution of correlation between methods, the shape of the point indicates the significance of correlation (triangle for significant correlations and inverted triangle for non-significant ones). (C) Top: The number of samples classified by each method as edited, unedited, or unassigned across multiple cancer types. Bottom: The proportion of each classification status by the four methods across each cancer type. Abbreviations: TCGA, The Cancer Genome Atlas; PDX, patient-derived xenograft; HLA, human leukocyte antigen. See also Table S1.

Figure 2

In-depth assessments of various methods based on colorectal cancer and lung adenocarcinoma (A and D) Heatmap shows the similarity in immunoediting scores calculated by different methods across TCGA CRC samples (A) and LUAD samples (D). (B and E) Bar chart and heatmap show the total number of samples classified by each method as edited, unedited, or unassigned in CRC (B) and LUAD (E). (C and F) UpSet plot shows the intersection of immunoedited samples identified by different methods in CRC (C) and LUAD (F). The left bar plot shows the total number of immunoedited samples per method. The top bar plot shows the shared samples among methods, with color indicating the number of methods involved. See also Figures S2 and S3.

Figure 3

Immunoediting analysis of PDX models and simulation datasets (A) Bar chart summarizes the total number of CRC samples classified by each method as edited, unedited, or unassigned. The bottom heatmap shows the detailed categorization of each sample. ES-CCF_without_p identified samples as edited when their ES-CCF was less than 0. (B) Heatmap shows the immunoediting scores calculated by different methods across PDX models and corresponding tumors, each row represents a patient. (C) Proportion of PDX models that exhibited reduced immune selection pressure compared to their original tumors. (D) CCF distribution of mutations in the simulated dataset under different immune selection. (E) Distribution of immunoediting scores obtained by different methods for simulated datasets under different immune selection, the two-sided Wilcoxon rank-sum tests are used to compare immunoediting scores across multiple groups and p values were corrected using the FDR method. Immunoediting scores were significantly decreased for all methods as immunoselection was increased (p < 0.05). See also Figure S3.

Figure 4

Immunoediting scores across methods under various neoantigen screening conditions (A and C) Distribution of immunoediting scores obtained by different methods for TCGA CRC samples after screening for neoantigens, distinguished by various peptide lengths (A) or binding affinities (C). FDR-adjusted p values from two-sided Wilcoxon rank-sum tests are indicated. (B and D) Density plot shows the number of neoantigens in CRC samples after screening for neoantigens based on different peptide lengths (B) or binding affinities (D). (E) Neoantigen sharing identified by NetMHCpan, MHCFlurry, and MixMHCpred. (F) Distribution of immunoediting scores obtained by different methods using neoantigens predicted by NetMHCpan, MHCFlurry, and MixMHCpred. FDR-adjusted p values from two-sided Wilcoxon rank-sum tests are indicated. Statistical significance is denoted by p values as ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, ns = not significant. See also Figures S4.

Figure 5

Impact of mutation number and HLA LOH on computational methods of immunoediting (A) Spearman correlation coefficients between immunoediting scores with mutation number, TMB, neoantigen number, and neoantigen burden. Correlation coefficients and significance are shown in each cell (hypermutated samples were excluded from the correlation analysis). (B) Distribution of immunoediting scores across low-, medium-, and high-TMB groups. FDR-adjusted p values from two-sided Wilcoxon rank-sum tests are indicated. (C) Top: Distribution of immunoediting scores obtained from different methods in each group of samples after grouping the samples according to the number of mutations. Bottom: Variance change in immunoediting scores among samples from different groups. The lines follow a fitted linear regression. Linear regression coefficients and p values are shown in red text. (D) Density plot shows the CCF of neoantigen and non-antigenic mutations in CRC samples. (E) Distribution of immunoediting scores obtained from different methods between samples with and without HLA LOH. p values from two-sided Wilcoxon rank-sum tests are indicated. (F) The proportion of samples identified as immunoedited by different methods in TCGA CRC samples with or without HLA LOH. p values from Fisher’s exact tests are indicated. Statistical significance is denoted by p values as ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, ns = not significant. See also Figure S5.

Figure 6

Comparison of immunoediting scores obtained with different HLA alleles selection strategies and different reference datasets (A) Distribution of immunoediting scores obtained from different methods using both common and sample-specific HLA alleles. (B) Scatterplots compare the immunoediting scores generated by common and sample-specific HLA alleles, the red dashed line represents a value of ${OEratio}_{refSample}$ = 1, ${OEratio}_{codon}$ = 1, HBMR = 1, or ES-CCF = 0. Spearman correlation coefficients and significance are shown. The lines follow a fitted linear regression. (C) Distribution of ${OEratio}_{refSample}$ and ${OEratio}_{codon}$ for TCGA CRC samples when using different reference datasets, the red dashed line represents a value of ${OEratio}_{refSample}$ = 1 or ${OEratio}_{codon}$ = 1. FDR-adjusted p values from two-sided Wilcoxon rank-sum tests among groups are demonstrated in Table S2. (D) Top: Distribution of TMB for different reference datasets from TCGA; bottom: the total number of samples for different reference datasets from TCGA. (E) Distribution of immunoediting scores for TCGA colorectal cancer samples obtained using OE-ratios with different reference datasets from other studies. Statistical significance is denoted by p values as ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001. See also Figure S6.

Figure 7

Influence of reference datasets and clinical variables on the quantification of immunoediting (A–D) Distribution of immunoediting scores by method in TCGA colorectal cancer samples, grouped by sex (A), age (B), stage (C), and race (D). Statistical significance is denoted by p values as ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001. See also Figure S7.