A Novel Tool for the Risk Assessment and Personalized Chemo-/Immunotherapy Response Prediction of Adenocarcinoma and Squamous Cell Carcinoma Lung Cancer

Introduction

The prevalence and cancer-specific death rate of lung cancer (LC) have risen in recent decades. According to the data from GLOBOCAN, approximately 2.09 million patients were newly diagnosed with LC in 2018, and approximately 1.76 million patients died of LC in 2018.¹ The prevalence of LC is increasing in developing regions, as expected, and is mostly associated with an increasing number of smokers.² The 5-year survival of LC patients is as low as 10% to 20% around the world, despite the development of clinical treatments, including surgery, chemotherapy, and immunotherapy.³ From the perspective of histology, non-small-cell LC (NSCLC) and small-cell LC (SCLC) are the two major subtypes. Moreover, two important typical subtypes are separated because of the heterogeneous background of NSCLC: Adenocarcinoma LC (LUAD) and squamous cell carcinoma LC (LUSC).

The molecular characteristics of LUAD and LUSC have been well described. Comprehensive molecular profiling of LUAD was performed by the Cancer Genome Atlas Research Network and revealed three main subtypes: Terminal respiratory unit (TRU), proximal proliferative (PP), and proximal inflammatory (PI). The TRU subtype of LUAD is characterized by the features of EGFR mutation and has a preferable prognosis. The PP subtype is characterized by enrichment of KRAS mutations, while PI subtype patients are characterized by commutation of both the TP53 and NF1 genes.⁴ Four molecular subtypes were also defined in LUSC patients: Primitive, classical, basal, and secretory. The primitive subtype has the activated function of cellular proliferation, the classical subtype exhibits the distinctive functional theme of xenobiotic metabolism, the secretory subtype contains the major distinctive function of the immune response, and the basal subtype is characterized by the adhesion function.⁵ Immune infiltration in the tumour microenvironment (TME) has been reported to be closely associated with tumorigenesis and prognosis in several tumours. Recently, a series of studies reported an immune molecular subtype among tumours, including LC. Patients with an activated immune environment showed the best clinical outcome, while the immunosuppressed subgroup presented the worst clinical outcome.^6–9

Several studies have also tried to construct the prognostic signature of LC. Li et al¹⁰ reported a 25 immune gene pair-based prognostic signature based on 2414 NSCLC patients, with a predicted C-index value of 0.64. Liu et al¹¹ generated two signatures based on autophagy-associated genes to separately predict the prognosis of LUAD and LUSC, but these signatures lacked universal applicability for both LUAD and LUSC patients. In the current study, we explored the consensus prognostic genes in both LUAD and LUSC and established a universally applicable prognostic signature for most LC patients.

Methods and Materials

Patients Summary and mRNA Expression Profile

A total of 453 patients from The Cancer Genome Atlas (TCGA)-LUAD and 452 patients from TCGA-LUSC cohorts were first enrolled for the analysis. All patients had gene expression profiles and matched clinical information. All these data were downloaded by the R package “TCGAbiolinks”.¹² Importantly, patients with an overall survival (OS) time of less than one month were excluded to avoid potential bias. For the gene expression profile, genes with zero values in more than 10% of samples were also excluded. The count data were transformed to transcripts per kilobase million values (TPM) and then with the final form after log2 (TPM+1) transformation for subsequent analysis.¹³ For the clinical information, we separated patients into nonsmokers who had never smoked, ever smokers, who had smoked in the past, and smokers, who currently smoked.

Construction the Prognostic Signature

Univariate Cox analysis was applied to explore the prognostic genes in the LUAD and LUSC cohorts. Then, a Venn plot was selected and displayed the consensus prognostic genes for subsequent analysis. The least absolute shrinkage and selection operator (LASSO) regression analysis was conducted by the “glmnet” package¹⁴ and was employed to select and regularize all the input factors to generate a statistical model with high prediction accuracy and interpretability. As a kind of regularization method, LASSO can avoid overfitting, can more effectively perform feature selection and optimize the form of the regression model than other common ordinary least-squares subset selection methods and is especially suited for situations where there are numerous potential predictors but where we think only a few are likely relevant. These genes selected via LASSO analysis were used to calculate the risk score of each patient with the gene expression value and the coefficient.

where C_i is the coefficient of the gene, and E_i is the expression value of the corresponding gene.

Validation of the Prognostic Value of the Signature

To validate the prognostic value of the established signature in external lung cancer cohorts, we used online tools. In SurvExpress,¹⁵ the gene symbols of 23 genes pointed out from the LASSO analysis were input for the analysis, the risk score was calculated automatically, and the deficiency of no more than five genes was allowed. Another online website, Kaplan–Meier Plotter,¹⁶ was also used. A total of 1144 lung cancer patients, including both LUAD and LUSC patients, and the clinical information on sex, smoking status and tumour stage were also recorded. The risk score for the analysis in Kaplan–Meier Plotter was calculated along with the index of the 23 selected genes.

Gene Set Enrichment Analysis

To investigate the mechanism of tumorigenesis of lung cancer, we used gene set enrichment analysis (GSEA) to explore the activated signalling pathways in high-risk and low-risk patients in both the LUAD and LUSC cohorts.¹⁷ The 186 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were downloaded from MSigDB (https://www.gsea-msigdb.org/gsea/msigdb/genesets.jsp?collection=CP:KEGG). Enrichment analysis is a means to characterize biological attributes in the given KEGG gene sets. We calculated an enrichment score (ES) for each gene set by walking down the rank list and finding the maximum deviation from zero of a running-sum, weighted Kolmogorov–Smirnov-like statistic. Next, the normalized enrichment score (NES) was calculated by dividing positive and negative ES by the mean of positive or negative pES, respectively. Finally, to estimate the false discovery rate (FDR), a null distribution of NES values is generated using a list of background gene sets.¹⁸

Immune Infiltration Evaluation

The infiltration of 28 types of immunocytes into tumours was evaluated via single-sample gene set enrichment analysis (ssGSEA).¹⁹ For each patient, the ssGSEA gives a score representing the infiltration of each type of immunocyte. Then, the correlation between the risk score and infiltration of 28 types of immunocytes was evaluated and displayed by a lollipop plot, and only the cell type with a P value less than 0.05 was displayed.

Personalized Treatment Prediction

Based on the drug sensitivity and phenotype data from GDSC 2016 (https://www.cancerrxgene.org/), we employed the R package “pRRophetic” to predict the chemotherapeutic sensitivity for each LC patient.²⁰ The estimated IC50 (lower IC50 indicates increased sensitivity to treatment) of each sample treated with a specific chemotherapy drug was obtained by ridge regression,²¹ and prediction accuracy was measured through 10-fold cross-validation.²² For immunotherapy, we harnessed subclass mapping to predict the clinical response to immune checkpoint blockade.¹³ In this manner, we retrieved a published dataset containing 47 patients with melanoma who responded to immunotherapies.²³

Results

Selection of the Prognostic Genes

First, we used univariate Cox regression analysis to select prognostic genes. In LUAD, we identified 1822 genes associated with poor prognosis and 752 genes linked with favourable prognosis (Figure 1A). For LUSC, a total of 308 genes were reported to promote tumorigenesis, while another 1039 genes acted as protectors (Figure 1B). To construct a prognostic model that was suitable for both LUAD and LUSC, we combined oncogenes and suppressor genes in both groups. Finally, 18 oncogenes and 29 suppressor genes were chosen for LASSO analysis (Figure 1C). The expression of these 47 genes and their association with clinical features are shown in Figure 1D–E.

Figure 1 Construct the prognostic signature PLCRC. (A) Volcano plot shows the prognostic oncogene and suppress gene among LUAD cohort; (B) Volcano plot shows the prognostic oncogene and suppress gene among LUSC cohort; (C) Venn diagram shows the intersect oncogene and suppress gene in both LUAD and LUSC cohorts; (D) Heatmap shows the 47 intersected genes distribution in LUAD cohort; (E) Heatmap shows the 47 intersected genes distribution in LUSC cohort.

Construction of the Prognostic Signature

To construct the prognostic signature based on the 47 genes, we conducted LASSO analysis. Twenty-three genes were retained in the formula, and the index was also calculated by LASSO (Figure S1A and B). The pan-lung cancer risk score (PLCRC) of each patient was calculated via the following formula: PLCRC = 0.00542 * ELK3 + 0.26072 * EMC6 + 0.03244 * BCAR1 + 0.06126 * TBC1D1 + 0.00671 * SNAI1 + 0.09366 * FLNC + 0.07842 * RAB27B + 0.0395 * FSTL3 + 0.05682 * ID1 + 0.02545 * FAM83A - 0.03249 * ANKRD65 - 0.03771 * MYLIP - 0.12505 * NRTN - 0.01225 * SNHG10 - 0.03846 * RFXAP - 0.02943 * ZNF57 - 0.20389 * MRPL54 - 0.20496 * INAFM2 - 0.01993 * TMEM60 - 0.06058 * PHKG2 - 0.02463 * C3ORF62 - 0.19104 * ZNF394 - 0.03237 * NPRL2. In the entire cohort, patients were separated into high-risk and low-risk groups via the median PLCRC value. We observed that patients in the high-risk subgroup suffered from extremely pejorative OS compared with the low-risk subgroup (p < 0.001, HR: 2.75, 95% CI: 2.131–3.546, Figure 2A), and the prognostic value of the PLCRC in all 905 patients was assessed by ROC curve analysis and showed a moderate prognostic value (AUC: 0.643, 95% CI: 0.603–0.684, Figure 2B). In addition, we evaluated the prognostic value of the PLCRC in LUAD and LUSC patients separately. In the LUAD patients, the PLCRC also presented a relatively good application (AUC: 0.636, 95% CI: 0.573–0.699, Figure 2D), and patients in the high-risk subgroup showed a 2.51-fold-change risk of poor prognosis compared with the low-risk subgroup (p < 0.001, HR: 2.51, 95% CI: 1.719–3.676, Figure 2C). For LUSC patients, the OS outcome in the two risk groups also separated significantly (p < 0.001, HR: 2.89, 95% CI: 2.038–4.088, Figure 2E), with a preferred prognostic value (AUC: 0.656, 95% CI: 0.602–0.710, Figure 2F).

Figure 2 Prognostic value of PLCRC in entire, LUAD and LUSC cohorts. (A) Kaplan-Meier (K-M) plot shows the different OS dichotomized by PLCRC in entire cohort; (B) The receiver operating characteristic (ROC) curve shows the prognostic value of PLCRC in entire cohort; (C), K-M plot shows the different OS dichotomized by PLCRC in LUAD cohort; (D) ROC curve shows the prognostic value of PLCRC in LUAD cohort; (E) K-M plot shows the different OS dichotomized by PLCRC in LUSC cohort; (F) ROC curve shows the prognostic value of PLCRC in LUSC cohort.

PLCRC is Independent Prognostic Factor for LC Patients

To further evaluate the clinical application value of the PLCRC, we performed multivariate Cox regression analysis of the PLCRC and clinical features. In the entire cohort, we revealed that the PLCRC value was an independent prognostic factor (p < 0.001, HR: 2.281, 95% CI: 2.161–3.684, Figure 3A, Table 1) after adjusting for the impact of age, sex, smoking status, and tumour stage. We also found that age higher than 70 years old and advanced tumour stage were factors impacting OS in the entire cohort. Moreover, to improve the prognostic value, we newly defined a combined model that enrolled the PLCRC and all four other clinical features. The combined model showed the best prognostic value assessed via ROC curve analysis (combined, AUC: 0.724, 95% CI: 0.670–0.778; PLCRC, AUC: 0.712, 95% CI: 0.657–0.767; sex: AUC: 0.484, 95% CI: 0.432–0.535; age, AUC: 0.527, 95% CI: 0.474–0.580; smoking, AUC: 0.490, 95% CI: 0.442–0.538; stage, AUC: 0.619, 95% CI: 0.563–0.674, Figure 3B). Subsequently, we conducted multivariate Cox regression analysis on the LUAD and LUSC subgroups. The PLCRC acted as an independent prognostic factor in LUAD patients as well (p < 0.001, HR: 2.06, 95% CI: 1.35–3.10, Figure 3C, Table 1). The combined model of the PLCRC plus clinical features showed a prognostic AUC value as high as 0.793 (95% CI: 0.708–0.877, Figure 3D), which was better than that in the entire cohort. For LUSC patients, the PLCRC identified the patients with high risk better, and patients in the high-risk group had an increased risk fold-change of 3.22 compared with those in the low-risk group (p < 0.001, HR: 3.22, 95% CI: 2.23–4.60, Figure 3E, Table 1); however, the overall combined model did not increase the prognostic value compared with the PLCRC alone (combined, AUC: 0.660, 95% CI: 0.585–734; PLCRC, AUC: 0.671, 95% CI: 0.597–0.746; Figure 3F).

Table 1 The Independent Prognostic Value of Risk Score in Entire, LUAD and LUSC Cohort

Figure 3 Independent and combined prognostic value of PLCRC. (A) Forest plot shows the independent prognostic value of PLCRC after adjusting clinical features in entire cohort; (B) The receiver operating characteristic (ROC) curve shows the comparison and combined prognostic value of PLCRC and clinical features in entire cohort; (C) Forest plot shows the independent prognostic value of PLCRC after adjusting clinical features in LUAD cohort; (D) ROC curve shows the comparison and combined prognostic value of PLCRC and clinical features in LUAD cohort; (E) Forest plot shows the independent prognostic value of PLCRC after adjusting clinical features in LUSC cohort; (F) ROC curve shows the comparison and combined prognostic value of PLCRC and clinical features in LUSC cohort; *P < 0.05; **P < 0.01; ***P < 0.001.

Validation the Prognostic Value of the PLCRC in SurvExpress

To validate the prognostic value of the PLCRC, we input the 23 genes in SurvExpress and obtained the K-M plot from high- and low-risk score groups in different databases. The modified PLCRC in nine GEO databases was calculated automatically by SurvExpress, and we separated the patients into high- and low-risk groups by the median value of the modified PLCRC to confirm that the patient number in each group was consistent. Notably, the deficiency of no more than five genes among the 23 input genes was allowed. The stable prognostic value of the PLCRC was confirmed, and patients in the modified PLCRC-defined high-risk group all showed significantly poorer clinical outcomes than those in the low-risk group in all nine GEO databases (all p < 0.001, Figure 4). The HR value in different GEO databases ranged from 2.41 to 8.91, the highest HR belonged to GSE31210, there are only 4 patients dead among the 113 patients in the low-risk group, while another 31 patients among the 113 patients in the high-risk group met the dead end during the follow-up (p < 0.001, HR: 8.91, 95% CI: 3.14–25.25, Figure 4).

Figure 4 Validation the prognostic value of PLCRC via SurvExpress. Nine cohorts collected from GEO database, the modified PLCRC was automatically calculated by SurvExpress based on the expression of 23 selected genes.

Validation of the Prognostic Value of the PLCRC in Kaplan–Meier Plotter

We also evaluated the prognostic value of the PLCRC in clinical subgroups via Kaplan–Meier Plotter. A total of 1144 patients were recorded in the Kaplan–Meier Plotter, with gene expression of 23 genes and clinical information. The PLCRC of each patient was calculated along with the index and expression of the 23 LASSO analysis identified genes. We found that the PLCRC-defined high-risk group showed a significantly severe clinical outcome (p < 0.001, HR: 1.96, 95% CI: 1.65–2.32, Figure 5A). The opposite two risk groups were identified by the PLCRC in both the male (p < 0.001, HR: 1.57, 95% CI: 1.28–1.92, Figure 5B) and female (p < 0.001, HR: 2.96, 95% CI: 2.05–4.28, Figure 5C) subgroups, demonstrating its stable prognostic value. For smoking status, we revealed that the prognostic value was useful for smokers but not nonsmokers (n=300, p < 0.001, HR: 2.21, 95% CI: 1.44–3.39, Figure 5D). Regarding tumour stage, PLCRC showed a better prognostic value in advanced stage (stage I: p = 0.094, HR: 1.48, 95% CI: 0.93–2.34, Figure 5E; stage II: p < 0.001, HR: 3.47, 95% CI: 2.43–4.95, Figure 5E).

Figure 5 Validation the prognostic value of PLCRC via Kaplan–Meier Plotter. (A) Kaplan-Meier (K-M) plot shows the different overall survival (OS) dichotomized by PLCRC in entire 1144 patients; (B) K-M plot shows the different OS dichotomized by PLCRC in male subgroup; (C) K-M plot shows the different OS dichotomized by PLCRC in female subgroup; (D) K-M plot shows the different OS dichotomized by PLCRC in smoker subgroup; (E) K-M plot shows the different OS dichotomized by PLCRC in Stage I subgroup; (E) K-M plot shows the different OS dichotomized by PLCRC in Stage II subgroup.

Pathways Derives the Prognosis of LC

To identify the pivotal signalling pathways impacting and promoting the tumorigenesis and clinical outcome of LC patients, we employed GSEA to identify the most significant pathways in LUAD and LUSC patients. We revealed that the extracellular matrix pathway and focal adhesion pathways activated in the PLCRC defined high-risk groups in LUAD and LUSC (Figure 6A). In addition, we also observed activated immune-associated signalling pathways in high-risk patients, and the cytokine and receptor interaction pathway and leukocyte transendothelial migration pathway were all significantly activated in both LUAD and LUSC high-risk patients (Figure 6B). More details of the activated signatures are listed in Table 2.

Table 2 Activated Signaling Pathways in High-Risk Patients of LUAD and LUSC

Figure 6 Signaling pathways and immune infiltration involved in LC. (A) Activated signaling pathways in the PLCRC defined high-risk LUAD patients; (B) Activated signaling pathways in the PLCRC defined high-risk LUSC patients; (C) Correlation between PLCRC value and infiltration of immunocytes in LUAD; (D) Correlation between PLCRC value and infiltration of immunocytes in LUSC.

To further understand the mechanism of how the tumour immune microenvironment impacts the tumorigenesis of LC, we compared the correlation between 28 types of immunocytes and the PLCRC value. Interestingly, we found that the PLCRC was positively associated with the infiltration of numerous immunocytes, especially neutrophils, central memory CD8 T cells, and natural killer T cells (Figure 6C and D).

Identify the Personalized Treatment for LC Patients

Because we revealed the different activated signalling pathways and discrepant infiltration of immunocytes in the high-PLCRC and low-PLCRC subgroups, we tried to search for personalized therapy for LC patients. With the help of 138 recorded drugs in the R package “pRRophetic”,²⁰ we revealed that ten drugs were more suitable for patients with a low PLCRC than for those with a high PLCRC in both LUAD (all P < 0.05, Figure 7A) and LUSC (all P < 0.05, Figure 7B) patients. We were excited to find that patients with a high PLCRC could benefit more from anti-cytotoxic T lymphocyte-associated protein 4 (CTLA4) immunotherapy, and that this was true for both LUAD (Figure 8A) and LUSC (Figure 8B) patients.

Figure 7 Precision chemotherapy drugs for LC. (A) Personalized therapeutic drugs for LUAD patients with low PLCRC; (B) Personalized therapeutic drugs for LUSC patients with low PLCRC.

Figure 8 Patients with high PLCRC are more suitable to anti-CTLA4 therapy. (A) Response results to anti-CTLA4 and anti-PD1 therapy of LUAD patients; (B) Response results to anti-CTLA4 and anti-PD1 therapy of LUAD patients.

Discussion

LC is still the most malignant tumour, with a poor 5-year survival rate of only 19%.²⁴ For the early stage of NSCLC, in which the tumour stage remains at stage I or stage II, surgery is the most frequently recommended clinical treatment.²⁵ The 5-year survival rate of stage IA patients after surgery rises to 77–92%, the survival rate of stage IB patients is 68%, and the survival rate of stage II NSCLC patients is almost 60%.²⁶ For advanced NSCLC, the clinical treatment paradigm has improved, mostly owing to the new findings of molecular targeted therapies.²⁷ Epidermal growth factor receptor (EGFR) mutations are present in approximately one-fifth of NSCLC patients, mostly in LUAD patients.²⁸ Tyrosine kinase inhibitors (TKIs) are effective for the clinical treatment of LUAD patients with EGFR mutations, with an improvement in progression-free survival compared to chemotherapy; however, the toxic effects are nonnegligible.²⁹ Moreover, secondary mutations, activation of bypass pathways, alternative activation of other pathways and histological transformation can lead to resistance to TKI treatment in LUAD patients.^30,31 In addition to EGFR mutations, researchers have also revealed genetic alterations and developed targeted drugs, including ALK gene rearrangements,³² ROS1 gene rearrangements,³³ RET fusions,³⁴ MET mutations,³⁵ BRAF mutations,³⁶ and KRAS mutations.³⁷ Although anti-PD1/PD-L1 therapy is applied for the treatment of advanced NSCLC, only 14–20% of patients have rapid and durable responses.^38–40 Therefore, it is urgent to construct a prognostic signature to identify LC patients with a high risk of poor OS survival at the early stage and choose the appropriate therapy.

For the early diagnosis of LC, the National Lung Screening Trial (NLST) provides the clinical recommendation for low-dose computed tomography of the chest (LDCT) screening. Participants were randomly assigned to radiography or LDCT and screened at baseline with two annual follow-up scans. Within the maximum follow-up time of 7 years, the LDCT group had a 20% decrease in mortality and a 6.7% reduction in all-cause mortality.⁴¹ Weiss et al⁴² reported that methylation of SHOX2 and PTGER4 in plasma DNA is more effective in identifying LC than in identifying nonmalignant lung diseases. Chen et al⁴³ constructed a five-gene signature to calculate the risk score based on the expression of DUSP6, MMD, STAT1, ERBB3 and LCK detected by RT-PCR.

In the current study, we first collected data from 453 LUAD patients and 452 LUSC patients from the TCGA database after removing patients lacking clinical records or those with OS times less than one month. The oncogenes and suppressive genes of both LUAD and LUSC were selected as the input genes for LASSO Cox regression analysis. A total of 23 genes remained for the construction of the PLCRC signature, which was calculated by the index generated from LASSO analysis and the log2 transformed TMP value of each gene. The moderate prognostic value of the PLCRC was observed in the entire cohort, LUAD cohort and LUSC cohort, and the PLCRC can dichotomize patients with significantly different OS. To assess the application of the PLCRC in patients with different backgrounds, we further performed multivariate Cox regression analysis. The most exciting finding was that the PLCRC was an independent prognostic factor of both LUAD and LUSC patients after adjusting for age, sex, smoking status, and tumour stage. Moreover, personalized therapy is necessary for patients. We revealed that the PLCRC can also act as a guideline for clinical therapy for both LUAD and LUSC. Chemotherapy was more suitable for patients with a low PLCRC, while patients with a high PLCRC can benefit more from anti-CTAL4 therapy.

The clinical application of a prognostic signature is validated in multiple external cohorts. With the help of SurvExpress, we successfully validated the prognostic value of the PLCRC in 1342 LC patients, recorded in GSE19118, GSE30219, GSE31210, GSE37745, GSE42127, GSE3141, GSE8894, GSE13213, and GSE17710. The prognostic value of the PLCRC in clinical subgroups was also validated by Kaplan–Meier Plotter. The PLCRC could dichotomize patients with significantly different OS in males, females, smokers, and especially stage II patients. We also studied the potential mechanism of tumorigenesis in LUAD and LUSC patients and linked the poor prognosis of LC patients to the abundant infiltration of immunocytes.

Conclusion

We established and validated a newly defined prognostic signature, the PLCRC, for LC patients, and this signature is universally useful and stable for the prognosis of both LUAD and LUSC patients. Patients with a low PLCRC benefit more from chemotherapy, while patients with a high PLCRC benefit more from anti-CTLA4 immunotherapy.

Consent for Publication

All authors have read the final manuscript and approve of the publication of its content in entirety.