Construction and validation of an immunoediting-based optimized neoantigen load (ioTNL) model to predict the response and prognosis of immune checkpoint therapy in various cancers

Constructing the neoantigen quantitative model with immune elimination capability

We developed ioTNL, an algorithm model that evaluates the concentration of neoantigens in tumor clones with immune elimination capability (Supplementary Figure 1). If a tumor clone acquired immune escape and tolerance to immunoediting neoantigens, it would not contribute to the tumor immunogenicity and would be discarded while calculating the neoantigen load. Only the tumor clones with immune clearance ability would be retained (Figure 1A). The calculation of the ioTNL comprises several steps. First, the cluster of each tumor clone is identified based on the allele frequency and corresponding ploidy of the detected mutations. Second, the immune editing score of each clone is calculated. The tumor was heterogeneous with different clones, and the immune stages of clones were various. We used the ratio of the neoantigen to non-silent mutation to calculate the immune editing score. A tumor clone with a high immune editing score was considered to be in the immune escape stage, otherwise, in the immune elimination stage. Subsequently, the neoantigen load indexes of the immune elimination clones were estimated. Neoantigens are identified as 8 to 11 amino acids for class I MHC-binding peptides, which arise from non-silent mutations. A previous study reported that insertions and deletions (InDels) generate a higher number of neoantigens than missense mutations. Based on the hypothesis that InDel mutations may be more immunogenic than missense mutations, neoantigens derived from InDels were also included in the ioTNL model [9]. Therefore, the design of the ioTNL model is consistent with the idea that truncal neoantigens in clonal have higher immunogenicity and correlated more with the response than subclonal neoantigens [2]. Finally, the sum of the cancer cell fraction of each neoantigen in the immune elimination clones was calculated as the ioTNL score.

Figure 1. The model of immunoediting-based optimized neoantigen load (ioTNL). (A) Illustration showing the concept of ioTNL model. Left, due to the instability of the genome, the tumor produces different clones. If tumor clones acquired immune escape and immune tolerance to neoantigens, they would not contribute to the tumor immunogenicity and would be discarded. Only the tumor clones with immune clearance ability would be retained. Right, we demonstrated a hypothetical tumor of three clones. Tumor cells from clone 1 and clone 2 were immune-eliminated cells while tumor cells from clone 3 were immunoedited cells. We assumed that this tumor had five neoantigens so that the tumor neoantigen load was equal to 5. However, neoantigens from immunoedited cells were excluded so that the ioTNL score of this tumor was equal to 1.1. (B) Distribution of ioTNL score in the NSCLC cohort. Patients with objective response (ORR) were labeled in blue while patients with non-objective response (NOR) were labeled in cyan. The scores of ioTNL were transformed into log10 format. (C) Boxplot of the distribution of ioTNL scores between patients with ORR and NOR. (D) Sensitivity and specificity of ioTNL in predicting ORR and DCB in NSCLC cohort. (E) Barplots of ORR rate (left) and DCB rate (right) between ioTNL-H group and ioTNL-L group.

Better response to checkpoint inhibitors is associated with higher ioTNL in multiple cancers

To investigate the implications of the ioTNL in the response to checkpoint inhibitor treatment, we considered a non-small cell lung cancer (NSCLC) cohort from the study by Fang et al. [10]. The cohort included 65 patients who were treated with second-line anti-PD-1/PD-L1 therapy. Previous studies revealed that the tumor mutation burden (TMB) was associated by whole-exome sequencing (WES), and the ioTNL score of each patient was calculated (Supplementary Table 2, Figure 1B). We found that patients with an objective response rate (ORR; [CR/PR]) had significantly higher ioTNL than those with a no objective response (NOR, SD/PD; p = 0.045, Figure 1C), as well as in a durable clinical benefit (DCB, p = 0.019, Supplementary Figure 2A). The ROC curve showed that the AUCs were 0.676 and 0.689 in predicting the ORR and DCB, respectively (Figure 1D). Patients were divided into ioTNL high and low groups based on the Youden index (the maximum sensitivity and the best specificity under the AUC curve). Patients with high ioTNL had a higher ORR (45% vs. 11%, p = 0.007) and DCB (55% vs. 15.6%, p = 0.002) than those with low ioTNL (Figure 1E).

We analysed three more clinical cohorts of different cancer types treated with immune checkpoint inhibitors (ICIs) to further validate the prediction efficacy of the ioTNL. The three cohorts were a skin cutaneous melanoma (SKCM) cohort of 64 patients who underwent nivolumab treatment (Riaz et al.), [11] an intrahepatic cholangiocarcinoma (ICC) cohort of 17 patients (Feng et al.) [12] and a nasopharyngeal carcinoma (NPC) cohort of 61 patients (Fang et al.) [13] The ioTNL scores of 33 patients who received nivolumab as the second-line therapy were evaluated in the SKCM cohort (Supplementary Table 3). A relatively higher median ioTNL score was observed in patients who achieved an ORR (p = 0.058, Figure 2A, left). Moreover, patients in the ioTNL-high group, determined by the Youden index described above, had a significant ORR (41.7% vs. 5.6%, p = 0.0256, Figure 2B, left) compared with those in the ioTNL-low group. A similar result was achieved in the second validation cohort of 61 patients with advanced NPC who received anti-PD-1 therapy alone or in combination with chemotherapy. In contrast to most anti-PD-1-responsive solid tumors, NPC has a modest mutation burden (Supplementary Table 4). Patients with a high ioTNL had a significantly better ORR (p = 0.0025, Figure 2B, middle) and DCB (p = 0.008, Supplementary Figure 2B). Furthermore, we collected the data of 17 patients with ICC who received the combination of nivolumab and chemotherapy (Supplementary Table 5). Patients who achieved CR or PR showed a tendency of higher ioTNL than those who did not (p = 0.139, Figure 2A, right); however, the tendency was not significant due to the small cohort size. In addition, patients in the ioTNL-high group demonstrated a higher ORR than those in the ioTNL-low group (80% vs. 33.3%, p = 0.131, Figure 2B, right).

Figure 2. Validation of ioTNL in multiple datasets. (A) Boxplots of the distribution of ioTNL scores between patients with ORR and NOR in the SKCM cohort (left), the NPC cohort (middle) and the ICC cohort (right). (B) Barplots of ORR rate between ioTNL-H group and ioTNL-L group in the SKCM cohort (left), the NPC cohort (middle) and the ICC cohort (right).

Higher ioTNL predicts a more favourable prognosis in multiple cancers

We examined the association between the ioTNL and prognosis in the NSCLC, SKCM, NPC and ICC cohorts to understand the implications of immunoediting heterogeneity on patient survival. We speculated that patients with a low immune editing score, the neoantigens in which were considered to be in the immune elimination stage, are more likely to benefit from checkpoint inhibitor therapy. As expected, we observed that higher ioTNL scores had a significantly positive association with better survival in these patients. In the NSCLC cohort, progression-free survival (PFS) in patients with a high ioTNL was longer than in those with a low ioTNL (median 161 days vs. 61 days, P < 0.001, Figure 3A). In the melanoma cohort, the overall survival (OS) was prolonged (p = 0.011, Figure 3B) in patients with a high ioTNL than in those with a low ioTNL. Analysis of the NPC cohort also indicated that the OS was significantly longer in the ioTNL-high than in the ioTNL-low group (median 105 days vs. 52.5 days, p = 0.047, Figure 3C). Analysis of the ICC cohort revealed a longer PFS among patients in the ioTNL-high group (p = 0.025, Figure 3D). The hazard ratios of above four datasets revealed high ioTNL level was a positive factor for prognosis (Figure 3E).

Figure 3. Prognosis analysis of ioTNL in multiple datasets. Kaplan-Meier estimator was used to visualize the survival of ioTNL-H group and ioTNL-L group in the NSCLC cohort (A), the SKCM cohort (B), the NPC cohort (C) and the ICC cohort (D). (E) Forest plot for the hazard ratios of ioTNL in four cohorts. Larger boxes indicate statistical significance.

Association between the ioTNL score and other immunotherapy biomarkers

We analysed the association between TMB as well as TNL and clinical outcomes to compare the predictive performance of ioTNL with known genomic biomarkers of immunotherapy (Supplementary Figure 3). The TMB could significantly predict response only in the NSCLC cohort (p = 0.01, Supplementary Figure 3A, left). Similarly, TNL could also predict response only in the NSCLC cohort (p = 0.0475, Supplementary Figure 3A, right). In the other three cohorts, patients with a high TMB or a high TNL had a tendency of better response; however, the outcomes were not significant (Supplementary Figure 3B–3D). We used the ROC curve to measure the predictive efficiency in the four cohorts. We found that the AUC of the ioTNL was greater than that of the TMB and TNL in all four cohorts (Figure 4A–4D, left). Following this, the association between the TMB as well as TNL and survival were analyzed. TMB could predict survival in two cohorts, namely, NSCLC (P < 0.001, Figure 4A, middle) and SKCM (p = 0.027, Figure 4B, middle). Similarly, TNL could predict survival in the NSCLC (p = 0.0037, Figure 4A, right) and the SKCM cohorts (p = 0.03, Figure 4B, right). However, the ioTNL could predict survival in all four cohorts and had more significant p values (Figure 3A–3D). The association between biomarkers of the microenvironment and the ioTNL was compared in the SKCM cohort. Immune checkpoint genes (PD-1, PD-L1 and CTLA4), immune cell abundance (CD8+ T cell and CD4+ T cell) as well as immune signatures (cytolytic, IFN-gamma and T cell exhaustion signatures) were identified from RNA transcriptome data. None of these microenvironment biomarkers correlated with the ioTNL (Spearman correlation <0.2, Figure 4E). Taken together, our data demonstrated that the ioTNL is a better genomic biomarker than the TMB and TNL in predicting the efficacy of immunotherapy. In addition, the ioTNL is independent of the microenvironment biomarkers.

Figure 4. Association between ioTNL score and other immunotherapy biomarkers. Comparison of the response prediction among TMB, TNL and ioTNL in the NSCLC cohort (A, left), the SKCM cohort (B, left), the NPC cohort (C, left) and the ICC cohort (D, left). Kaplan-Meier analysis of patient survival of TMB and TNL in the NSCLC cohort (A, middle and right), the SKCM cohort (B, middle and right), the NPC cohort (C, middle and right) and the ICC cohort (D, middle and right). (E) Spearman correlation of ioTNL with genomic biomarkers and microenvironment-related biomarkers in SKCM cohort. Positive correlations are labeled in blue and negative correlation are labeled in red. Size of circle represented the value of correlation.

Panel-based immune editing score correlated with the TMB and ORR of immunotherapy in different cancers

Compared with WES, target sequencing on an appropriate region size of the panel would be a cost-effective alternative for clinical detection. To facilitate the application of immunoediting status to clinical detection and the assessment of the ioTNL, we designed YuceOne, a 1.4M panel assembled by screening genomic regions and HLA genes that are susceptible to neoantigens. We collected the data of 2,469 patients across 20 cancer types and applied them on YuceOne to detect mutations and predict neoantigens (Supplementary Table 6). The same approach was followed to calculate the immune editing scores in each clone cluster for each patient. We demonstrated the landscape of immune editing scores across 20 cancer types and found that the majority of the clones had scores ranging from 0 to 10, with a median value of approximately 1 (Figure 5A). Carcinomas such as LUSC, SKCM, PRAD, BLCA and STAD had lower median immune editing scores than 1, implying stronger immune clearance in the clones. On the other hand, carcinomas such as SARC, LIHC and THYM had higher median immune editing scores than 1, implying more severe immune escape in the clones (Figure 5A).

Figure 5. Landscape of immune editing score and application of ioTNL on panel-based Yuce cohort. (A) Landscape of immune editing score in 20 cancer types. The immune editing scores were scale in log10 on y-axis. 20 cancer types were label on x-axis and arranged in ascending order by their median immune editing score. (B) Boxplots of the distribution of ioTNL scores between patients with DCB and NDB (left), also with ORR and NOR (right). (C) Barplots of DCB rate (left) and ORR rate (right) between ioTNL-H group and ioTNL-L group. (D) Kaplan-Meier analysis of patient progression-free survival between ioTNL-H group and ioTNL-L group. (E) Comparison of sensitivity and specificity between ioTNL and TMB in predicting patient DCB (left) and ORR (right).

A previous study reported a significant correlation between the median TMB and the ORR following anti–PD-(L)1 therapy [1]. To explore whether the immune editing score has a similar correlation with the TMB and ORR, we collected the median TMB and response data of ICI monotherapy for each cancer type. We plotted the median TMB and the ORR for anti–PD-(L)1 therapy against the corresponding median immune editing score across multiple cancer types (Supplementary Table 7). There was a significant negative correlation between the immune editing score and the TMB (R = −0.54, p = 0.014, Supplementary Figure 4A). Carcinomas with a high TMB such as SKCM, LUSC and BLCA had lower median immune editing scores, indicating stronger immune clearance in the tumor clones. We also observed a similar tendency of negative correlation between the immune editing score and clinical response; however, the outcomes were not significant (R = −0.45, p = 0.12, Supplementary Figure 4B). The lower the median immune editing score in one cancer type, the higher the ORR for its anti–PD-(L)1 therapy.

ioTNL identified by the target panel predicted the clinical outcome of immunotherapy

The datasets used for evaluating and verifying the ioTNL were derived from WES. To expand the application of the ioTNL model and explore whether it can predict the response to immunotherapy on a panel level, we collected the data of 21 patients with NSCLC with complete clinical outcomes who underwent anti-PD-1/PD-L1 therapy and identified the ioTNL using YuceOne. Consistent with the previously reported findings, the ioTNL scores could successfully predict the response to immunotherapy (Supplementary Table 8). Patients with a DCB and ORR showed a higher ioTNL (Figure 5B). Patients in the ioTNL-high group, had a significant DCB (100% vs. 21%, p = 0.001, Figure 5C) and ORR (57% vs. 7%, p = 0.025, Figure 5D) compared with those in the ioTNL-low group. Meanwhile, PFS analysis revealed that patients with higher ioTNL had a longer survival time (P < 0.001, Figure 5D). For comparison, the TMB and TNL were also evaluated in our cohort (Supplementary Table 2), and the ROC curves showed that the ioTNL performed better in predicting the response to immunotherapy (Figure 5E).

Mutation calling

After removing the reads containing sequencing adapters and low-quality reads with more than five ambiguous bases, high-quality reads were aligned to the NCBI human reference genome (hg19) using BWA (v0.5.9) [20] with default parameters. Picard (v1.54) [21] (http://picard.sourceforge.net/) was used to mark duplicates, followed by the use of the Genome Analysis Toolkit [22] (v1.0.6076, GATK IndelRealigner) to improve alignment accuracy.

Somatic SNVs were detected by VarScan2.2.5 [23] based on the BWA align algorithm and high-confidence somatic SNVs were called if the following criteria were met: (I) both the tumor and normal samples were covered sufficiently (≥ 10×) at the genomic position; (II) the variants were supported by at least 5% of the total reads in the tumor in contrast to less than 2% in the normal; (III) the variants were supported by at least three reads in the tumor; (IV) the distance between adjacent somatic SNVs were over 10 bp; (V) the mapping qualities of reads supporting the mutant allele in the tumor were significantly higher than 30 (Wilcoxon rank sum test, P < 0.2); (VI) the base qualities of reads supporting the mutant allele in the tumor were significantly higher than 20 (Wilcoxon rank sum test, P < 0.05); (VII) the mutations were not enriched within the 5 bp 5’ or 3’ of the read end (Wilcoxon rank sum test, P < 0.1); (VIII) the mutant allele frequency changes between the tumor and blood were statistically significant (Fisher’s exact test, P < 0.05).

High-confidence somatic InDels were called based on the following steps: (I) candidate somatic InDels were predicted with the GATK somatic InDel detector with default parameters; (II) for each predicted somatic InDel, local realignment was performed with combined normal and tumor bam files; (III) frequent of variant reads less than 10% were filtered. (IV) high-confidence somatic InDels were defined after filtering germline events.

Finally, all mutations were re-annotated using in-house annotation software based on SnpEff [24].

Copy number analysis and tumor purity assessment

Copy number variants (CNV) were called by exome-wide profile comparisons between tumors and matched peripheral blood using CNVkit (v0.8.1). Allele-specific copy number and tumor purity of the tumor genome were assessed using ascatNgs (v3.1.0).

HLA class I neoantigen predictions, prioritization, and selection

HLA genotype is identified with combination use of Polysover [25] and BWA-HLA. 9-mer polypeptides centered on mutated residues were scanned to identify candidate peptides binding to HLA class I, i.e., peptide sequences surrounding mutated amino acids resulting from missense mutations, frame-shift or non-frame-shift indels. The affinity of 8–11 mer peptides binding to HLA class I was predicted using the NetMHCPan4.0 [26] binding algorithm.

Tumor clone determination

PyClone (v0.13.0) was used to infer the cancer cell fraction (CCF) of mutations in tumors. Major copy number and minor copy number of each mutation were acquired from the result of ascatNgs software. Moreover, set prior to major_copy_number. For some samples whose copy number of each SNV and tumor purity information were not accessible, those parameters were set to the default value and set prior to total_copy_number. PyClone was deployed with 10,000 iterations and a burn-in of 1000 for all samples.

Calculation of immune editing score and ioTNL score

ioTNL was calculated by counting neoantigen concentration selected from tumor clones which could induce immune elimination. We calculated the ratio of neoantigen to non-silent mutation as immune editing score to quantify the immune editing stages for tumor clones.

$$Immuneeditingscore=\frac{n_{neoantigen}}{n_{nonsilentmutation}}$$

Tumor clonal architecture was calculated by PyClone. Cancer cell fraction of each neoantigen in the immune elimination clone were sum up as final ioTNL.

$$ioTNL={\displaystyle \sum _{i=1}^{n}loa{d}_i\times CC{F}_i}$$

Where n represents the total number of immune elimination clone, load_i is the number of neoantigens in the immune elimination clone i and CCF_i is the cancer cell fraction in the immune elimination clone i. We found the best cutoff of immune editing score from 0.5 to 1.5 in increment of 0.1 when the best prediction performance of ioTNL was found. The best prediction performance of ioTNL was found as maximum area of ROC curve. The cutoff of ioTNL was identified by Youden index (the maximum sensitivity and the best specificity under the AUC curve).

Statistical analysis

Categorical variables were evaluated with Fisher’s exact tests. Unpaired Mann–Whitney U test was used to compare differences for continuous variables between groups. Correlation analysis was assessed by Pearson coefficient. ROC analysis was done using the ROCR package in R. Significance of overall survival (OS) and progression-free survival (PFS) was determined via Kaplan-Meier analysis with log-rank analysis. The hazard ratio was calculated by the coxph function of the survival package in R. All statistical analysis was performed in the R statistical environment version 3.6.1. All tests were two-tailed and p-value < 0.05 was considered significant.

Ethics approval and consent to participate

A 68-year-old squamous cell carcinoma of the lung cancer patient01 was enrolled in Shanghai Tenth People’s Hospital on 2018. The study was approved by Shanghai Tenth People’s Hospital ethics committee. The patient provided written informed consent.

Availability of data and material

All of the published data used in the article can be found on the NCBI website and supplementary material, all the material used in the experiment can be bought. Codes and YuceOne (panel-based) data are available upon reasonable request.