Accurate Prediction of Immune Therapy Response in Non-Small Cell Lung Cancer

Twenty patients with non-small cell lung cancer treated with PD-1/PD-L1 inhibitors in The Second Xiangya Hospital, Central South University who had genomic profiling of whole exome sequencing (WES) before treatment were included in our GloriousMed cohort (). 2

TMB was defined as the total number of somatic mutations per exome in megabases. PD-L1 staining was evaluated centrally by IHC using 22C3 antibody and an automated staining procedure developed by Dako. The percentage of PD-L1 expression was scored by a qualified pathologist in samples with a minimum of 100 viable tumor cells.

Objective response was assessed by investigator-assessed RECIST 1.1 criteria every 6 weeks (two cycles of ICB administration). The complete response (CR) or partial response (PR) was considered as responders, whereas patients with stable disease (SD) or progressive disease (PD) were considered as non-responders.

All patients collection and usage were in accordance with the principles of the Declaration of Helsinki and approved by the Institution Review Board of The Second Xiangya Hospital, Central South University. The written informed consent for sample acquisition was obtained from all patients. All data were deidentified.

Three independent public cohorts including Hellmann cohort (25), Rizvi cohort (26), and Samstein cohort (27) were also used in this study. The data for the three independent cohorts were retrieved from published articles (Supplementary Table S2). Hellmann cohort included 75 NSCLC patients treated with combined PD-1 and CTLA-4 blockade. Rizvi cohort included 34 NSCLC patients that treated with pembrolizumab. The Samstein cohort contained 1,662 patients received immunotherapy from 11 different cancers.

DNA was extracted from FFPE-fixed tumor tissue using QIAamp DNA FFPE Tissue Kit (Qiagen), and Genomic DNA (gDNA) was extracted from white blood cells using the Blood Genomic DNA Mini Kit (Cwbiotech). Integrated DNA Technologies's xGen Exome Research Panel v1.0 according to the standard procedures (IDT) were used to capture whole exome. For each sample, 200 to 500 ng FFPE DNA or 500 ng gDNA was then used for library preparation and quantification guided by KAPA Hyper Prep protocols (KAPA). Libraries were then purified by AMPure XP (Beckman) and quantified by Qubit™ dsDNA HS Assay Kit (Thermo Fisher). Final library was sequenced on the Illumina Novoseq6000 (PE150). Sequencing adapters were trimmed by Trimmomatic from the raw data (28). The reads after adapter trimming were then aligned with the human reference genome hg19 by BWA (29). Duplicated reads were removed by Picard. Mapped reads were also realigned to the genome by Genome Analysis Tool Kit. Somatic mutations were called by Mutect2 with a paired workflow. Variants were then annotated by ANNOVAR and self-development code (30). An in-house script was used to verify the human identity concordance of paired samples. Somatic mutations were filtered with the following rules: (1) base quality value ≥20; (2) mutation reads depth ≥10; (3) variant allele frequency ≥5%; (4) reads supporting variation <4 and frequency <2% in normal, tumor abundance/normal abundance ≥8; (5) no strand bias (GATK parameter FS > 60 for SNP and FS >200 for indel); (6) discard synonymous mutations.

TMB and PD-L1 expression of Hellmann cohort and Rizvi cohort were retrieved from published articles. MATH was calculated through R package maftools for GloriousMed, Hellmann and Rizvi cohorts (31). Correlation among TMB, MATH, and PD-L1 expression (%) were examined by the Pearson rank correlation method. Correlation between TMB or MATH and grouped PD-L1 expression were examined by the Wilcoxon signed-rank test.

The overall survival (OS) was defined from the start of ICIs treatment until death due to any cause. And the progression-free survival (PFS) was defined as the time from the start of ICIs treatment until disease progression. Of notes, the Samstein cohort merely published OS data and Rizvi cohort provided PFS data. The Kaplan-Meier method was used to estimate OS or PFS, and the log-rank test was used to compare the survival curves. All tests with a p value ≤ 0.05 were considered statistically significant.

The detailed profiles of genes involved in HRR, MMR, BER, JAK, MAPK, PI3K, NF-κB, and Wnt pathways were listed in. At first, mutational status of aforementioned six immune-related pathways in every sample was classified into two categories: the first one assigned with 0 (no non-synonymous mutation) and the second with 1 (at least one non-synonymous mutation). Then, DDR pathway mutation status of each sample was classified into three groups based on the mutational status of HRR, MMR, and BER. “N” represented no mutation in HRR, MMR, or BER, “C” was stood for co-mutation between HRR and MMR or BER, and “S” was other cases. In addition, the mutational status of PI3K, JAK, and NF-κB were integrated as one variable by summing the mutational status. 4

Three models were constructed, one model with TMB, PD-L1 expression, MATH, and immune-related pathways, called “pathway-model”; a second with TMB, PD-L1 expression, and MATH, called “tri-model”; the last one, called “bivariate-model”, with TMB and PD-L1 expression (Table 1). Both TMB and MATH were z-score normalized. PD-L1 expression was stratified as 0% (Z), 1%-49% (L), ≥50% (H), or unknown (N). And immune-related pathways were processed according to Immune-Related Pathway Selection. All of the models were trained via Bayesian Regularized Neural Networks (BRNN) algorithm using corresponding variables with 2 layers and default hyperparameters from R package caret (32), and the resampling method “boot” was used to choose the optimal model. The cutoff value of single-factor variable, TMB, PD-L1 expression and MATH was estimated by BRNN algorithm as well. Fifty-three patients of the Hellmann cohort were used as the training set, and remaining 22 patients were validation set. Rizvi cohort and GloriousMed cohort were processed as above description and were used as testing cohort.

Receiver operating characteristic (ROC) curves were constructed with the predictor estimated from each of the previous models and single-factor variables with roc function of R package pROC (33). Benefit probability of each patient was extracted from prediction results, and DCB/NDB information was provided by the cohorts. Differences between DCB and NDB with benefit probability were examined by the Wilcoxon signed-rank test.

The clinical information, RNA expression, mutational status and protein array of The Cancer Genome Atlas Lung Adenocarcinoma (TCGA LUAD) and Lung Squamous Cell Carcinoma (LUSC) patients were retrieved from TCGA database. The patients with EGFR exon 18–21 mutations and ALK gene fusions were filtered to avoid make a disturbance for the analysis. In the signature score analysis, the expression of genes in a signature was normalized in the form of fragments per kilobase of exon model per million mapped fragments (FPKM). Then, a principal component analysis (PCA) was performed, and PC1 was extracted to serve as gene signature score (34). The 18 signatures and their gene sets were summarized from published papers (34–38). The significantly differential expression analysis was based on DESeq2 (39). The row counts of LUAD and LUSC patients were used as input for DESeq2. The differential expression genes were defined as the genes with absolutely log2Foldchange ≥ 1 and p-value ≤ 0.05. The oncoplot of top 30 mutated genes were drawn by using R package maftools (31).

The previous studies documented that higher TMB or PD-L1 expression correlated with better outcomes as compared with lower TMB or PD-L1 expression (11, 12, 25, 40). However, in 70 of 75 patients from Hellmann cohort who had all three biomarkers data, correlation between TMB and PD-L1 expression was not significant (R=-0.14, p-value=0.24). TMB of some patients was more than 10 but PD-L1 expression was less than 25% (Figure 1A). The results might reveal the biomarkers were not consistent in response prediction of ICIs treatment. In the meantime, the novel biomarker MATH was not significantly correlated with PD-L1 expression (R = −0.2, p-value = 0.099) or TMB (R = 0.14, p-value = 0.24) as well (Figures 1B, C). We further explored the correlation between stratified PD-L1 expression and TMB or MATH by stratifying PD-L1 expression as 0% (Z), 1% to 49% (L), ≥50% (H), and unknown (N). Neither MATH nor TMB showed a significant difference with any PD-L1 expression groups (Figures 1D, E). The Rizvi and GloriousMed cohort showed the consistent correlation results as well (Supplementary Figure 1). This lack of correlation suggested that TMB, PD-L1 expression, and MATH are independent predictive measures of response to ICIs treatment, and a robust model should be constructed to unify these variables.

A prior study has shown that co-mutation information of DNA damage response (DDR) pathway can be used as a predictor of response to immune checkpoint blockade, and the mutation of the DDR solved the problem of difficulty in determining an optimal TMB threshold (22). This finding provided a new way to predict the response of immunotherapy. Besides DDR pathway, we selected six pathways, homologous recombination repair (HRR), Janus kinase (JAK), mitogen-activated protein kinase (MAPK), phosphoinositide 3-kinase (PI3K), and nuclear factor kappa-light-chain-enhancer of activated B cells (NF-κB), Wnt, through literature survey, which are associated with tumor immunity or immunotherapy escape (41, 42). We also collected the mutational status of these pathways from Samstein cohort treated with ICIs (27) and explored its correlation with the overall survival (OS). The results showed that patients with mutations in any of six pathways had better survival than those without mutation (Figure 2). Furthermore, the results also revealed the selected pathways could be used as biomarkers to distinguish the prognosis for ICIs treatment.

We extracted 70% patients from Hellmann cohort, which totally included 75 NSCLC patients, as training data set (25) and the rest 30% patients were used to validate the models. Three different models were trained by using the training data set with different variables and were adjusted with clinical benefit as outcomes (Table 1). The pathway-model contains seven variables, including TMB, PD-L1 expression, MATH and the mutational status of six immune-related pathways (Figure 3). The mutational status of JAK, MAPK, and PI3K was integrated into one variable to improve the prediction accuracy. ROC curves based on the predictor for each of the three models estimated on Hellmann cohort (22 patients) were available and the results showed that the pathway-model was more predictive than other two models (AUC is 0.87, 0.83, and 0.59 for pathway-, tri-, and bivariate-model). The AUC of pathway-model was higher than single-factor variables containing TMB, PD-L1 expression, and MATH as well (AUC is 0.56, 0.49, and 0.69 for TMB, PD-L1 expression, and MATH) (Figure 4A and Table 2). We also checked the prediction benefit probability, a quantitative output generated from the model which represents the likelihood of immunotherapy response, of each patient compared with real clinical benefit information among three models. The benefit probability generated from pathway-model and tri-model are significantly higher in DCB group than in NDB group (p-value is 0.0024 for pathway-model and 0.0066 for tri-model), however, the median benefit probability of pathway-model (0.70) was higher than tri-model (0.46). The difference of benefit probability was not significant in other models and single factors (Figure 4B).

We further tested the predictive ability of pathway-model in Rizvi cohort (26), consisting of 34 NSCLC patients treated with pembrolizumab, with all predictive variables and clinical benefit information available. The results showed that pathway-model could more accurately predict the clinical benefit of ICIs than other two models and single-factor variables (AUC is 0.74 for pathway-model, 0.67 for tri-model, 0.68 for bivariate-model, 0.63 for TMB, 0.72 for PD-L1 expression, and 0.55 for MATH) (Figure 4C and Table 2). The benefit probability of patients in DCB and NDB groups was significantly different as well (p-value is 0.0017, Figure 4D). The survival analysis indicated that the high benefit probability group also showed a better PFS (Figure 4E).

Finally, we tested pathway-model in GloriousMed cohort with 20 NSCLC patients, who were treated by PD-1/PD-L1 inhibitors (Supplementary Table S1). The accuracy of pathway-model was much higher than tri-model and bivariate-model (AUC is 0.80 for pathway-model, 0.47 for tri-model and 0.64 for bivariate-model) (Figure 5A and Table 2). Even though, the benefit probability was not significantly different between DCB and NDB group (p-value is 0.08 for pathway-model), all DCB patients have a predictive benefit probability higher than 0.5 (Figure 5B). Thus, pathway-model can be generalized in clinical to improve the prediction accuracy of the response to immunotherapy.

We predicted the benefit probability of TCGA LUAD and TCGA LUSC cohorts without EGFR exon 18-21 mutations and ALK gene fusions patients in immunotherapy with pathway-model and classified patients to two groups at the median cut-point. Then, we calculated signature scores of 18 gene sets with principle component analysis (PCA) method. In TCGA LUAD cohort, thirteen signatures are significantly different between high benefit probability group and low probability group (Figure 6A). In consideration of TMB, and mutational status of DDR and Wnt pathways are included in prediction model, the benefit probability difference in DDR, WNT target, DNA repair–related signatures and cell cycle were expected. The signature score of CD 8 T effector and Immune Checkpoint were higher in high probability group than in that of low group, while the signature score of EMT3 and FGFR3 related was lower in high probability group (Figure 6A). However, in LUSC cohort, we did not find significant difference between high and low benefit probability groups as LUAD cohort (Figure 6D).

Furthermore, we analyzed the differential expression genes between high benefit probability groups and low group in LUAD and LUSC respectively (Figures 6B, E, Supplementary Table S4). There are 153 differential expression genes (106 up-regulated) in LUAD, including AFP and G6PC, which related to P53 downstream pathway and FOXO pathway. In LUSC, there are 120 differential expression genes (50 up-regulated) including FGF3 and DLK1, which related to FGFR pathway and NOTCH pathway. Apart from that, part of the top 30 mutated genes, such as KRAS and PTPRD, have different mutation pattern between high benefit probability group and low group, as well as between LUAD and LUSC (Figures 6C, F).

Above all, the comprehensive analysis of TCGA LUAD and LUSC cohorts imply that high benefit probability patients from pathway-model is associated with immune response.

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.