Wenhao Ouyang1,2†
Yupeng Jiang1†
Shiyi Bu1,2
Tiantian Tang1,2Linjie Huang1,2Ming Chen1,2Yujie Tan1Qiyun Ou1,3Luhui Mao1Yingjie Mai1
Herui Yao1*
Yunfang Yu1,4*Xiaoling Lin1,2*- 1Guangdong Provincial Key Laboratory of Malignant Tumor Epigenetics and Gene Regulation, Department of Medical Oncology, Breast Tumor Centre, Phase I Clinical Trial Centre, Sun Yat-Sen Memorial Hospital, Sun Yat-Sen University, Guangzhou, China
- 2Department of Pulmonary and Critical Care Medicine, Sun Yat-Sen Memorial Hospital, Sun Yat-Sen University, Guangzhou, China
- 3Department of Ultrasound in Medicine, Sun Yat-Sen Memorial Hospital, Sun Yat-Sen University, Guangzhou, China
- 4Artificial Intelligence and Digital Media Programme, Division of Science and Technology, Beijing Normal University-Hong Kong Baptist University United International College, Hong Kong Baptist University, Zhuhai, China
Background: Lung adenocarcinoma (LUAD), the most common subtype of non-small cell lung cancer (NSCLC), is associated with poor prognosis. However, current stage-based clinical methods are insufficient for survival prediction and decision-making. This study aimed to establish a novel model for evaluating the risk of LUAD based on hypoxia, immunity, and epithelial-mesenchymal transition (EMT) gene signatures.
Methods: In this study, we used data from TCGA-LUAD for the training cohort and GSE68465 and GSE72094 for the validation cohorts. Immunotherapy datasets GSE135222, GSE126044, and IMvigor210 were obtained from a previous study. Using bioinformatic and machine algorithms, we established a risk model based on hypoxia, immune, and EMT gene signatures, which was then used to divide patients into the high and low risk groups. We analyzed differences in enriched pathways between the two groups, following which we investigated whether the risk score was correlated with stemness scores, genes related to m6A, m5C, m1A and m7G modification, the immune microenvironment, immunotherapy response, and multiple anti-cancer drug sensitivity.
Results: Overall survival differed significantly between the high-risk and low-risk groups (HR = 4.26). The AUCs for predicting 1-, 3-, and 5-year survival were 0.763, 0.766, and 0.728, respectively. In the GSE68465 dataset, the HR was 2.03, while the AUCs for predicting 1-, 3-, and 5-year survival were 0.69, 0.651, and 0.618, respectively. The corresponding values in the GSE72094 dataset were an HR of 2.36 and AUCs of 0.653, 0.662, and 0.749, respectively. The risk score model could independently predict OS in patients with LUAD, and highly correlated with stemness scores and numerous m6A, m5C, m1A and m7G modification-related genes. Furthermore, the risk model was significantly correlated with multiple immune microenvironment characteristics. In the GSE135222 dataset, the HR was 4.26 and the AUC was 0.702. Evaluation of the GSE126044 and IMvigor210 cohorts indicated that PD-1/PD-LI inhibitor treatment may be indicated in patients with low risk scores, while anti-cancer therapy with various drugs may be indicated in patients with high risk scores.
Conclusion: Our novel risk model developed based on hypoxia, immune, and EMT gene signatures can aid in predicting clinical prognosis and guiding treatment in patients with LUAD.
1 Introduction
Lung adenocarcinoma (LUAD) is the most common subtype of non-small cell lung cancer (NSCLC) (Hyuna et al., 2021). Despite advances in standard treatment strategies based on clinical stage, the survival rate remains poor among patients with LUAD (Bi et al., 2020; Siegel et al., 2021), and the associated tumors are highly heterogeneous. Thus, developing a method for accurately stratifying risk and guiding treatment is essential.
Hypoxic conditions in the tumor microenvironment (TME) and immune microenvironment play a crucial role and are regarded as the major drivers of malignancy in LUAD. Further, both environments are strongly associated with malignant progression, therapeutic resistance, and poor prognosis (Wang D. D. et al., 2021; Wu et al., 2021; Zhang Y. et al., 2021). Several studies have recently shown that a hypoxic stimulus can alter the TME, decreasing the proportion of immune cells and increasing the expression of immunosuppressive cytokines (Zeng et al., 2021). Thus, hypoxia is considered the major immunosuppressive mechanism during cancer development (Labiano et al., 2015). Moreover, hypoxic stimulation can activate epithelial-mesenchymal transition (EMT), a key link in cancer progression (Jiang et al., 2011). Despite these findings, reliable prognostic signatures based on the fundamental combination of hypoxia, immunity, and EMT gene signatures have yet to be established.
Hence, to aid in improving clinical management strategies, the present study aimed to establish a novel model for evaluating LUAD risk based on genes related to hypoxia, immunity, and EMT.
2 Materials and Methods
2.1 Data Acquisition
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.
2.2 Screening of Differentially Expressed Hypoxia-, Immunity-, and EMT-Related Genes
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.
2.3 Functional Exploration of DEGs
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.
2.4 Construction and Verification of the Risk Model
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.
2.5 Selection of m6A, m5C, m1A and m7G Genes
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).
2.6 Differential Analysis of Immune Cell Infiltration, Immune Function, and Immune Checkpoint Function and the Validation of Immunotherapeutic Responses
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).
2.7 Predicting Anti-Cancer Drug Response
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.
2.8 Statistical Analysis
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.
3 Results
Totals of 500 and 840 patients with LUAD were selected from the training and validation sets, respectively. The detailed clinical features of these patients are summarized in Supplementary Table S1. The flowchart of the study is shown in Figure 1.
FIGURE 1. Flow chart of this study.
3.1 Differentially Expressed Hypoxia-, Immune-, and EMT-Related Genes
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).
FIGURE 2. Separativelly screening of differentially expressed hypoxia, immunity and EMT related genes. (A) Volcano plots showing the hypoxia-related DEGs. (B) Heatmaps of differentially expressed hypoxia-related mRNAs. (C) GO enrichment of hypoxia-related DEGs. (D) KEGG pathways of hypoxia-related DEGs. (E) Volcano plots showing the immune-related DEGs. (F) Heatmaps of differentially expressed immune-related mRNAs. (G) GO enrichment of immune-related DEGs. (H) KEGG pathways of immune-related DEGs. (I) Volcano plots showing the EMT-related DEGs. (J) Heatmaps of differentially expressed EMT-related mRNAs. (K) GO enrichment of EMT-related DEGs. (L) KEGG pathways of EMT-related DEGs.
FIGURE 3. Integrally screening of differentially expressed hypoxia, immunity and EMT related genes. (A) Volcano plots of the integrated DEGs. (B) Heatmaps of the integrated mRNAs. (C) GO enrichment based of integrated DEGs. (D) KEGG pathways of integrated DEGs.
3.2 Functional Analysis of Hypoxia-, Immune-, and EMT-Related DEGs
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).
3.3 Predictive Ability of the Risk Score
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).
FIGURE 4. Venn diagram of the intersected genes and DEGs. (A) Venn diagram of the intersected genes among the cohorts of TCGA, GSE68465 and GSE72094. (B) Venn diagram of the intersected hypoxia, immune and EMT related DEGs.
FIGURE 5. Construction of risk score formula. (A) Forest plots showing the results of survival related gene via univariate Cox regression analysis between interacted genes and OS. (B) LASSO coefficient profile plots of the 57 prognostic related genes showing that the variations in the size of the coefficients of parameters shrink with an increasing value of the k penalty. (C) Penalty plot for the LASSO model for the 57 prognostic genes with error bars denoting the standard errors.
FIGURE 6. Prognostic analysis of the risk score formula in the training set. (A) Distribution of risk score for the training set. (B) Patterns of the survival time and survival status between the high-risk and low-risk groups for training set. (C) Heatmaps of the 27 prognostic genes for each patient in training set. (D) Time-related ROC analysis proved the prognostic performance of the risk score in the training set. (E) Kaplan-Meier survival curve of the patients in the high-risk and low-risk groups for OS in the training set.
3.4 Stability of the Risk Score Formula Constructed Using Hypoxia-Related Genes
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).
FIGURE 7. Prognostic analysis of the 27-gene signature model in the validation sets GSE68465 and GSE72094. (A) Distribution of risk score for the GSE68465. (B) Patterns of the survival time and survival status between the high-risk and low-risk groups for GSE68465. (C) Heatmaps of the 27 prognostic genes for each patient in GSE68465. (D) Time-related ROC analysis proved the prognostic performance of the risk score in the GSE68465. (E) Kaplan-Meier survival curve of the patients in the high risk score and low risk score groups for OS in the GSE68465. (F) Distribution of risk score for the GSE72094. (G) Patterns of the survival time and survival status between the high-risk and low-risk groups for GSE72094. (H) Heatmaps of the 27 prognostic genes for each patient in GSE72094. (I) Time-related ROC analysis proved the prognostic performance of the risk score in the GSE72094. (J) Kaplan-Meier survival curve of the patients in the high-risk and low-risk groups for OS in the GSE72094.
3.5 Subgroup Analysis Using the Risk Score Formula
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).
FIGURE 8. Clinical subgroups analysis between high risk score group and low risk score group. (A) Time-related ROC analysis and Kaplan-Meier survival curve of the patients in stage I and stage II. (B) Time-related ROC analysis and Kaplan-Meier survival curve of the patients in stage III and stage IV. (C) Time-related ROC analysis and Kaplan-Meier survival curve of the male patients. (D) Time-related ROC analysis and Kaplan-Meier survival curve of the female patients. (E) Time-related ROC analysis and Kaplan-Meier survival curve of the patients more than 70 years old. (F) Time-related ROC analysis and Kaplan-Meier survival curve of the patients less than or equal to 70 years old.
3.6 Functional Analysis
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.
FIGURE 9. GO and KEGG enrichment analysis between high risk score group and low risk score group. (A) GO enrichment between the high-risk patients and low-risk patients in TCGA cohort. (B) KEGG pathways between the high-risk patients and low-risk patients in TCGA cohort.
3.7 Tumor Stemness Analysis and Gene Mutation Landscape
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).
FIGURE 10. Analysis of risk score with tumor stem cell score and Differential expression of m6A related gene between high risk score group and low risk score group. (A) The relationship between risk score and RNAss. (B) The relationship between risk score and DNAss. (C) Differential expression analysis of m6A related gene between high risk score group and low risk score group. (D) Differential expression analysis of m5C related gene between high risk score group and low risk score group. (E) Differential expression analysis of m1A related gene between high risk score group and low risk score group. (F) Differential expression analysis of m7G related gene between high risk score group and low risk score group.
3.8 Expression of m6A, m5C, m1A and m7G Modification-Related Genes
Previous research has indicated that m6A, m5C, m1A and m7G 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 m6A modification genes WTAP, HNRNPA2B1, IGF2BP2, HNRNPC, CBLL1, ELAVL1, RBM15B, LRPPRC, and ELAVL1, the expression of m5C 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 m6A modification gene METTL3, the expression of m5C modification gene NSUN7 and NSUN6 were significantly higher in the low-risk group (Figures 10C–F).
3.9 Analysis of Immune Status
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.
FIGURE 11. Exploration of tumor immune microenvironment between high risk score group and low risk score group. (A) Heatmap for immune responses based on CIBERSORT, TIMER, CIBERSORT-ABS, QUANTISEQ, MCPCOUNTER, XCELL, EPIC algorithms among high risk score group and low risk score group. (B) ssGSEA for the association between immune cell subpopulations and related functions. (C) ESTIMATE score. (D) Immune score. (E) Stromal score. (F) The relationship between risk score and the ESTIMATE score. (G) The relationship between risk score and the Immune score. (H)The relationship between risk score and the Stromal score. *p < 0.05; **p < 0.01; ***p < 0.001, ns, not significant.
3.10 Analysis of Anti-Cancer Treatment Sensitivity
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.
FIGURE 12. Validation of the risk score formula for immunotherapy. (A) Kaplan-Meier survival curve of GSE135222 cohort for PFS. (B) Time-related ROC analysis proved the prognostic performance of the risk score in the GSE135222. (C) The difference of risk score in the subgroup of PD-1 treatment response in GSE126044. (D) Kaplan-Meier survival curve of IMvigor210 cohort for OS. (E) The difference of Risk score in the subgroup of PD-L1 treatment response. (F) The expression of immue-related checkpoints among high and low risk groups in IMvigor210 cohort. *p < 0.1; *p < 0.05; **p < 0.01; ***p < 0.001.
FIGURE 13. The difference of multiple anti-cancer drugs sensitivity between high risk score group and low risk score group.
4 Discussion
The current era of precision medicine highlights an urgent need to establish a more precise method for evaluating prognosis and guiding treatment in patients with LUAD. Hypoxia, the immune microenvironment, and EMT play crucial roles in tumorigenesis, progression, and drug resistance in LUAD (Isomura et al., 2021; Wu et al., 2021; Yoshida et al., 2021).
Numerous studies have identified the existence of a hypoxic area as one of the key characteristics of cancer growth (Shen et al., 2015). Indeed, hypoxia promotes cancer metastasis and reduces the survival rate in patients with cancer (Nobre et al., 2018), and the expression of hypoxia genes has been shown to increase metabolism in lung adenocarcinoma cells (Smolle et al., 2020). At the same time, the hypoxic microenvironment of the tumor suppresses the ability of immune cells to detect and kill tumor cells (Vito et al., 2020). Additional studies have reported that the hypoxic microenvironment influences the effect of tumor chemotherapy and immune checkpoint inhibitors, which explain the increasing mortality rate among patients with LUAD (Ando et al., 2019; Daniel et al., 2019; Gao et al., 2019). Hypoxia can also contribute to EMT in patients with LUAD (Ando et al., 2019). Moreover, the immune system plays a vital role in the development and progression of malignant tumors (Lin et al., 2021). Immunotherapy is a novel treatment for LUAD that has achieved multiple satisfactory results (Li et al., 2020). EMT has also been shown to play a critical role in tumor development from initiation to metastasis (Taki et al., 2021). After EMT, LUAD cells can produce more extracellular matrix, which can hasten tumor metastasis, aggravate immune evasion, and induce drug resistance (Pastushenko and Cédric, 2019; Cui et al., 2020; Taki et al., 2021). Based on these findings, the present study combined hypoxia-, immunity-, and EMT-related gene signatures to construct a prognostic model for LUAD risk. To our best knowledge, this is the first study to combine these gene signatures to predict prognosis in patients with LUAD.
Using the LASSO Cox algorithm, 200 hypoxia-related genes, 2,498 immune-related genes, and 200 EMT-related genes were used to identify the most robust biomarkers and establish a novel risk score. In total, 27 related genes were included in the risk formula for LUAD prognosis. Using this formula, we classified patients with LUAD into the high- and low-risk groups. The formula had AUCs of 0.763, 0.766, and 0.728 for predicting 1-, 3-, and 5-year OS, respectively, indicating that it has high accuracy and reliability. Further, OS was significantly lower in the high-risk group than in the low-risk group, and the risk score exhibited high predictive capability in the GSE68465 and GSE72094 validation sets.
Subsequent subgroup analysis by sex, age, and stage indicated that the formula exhibited good predictive capability across all categories. Prognosis was accurately predicted in male and female patients and in patients aged >70 years or <70 years. Importantly, patients in the high-risk group also had significantly lower OS than patients in the low-risk group, regardless of disease stage, indicating the need for a gene-based classification for clinical use. Functional analysis between each group revealed strong associations between a high risk score and genes related to the humoral immune response, collagen-containing extracellular matrix, and focal adhesion. All of these are highly correlated with the anti-tumor response, tumor metastasis, drug resistance, and tumor progression (Murphy et al., 2012; Lovitt et al., 2018; Kosibaty et al., 2021).
Stemness-related biomarkers in tumor cells are highly correlated with drug resistance, cancer recurrence, and tumor proliferation (Liu et al., 2021). In this study, stemness-related biomarkers were positively associated with the risk score, demonstrating the prognostic value of the formula. Modification of m6A, m5C, m1A, and m7G are the common type of modification in RNA and plays critical roles in cancer development (Teng et al., 2021; Xu F. et al., 2021; Xu R. et al., 2021). And RNA methylation highly interconnected with hypoxia, immune response and EMT(Lin et al., 2019; Chao et al., 2020; Wang E. et al., 2021). In previous study, researcher identified hypoxia can induced the sumolytion of m6A enzyme(Hou et al., 2021), hypoxia-inducible factor-1 alpha (HIF-1α) can drive m5C modification to promote tumorigenesis(Wang J. Z. et al., 2021), m1A and m6A modification can significantly affect the infiltration of immune cells(Cai et al., 2021; Gao et al., 2021), m7G modification can drives immune evasion (Devarkar et al., 2016), and m6A, m7G modification can induce EMT in cancer development (Yu X. et al., 2020; Xia et al., 2021). Hence we investigate the correlation of RNA methylation with this risk score. In this study, the expression of WTAP, HNRNPA2B1, IGF2BP2, HNRNPC, CBLL1, ELAVL1, RBM15B, LRPPRC, ELAVL1, ALYREF, NSUN1, NSUN2, METTL1, BUD23, RNMT, METTL3, NSUN7 and NSUN6 differed significantly between the high-risk and low-risk groups. These results further support the value of our risk score formula.
Immune cells in the TME are associated with the development of cancer (Bruni et al., 2020). Our risk formula was highly correlated with markers of the immune microenvironment. Previous studies have reported that the characteristics of the immune microenvironment can predict sensitivity to immune checkpoint inhibitor treatment in patients with LUAD (Yu et al., 2020b; Park et al., 2020). In this study, patients in the low-risk group had higher immune activity. To validate the prognostic value of the risk score for immunotherapy sensitivity, we used two external datasets (GSE135222, GSE126044) containing information from patients with NSCLC treated with anti-PD-1/PD-L1 therapy. Our findings indicated that the risk score formula was associated with PFS and treatment responses in patients with NSCLC undergoing anti-PD-1/PD-L1 therapy. The IMvigor210 study examined the effect of PD-L1 inhibitor treatment in patients with metastatic urothelial cancer (Mariathasan et al., 2018; Wang S. et al., 2021). Because relatively few patients were enrolled in the GSE135222 and GSE126044 cohorts, we further evaluated the formula in the IMvigor210 cohort. Our results suggest that the risk score formula is not only suitable for predicting the effect of anti-PD-1/PD-L1 treatment in lung cancer, but that it may also be applicable in patients with other cancer types. Hence, a low-risk score may be an indicator for immunotherapy.
Previously, Sun et al. (2020) reported that the hypoxia-related signature could aid in predicting OS in patients with early-stage LUAD. However, whether hypoxia-related gene signatures can be used to develop a simple predictive formula for late-stage LUAD outcomes or immunotherapy sensitivity remains unknown. Other studies have also attempted to use immune- or EMT-related gene signatures to establish a prognostic model for LUAD (Tang et al., 2020; Wang et al., 2020). However, few of these studies used an authentic immune therapy cohort for validation. Given that multiple factors can have a substantial effect on the prognosis of LUAD and the intimate interconnections among hypoxia, immune responses, and EMT function, we aimed to establish a novel prognostic model based on the integration of multiple gene signatures. Our analysis indicated that the prognostic formula developed in this study exhibits precision for both early- and late-stage LUAD and is valid for predicting sensitivity to immunotherapy based on findings from a clinical cohort.
Currently, there are only a few methods for evaluating tumor sensitivity to chemotherapy (Rochigneux et al., 2020). In this study, the risk score was positively associated with drug sensitivity to cisplatin, docetaxel, paclitaxel and gemcitabine. This indicates that the risk score can be used to determine the appropriateness and benefit of chemotherapy in patients with LUAD, which may aid in the development of individualized treatment strategies.
This study also had some limitations. As the study was based on information within public databases, real-world prospective cohort studies are required to validate the risk score formula. The sequencing methods of the cohort included in this study were different, this may also affect the accuracy of this formula. Furthermore, most patients were Caucasian, highlighting the need to evaluate the predictive ability of the risk score in patients of other races.
5 Conclusion
In summary, this study established a novel 27-gene prognostic risk score for LUAD. The risk score was independently associated with OS in patients with LUAD and with functional analysis, tumor stemness, RNA methylation analysis, the immune microenvironment, and treatment response. Further, it accurately predicted prognosis in subgroups according to age, sex, and disease stage. These findings indicate that molecular risk stratification may be useful for predicting prognosis and guiding treatment in patients with LUAD.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Ethics Statement
The studies involving human participants were reviewed and approved the by Ethics and Research Committees of Sun Yat-Sen Memorial Hospital. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.
Author Contributions
All authors contributed to study conception and design. WO, YJ, and YY analyzed and interpreted the data and performed the statistical analysis; WO, YJ, YY, and HY drafted and revised the manuscript. All authors provided administrative, technical, and material support. YY, HY, and XL supervised the study. All authors reviewed the manuscript and approved the final version of the manuscript.
Funding
This study was supported by grants 2017A030313822 and 2017A030313828 from the Natural Science Foundation of Guangdong Province, grant 2020ZX09201021 from the National Science and Technology Major Project, grant YXRGZN201902 from the Medical Artificial Intelligence Project of Sun Yat-Sen Memorial Hospital; Grant No. 2021M703771 from China Postdoctoral Science Foundation Funded Project; Grants 81572596, 81972471, and 82073408 from the National Natural Science Foundation of China, grant 201704020131 from the Guangzhou Science and Technology Major Program, grant 2017B030314026 from the Guangdong Science and Technology Department, grant 2018007 from the Sun Yat-Sen University Clinical Research 5010 Program, grant SYS-C-201801 from the Sun Yat-Sen Clinical Research Cultivating Program, grant A2020558 from the Guangdong Medical Science and Technology Program, grant SYSU-81000-20200311-0001 and SYSU-05160-20200506-0001 from Tencent Charity Foundation.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.758777/full#supplementary-material
References
Ando, A., Hashimoto, N., Sakamoto, K., Omote, N., Miyazaki, S., Nakahara, Y., et al. (2019). Repressive Role of Stabilized Hypoxia Inducible Factor 1α Expression on Transforming Growth Factor β‐induced Extracellular Matrix Production in Lung Cancer Cells. Cancer Sci. 110 (6), 1959–1973. doi:10.1111/cas.14027
Barbieri, I., and Kouzarides, T. (2020). Role of RNA Modifications in Cancer. Nat. Rev. Cancer 20 (6), 303–322. doi:10.1038/s41568-020-0253-2
Bi, K.-W., Wei, X.-G., Qin, X.-X., and Li, B. (2020). BTK Has Potential to Be a Prognostic Factor for Lung Adenocarcinoma and an Indicator for Tumor Microenvironment Remodeling: A Study Based on TCGA Data Mining. Front. Oncol. 10, 424. doi:10.3389/fonc.2020.00424
Bruni, D., Angell, H. K., and Galon, J. (2020). The Immune Contexture and Immunoscore in Cancer Prognosis and Therapeutic Efficacy. Nat. Rev. Cancer 20 (11), 662–680. doi:10.1038/s41568-020-0285-7
Cai, C., Long, J., Huang, Q., Han, Y., Peng, Y., Guo, C., et al. (2021). M6A "Writer" Gene METTL14: A Favorable Prognostic Biomarker and Correlated with Immune Infiltrates in Rectal Cancer. Front. Oncol. 11, 615296. doi:10.3389/fonc.2021.615296
Chao, Y., Shang, J., and Ji, W. (2020). ALKBH5-m6A-FOXM1 Signaling axis Promotes Proliferation and Invasion of Lung Adenocarcinoma Cells under Intermittent Hypoxia. Biochem. Biophysical Res. Commun. 521 (2), 499–506. doi:10.1016/j.bbrc.2019.10.145
Charoentong, P., Finotello, F., Angelova, M., Mayer, C., Efremova, M., Rieder, D., et al. (2017). Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep 18 (1), 248–262. doi:10.1016/j.celrep.2016.12.019
Cho, J.-W., Hong, M. H., Ha, S.-J., Kim, Y.-J., Cho, B. C., Lee, I., et al. (2020). Genome-wide Identification of Differentially Methylated Promoters and Enhancers Associated with Response to Anti-PD-1 Therapy in Non-small Cell Lung Cancer. Exp. Mol. Med. 52 (9), 1550–1563. doi:10.1038/s12276-020-00493-8
Cui, J., Song, Y., Han, X., Hu, J., Chen, Y., Chen, X., et al. (2020). Targeting 14-3-3ζ Overcomes Resistance to Epidermal Growth Factor Receptor-Tyrosine Kinase Inhibitors in Lung Adenocarcinoma via BMP2/Smad/ID1 Signaling. Front. Oncol. 10, 542007. doi:10.3389/fonc.2020.542007
Daniel, S. K., Sullivan, K. M., Labadie, K. P., and Pillarisetty, V. G. (2019). Hypoxia as a Barrier to Immunotherapy in Pancreatic Adenocarcinoma. Clin. Transl Med. 8 (1), 10. doi:10.1186/s40169-019-0226-9
Devarkar, S. C., Wang, C., Miller, M. T., Ramanathan, A., Jiang, F., Khan, A. G., et al. (2016). Structural Basis for m7G Recognition and 2′-O-Methyl Discrimination in Capped RNAs by the Innate Immune Receptor RIG-I. Proc. Natl. Acad. Sci. USA 113 (3), 596–601. doi:10.1073/pnas.1515152113
Dib, L., San-Jose, L., Ducrest, A.-L., Salamin, N., and Roulin, A. (2017). Selection on the Major Color Gene Melanocortin-1-Receptor Shaped the Evolution of the Melanocortin System Genes. Ijms 18, 2618. doi:10.3390/ijms18122618
Gao, X. Z., Wang, G. N., Zhao, W. G., Han, J., Diao, C. Y., Wang, X. H., et al. (2019). Blocking OLFM4/HIF‐1α axis Alleviates Hypoxia‐induced Invasion, Epithelial-Mesenchymal Transition, and Chemotherapy Resistance in Non‐small‐cell Lung Cancer. J. Cel Physiol 234, 15035–15043. doi:10.1002/jcp.28144
Gao, Y., Wang, H., Li, H., Ye, X., Xia, Y., Yuan, S., et al. (2021). Integrated Analyses of m1A Regulator-Mediated Modification Patterns in Tumor Microenvironment-Infiltrating Immune Cells in colon Cancer. Oncoimmunology 10 (1), 1936758. doi:10.1080/2162402x.2021.1936758
Hou, G., Zhao, X., Li, L., Yang, Q., Liu, X., Huang, C., et al. (2021). SUMOylation of YTHDF2 Promotes mRNA Degradation and Cancer Progression by Increasing its Binding Affinity with m6A-Modified mRNAs. Nucleic Acids Res. 49 (5), 2859–2877. doi:10.1093/nar/gkab065
Hyuna, S., Jacques., F., Siegel, R. L., Mathieu., L., Isabelle, S., Jemal, A., et al. (2021). Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J. Clin. 71 (3), 209–249. doi:10.3322/caac.21660
Isomura, H., Taguchi, A., Kajino, T., Asai, N., Nakatochi, M., Kato, S., et al. (2021). Conditional Ror1 Knockout Reveals Crucial Involvement in Lung Adenocarcinoma Development and Identifies Novel HIF-1α Regulator. Cancer Sci. 112 (4), 1614–1623. doi:10.1111/cas.14825
Jiang, J., Tang, Y.-l., and Liang, X.-h. (2011). EMT: a New Vision of Hypoxia Promoting Cancer Progression. Cancer Biol. Ther. 11 (8), 714–723. doi:10.4161/cbt.11.8.15274
Jung, H., Kim, H. S., Kim, J. Y., Sun, J.-M., Ahn, J. S., Ahn, M.-J., et al. (2019). DNA Methylation Loss Promotes Immune Evasion of Tumours with High Mutation and Copy Number Load. Nat. Commun. 10 (1), 4278. doi:10.1038/s41467-019-12159-9
Kosibaty, Z., Murata, Y., Minami, Y., Noguchi, M., and Sakamoto, N. (2021). ECT2 Promotes Lung Adenocarcinoma Progression through Extracellular Matrix Dynamics and Focal Adhesion Signaling. Cancer Sci. 112 (2), 703–714. doi:10.1111/cas.14743
Labiano, S., Palazon, A., and Melero, I. (2015). Immune Response Regulation in the Tumor Microenvironment by Hypoxia. Semin. Oncol. 42 (3), 378–386. doi:10.1053/j.seminoncol.2015.02.009
Li, F., Huang, Q., Luster, T. A., Hu, H., Zhang, H., Ng, W.-L., et al. (2020). In Vivo Epigenetic CRISPR Screen Identifies Asf1a as an Immunotherapeutic Target in Kras-Mutant Lung Adenocarcinoma. Cancer Discov. 10 (2), 270–287. doi:10.1158/2159-8290.cd-19-0780
Lin, J., Wu, C., Ma, D., and Hu, Q. (2021). Identification of P2RY13 as an Immune-Related Prognostic Biomarker in Lung Adenocarcinoma: A Public Database-Based Retrospective Study. PeerJ 9, e11319. doi:10.7717/peerj.11319
Lin, X., Chai, G., Wu, Y., Li, J., Chen, F., Liu, J., et al. (2019). RNA m6A Methylation Regulates the Epithelial Mesenchymal Transition of Cancer Cells and Translation of Snail. Nat. Commun. 10 (1), 2065. doi:10.1038/s41467-019-09865-9
Liu, X., He, M., Li, L., Wang, X., Han, S., Zhao, J., et al. (2021). EMT and Cancer Cell Stemness Associated with Chemotherapeutic Resistance in Esophageal Cancer. Front. Oncol. 11, 672222. doi:10.3389/fonc.2021.672222
Lovitt, C. J., Shelper, T. B., and Avery, V. M. (2018). Doxorubicin Resistance in Breast Cancer Cells Is Mediated by Extracellular Matrix Proteins. BMC Cancer 18 (1), 41. doi:10.1186/s12885-017-3953-6
Luo, Q., and Vögeli, T.-A. (2020). A Methylation-Based Reclassification of Bladder Cancer Based on Immune Cell Genes. Cancers 12, 3054. doi:10.3390/cancers12103054
Mariathasan, S., Turley, S. J., Nickles, D., Castiglioni, A., Yuen, K., Wang, Y., et al. (2018). TGFβ Attenuates Tumour Response to PD-L1 Blockade by Contributing to Exclusion of T Cells. Nature 554 (7693), 544–548. doi:10.1038/nature25501
Murphy, M. A., O'Leary, J. J., and Cahill, D. J. (2012). Assessment of the Humoral Immune Response to Cancer. J. Proteomics 75 (15), 4573–4579. doi:10.1016/j.jprot.2012.01.021
Nobre, A. R., Entenberg, D., Wang, Y., Condeelis, J., and Aguirre-Ghiso, J. A. (2018). The Different Routes to Metastasis via Hypoxia-Regulated Programs. Trends Cel Biol. 28 (11), 941–956. doi:10.1016/j.tcb.2018.06.008
Park, C., Na, K. J., Choi, H., Ock, C.-Y., Ha, S., Kim, M., et al. (2020). Tumor Immune Profiles Noninvasively Estimated by FDG PET with Deep Learning Correlate with Immunotherapy Response in Lung Adenocarcinoma. Theranostics 10 (23), 10838–10848. doi:10.7150/thno.50283
Pastushenko, I., and Blanpain, C. (2019). EMT Transition States during Tumor Progression and Metastasis. Trends Cel Biol. 29 (3), 212–226. doi:10.1016/j.tcb.2018.12.001
Rochigneux, P., Garcia, A. J., Chanez, B., Madroszyk, A., Olive, D., and Garon, E. B. (2020). Medical Treatment of Lung Cancer: Can Immune Cells Predict the Response? A Systematic Review. Front. Immunol. 11, 1036. doi:10.3389/fimmu.2020.01036
Schabath, M. B., Welsh, E. A., Fulp, W. J., Chen, L., Teer, J. K., Thompson, Z. J., et al. (2016). Differential Association of STK11 and TP53 with KRAS Mutation-Associated Gene Expression, Proliferation and Immune Surveillance in Lung Adenocarcinoma. Oncogene 35 (24), 3209–3216. doi:10.1038/onc.2015.375
Shen, H., Feng, G., Cui, J., Du, Q., Qin, Y., Cai, J., et al. (2015). Clinical Implications of Serum Hypoxia Inducible Factor-1α and Vascular Endothelial Growth Factor in Lung Cancer. Tumori 101 (4), 404–411. doi:10.5301/tj.5000320
Siegel, R. L., Miller, K. D., Fuchs, H. E., and Jemal, A. (2021). Cancer Statistics, 2021. CA A. Cancer J. Clin. 71 (1), 7–33. doi:10.3322/caac.21654
Smolle, E., Leko, P., Stacher-Priehse, E., Brcic, L., El-Heliebi, A., Hofmann, L., et al. (2020). Distribution and Prognostic Significance of Gluconeogenesis and Glycolysis in Lung Cancer. Mol. Oncol. 14 (11), 2853–2867. doi:10.1002/1878-0261.12780
Sun, J., Zhao, T., Zhao, D., Qi, X., Bao, X., Shi, R., et al. (2020). Development and Validation of a Hypoxia-Related Gene Signature to Predict Overall Survival in Early-Stage Lung Adenocarcinoma Patients. Ther. Adv. Med. Oncol. 12, 1758835920937904. doi:10.1177/1758835920937904
Taki, M., Abiko, K., Ukita, M., Murakami, R., Yamanoi, K., Yamaguchi, K., et al. (2021). Tumor Immune Microenvironment during Epithelial-Mesenchymal Transition. Clin. Cancer Res. 27, 4669–4679. doi:10.1158/1078-0432.CCR-20-4459
Tang, Y., Jiang, Y., Qing, C., Wang, J., and Zeng, Z. (2020). Systematic Construction and Validation of an Epithelial-Mesenchymal Transition Risk Model to Predict Prognosis of Lung Adenocarcinoma. Aging 13 (1), 794–812. doi:10.18632/aging.202186
Teng, P. C., Liang, Y., Yarmishyn Aliaksandr, A., Hsiao, Y. J., Lin, T. Y., Lin, T. W., et al. (2021). RNA Modifications and Epigenetics in Modulation of Lung Cancer and Pulmonary Diseases. Int. J. Mol. Sci. 22, 10592. doi:10.3390/ijms221910592
Vito, A., El-Sayes, N., and Mossman, K. (2020). Hypoxia-Driven Immune Escape in the Tumor Microenvironment. Cells 9, 992. doi:10.3390/cells9040992
Wang, D. D., Shaver, L. G., Shi, F. Y., Wei, J. J., Qin, T. Z., Wang, S. Z., et al. (2021). Comparative Efficacy and Safety of PD-1/pd-L1 Immunotherapies for Non-small Cell Lung Cancer: a Network Meta-Analysis. Eur. Rev. Med. Pharmacol. Sci. 25 (7), 2866–2884. doi:10.26355/eurrev_202104_25541
Wang, E., Li, Y., Ming, R., Wei, J., Du, P., Zhou, P., et al. (2021). The Prognostic Value and Immune Landscapes of a m6A/m5C/m1A-Related LncRNAs Signature in Head and Neck Squamous Cell Carcinoma. Front. Cel Dev. Biol. 9, 718974. doi:10.3389/fcell.2021.718974
Wang, J. Z., Zhu, W., Han, J., Yang, X., Zhou, R., Lu, H. C., et al. (2021). The Role of the HIF‐1α/ALYREF/PKM2 axis in Glycolysis and Tumorigenesis of Bladder Cancer. Cancer Commun. 41 (7), 560–575. doi:10.1002/cac2.12158
Wang, Q., Li, M., Yang, M., Yang, Y., Song, F., Zhang, W., et al. (2020). Analysis of Immune-Related Signatures of Lung Adenocarcinoma Identified Two Distinct Subtypes: Implications for Immune Checkpoint Blockade Therapy. Aging 12 (4), 3312–3339. doi:10.18632/aging.102814
Wang, S., Ma, H., Li, H., Liu, Q., Huang, S., Huang, L., et al. (2021). Alternatively Expressed Transcripts Analysis of Non-small Cell Lung Cancer Cells under Different Hypoxic Microenvironment. J. Oncol. 2021, 5558304. doi:10.1155/2021/5558304
Wu, J., Li, L., Zhang, H., Zhao, Y., Zhang, H., Wu, S., et al. (2021). A Risk Model Developed Based on Tumor Microenvironment Predicts Overall Survival and Associates with Tumor Immunity of Patients with Lung Adenocarcinoma. Oncogene 40, 4413–4424. doi:10.1038/s41388-021-01853-y
Xia, P., Zhang, H., Xu, K., Jiang, X., Gao, M., Wang, G., et al. (2021). MYC-targeted WDR4 Promotes Proliferation, Metastasis, and Sorafenib Resistance by Inducing CCNB1 Translation in Hepatocellular Carcinoma. Cell Death Dis 12 (7), 691. doi:10.1038/s41419-021-03973-5
Xu, F., Huang, X., Li, Y., Chen, Y., and Lin, L. (2021). m6A-related lncRNAs Are Potential Biomarkers for Predicting Prognoses and Immune Responses in Patients with LUAD. Mol. Ther. - Nucleic Acids 24, 780–791. doi:10.1016/j.omtn.2021.04.003
Xu, R., Pang, G., Zhao, Q., Yang, L., Chen, S., Jiang, L., et al. (2021). The Momentous Role of N6‐methyladenosine in Lung Cancer. J. Cel Physiol 236 (5), 3244–3256. doi:10.1002/jcp.30136
Yoshida, C., Kadota, K., Ikeda, T., Ibuki, E., Go, T., Haba, R., et al. (2021). Tumor-associated Macrophage Infiltration Is Associated with a Higher Rate of Tumor Spread through Air Spaces in Resected Lung Adenocarcinomas. Lung Cancer 158, 91–96. doi:10.1016/j.lungcan.2021.06.009
Yu, X., Zhao, H., and Cao, Z. (2020). The m6A Methyltransferase METTL3 Aggravates the Progression of Nasopharyngeal Carcinoma through Inducing EMT by m6A-Modified Snail mRNA. Minerva Med. doi:10.23736/S0026-4806.20.06653-7
Yu, Y., Lin, D., Li, A., Chen, Y., Ou, Q., Hu, H., et al. (2020a). Association of Immune Checkpoint Inhibitor Therapy with Survival in Patients with Cancers with MUC16 Variants. JAMA Netw. Open 3 (6), e205837. doi:10.1001/jamanetworkopen.2020.5837
Yu, Y., Zeng, D., Ou, Q., Liu, S., Li, A., Chen, Y., et al. (2019). Association of Survival and Immune-Related Biomarkers with Immunotherapy in Patients with Non-small Cell Lung Cancer. JAMA Netw. Open 2 (7), e196879. doi:10.1001/jamanetworkopen.2019.6879
Yu, Y., Zhang, W., Li, A., Chen, Y., Ou, Q., He, Z., et al. (2020b). Association of Long Noncoding RNA Biomarkers with Clinical Immune Subtype and Prediction of Immunotherapy Response in Patients with Cancer. JAMA Netw. Open 3 (4), e202149. doi:10.1001/jamanetworkopen.2020.2149
Zeng, F., Zhang, Y., Han, X., Zeng, M., Gao, Y., and Weng, J. (2021). Employing Hypoxia Characterization to Predict Tumour Immune Microenvironment, Treatment Sensitivity and Prognosis in Hepatocellular Carcinoma. Comput. Struct. Biotechnol. J. 19, 2775–2789. doi:10.1016/j.csbj.2021.03.033
Zhang, A., Yang, J., Ma, C., Li, F., and Luo, H. (2021). Development and Validation of a Robust Ferroptosis-Related Prognostic Signature in Lung Adenocarcinoma. Front. Cel Dev. Biol. 9, 616271. doi:10.3389/fcell.2021.616271
Keywords: lung adenocarcinoma, hypoxia, immune, EMT, gene signature, immunotherapy response
Citation: Ouyang W, Jiang Y, Bu S, Tang T, Huang L, Chen M, Tan Y, Ou Q, Mao L, Mai Y, Yao H, Yu Y and Lin X (2022) A Prognostic Risk Score Based on Hypoxia-, Immunity-, and Epithelialto-Mesenchymal Transition-Related Genes for the Prognosis and Immunotherapy Response of Lung Adenocarcinoma. Front. Cell Dev. Biol. 9:758777. doi: 10.3389/fcell.2021.758777
Received: 15 August 2021; Accepted: 28 December 2021;
Published: 24 January 2022.
Edited by:
Daniele Vergara, University of Salento, ItalyReviewed by:
Jie Meng, Central South University, ChinaWangjun Liao, Southern Medical University, China
Copyright © 2022 Ouyang, Jiang, Bu, Tang, Huang, Chen, Tan, Ou, Mao, Mai, Yao, Yu and Lin. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Herui Yao, yaoherui@mail.sysu.edu.cn; Yunfang Yu, yuyf9@mail.sysu.edu.cn; Xiaoling Lin, LinXL9@mail.sysu.edu.cn
†These authors have contributed equally to this work and share first authorship