Neoantigen immunogenicity landscapes and evolution of tumor ecosystems during immunotherapy with nivolumab

Abstract

Neoantigen immunoediting drives immune checkpoint blockade efficacy, yet the molecular features of neoantigens and how neoantigen immunogenicity shapes treatment response remain poorly understood. To address these questions, 80 patients with non-small cell lung cancer were enrolled in the biomarker cohort of CheckMate 153 (CA209-153), which collected radiographic guided biopsy samples before treatment and during treatment with nivolumab. Early loss of mutations and neoantigens during therapy are both associated with clinical benefit. We examined 1,453 candidate neoantigens, including many of which that had reduced cancer cell fraction after treatment with nivolumab, and identified 196 neopeptides that were recognized by T cells. Mapping these neoantigens to clonal dynamics, evolutionary trajectories and clinical response revealed a strong selection against immunogenic neoantigen-harboring clones. We identified position-specific amino acid and physiochemical features related to immunogenicity and developed an immunogenicity score. Nivolumab-induced microenvironmental evolution in non-small cell lung cancer shared some similarities with melanoma, yet critical differences were apparent. This study provides unprecedented molecular portraits of neoantigen landscapes underlying nivolumab's mechanism of action.

Extended Data Fig. 1 ∣. Mutational landscape is altered between responders and non-responders pre-therapy and on-therapy.

(a) Gene alteration analysis pre-therapy vs on-therapy shows non-responders have enrichment in SKT11 and KEAP1 mutations. Using selected genes related to NSCLC as highlighted by TCGA, a comparison of mutation frequency pre-therapy within responders and non-responders is shown. (b) On-therapy non-responders are enriched in STK11 mutations. Using selected genes related to NSCLC, a comparison of mutation frequency on-therapy within responders and non-responders is shown. (c) Quantification of overall mutational differences pre-therapy between responders and non-responders with odds ratio, confidence intervals, and p-value (Fisher’s exact test) shown for the top genes. (d) Quantification of overall mutational differences on-therapy between responders and non-responders with odds ratio, confidence intervals, and p-value (Fisher’s exact test) shown for the top genes. (e) Antigen presentation machinery genes show broad deletions by a GISTIC threshold of less than or equal to −1, which were not specific to any response group. (f) HLA evolutionary divergence (HED) was analyzed based on HLA-types from Polysolver using previously described methods. Splitting HED by optimal threshold test in R package Survival, patients with High HED > 5.8 were more likely to benefit from therapy (log-rank test, p = 0.02). (g) HED corrected for LOH using LOHLA was performed by removing those HLA alleles which had LOHHLA unique Pval <0.05 (log-rank test, p = 0.043). All boxplots represent the median, two hinges, and two whiskers, with the lower and upper hinges corresponding to the first and third quartile and the whiskers extending to the 1.5*IQR.

Extended Data Fig. 2 ∣. Somatic mutations occurring in pathways associated with response.

(a) Pathway enrichment analysis identifies dissolution of fibrin clot pathway enriched in non-responders pre-therapy. Enrichment analysis of genes mutated pre-therapy in responders and non-responders with a minimum gene set size of n = 5 highlights pathways with enriched mutations (Fisher’s exact test). (b) Pathway enrichment analysis identified activation of HOX genes pathway enriched in non-responders on-therapy. Enrichment analysis of genes mutated on-therapy in responders and non-responders with a minimum gene set size of n = 5 highlights pathways with enriched mutations (Fisher’s exact test). For all plots (p < 0.05 = *, p < 0.01 = **, p < 0.001 = ***). All boxplots represent the median, two hinges, and two whiskers, with the lower and upper hinges corresponding to the first and third quartile with and whiskers extending to the 1.5*IQR.

Extended Data Fig. 3 ∣. Associations of tumor ecosystem changes and clinical outcome.

(a) Deconvolution analysis of pre-therapy samples. Associations of tumor ecosystem changes and clinical outcome. (a) Deconvolution analysis of pre-therapy samples. ssGSEA deconvolution analysis of immune signatures for pathways and cell types related to immunotherapy were used to score pre-treatment RNAseq samples (n = 24). (b) Deconvolution of on-therapy samples. On-therapy RNAseq data were used for the deconvolution of immunotherapy-related cell types and pathways with unsupervised clustering. (c) Nivolumab-activated gene set common between melanoma and NSCLC patients. Paired differential expression analysis during treatment for CHECKMATE-038 and CHECKMATE-153 were compared, highlighting the genes with log FC > 0 and <0 which have a p < 0.05 and base mean expression >200 (black = commonly upregulated between CHECKMATE-038 (melanoma) and CHECKMATE-153 (NSCLC), purple = unique to melanoma, orange = unique to NSCLC, red = commonly downregulated between melanoma and NSCLC). (d) Commonly enriched pathways activated by nivolumab treatment in melanoma and NSCLC. Ranked gene lists were generated by long fold change on-therapy for CHECKMATE-038 and CHECKMATE-153 and then used for GSEA analysis of the REACTOME database. Pathways with positive NES scores for both CHECKMATE-038 and CHECKMATE-153 demonstrate common pathways increased on-therapy (d), while those with common negative NES scores demonstrate shared pathways that are downregulated on-therapy (e). (f) Overall survival analysis of n = 58 patients split by median SNV load pre-treatment and on-treatment, and indel-derived neoantigen load pre-treatment and on-treatment. (log rank p-value = 0.85, p-value = 0.082, p-value = 0.24, p-value = 0.044 respectively). Change in SNV-derived neoantigens as predicted by NetMHC rank<2 was not significantly different (log rank p-value = 0.18). (g) Tumor purity estimates from average VAF for paired n = 36 patients separated by response (Wilcoxon test). All boxplots represent the median, two hinges, and two whiskers, with the lower and upper hinges corresponding to the first and third quartile and the whiskers extending to the 1.5*IQR.

Extended Data Fig. 4 ∣. Transcriptome analysis reveals differential gene expression programs and immune processes in nivolumab-treated NSCLC versus melanoma.

(a) Reactome pathway analysis identifies four distinct groups of pathways for CHECKMATE-153 and CHECKMATE-038 studies combined. Full Reactome database GSEA analysis based on differential expression of on- versus pre-treatment comparing CHECKMATE-038 (melanoma) and CHECKMATE-153 (NSCLC). Differential expression of RNAseq data was performed on CHECKMATE-153 and CHECKMATE-038 independently to identify changes in gene expression on-therapy. NES scores are plotted in a heatmap. (b) The top region of the heatmap represents commonly downregulated pathways, and those common pathways with the most negative NES scores are shown. (c, d) CHECKMATE-153-specific pathways and then CHECKMATE-038-specific pathways are shown with divergent NES scores. (e) Gene expression programs commonly enriched.

Extended Data Fig. 5 ∣. Clonal architecture, tumor evolution, and nivolumab treatment.

(a) Clonal and sub-clonal abundance analysis of overall survival. Categorizing variants into clonal (CCF > = 0.95) and sub-clonal (CCF < 0.95) was used to assess survival difference pre- and on-therapy (log-rank test), where we note that on-therapy clonal abundance was associated with poor outcomes (log-rank test, p-value = 0.033). (b) Clonal and sub-clonal variants were compared across RECIST criteria – CR/PR (n = 8), SD (n = 13), PD (n = 15) – by the percent of total variants per sample both pre-therapy and on-therapy (Wilcoxon test). On-therapy changes in clonal and sub-clonal abundance were noted with reduced clonal and increased sub-clonal abundance in responders on-therapy. (c, d) Clonal CCF is increased in non-responders. Relative clonal CCF measured by PhylogicNDT was graphed for pt1643 and pt1683 demonstrating reduced tumor clonal abundance of clonal populations in pt1683. (e) There are no differences in the number of clones with CCF > = 50 across response categories. Using PhylogicNDT clusters, the average number of clones with CCF > = 50 is shown across response groups CR/PR (n = 8), SD (n = 13), PD (n = 15) (p = 0.82, Kruskal-Wallis test). (f) Peptides selection criteria for screening. Outline of peptide selection for study, where neoantigens were derived from those variants either reduced, gained, or maintained on-therapy. (g) Neoantigens for screens include those in gain, loss, and kept categories (gain >10% CCF increase, loss >10% CCF decreased, kept = all others). Scatterplot of those variants selected for immunogenicity screening with associated variant allele frequency (VAF) pre-therapy and on-therapy. (h) Table showing number of variants per patient that were kept, gained, or lost. (i) Outline of the INDEL peptides per patient that were selected for immunogenicity screening. (j) Breakdown of the 8 HLAs tested and which HLA was identified as a potential binder for those peptides. All boxplots represent the median, two hinges, and two whiskers, with the lower and upper hinges corresponding to the first and third quartile and the whiskers extending to the 1.5*IQR.

Extended Data Fig. 6 ∣. Large-scale neopeptide screening process in nivolumab-treated NSCLC patients, related to Fig. 3.

(a) There is no significant difference in binders across response groups. Neopeptides screened for immunogenicity that bound MHC (binders) shown as percentage of peptides tested per patient. Results shown by response category CP/PR (n = 6), SD (n = 3), PD (n = 5). (b) Delta CCF is not different between identified non-binders and binders. Neopeptides tested for immunogenicity that did not bind MHC (non-binder, n = 951) and those that did bind MHC (binders, n = 502) plotted by their associated Delta CCF on-therapy p = 0.18. (c) Delta CCF is not changed between binders and non-binders based on clonality. Clonal and sub-clonal Delta CCF split by MHC binding status. (d) binders are more strongly reduced in responders compared to non-responders. Reduced CCF on-therapy in responders’ neoantigens that were binders (p = 0.024, two-sided T-test). Data from non-responders (PD) and responders (PR/SD) plotted by binding status. (e) The percentage of Tetramer+ in response categories CR/PR (n = 6), SD (n = 3), PD (n = 5). Increase in Tetramer+ frequency in CR/PR compared to SD. Plot of Tetramer+ variants showing the percentage of total binders per patient in each response category. (f) Clonal variants with increased Tetramer+ T cells on-therapy are more strongly selected against on-therapy (p = 0.061, two-sided T-test). (g) Variants associated with clusters containing Tetramer+ T cells are decreased on-therapy, preferentially in responders (PR/SD). Plots showing change in CCF on-therapy. Clusters with mean CCF < 0.5 were evaluated for how variants on-therapy changed depending on whether they belonged to a cluster with or without a Tetramer+ variant in the cluster. Representative flow cytometry analysis of Tetramer screen showing positive control peptides from CMV and a peptide identified as Tetramer+ that was increased on-therapy. (h) Example representative flow cytometry Tetramer staining pre-therapy and on-therapy for patient 1695 showing two positive neoantigens with low frequency. (i) Example representative flow Tetramer staining of positive control peptides for CMV FLU peptides for patients 1615 and 1668 pre- and on-therapy. (j). Control peptides (n = 35) summary showing control peptide frequency does not change on-therapy for those which stained positive pre-therapy. The log2 (percentage of Tetramer+ T cells) for CMV, EBV, and Flu peptides which were used as positive controls separated by pre-therapy and on-therapy for all patients (p = 0.3, Wilcoxon test). (k). Change in the log2 (percentage) of T cells that are Tetramer+ for control peptides related to CMV, EBV and Flu, separated by response groups PR (n = 19), SD (n = 4), PD (n = 12) (p = 0.16, 0.88, 0.79, Wilcoxon test). All boxplots represent the median, two hinges, and two whiskers, with the lower and upper hinges corresponding to the first and third quartile and the whiskers extending to the 1.5*IQR.

Extended Data Fig. 7 ∣. Somatic mutations enriched in immunoreactive clones that contain experimentally validated immunogenic peptides.

Here, we utilized the clonal reconstructions of each patient with matched pre- and on-therapy data. Clones from patients with an experimentally validated peptide were quantified for gene overlap – that is, how many clones a non-synonymous mutation in a particular gene occurred. (a) Using PhylogicNDT clonal reconstructions, we mapped immunogenic neoantigens back to each clone and quantified the total number of times a gene showed up an immunogenic neoantigen-containing clone. This data was used for a waterfall plot ranked by the mean CCF of the clones from which somatic mutations were derived. (b) Based on results above in panel (a), we focused only on those immunoreactive clones that were reduced on-therapy and show that many somatic mutations only appear once, but there are a handful of mutations – including TTN, RYR2, and KRAS – which occur at higher frequencies in these clones that are reduced on-therapy and are known to be targeted by the immune system.

Extended Data Fig. 8 ∣. Diversity of physiochemical properties of neopeptides.

(a) Known factors that influence MHC I binding are significantly different in univariate comparisons on binders and non-binders. The comparison of factors that influence MHC I binding shows significant differences in univariate comparisons between binders (n = 437) and non-binders (n = 816). non-binders (pink) to binders (green) from left to right – MT rank (p < 2.2e-16, two-sided Wilcoxon), hydrophobicity (p = 0.0011, two-sided Wilcoxon), agretopicity (p = 0.042, two-sided Wilcoxon), aromaticity (p = 0.0015, two-sided Wilcoxon), and Grantham’s distance (p = 0.38, two-sided Wilcoxon). (b) Factors that influence HLA binding are not significantly different in univariate comparisons on Tetramer+ (n = 171) and Tetramer− (266) peptides. Comparing Tetramer+ (Blue) and Tetramer− (green) from left to right – MT rank (p = 0.00041, two-sided Wilcoxon), hydrophobicity (p = 0.27, two-sided Wilcoxon), agretopicity (p = 0.9, two-sided Wilcoxon), aromaticity (p = 0.56, two-sided Wilcoxon), and Grantham’s distance (p = 0.23, two-sided Wilcoxon). (b) Biochemical properties of Tetramer− and Tetramer+ peptides show differences in positions 3–8. Logos plots for 9mers from Tetramer+ and Tetramer− neopeptides highlighting biochemical properties of amino acid residues and their abundance in the peptide pools. (c) Amino acid composition of all Tetramer+ peptides across all HLA types screened compared to a background amino acid composition of Tetramer+ peptides with enrichment shown in purple. (d) Amino acid composition of all Tetramer− peptides across all HLA types screened compared to a background amino acid composition of Tetramer− peptides with enrichment shown in purple. (e) Ratio of amino acid composition of (Tetramer+ peptides)/(Tetramer− peptides) (purple= enriched in Tetramer+, yellow= enriched in Tetramer−). All boxplots represent the median, two hinges, and two whiskers, with the lower and upper hinges corresponding to the first and third quartile and the whiskers extending to the 1.5*IQR.

Extended Data Fig. 9 ∣. HLA- A*02:01-specific features of immunogenicity.

(a) HLA -A*02:01-screened peptides were derived from n = 4 patients; n = 109 non-binders, n = 108 Tetramer−, n = 44 Tetramer+. (b) Amino acid composition of HLA -A*02:01 Tetramer+ neopeptides compared to a background amino acid composition of HLA -A*02:01 Tetramer− peptides with enrichment shown in purple. (c) HLA -A*02:01 peptides were analyzed via iFeature python package and ranked for their relationship immunogenicity. Stacked bar charts to the (right) of the ranked waterfall plot show the distribution of features and which sequence position they are linked to for the top, middle, and bottom features. The lowest ranked features are those with the highest association to immunogenicity, primarily containing features of peptides at positions 4 and 5. (d) Univariate analysis of iFeatures results comparing Tetramer+ versus Tetramer− peptides was performed separating features into two groups – those which had significant difference of p < 0.01 and those with no significant difference p > 0.01. Significant features were primarily associated with position 4 of peptides. (e) Using HLA -A*02:01 peptides and a Tetramer+ cutoff of > 20 cells, univariate testing of all features from the amino acid index identified sequence position 4 JOND920102 (relative mutability) (p = 0.018, Wilcoxon test). (f) The relative mutability index as detailed in iFeatures’ use of the amino acid index, along with a diagram of features we have highlighted for immunogenic peptides of HLA -A*02:01 (red=anchor residues, green=high mutability at position 4, dashed lines= TCR contact positions. (g) Using Tetramer definition of > 10 cells, HLA -A*02:01 screened peptides are plotted for their mutability scores at position 4 (non-binder vs Tetramer− p = 0.28, non-binder vs Tetramer+ p = 0.026, Tetramer− vs Tetramer+ p = 0.093, Wilcoxon tests). (h) MT ranks from NetMHCpan analysis are plotted for HLA -A*02:01-screened peptides (non-binder vs Tetramer− p = 0.49, non-binder vs Tetramer+ p = 9.4e-5, Tetramer− vs Tetramer+ p = 5.9e-12, Wilcoxon tests). (i) Development of a logistic regression model identifying Tetramer+ T cells as the response variable, using the predictors sequence position 4 mutability and MT rank (intercept = −1.35484, SeqPos4JOND920102 coefficient=0.01575, MT rank coefficient = −0.22194). (J) HLA-A*02:01 combined model improves identification of Tetramer+ peptides. Model analysis on test data of HLA -A*02:01 screened peptides with a Tetramer cutoff of >10 cells (non-binder vs Tetramer− p = 0.000084, non-binder vs Tetramer+ p = 0.0035, Tetramer− vs Tetramer+ p = 0.05, Wilcoxon tests). All boxplots represent the median, two hinges, and two whiskers, with the lower and upper hinges corresponding to the first and third quartile and the whiskers extending to the 1.5*IQR.

Fig. 1 ∣. Neoantigen screening design and genomic features of NSCLC cohort.

a, WES was performed on n = 58 pre-therapy samples and n = 42 on-therapy samples, which includes n = 36 cases with matched pre-/on-therapy data. After neoantigen prediction n = 14 patients underwent immunogenicity screening, and the subsequent data were used to model features driving immunogenicity. b, Mutational landscape pre-therapy of common somatic mutations in NSCLC (n = 58). Top 3 mutational signature scores by DeconstructSig (SBS1, SBS4, SBS87). Frequency of clonal/sub-clonal variants (defined by ≥0.95 blue bar CCF or <0.95 CCF red bar). Fraction copy number change per samples represented as bar plot (bottom panel). del, deletion; ins, insertion; NSQ, nonsquamous;SQ, squamous.

Fig. 2 ∣. Genomic and microenvironmental sculpting in NSCLC tumors during nivolumab immunotherapy.

a, Whole transcriptome RNAseq of n = 24 pre-treatment biopsies were deconvoluted using ssGSEA and then compared pre-treatment between responders n = 8 and non-responders n = 16 (two-sided t-test, only significant P values shown on plot above comparison). b, GSEA Gene Oncology (GO) analysis identifies responders as having increased lymphocyte activation on-therapy. Differential expression of all RNAseq samples on-therapy versus pre-therapy was performed (total differentially expressed genes (DEGs) n = 145, non-responders DEGs n = 52, responders DEGs n = 247 and cutoff of P < 0.01 were used for functional gene set enrichment via GO pathway analysis). c, Paired sample analysis of WES data from n = 36 samples was used to calculate the change in SNVs. d, OS split by median change in SNVs and INDELS from pre-treatment, on treatment. The high group (red) represents patients that had higher than median and blue lower than median for the given plot (n = 36, log-rank test). PFS presented with the high group (red) represents patients that had an increase in SNV abundance (n = 36, P = 0.046, log-rank test). e, Mutational expression analysis comparing on-/pre-therapy expression categorized by their DNA-based clonal changes (contraction, persistence, expansion) (Wilcoxon test). All box plots represent the median, two hinges and two whiskers, with the lower and upper hinges corresponding to the first and third quartile and the whiskers extending to 1.5 times the interquartile range (IQR). aDC, activated dendritic cell; APM1, antigen processing machinery 1; APM2, antigen processing machinery 2; DC, dendritic cell; iDC, interstitial dendritic cell; pDC, plasmacytoid dendritic cell; IFNG, interferon gamma; T_CM, central memory T cell; T_EM, effector memory T cells; T_FH, T follicular helper cell; Tgd, gamma delta T cell; T_H1, T helper 1 cell; T_H17, T helper 17 cell; T_H2, T helper 2 cell; T_reg, regulatory T cell; P.adjust, adjusted P value; CPM, counts per million.

Fig. 3 ∣. Nivolumab-dependent clonal evolution and response.

a, Genomic contraction is increased in tumors with clinical benefit (PR and SD) from nivolumab. Mutation categories were defined as follows: mutation with similar CCFs in both pre-therapy and on-therapy samples (genomic persistence), increased CCF or novel in on-therapy samples (genomic expansion) and decreased CCF/lost in on-therapy samples (genomic contraction). The net change was calculated as the fraction of mutations undergoing genomic contraction subtracted from the fraction of mutations undergoing genomic persistence. Further genomic characteristics are represented in the heat map under the waterfall plot scaled for each data field (PFS months, OS months, SNV load, number of SNVs undergoing contraction, number of SNVs undergoing expansion, number of SNVs undergoing persistence). b, Genomic contraction is associated with improved OS (P = 0.059, log-rank test). c,d, Responders have increased clonal contraction on therapy (c) and nonresponders have increased clonal expansion on therapy (d). Clonal reconstruction of pre-therapy and on-therapy paired samples of a responding patient (1683) and patient with PD (1643). From left to right, clonal hierarchy shown via phylogenic tree, changes in CCF over time per clone or subclone and changes in CCF of variants per clone from pre-therapy to on-therapy samples. e, Density histograms of CCF values from pretreatment and on-treatment samples from partial responders, patients with PD and patients with SD. f, The delta CCF of variants that belonged to clusters with a mean CCF ≥ 50 and < 50 was graphed for each RECIST response group PR (n = 8), SD (n = 13) and PD (n = 15; two-sided Wilcoxon test). All box plots represent the median, two hinges and two whiskers, with the lower and upper hinges corresponding to the first and third quartiles and the whiskers extending to 1.5 times the IQR.

Fig. 4 ∣. Neoantigen immunogenicity in patients with NSCLC treated with nivolumab.

a, Peptides were screened against patient-specific MHC class I and, if they bind (binders), tetramers were made with fluorescent tags to stain peripheral blood mononuclear cells and were analyzed via flow cytometry with a dual tetramer-positivity assay. b, Representative flow cytometry with dual tetramer staining of a patient at multiple time points (before therapy, 3 weeks on therapy and 4 weeks on therapy). Drawn circle highlights positive cells for visualization. The percentage of total T cells is shown in the upper corner of each plot. c, Shown are data per patient within the neoantigen screen and the percentage of tetramer+ peptides of all predicted candidate neoantigens. d, Tetramer+ T cells are increased on therapy. Tetramer+ peptides were grouped by those that had increasing T cells from pre-therapy to on-therapy samples versus those that had decreasing peptides from pre-therapy to on-therapy. e, The log₂ fold change of tetramer+ T cells (on-/pre-therapy) separated by RECIST response PR (n = 106), SD (n = 27) and PD (n = 63; two-sided Wilcoxon test). f, Tetramer+ variants with increased T cells (n = 131) on therapy had a larger reduction in CCF compared to tetramer+ variants with decreased T cells (n = 65) on therapy (P = 0.049, t-test). g, Heat map showing all neopeptides screened for T cell immunogenicity. Labeled color bars indicate clinical response and T cell change, as well as the percentage of non-binders, tetramer− peptides and Tetramer+ peptides. h, Analysis of shared public immunogenic neopeptides identified KRAS-derived peptides that were common across patients. i, Variants that result in the KRAS-identified immunogenic peptides were analyzed for their CCF before therapy and on therapy separated by RECIST response group PD (n = 4), PR (n = 4) and SD (n = 4; paired t-test). All box plots represent the median, two hinges and two whiskers, with the lower and upper hinges corresponding to the first and third quartiles and the whiskers extending to 1.5 times the IQR.

Fig. 5 ∣. T cell recognition of neoantigens, clonal evolution and immunotherapeutic response.

a, Frequencies of binders and non-binders were grouped by response and compared by Fisher’s test (P < 0.0001). b, Tetramer+ and tetramer− peptides categorized by response groups and tested by Fisher’s test (P = 0.0245). Abs, absolute; Pos., position; coef, coefficient. c, Distribution of tetramer+ variants CCF before therapy and on therapy. d, Change in CCF (delta CCF) from pre-therapy to on-therapy samples was analyzed for tetramer− (n = 406) and tetramer+ (n = 196; P = 0.00077, two-sided t-test) variants. e, Comparing the delta CCF for tetramer+ and tetramer− peptides in clonal (CCF ≥ 0.95) and sub-clonal (CCF < 0.95; P = 0.000043, two-sided t-test) clusters. f, Comparing delta CCF of tetramer+ and tetramer− variants between nonresponders (PD) and responders (PR and SD) identified a significant reduction in CCF of tetramer+ variants in responders (P = 0.00078, two-sided t-test). g,h, Clonal reconstruction of a patient with progressive disease (g) and with responsive disease (h) using PhylogicNDT for clustering variants with the tetramer+ peptides shown as (blue) dots. i, The delta CCF for those variants belonging to tetramer+ clusters and those clusters without tetramer+ variants are graphed for clusters with mean CCF < 0.5 and mean CCF > 0.5 (P < 2.2 × 10⁻¹⁶, two-sided t-test). j, Delta CCF was significantly different in variants that belong to clusters that have a tetramer+ variant (PR P < 2.2 × 10⁻¹⁶, two-sided t-test). All box plots represent the median, two hinges and two whiskers, with the lower and upper hinges corresponding to the first and third quartiles and the whiskers extending to 1.5 times the IQR.

Fig. 6 ∣. Physiochemical properties in experimentally validated immunogenic neoantigens.

Immunogenic peptide feature extraction was performed for all n = 1,253 9mers in our experimentally validated neoantigen cohort and subsequently combined with known features of immunogenicity in a lasso regression model. a, Importance scores based on beta coefficients of the cross-validated lasso model are shown. b, Summary of amino acid index features from the lasso model with the absolute value of beta coefficients >0 summarized by frequency. c, Test set validation on n = 252 held-out peptides scored using our feature identification model, showing higher scores in immunogenic peptides (Wilcoxon, P = 4.1 × 10⁻⁸). d, Secondary validation on n = 258 peptides from the TESLA consortium dataset with matched HLA types, and 9mer peptides were scored using our feature identification model (Wilcoxon, P = 0.049). e, Validation cohort of n = 593 peptides curated from literature with matched HLA types and 9mer peptides scored using our feature identification model (Wilcoxon, P = 0.023). Utilizing all n = 58 pre-therapy samples from our study, we scored all possible 9mer peptides generated from nonsynonymous mutations, >140,000 peptides, identifying those most likely to be immunogenic. f, Variants associated with highly ranked peptides were more likely to be reduced on therapy graphed by RECIST response groups PR (n = 8), SD (n = 13) and PD (n = 15). All box plots represent the median, two hinges and two whiskers, with the lower and upper hinges corresponding to the first and third quartiles and the whiskers extending to 1.5 times the IQR.