A total of 29 TME signatures were obtained from the previous research (Bagaev et al., 2021), including four groups of signatures, anti-tumor microenvironment (MHCⅠ, MHCⅡ, coactivation molecules, cytotoxic cells, T cells, T cells trafficking, B cells, M1 signature, NK cells, Th1 signature, and anti-tumor cytokines), pro-tumor microenvironment (Treg, Treg traffic, MDSC, MDSC traffic, neutrophil signature, granulocyte traffic, macrophages, Th2 signature, macrophages/DC traffic, and pro-tumor cytokines), angiogenesis fibrosis (angiogenesis, endothelium, cancer-associated fibroblasts (CAFs), matrix, matrix remodeling), and malignant properties (proliferation rate signature and EMT).

The enrichment score of 29 TME signatures was measured by single sample gene set enrichment analysis (ssGSEA) (Hänzelmann et al., 2013). ConsensusClusterPlus R package (Wilkerson and Hayes, 2010) was utilized to construct consensus matrix based on the ssGSEA score of TME signatures in TCGA dataset. 1 - Pearson correlation was selected as distance and KM algorithm was used for repeating 500 times of bootstraps with each bootstrap having 80% samples of TCGA dataset. The optimal cluster number (k) was selected according to the cumulative distribution function (CDF) curves and consensus matrix.

For the validation of the TME-based subtyping in the external datasets (GSE72094 and GSE50081), a support vector machine (SVM) model was used (Huang et al., 2018) (LUAD samples within TCGA dataset were randomly grouped into training set and testing set with a ratio of 7: 3).

In order to observe the methylation difference of different subtypes, we obtained the methylation data set of HumanMethylation450 from the TCGA database, extracted the methylation signals of each sample, and completed the missing values using the KNN method. Limma was used to analyze the methylation difference of each subtype (P.val< 0.05 and |logFC|>log2 (1.1)). In addition, we annotated the methylation sites to the gene promoter region to obtain genes regulated by methylation, The Biological Pathway to Obtain Methylation Disorder of Each Subtype by Function Enrichment Analysis.

Given that TME signatures were significantly related to LUAD prognosis, we attempted to identify molecular subtypes based on their enrichment scores. By using consensus clustering, we determined four molecular subtypes (C1, C2, C3, and C4) according to CDF and consensus matrix (Supplementary Figure S1). Four subtypes had distinct enrichment patterns of 29 TME signature as shown in the heatmap (Figure 2A). C3 and C4 subtypes had higher enrichment of anti-tumor signatures compared to C1 and C2. C4 subtype had evidently lower enrichment of EMT signature and tumor proliferation rate compared to C1, C2 and C3 subtypes. Principle component analysis (PCA) displayed the different distribution of four subtypes based on the TME signatures (Figure 2B). Significant differences were shown among four subtypes on the enrichment of all TME signatures (Figure 2C). In addition, we analyzed the oncogenic activity of four subtypes, and the result showed different activation of these oncogenic pathways (Figures 2D,E). PI3K and hypoxia were activated in C1, EGFR and TGF-β were activated in C2, JAK-STAT and NFκB were activated in C3, and P53 signaling was activated in C4. Moreover, survival analysis revealed that C1 subtype had the worst prognosis while C4 had the longest overall survival (p < 0.0001, Figure 2F). Different activation of these pathways may indicate different TME-based molecular mechanisms of tumor progression. Furthermore, we evaluated the distribution of different clinical characteristics in four subtypes. The result exhibited an evident trend that the samples with advanced stages were more distributed in C1 (Supplementary Figure S2). Female patients and the patients with age >60 had a higher proportion in C4 compared with that in C1-C3 (Supplementary Figure S2). Not surprisingly, C1 had the largest number of the samples with dead status than other subtypes (Supplementary Figure S1).

Genomic instability has been demonstrated to be associated with tumor development. We obtained a series of genomic characteristics of TCGA-LUAD data from a pan-cancer research (pan-cancer) (Thorsson et al., 2018). C1 subtype had relatively high scores of tumor mutation burden (TMB), aneuploidy, homologous recombination deficiency, loss of heterozygosity, purity, and ploidy, while C3 showed relatively high score of intratumor heterozygosity (Figure 3A). In the previous pan-cancer research, they identified six immune subtypes of LUAD (C1, C2, C3, C4, and C6 immune subtypes). We analyzed the distribution of previous immune subtypes in our TME-based subtypes (Figure 3B). C1 immune subtype (also known as wound healing) mostly distributed in C1 and C2 TME-based subtypes. C2 immune subtype (also known as IFN-γ dominant) mostly accumulated in C1 and C3 TME-based subtypes. C4 TME-based subtype had the highest proportion of C3 immune subtype (also known as inflammatory). The different distribution of previous immune subtypes in our TME-based subtypes also supported the distinct TME characteristics of four subtypes. In addition, we evaluated the gene mutations in four subtypes, and observed that TP53, LRP1B, and SPTA1 contributed high somatic mutation frequencies (Figure 3C). In addition, we analyzed the different DNA methylation sites of each subtype in the genome. Among them, C1 has the most differential methylation sites, C2 has only a small amount of DNA methylation differences, and the different methylation sites of each subtype overlap less (Figure 3D). Further functional analysis showed that the methylation site of C1 imbalance was mainly related to LUNG CANCERALVEOLAR CELL CARCINOMA, INCLUDED, the methylation site of C2 imbalance was mainly related to Type II interaction signaling (IFNG), the methylation site of C3 imbalance was mainly related to ER Phagosome pathway, and the methylation site of C4 imbalance was mainly related to detection of chemical stimulus involved in sensor performance of smart (Figure 3E), These results indicate that different molecular subtypes may have different apparent disorder patterns.

To further understand the different molecular mechanism of tumor development in four subtypes, we analyzed the biological pathways using GSEA. Different pathways were enriched in four subtypes. In C1 subtype, cell cycle and DNA repair-related pathways were relatively activated such as mismatch repair, base excision repair, homologous recombination, and DNA replication (Figure 4A). In C2 subtype, EMT-related pathways were significantly enriched such as ECM receptor interaction, tight junction, TGF-β signaling pathway, and focal adhesion (Figure 4B). In C3 subtype, cell cycle and immune-related pathways were activated such as DNA replication, cell cycle, cytokine-cytokine receptor interaction, antigen processing and presentation, chemokine signaling pathway, and Toll-like receptor signaling pathway (Figure 4C). In C4 subtype, immune-related pathways were also evidently activated such as chemokine signaling pathway, antigen processing and presentation, cytokine-cytokine receptor interaction, and complement and coagulation cascades (Figure 4D).

Additionally, similar results were carried out in hallmark pathways (Figure 4E). Cell cycle-related pathways were much enriched in C1 subtype. C2 subtype showed activated EMT, angiogenesis, hypoxia, Notch signaling, and TGF-β signaling. Immune-related pathways were significantly enriched in both C3 and C4 subtypes. Besides, C4 subtype also had relatively high enrichment of metabolic pathways such as heme metabolism, fatty acid metabolism, adipogenesis, xenobiotic metabolism, and bile acid metabolism. The above results implied that these differently enriched pathways may result in different TME characteristics in four subtypes.

As four TME-based subtypes showed different TME scores, prognosis and activated pathways, we then identified the DEGs among four subtypes. A total of 353 DEGs (135 upregulated and 218 downregulated) were screened in C1 vs other, 91 DEGs (9 upregulated and 82 downregulated) were screened in C2 vs other, 171 DEGs (161 upregulated and 10 downregulated) were screened in C3 vs other, and 396 DEGs (223 upregulated and 173 downregulated) were screened in C4 vs other. Functional analysis on all upregulated DEGs unveiled different biological function of upregulated DEGs in different subtypes. Cell cycle-related pathways and processes were enriched in C1; stromal-related processes were enriched in C2; immune-related pathways and processes were enriched in C3 (Supplementary Figure S3).

We screened a total of 648 DEGs among four subtypes in TCGA dataset after removing the duplicate DEGs. Univariate Cox regression analysis detected a total of 164 risk genes and 254 protective genes within 648 DEGs (p < 0.05, Supplementary Table S1). Next, we applied Lasso regression and stepAIC to dig out key prognostic genes from the above risk and protective genes. Lasso analysis compressed the coefficients to zero and remained nine prognostic genes when the lambda value = 0.0689 (Figures 7A,B). Subsequently, stepAIC was performed on the nine prognostic genes and further compress the number of genes. Consequently, five prognostic genes were remained, including PTTG1, MS4A1, ZNF750, RHOV, and KRT6A. The 5-gene prognostic model was determined as: Risk Score = 0.206*PTTG1—0.155*MS4A1—0.12*ZNF750 + 0.136*RHOV + 0.05*KRT6A.

Each LUAD sample obtained a risk score, and was stratified into high-risk and low-risk groups referring z-score = 0 as a cut-off (Figure 7D). Two risk groups exhibited different enrichment pattern of survival status, with a higher density of dead samples in high-risk groups. Five prognostic genes showed distinct expression patterns in two risk groups, where ZNF570 and MS4A1 were relatively upregulated in low-risk group, while RHOV, PTTG1, and KRT6A were upregulated in high-risk group. There were 242 and 245 LUAD samples in high- and low-risk groups respectively, and two groups exhibited distinct prognosis (p < 0.0001, Figure 7E). ROC curve analysis revealed that the risk model had favorable AUC in predicting survival at 1, 3, and 5 years with the scores of 0.75, 0.72, and 0.67 respectively (Figure 7F). Furthermore, we validated the 5-gene risk model in two independent datasets (GSE72094 and GSE50081). The validation results were consistent with the TCGA dataset (Figures 7G–J), which demonstrated the effectiveness and reliability of the risk model.

The relation of risk score to different clinical characteristics was assessed and a trend showed that the risk score was higher in the advanced stages compared with early stages (Supplementary Figure S4). A significant difference was also observed in different ages and genders. Moreover, we compared the risk score in four TME-based subtypes, and the results showed that the prognosis of subtypes was consistent with their risk levels. C1 with the worst prognosis showed the highest risk score, which was consistent with the previous results (Supplementary Figure S4; Supplementary Figure S2). In different clinical characteristics, the risk model also showed a favorable performance in dividing samples into high-risk and low-risk groups (Supplementary Figure S4).

Next, we compared the difference of two risk groups in biological pathways and immune characteristics. GSEA on KEGG pathways revealed that high-risk group had relatively activated pathways of purine metabolism, pyrimidine metabolism, citrate cycle TCA cycle, oxidative phosphorylation, and DNA replication (Figure 8A). Consider that immunity is related to the tumor, thus we determine immune characteristics. In terms of tumor microenvironment, two risk groups had distinguished infiltration levels that low-risk group had higher infiltration of both immune cells and stromal cells than high-risk group (Figure 8B). Of 22 immune cells, 11 immune cells were differentially distributed in two risk groups (Figure 8C). Low-risk group had higher infiltration of resting dendritic cells, memory B cells, resting memory CD4 T cells, and resting mast cells than high-risk group, while M0 macrophages and M1 macrophages were lower enriched in low-risk group. Supportively, risk score was significantly correlated with resting memory CD4 T cells and M0 macrophages (Figure 8D). In the relation of risk score with oncogenic pathways, EGFR, hypoxia, PI3K, and VEGF pathways were positively correlated with risk score (R = 0.42, 0.42, 0.36, and 0.28, respectively) (Figure 8E), suggesting that these pathways may be highly involved in the TME modulation and tumor progression. Moreover, risk score was strongly associated to 29 TME pathways (Figure 8F).