Development of cancer-associated fibroblast-related gene signature for predicting the survival and immunotherapy response in lung adenocarcinoma

Data acquisition and preprocessing

The transcriptome and somatic mutational profiles and clinical data of LUAD were extracted from TCGA-GDC (https://portal.gdc.cancer.gov/). GSE41271, with clinical information obtained from GEO, was treated as the external testing group. The immune-related functional gene sets were downloaded from GSEA (http://www.gsea-msigdb.org/gsea/index.jsp). The transcriptome profiling and clinical variables of the IMvigor210 cohort were curated from a freely available software and data package that can be downloaded from http://research-pub.gene.com/IMvigor210CoreBiologies. Raw count data were first transformed into the TPM value. Our study did not require approval from the ethics committee, as it used open-access data retrieved from the TCGA or GEO database.

Abundance of tumor-infiltrating immune cells

The MCPcounter algorithm was utilized to quantify 11 immune/stromal cells (T cells, CD8⁺ T cells, cytotoxic lymphocytes, NK cells, B lymphocytes, monocytes, macrophage, myeloid dendritic cells, neutrophils, endothelial cells, and CAFs) [22, 23]. The fractions of intratumoral immune/stromal cells using distinct algorithms, including MCPcounter, were obtained from TIMER online (http://timer.cistrome.org/) [24]. In addition, the proportions of the 16 immune cell types were evaluated via the CIBERSORT algorithm [25].

Selection of hub modules related to immune or stromal cells

The R package of WGCNA was employed to construct the weight co-expression network based on the gene expression matrix of TCGA-LUAD. Briefly, a similarity matrix containing the expression levels of paired genes was generated and then converted into the adjacency matrix based on the adjacency between the paired genes. Parameter β was then defined and used to construct a weighted proximity matrix that matched the gene distribution. We performed hierarchical clustering with the minimum size to build dynamic trees. Module eigengenes (MEs) and clustering heatmaps were used to characterize the hub modules. The correlation analysis between the module genes and fraction of immune/stromal cells was calculated using Pearson’s test [26, 27]. Through this, the hub gene modules related to the infiltrating immune/stromal cells were identified.

Construction and validation of CAFRS

A total of 457 LUAD patients from TCGA were randomly partitioned into the training and testing cohorts using the R “caret” package, according to the 1:1 ratio. Then, the candidate genes related to the immune/stromal cells were incorporated into the univariate Cox regression to identify the survival-related genes, followed by a LASSO regression analysis to exclude the overfitting genes. Based on the screened survival-related genes, a multivariate Cox regression analysis using forward and backward regression analyses was employed to construct the predictive model in the training cohort. The risk score of each patient was calculated using the following formula: Σgene expression level × regression coefficient.

LUAD patients from TCGA training, TCGA testing, and GEO cohorts were divided into high- and low-risk groups using the median value of the TCGA training cohort as the cutoff. Kaplan–Meier survival curves between high- and low-risk subgroups and time-dependent receiver operating characteristic (ROC) curves in three independent cohorts were plotted using R packages of “survival,” “survminer,” and “timeROC”. Moreover, the risk score and survival status of each individual were shown by the risk curves, patients’ survival scatter plots, and heatmaps.

Mutation landscape, gene set variation analysis (GSVA), and single sample gene set enrichment analysis (ssGSEA)

The R “maftools” package was used to analyze the mutation landscape of the top 15 genes and Epidermal growth factor receptor (EGFR) between the high- and low-risk subgroups. The survival difference of the subgroups stratified by the TP53 status and risk scores were computed. Referring to the hallmark gene set, the expression matrix of all genes was converted into the enrichment score of key oncogenic pathways/phenotypes using the R “gsva” package (method = “gsva”). Then, differentially expressed scores of the pathways/phenotypes were determined between the high- and low-risk subgroups using the R “limma” package. Similarly, the gene expression matrix was respectively converted into the proportions of immune cells and scores of immune functions using the R “gsva” package (method = “ssgsea”) and analyzed between the low- and high-risk subgroups using the Wilcox.test and presented as bar plots.

Prediction of the clinical response to immunotherapy or chemotherapy

Charoentong et al. [28] created a scoring scheme based on the gene sets of immune effector/suppressor cells, immune checkpoints, and the major histocompatibility complex to quantify the clinical response to immunotherapy, which was termed the immunophenoscore (IPS). In our study, the IPS of each patient was calculated. Furthermore, the following independent dataset was analyzed in our study: advanced urothelial carcinoma treated with anti-PD-L1 (IMvigor210) [29]. The proportion of different immunotherapy responses, including the complete response (CR), partial response (PR), stable disease (SD), and progressive disease (PD), were compared between high- and low-risk subgroups. Meanwhile, Kaplan–Meier curves of the overall survival (OS) stratified by CAFRS and ROC curves were delineated to present the predictive performance in OS and response to immunotherapy.

To identify the potential compounds that benefit from CAFRS, we calculated the half maximal inhibitory concentration (IC50) of compounds obtained from the Genomics of Drug Sensitivity Cancer website (https://www.cancerrxgene.org/) in LUAD patients using the R “pRRophetic” package and compared the IC50 values of compounds in the high- and low-risk subgroups. A low IC50 indicated a sensitive drug response.

Construction and validation of a nomogram for OS

Based on the independent prognostic variables, including risk signature and other clinical variables (P < 0.05), a novel nomogram was constructed for LUAD patients in the TCGA group using the R package “rms.” Calibration curves for the survival probability of 1, 3, and 5 years were plotted to ascertain the precision and accuracy of the synthesis model. The nomogram’s diagnostic capacity compared to other clinical variables was tested via the area under the curve (AUC) value of the ROC curves.

Expression verification of hub genes

A total of eight paired LUAD tumor and adjacent nontumor specimens were collected from the lung cancer center in West China Hospital. All patients provided their informed consent, and this work was approved by the Institutional Research Ethics Committee of West China Hospital. A quantitative real-time PCR (qRT-PCR) was employed to detect the mRNA expression of four candidate genes (SERPINE2, FGF10, LSAMP, PDE5A). Briefly, the total RNA was extracted from paired lung cancer tissues using TRIzol reagent (Invitrogen, USA). cDNA was synthesized from the total RNA using reverse transcription kits (TaKaRa, Japan). qRT-PCR was performed according to the manufacturer’s protocol using the TB Green Premix Ex Taq (TaKaRa, Japan). The primer sequences were obtained from Primer Bank (https://pga.mgh.harvard.edu/primerbank). The forward primer sequences of SERPINE2, FGF10, LSAMP, and PDE5A were as follows: 5`-TGGTGATGAGATACGGCGTAA-3`, 5`-CATGTGCGGAGCTACAATCAC-3`, 5`-AGAGTTCAGCCGGATCGGAA-3`, and 5`-GCAGAGTCCTCGTGCAGATAA-3`. The corresponding reverse primer sequences were as follows: 5`-GTTA GCCACTGTCACAATGTCTTT-3`, 5`-CAGGATGCTGTACGGGCAG-3`, 5`-CGTGCCTCGGTTAAAATCCAC-3`, and 5`-GTCTAAGAGGCCGGTCAAATTC-3`. The product length of these primers was 101 bp, 138 bp, 105 bp, and 83 bp, respectively.

Statistical analysis

All data were processed under R language (Version 3.6.1), and all statistical methods were conducted using the corresponding R packages. The Wilcox test was applied to compare the nonnormally distributed numerical variables, such as the gene expression level, proportion of immune/stromal cells, and scores of calculated pathways/phenotypes. When P < 0.05, the results were considered statistically different.

Construction of the weighted co-expression network and identification of the hub module

Through data preprocessing, a total of 16, 306 genes with large fluctuation (SD > 0.8) and 425 out of 477 patients were selected for further analysis. The fraction of intratumoral immune/stromal cells were quantified using the MCPcounter algorithm. Using Pearson’s correlation coefficient, we performed a co-expression analysis between the above candidate genes and infiltration of the immune/stromal cells. The threshold soft power was calculated, and value 5 was regarded as the optimal value to construct the hierarchical clustering tree (Figure 1A, 1B). A total of 17 modules were merged into 11 modules based on the module-separating threshold of 0.25 (Figure 1C, 1D). Finally, we chose the pink module (closely related to CAFs) containing 580 genes for the subsequent predictive model construction due to the highest correlation coefficient (R = 0.82, P <0.001) with CAFs (Figure 1E).

Determination of the CAFRS and validation of the predictive capability

The TCGA cohorts were randomized into training (n = 229) and internal validation (n = 228) subgroups. Clinical features of these two cohorts and GEO-GSE41271 cohorts were listed in Table 1. The predictive model was first constructed in the training cohort. A total of 580 CAFs-related genes were first incorporated into the univariate Cox regression analysis, and 31 survival-related genes were determined. To avoid overfitting of the predictive model, the survival-related genes were assessed using a LASSO regression analysis (Figure 2A, 2B), and 19 genes were screened out for the subsequent multivariate Cox regression analysis. Finally, seven genes (FGF10, SERINE2, LSAMP, STXBP5, PDE5A, GLI2, FRMD6) were identified and used to establish a predictive model (Figure 2C). The risk scores of CAFRS were calculated as follows: risk score = -0.420 × FGF10 + 0.199 × SERINE2 + -0.441 × LSAMP + 0.423 × STXBP5 + -0.527 × PDE5A + 0.466 × GLI2 + 0.331 × FRMD6. In addition, Kaplan–Meier curves showed that patients with a high expression of SERINE2, STXBP5, GLI2, and FRMD6 had a poor prognosis, while those with a low expression of FGF10, LSAMP, and PDE5A had a favorable prognosis (Supplementary Figure 1A–1G). The ROC curves illustrated that the three-year AUC of CAFRS reached 0.694, which is greater than that of each gene (Supplementary Figure 1H).

To validate the accuracy and reliability of the predictive model, patients from three independent cohorts were divided into low- and high-risk subgroup using the median value as the cutoff. The Kaplan–Meier survival curves found that patients with high-risk scores had a poorer outcome than those with low-risk scores, regardless of whether the TCGA training, TCGA testing, or GEO cohort was used (Figure 2D–2F). The ROC analysis revealed that the AUC values of 1-, 3-, and 5-year OS in the TCGA training cohort were 0.708, 0.785, and 0.720, respectively (Figure 2G), while the AUC values were 0.625, 0.607, and 0.642 for the 1-, 3-, and 5-year OS, respectively, in the TCGA testing cohort (Figure 2H). In the GEO cohort, the AUC values of the 1-, 3-, and 5-year OS were 0.598, 0.604, and 0.643, respectively (Figure 2I). The distribution of the patients` risk scores, survival status, and expression heatmaps of the prognostic genes in these three cohorts are shown in Supplementary Figure 1I–1K.

In addition, the riskscore was positively correlated with the expressions of poor prognostic genes-FRMD6, GLI2, SERPINE2 and STXBP5, but negatively correlated with the expressions of favorable prognostic genes-LSAMP and PDEA5 (Figure 3A). The expression profiling of the seven genes composed of the predictive model demonstrated distinct clustering which corresponds to low- and high-risk subgroup (Figure 3B).

Mutation landscape associated with CAFRS

The top 15 genes with high mutation rates (TP53, TTN, MUC16, RYR2, CSMD3, LRP1B, ZFHX4, USH2A, KRAS, XIRP2, FLG, SPTA1, NAV3, ZNF536, COL11A1) combined with EGFR were displayed as waterfall plots in both the high- and low-risk subgroups. Patients with mutations of these genes accounted for 94.32% and 88.24% in the high- and low-risk subgroups, respectively. Of note, the mutation rate of TP53 was 59% and 38% in the high- and low-risk subgroups, respectively, whereas the mutation rate of EGFR was 10% and 15% in the high- and low-risk subgroups, respectively (Figure 3C, 3D). As illustrated by the Kaplan–Meier curves, patients with a TP53 mutation had a shorter OS than those with the TP53 wild type (Figure 3E). Moreover, a combination of the risk score and TP53 mutation status could improve the predictive performance for OS because patients with the TP53 mutation type and high-risk scores had the worst outcome, while patients with the TP53 wild type and low-risk scores had a survival advantage (Figure 3F).

Clinical characteristics associated with CAFRS

As illustrated in Figure 4A, the high- and low-risk clustering was remarkably associated with the pathologic TNM stage and survival status. In particular, the proportion of patients with stage III+IV was higher. in the high-risk group than that in the low-risk group, while the reverse was true for patients with stage I+II Patients with stage III+IV had a remarkably higher riskscore than those with stage I+II (P<0.01, Figure 4B, 4C). Similarly, the proportion of patients dead was higher in the high-risk group than that in the low-risk group, while the reverse was true for patients alive. Patients dead had a remarkably higher riskscore than those alive (P<0.001, Figure 4D, 4E).

Phenotype characteristics associated with CAFRS

To explore the underlying mechanism by which CAFRS affects carcinogenesis and progression of LUAD, we analyzed the association of the risk score with the pathways/phenotypes. Through a differential GSVA score analysis between the low- and high-risk subgroups, 16 of 50 hallmark gene sets were observed to be remarkably altered (P < 0.01). In particular, the G2M checkpoint, epithelial-mesenchymal transition (EMT), hypoxia, glycolysis, PI3K-Akt signaling pathways were greatly enriched in the high-risk subgroup, while the fatty acid and bile acid metabolisms were enriched in the low-risk subgroup (Figure 5A). In addition, “Toll like receptor signaling pathway”, “Small cell lung cancer”, “pathways in cancer” and “Cytokine-cytokine receptor interaction” were significantly activated in the high-risk group when Kyoto Encyclopedia of Genes and Genomes (KEGG) was taken as the reference gene set (Figure 5B).

We quantified the proportion of intratumoral immune cells using the CIBERSORT algorithm and scored the immune functions using the ssGSEA method. Based on the matrix of immune cells and immune functions, the differential analysis between the high- and low-risk subgroups was performed using the Wilcox.test, and the results were presented as bar plots. As shown in Figure 5C, 5D, the high-risk subgroup was characterized by decreased infiltration of the CD4⁺ T cell memory resting and enrichment score of type II IFN response (P < 0.001), which indicated an immunosuppressive TME.

The role of CAFRS in predicting the benefits of immunotherapy and chemotherapy

Immunotherapy, represented by anti-PD-L1/anti-CTLA-4, has caused a breakthrough in antitumor treatment. In addition to PD-L1, TMB, and MSI-H/dMMR, the IPS was strongly recommended to predict the response to immunotherapy. Herein, the IPS was significantly higher in the low-risk subgroup compared with the high-risk subgroup (P < 0.001), indicating that patients in the low-risk subgroup were more susceptible to immunotherapy (Figure 6A). Furthermore, patients in the IMvigor210 cohort were divided into high- and low-risk subgroups using the median value of the CAFRS risk score. Patients in the low-risk subgroup exhibited a marked clinical response advantage and prolonged survival rate. The proportions of CR, PR, SD, and PD were 10%, 19%, 25%, and 46% in the low-risk subgroup, respectively, compared with the proportions of 7%, 10%, 18%, and 65% in the high-risk subgroup (Figure 6C). Similarly, the risk score of patients with CR was significantly lower than that of patients with PR, SD, and PD (Figure 6D). When the clinical response was dichotomized into response (CR + PR) and non-response (SD + PD) groups, the ratio of response to non-response was 29/71 in the low-risk subgroup. However, it fell to 17/83 in the high-risk subgroup. The risk score of the non-response group was significantly higher than that of the response group (Supplementary Figure 2A, 2B). In addition, the Kaplan–Meier curves revealed that patients in the low-risk subgroup survived longer than those in the high-risk subgroup (Figure 6E), which is in agreement with the results obtained from three independent cohorts. The AUC values of the ROC curves at 6 months, 12 months, and 18 months were 0.588, 0.642, and 0.624, respectively (Figure 6F), suggesting that this CAFRS can predict the prognosis of LUAD. Moreover, the established risk score is a poor prognostic factor for OS, independent of TMB, TNM stage, sex, chemotherapy, and immune subtype in the IMvigor210 cohort (Supplementary Figure 2C).

We investigated the differences between the IC50 values for common drugs between the low- and high-risk subgroups. As shown in Figure 6G–6I, the CAFRS score was markedly negatively associated with the IC50 values for paclitaxel and 5-fluorouracil (common chemotherapeutic agents) and positively associated with the IC50 values for erlotinib (common EGFR tyrosine kinase inhibitor [TKI]). Moreover, the low-risk subgroup presented remarkably higher IC50 values for paclitaxel and 5-fluorouraci but had lower IC50 values for erlotinib (Supplementary Figure 2D–2F). This indicated that LUAD patients with low-risk scores were more sensitive to erlotinib therapy but benefited little from conventional chemotherapy, which is likely due to drug resistance. In summary, CAFRS could predict anticancer drug responses in patients with LUAD.

Nomogram construction combined with clinical characteristics and CAFRS

The univariate and multivariable Cox regression analyses indicated that the CAFRS and TNM stages were independent prognostic factors for LUAD (Figure 7A). To explore a more accurate and reliable prediction tool, we established a nomogram to predict the survival probability of 1-, 3- and 5-year OS for LUAD (Figure 7B). During the 1-, 3- and 5-year OS, the calibration curves of the nomogram exhibited excellent concordance with the actual survival rate in two independent cohorts (TCGA-LUAD, GSE41271) (Figure 7C, 7D). The 5-year survival rate of stages I+II and the low-risk subgroup in the TCGA and GEO cohorts was up to 58.8% and 82.9%. However, the number of stages III+IV and the high-risk subgroup decreased, reaching 30% and 21.3%, respectively (Figure 7E, 7F). Moreover, compared with the single TNM stage or CRAFS, the merged nomogram had a better predictive performance in OS. The 5-year AUC value of the nomogram reached 0.742 and 0.702 in the TCGA and GEO cohorts, respectively, which is greater than that of the TNM stage and CAFRS (Figure 7G, 7H).

mRNA verification of differential expression of CAFRS genes

As shown in Figure 8A, four out of seven genes in the CAFRS were differentially expressed between the cancer and paratumor tissues from LUAD patients in TCGA. To prove the reliability of the CAFRS, the mRNA levels of these four genes were detected using a qRT-PCR. Consistent with the transcriptome data, SERPINE2 was more highly expressed in lung cancer tissues compared with that in paratumor tissues, while a low expression of FGF10, LSAMP, and PDE5A was found in lung cancer tissue (Figure 8B–8E).