Gene expression data, clinical survival information, and gene mutation information for patients with LUAD were downloaded from The Cancer Genome Atlas (TCGA) database (TCGA-LUAD) (Schabath et al., 2016) and the Gene Expression Omnibus (GEO) database (GSE68465, GSE72094) (Zhang A. et al., 2021). The TCGA-LUAD data were used for the training cohort, while those for GSE68465 and GSE72094 were used for the validation cohorts. The TCGA-LUAD dataset was delivered via an Illumina HiSeq 2000 microarray, the GSE68465 dataset was delivered via the Affymetrix Human Genome U133A Array, and the GSE72094 dataset was delivered via the Rosetta/Merck Human RSTA Custom Affymetrix 2.0 microarray. The “sva” package of R software was used to correct the batch effect between different datasets using the “combat” algorithm.

Hypoxia- and EMT-related genes were extracted from the hallmark gene set in the Molecular Signatures Database v7.0(MSigDB, www.gsea-msigdb.org), which includes 200 hypoxia genes and 200 EMT-related genes; 2,498 immune-related genes were acquired in the ImmPort (http://www.immport.org/). This study was approved by the Ethics and Research Committees of Sun Yat-Sen Memorial Hospital and Sun Yat-Sen University.

Information regarding the expression of 200 hypoxia-, 2,498 immune-, and 200 EMT-related genes was collected from the TCGA-LUAD database. Differentially expressed genes (DEGs) between LUAD and normal lung tissue were then identified using the Wilcoxon test according to |Log2FC| > 1 and p < 0.05 would considered as DEGs. Log2FC > 1 indicating upregulated genes and Log2FC < −1 indicating downregulated genes, respectively. Heat and volcano maps were then generated to show the expression of different genes.

An R software package (clusterprofiler, version 3.12) was used to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. Using Fisher’s exact test, those with false discovery rate (FDR)-corrected p values less than 0.05 were regarded as significant indicators.

First, RNA expression in the TCGA-LUAD, GSE68465, and GSE72094 datasets was cross-checked to identify co-expressed and differentially expressed hypoxia-, immunity-, and EMT-related genes. Consequently, univariate Cox analysis of overall survival (OS) was performed to screen for hypoxia-, immunity-, and EMT-related genes with prognostic values. Next, least absolute shrinkage and selection operator (LASSO) regression with 10-fold cross-validation was performed, and 1,000 cycles were run via the R software package “glmnet.” For each cycle, 1,000 random simulations were performed. Based on the optimal lambda value, the best possible gene was selected to construct the model, and a risk formula was established.

The risk scores were calculated according to the expression of each gene and its corresponding regression coefficients using the following equation: risk score = ∑genes Cox coefficient × gene expression. The patients were then categorized into high-risk and low-risk groups based on the optimal cutoff value, which was computed using the “surv_cutpoint” function in the “survminer” R package. Receiver operating characteristic curves were drawn via the R Package “survivalROC” to estimate the predictive sensitivity of the formula. Model effectiveness was evaluated in the validation set using the same coefficients and cutoff values used in the training set. We then evaluated whether the risk score formula exhibited independent prognostic value when combined with clinical variables via multiple regression analysis.

The expression matrices of m6A genes were including (METTL14, METTL3, RBM15, RBM15B, WTAP, CBLL1, ZC3H13, ALKBH5, FTO, YTHDC1, YTHDF1, YTHDC2, YTHDF2, IGF2BP1, YTHDF3, FMR1, HNRNPC, HNRNPA2B1, ELAVL1, and LRPPRC). The expression of m5C genes including (NSUN7, ALYREF, NSUN1, NSUN6, NSUN3, NSUN4, NSUN2 and NSUN5); The expression of m1A genes including (ALKBH3, ALKBH1 and YTHDF2); The expression of m7G genes including (METTL1, BUD23 and RNMT).

Immune cell infiltration was identified using timer 2.0 (cistrome.shinyapps.io/timer/) via the Timer, QUANTISEQ, CIBERSORT, CIBERSORT-ABS, XCELL, MCPCOUNTER, and EPIC algorithms. The “gsva” R package was used to process the single-sample gene set of the enrichment analysis (ssGSEA) to calculate the activity status of 13 immune-related pathways. The selection of immune-checkpoint genes was based on the findings of a previous study (Isomura et al., 2021). The ESTIMATE algorithm was used to calculate the stromal score, immune score, and ESTIMATE score of TCGA-LUAD samples.

Given the lack of information on immune therapy in the TCGA-LUAD cohort, the predictive capability of the risk score formula was evaluated using the GSE135222 (NSCLC), GSE126044 (NSCLC), and IMvigor210 (metastatic urothelial cancer) cohorts (Charoentong et al., 2017; Mariathasan et al., 2018; Jung et al., 2019; Yu et al., 2019; Yu et al., 2020a; Cho et al., 2020).

To evaluate the ability of the risk score to predict the chemotherapeutic response, the half-maximal inhibitory concentration (IC50) of common chemotherapeutic drugs was first calculated in the TCGA-LUAD training group, using the “pRRophetic” package in R software. The Wilcoxon rank test was then used to compare the difference in IC50 between the low- and high-risk groups. Finally, the R package “ggplot” was used to visualize the data.

DEGs were screened using the Wilcoxon test. Univariate Cox analysis of overall survival (OS) was performed to screen relevant genes with prognostic values. Kaplan–Meier survival curves were generated and compared between the two groups using the log-rank test. The associations between the risk score determined using the prognostic model and the stromal score, stemness score, and immune score were assessed using Spearman correlation analysis. All statistical analyses were performed using R version 4.0.0 (R-project.org) and its adequate packages. Statistical significance was set at p < 0.05.

In the training set, 66 of 169 hypoxia-related genes, 556 of 1,214 immune-related genes, and 81 of 177 EMT-related genes were differentially expressed between LUAD and adjacent normal tissues. Of these, 37 hypoxia-related genes, 345 immune-related genes, and 50 EMT-related genes were upregulated, while 29 hypoxia-related genes, 211 immune-related genes, and 31 EMT-related genes were downregulated (Figures 2A,B,E,F,I,J). In total, there were 703 of 1,560 DEGS, 432 and 271 of which were upregulated and downregulated, respectively (Figures 3A,B).

In the GO enrichment analysis, we identified the top 5 GO categories with significant enrichment of genes related to hypoxia, immunity, or EMT. The most significantly altered hypoxia-, immune-, and EMT-related genes were involved in the metabolic processing of ADP, in signaling receptor activator activity, and in extracellular matrix organization, respectively (Figures 2C,G,K; Supplementary Tables S2, S4, S6). We then performed KEGG analysis and identified the top 15 KEGG categories with significant enrichment of hypoxia-, immune-, and EMT-related genes. The altered hypoxia-related genes were mostly associated with glycolysis, while the altered immune-related genes were mostly associated with cytokine–cytokine receptor reactions. The EMT-related genes exhibiting the most significant alterations were involved in focal adhesion (Figures 2D,H,L) (detailed in Supplementary Tables S3, S5, S7). Further, when these gene signatures were combined, the most correlated GO and KEGG categories were signaling receptor activator activity cytokine–cytokine receptor reactions, respectively (Figures 3C,D) (Supplementary Tables S8, S9).

A total of 11,074 genes were co-expressed; among them, 430 of 668 hypoxia-, immune-, and EMT-related DEGs were selected (Figures 4A,B). Then, these 430 genes were used in the univariate Cox regression analysis. A total of 57 prognostic genes were identified (Figure 5A). To avoid overfitting the prognostic model, LASSO regression analysis was performed. Finally, 27 genes were selected and included in the risk score formula, as follows: Risk score = ADAM12 × 0.0537 + CCL20 × 0.1149 + LGR4 × 0.0481 − CTSG × 0.0435 + PDGFB × 0.2173 + INSL4×0.0526 + LIFR×0.0033 + LDHA × 0.1794 − FBP1 × 0.0417 − MAP3K8×0.3235 + SEMA3A × 0.0329 + MC1R × 0.1367 − CD79A × 0.1300 − WFDC2 × 0.0577 + PDYN × 0.2017 − GDF15 × 0.0710 + BCAN × 0.1043 + DDIT4 × 0.0715 − SPOCK1 × 0.0148 + TNFRSF11A × 0.1982 − CX3CR1×0.1238 − AKAP12 × 0.0061 + ANGPTL4 × 0.0227 + GPI × 0.1924 − CAT × 0.0789 + FURIN × 0.0187 + F2RL1 × 0.1408 (Figures 5B,C). Based on am optimistic cut off, 144 and 356 patients were categorized into the high-risk and low-risk groups, respectively (Figures 6A–C). Kaplan–Meier survival analysis revealed that OS was lower in the high-risk group than in the low-risk group (Figure 6E). The area under the curve (AUC) values for predicting 1-, 3-, and 5-year OS were 0.763, 0.766, and 0.728, respectively (HR = 4.26; 95% CI 3.15–5.75; p < 0.0001; Figure 6D). These results show that the risk model based on the 27 genes listed above had high accuracy in predicting the OS of patients with LUAD. Besides, we also proved the novel risk score independently predict the OS of LUAD (Supplementary Figure S1).

To check the stability of the model developed from the training set, patients in the validation sets (GSE68465 and GSE72094) were also divided into a high-risk group and a low-risk group according to the same cut-off value and risk formula as those in TCGA cohort (Figures 7A–C,F–H). The results indicated that OS was markedly lower in the high-risk group than in the low-risk group (Figures 7E,J). In the GSE68465 set, the AUCs for predicting 1-, 3-, and 5-year OS were 0.69, 0.651, and 0.618, respectively (HR = 2.03; 95% CI = 1.55–2.65; p < 0.0001). In the GSE72094 cohort, the corresponding AUCs were 0.653, 0.662, and 0.749, respectively (HR = 2.36; 95% CI = 1.63–3.43; p < 0.001; Figures 7D,I).

Next, we analyzed the association between clinical features (including stage, age, and sex) and the risk score in the TCGA-LUAD database. The risk score remained significantly effective across all subgroups based on tumor stage, sex, and age (Figure 8), supporting the reliability of the risk score formula. Moreover, in the univariate and multivariate Cox regression analysis, the risk score formula was identified as an independent prognostic indicator of poor outcomes in patients with LUAD (Supplementary Figures S1A,B).

Further analysis of the differences in enrichment pathways between the low-risk and high-risk groups showed that the most different pathways were related to the humoral immune response, collagen-containing extracellular matrix, and focal adhesion (Figures 9A,B; Supplementary Tables S10, S11). This may explain why OS was lower in the high-risk group than in the low-risk group.

Growing evidence indicates that increased expression of stemness-related biomarkers in tumor cells is highly correlated with drug resistance, cancer recurrence, and tumor proliferation (Luo and Vögeli, 2020). Hence, we assessed the correlations of the DNA stemness score (DNAss) and RNA stemness score (RNAss) with the risk score. The results indicate that the risk score was significantly positively correlated with the DNAss and RNAss (Figures 10A,B). Besides, this study also compared the gene mutation landscape between high and low risk score group. We found in high risk score group, the mutation frequency of TP53, TTN and KEAP were obviously higher than low risk score group (Supplementary Figure S2).

Previous research has indicated that m⁶A, m⁵C, m¹A and m⁷G modification, which were reversible epigenetic RNA process, significantly involved in the proliferation and migration of cancer cells (Dib et al., 2017; Barbieri and Kouzarides, 2020). In this study, the expression of m⁶A modification genes WTAP, HNRNPA2B1, IGF2BP2, HNRNPC, CBLL1, ELAVL1, RBM15B, LRPPRC, and ELAVL1, the expression of m⁵C modification gene ALYREF, NSUN1 and NSUN2 and the expression of m7G modification gene METTL1, BUD23 and RNMT were significantly higher in the high risk group, while the expression of m⁶A modification gene METTL3, the expression of m⁵C modification gene NSUN7 and NSUN6 were significantly higher in the low-risk group (Figures 10C–F).

The relationship between the risk score and the immune status of the patients in the TCGA cohort is shown in Figures 11A,B. There were significant alterations in immune checkpoint genes. Thus, we further compared the expression of immune checkpoint-related genes between the high-risk group and the low-risk group. The tumor immune microenvironment was also assessed using the immune score, ESTIMATE score, and stromal score (Yoshida et al., 2021). All three scores were negatively correlated with the risk score (Figures 11C–H), indicating stronger tumor immune activity in low-risk patients than in high-risk patients.

To verify the prognostic value of the risk score formula for immunotherapy sensitivity, we selected three immunotherapy datasets from patients with NSCLC and metastatic urothelial cancer. The risk score formula was associated with progression-free survival (PFS) in patients with NSCLC undergoing anti-PD-1/PD-L1 therapy in the GSE135222 cohort (HR = 4.26; 95% CI = 3.15–5.75; p = 0.04) (Figure 12A), and the AUC value for predicting the 12-month PFS was 0.702 (Figure 12B). In the GSE126044 cohort, the risk score was higher in patients with NSCLC who had experienced no benefit (disease progression [PD]) from nivolumab or pembrolizumab than in those who had experienced a benefit (partial response [PR] + stable disease [SD]) (p = 0.017) (Figure 12C). Furthermore, the risk score was associated with worse immunotherapy response in patients with metastatic urothelial cancer (Figures 12D,F). The risk score was also significantly correlated with several immune checkpoint-related genes: PD-1, CD8A, CTLA4, CXCL9, GZMA, HAVCR2, IDO1, PRF1, LAG3, IFNG, GZMB, and TBX2 (Figure 13G). Our analysis further revealed that a high risk score was associated with high sensitivity to common NCCN (National Comprehensive Cancer Network, https://www.nccn.org) recommended anti-LUAD drugs such as cisplatin, docetaxel, paclitaxel and gemcitabine (Figure 13). These results show that the risk score can be used as a potential predictor of chemosensitivity and that immunotherapy may be more appropriate for low-risk patients, while chemotherapy may be more appropriate for high-risk patients.