Identification of prognostic gene signature associated with microenvironment of lung adenocarcinoma

Data source and preprocessing

Fragments per Kilobase Million (FPKM) normalized expression profile data of LADC samples were downloaded from TCGA database using GDC data transfer tool and summarized into an expression matrix. The ensemble ids were converted into gene symbols according to the annotation file (Homo_sapiens.GRCh38.95.CRH.GTF). Then the xml formatted clinical information of LADC patients was downloaded and merged into a single matrix for further analysis. Three Gene Expression Omnibus (GEO) datasets, GSE3141, GSE30219 and GSE50081, which contained the microarray-based expression data of LADC patients and associated clinical information were downloaded from GEO website (https://www.ncbi.nlm.nih.gov/geo/) via GEOquery package in R software (Davis & Meltzer, 2007). All GEO datasets used GPL570 platform. The probe ids were converted into gene symbols according to related annotation file, and for multiple probes corresponding to the same gene, the average expression value was calculated.

Identification of differentially expressed genes related to LADC tumor microenvironment

ESTIMATE is an algorithm that infers the infiltration situation of immune cells and stromal cells in tumor tissue according to the transcriptome data of TME-related genes which contain a set of immune and stromal signature genes (Yoshihara et al., 2013). The expression data of these TME-related genes were extracted from TCGA LUAD dataset. Differentially expressed TME-related genes between LADC and normal samples were screened using Wilcoxon rank sum test. The False Discovery Rate (FDR) in multiple comparisons was controlled using Benjamini–Hochberg procedure (Benjamini & Hochberg, 1995). The screening criteria were |log2(Fold Change)| > 1 and FDR < 0.05.

Functional enrichment, KEGG and PPI network analysis

The clusterProfiler package was used for Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, and p. adjust (FDR) < 0.05 was considered statistically significant (Yu et al., 2012). Protein–protein (PPI) interactions network can visualize the patterns of molecular interactions and help to explain the mechanisms underlying phenotypes. To further explore the interactions among TME-related differentially expressed genes (DEGs), PPI network analysis was performed using the online database STRING with interaction score of 0.4 as the threshold (https://string-db.org/) (Szklarczyk et al., 2019).

The construction of LADC TME-related prognostic model

After removing the patients whose survival time was NA, the resulting 458 LADC patients in TCGA LUAD dataset were included in the Cox regression analysis. The univariate Cox model was used to determine the relationship between TME associated DEGs expression and OS. The p < 0.05 was considered statistically significant. The key genes were selected from the significant DEGs in the univariate analysis using least absolute shrinkage and selection operator (LASSO) regression analysis via glmnet package in R software (Friedman, Hastie & Tibshirani, 2010). The LASSO regression is a popular method for variable selection in fitting high-dimension generalized linear model, which can get a more refined model by constructing a penalty function to reduce the variable numbers and effectively avoid overfitting. Then, the selected key genes were included in multivariate Cox analysis, and risk score formula was constructed according the analysis results.

Validation of the gene signature for survival prediction in the testing dataset

The predictive performance of the TME-related gene signature was further validated in three GEO datasets. Samples in the testing datasets were divided into high- and low-risk group according to the formula of risk score derived from the training dataset, respectively. Kaplan–Meier (KM) survival analysis and receiver operating characteristic (ROC) curve were used to evaluate the predicting power of the gene signature, and the prognostic performance of other clinicopathological factors was also analyzed.

Function enrichment analysis of TME-related gene signature

We also identified pathways that were up- and down-regulated when the expression level of TME-related gene signature was changed by gene set enrichment analysis (GSEA 4.0) (Subramanian et al., 2005). The curated KEGG gene set was downloaded from the MSigDB database. Enrichment FDR values were based on 1,000 permutations, and FDR < 0.05 was considered to be statistical significance.

Evaluation of immune status between high-risk and low-risk groups stratified by prognostic model

To explore the potential relationship between immune system and TME-related gene signature, we analyzed the immune status of the high-risk and low-risk samples. First, using 29 immune signatures, we quantified the immune activities between high-risk and low-risk samples by single-sample gene-set enrichment analysis (ssGSEA) (Hänzelmann, Castelo & Guinney, 2013; He et al., 2018). Then, the ESTIMATE algorithm was used to calculate corresponding immune scores, stromal scores and tumor purities (Yoshihara et al., 2013), and the difference of tumor purities, expression of HLA genes between high-risk and low-risk samples was further analyzed.

Identification the TME-associated DEGs in LADC

Accumulating evidence suggests that tumor progression and patient’s prognosis were associated with the TME. To identify the prognostic LADC TME-related genes, differential expression analysis was conducted. First, the expression profiles of TME associated immune and stromal signature genes were extracted from TCGA LUAD dataset, which contained 141 immune and stromal signature genes respectively (Table S1). Then the DEGs between tumor and normal samples were identified by Wilcoxon rank sum test. Totally, 93 TME-related DEGs were identified, of which 23 DEGs were up-regulated and 70 DEGs were down-regulated (Fig. 1A), and the top 10 immune DEGs and stromal DEGs were presented in heatmap, respectively (Fig. 1B).

GO, KEGG and PPI analysis of TME-related DEGs

The GO enrichment and KEGG pathway analyses of the TME-related DEGs were conducted using clusterProfiler package in R environment (Yu et al., 2012). The DEGs were mainly associated with biological process (BP) of immune responses, such as negative regulation of immune system process (GO:0002683), regulation of inflammatory response (GO:0050727) and regulation of leukocyte activation (GO:0002694) (Fig. 2A). The KEGG analysis showed that the enriched pathways were Platelet activation (hsa04611), Staphylococcus aureus infection (hsa05150), Phagosome (hsa04145) and Leukocyte transendothelial migration (hsa04670) (Fig. 2B). The GO terms of the three categories, BP, CC, molecular function (MF) and KEGG results were presented in Table S2 and Table S3, respectively. Using STRING online tool, PPI network of LADC TME-related DEGs was constructed, which contained 93 nodes, 294 edges and average node degree was 6.32 (PPI enrichment p-value < 1.0E−16). PPI network also revealed that the stromal signature DEGs had extensively interactions with immune signature DEGs, which is consistent with previous reports that stromal cells and immune cells in the TME work together to form a microenvironment conducive to tumor growth. After removing the disconnected nodes, the PPI network was presented in Fig. 2C.

Construction and validation of TME-related gene Signature

To construct the TME-related gene signature, TCGA LUAD dataset (n = 458) were included in univariate Cox proportional hazard regression analysis and the resulting 23 significant DEGs (p < 0.05) were input in LASSO regression (Fig. 3A). Then, the seven key DEGs were selected to performed the multivariate Cox regression analysis (Fig. 3B). Finally, a prognostic model containing three genes (ADAM12, Bruton Tyrosine Kinase (BTK), ERG) was established to assess the prognosis of each patient as follows:

Risk score = 0.3619 × ADAM12 + (−0.4048) × BTK + (0.2396) × ERG. The detailed information of multivariate Cox regression was presented in Table 1.

The KM plot and ROC curve were used to evaluate the performance of three-gene signature in predicting the outcome of the LADC patients. In training dataset, the OS between the low- and high-risk groups classified by our prognostic model was significantly different (p = 2.359E−05) (Fig. 4A). The area under curve (AUC) of ROC were 0.738 (Fig. 4B). Then, the performance of our model in LADC patients were further assessed with other common prognostic factors by univariate and multivariate Cox regression analysis. Although univariate Cox analysis indicated that tumor stage, T, N stage and our model all had prognostic effect (Fig. 4C), only three-gene signature can be used as independent prognostic factor (p < 0.001, Fig. 4D). In line with the results in the training dataset, the TME-related gene signature can well stratify the samples in three GEO testing datasets as low-risk and high-risk group (Figs. 5A–5C). The AUC of ROC curves of 3 and 5 years in the testing dataset were 0.646, 0.635 (GSE30219), 0.718, 0.569 (GSE3141) and 0.643, 0.65 (GSE50081) (Figs. 5D–5F). These results indicated that the three-gene prognostic model was robust in predicting the outcome of LADC patients.

Function enrichment analysis of TME-related gene signature

To explore the underlying mechanisms of the prognostic effects of three-gene signature, GSEA enrichment analysis was performed. The results suggested that highly expressed ADAM12 correlated with the activation of pathways such as SCLC, pathway in cancer, transforming growth factor beta signaling pathway, while lowly expressed ADAM12 associated with the metabolism pathway such as butanoate metabolism, fatty acid metabolism, histidine metabolism (Fig. 6A). Highly expressed BTK and ERG mainly correlated with the immune associated pathways such as cytokine–cytokine receptor interaction, JAK-STAT signaling pathway, while lowly expressed BTK and ERG associated with DNA replication, spliceosome pathway (Figs. 6B and 6C).

Evaluation the immune status between low-risk and high-risk groups

To further explore the relationship between three-gene signature and immune system, ssGSEA method was used to assess the overall immune status of high-risk and low-risk groups by analyzing the expression profiles of the 29 immune signature genesets. The heatmap showed that in TCGA and three GEO datasets, the immune status of the low-risk and high-risk samples showed a certain degree of heterogeneity. Except GSE50081 (Fig.7D), the low-risk group in TCGA (Fig.7A), GSE30219 (Fig.7B) and GSE3141 (Fig.7C) showed more immune activities than that of high-risk group. Consistent with ssGSEA results, except GSE50081 (Fig. 8D), tumor purity of low-risk groups in TCGA (Fig.8A), GSE30219 (Fig.8B) and GSE3141 (Fig.8C) were significantly higher than that of high-risk groups, which suggested more infiltrated immune and stromal cells in the TME of low-risk samples. The HLA related genes play a key role in immune regulation. The analysis showed that the expression of key HLA genes in all four datasets was significantly higher in the low-risk group than in the high-risk group (Figs. 8E–8H). Interestingly, although there is no significant difference in the overall immune status and tumor purity between the high- and low-risk groups in the GSE50008 dataset, significant differences in the expression of HLA genes suggest that even if there is the relative same proportion of immune cells infiltration, the difference in the function of immune cells can also affect the prognosis of patients. However, the more specific mechanisms may need further research.