Gene patterns linked to lung adenocarcinoma outcomes and tumor environment

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,,and, which contained the microarray-based expression data of LADC patients and associated clinical information were downloaded from GEO website () via GEOquery package in R software (). All GEO datasets usedplatform. 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. GSE3141 GSE30219 GSE50081 https://www.ncbi.nlm.nih.gov/geo/↗ GPL570 [Davis & Meltzer, 2007]

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 (). 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 (). The screening criteria were |log2(Fold Change)| > 1 and FDR < 0.05. [Yoshihara et al., 2013] [Benjamini & Hochberg, 1995]

The clusterProfiler package was used for Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, and. adjust (FDR) < 0.05 was considered statistically significant (). 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/↗ p [Yu et al., 2012] [Szklarczyk et al., 2019]

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< 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 (). 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. p [Friedman, Hastie & Tibshirani, 2010]

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.

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) (). 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. [Subramanian et al., 2005]

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) (;). Then, the ESTIMATE algorithm was used to calculate corresponding immune scores, stromal scores and tumor purities (), and the difference of tumor purities, expression of HLA genes between high-risk and low-risk samples was further analyzed. [Hänzelmann, Castelo & Guinney, 2013] [He et al., 2018] [Yoshihara et al., 2013]

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 (). 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 (), and the top 10 immune DEGs and stromal DEGs were presented in heatmap, respectively (). Table S1 Fig. 1A Fig. 1B

The GO enrichment and KEGG pathway analyses of the TME-related DEGs were conducted using clusterProfiler package in R environment (). The DEGs were mainly associated with biological process (BP) of immune responses, such as negative regulation of immune system process (), regulation of inflammatory response () and regulation of leukocyte activation () (). The KEGG analysis showed that the enriched pathways were Platelet activation (hsa04611),infection (hsa05150), Phagosome (hsa04145) and Leukocyte transendothelial migration (hsa04670) (). The GO terms of the three categories, BP, CC, molecular function (MF) and KEGG results were presented inand, 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-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. GO:0002683 GO:0050727 GO:0002694 [Yu et al., 2012] Fig. 2A Fig. 2B Table S2 Table S3 Fig. 2C Staphylococcus aureus p

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

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 (= 2.359E−05) (). The area under curve (AUC) of ROC were 0.738 (). 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 (), only three-gene signature can be used as independent prognostic factor (< 0.001,). 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 (–). The AUC of ROC curves of 3 and 5 years in the testing dataset were 0.646, 0.635 (), 0.718, 0.569 () and 0.643, 0.65 () (–). These results indicated that the three-gene prognostic model was robust in predicting the outcome of LADC patients. GSE30219 GSE3141 GSE50081 p p Fig. 4A Fig. 4B Fig. 4C Fig. 4D Figs. 5A 5C Figs. 5D 5F

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 (). 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 (and). Fig. 6A Figs. 6B 6C

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(), the low-risk group in TCGA (),() and() showed more immune activities than that of high-risk group. Consistent with ssGSEA results, except(), tumor purity of low-risk groups in TCGA (),() and() 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 (–). Interestingly, although there is no significant difference in the overall immune status and tumor purity between the high- and low-risk groups in thedataset, 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. GSE50081 GSE30219 GSE3141 GSE50081 GSE30219 GSE3141 GSE50008 Fig.7D Fig.7A Fig.7B Fig.7C Fig. 8D Fig.8A Fig.8B Fig.8C Figs. 8E 8H