Identifying Molecular Subtypes and 6-Gene Prognostic Signature Based on Hypoxia for Optimizing Targeted Therapies in Non-Small Cell Lung Cancer

Introduction

Lung cancer, as the most frequently diagnosed cancer all over the world, shows the highest incidence in males.¹ A major histological subtype of lung cancer is non-small cell lung cancer (NSCLC), of which lung adenocarcinoma and lung squamous cell carcinoma are the main subtypes. Tobacco smoking is the most common risk factor for lung cancer, and it is even the leading cause of more than 80% of all lung cancer cases in countries with a high smoking rate.² But still, non-smokers who may be exposed to second-hand smoking and occupational carcinogens may also develop lung cancer. Tobacco prevention, for example, the use of e-cigarettes to eliminate carcinogenic chemicals in cigarettes, is an effective strategy for preventing lung cancer, but it is far from effective to fight against lung cancer. A number of studies have shown that genetic susceptibility is an important cause for lung cancer. Kirsten rat sarcoma (KRAS) and epidermal growth factor receptor (EGFR) are the most common genes detected in NSCLC patients.³

Oxygen deprivation (hypoxia) is a critical characteristic in solid tumors, and it can facilitate cancer progression and metastasis. Arterial oxygen tension is 5.73 kPa in the lung, and an oxygen tension below 2.03–3.04 kPa is defined as hypoxia. Evidence suggested that a significantly low oxygen below 0.02 kPa can induce oriented mutations and improve resistance to ionizing photon irradiation,^4,5 indicating that hypoxia could lead to unfavorable survival of cancer.⁶ Inspired by such an important role in cancers, researchers started to consider hypoxia as a therapeutic target in many cancers including in NSCLC, and have discovered encouraging results.⁷ Hypoxia-targeted therapies are effective to suppress drug resistance of hypoxia, and have shown a strong performance when combined with other therapies, such as radiotherapy and chemotherapy in clinical trials.⁸

To design a personalized treatment for optimized outcomes, selecting suitable NSCLC patients for receiving hypoxia-targeted therapies is important. To this end, many studies have explored a series of gene signatures to more accurately identify lung cancer patients with different hypoxia levels.^9–11 Also, a number of biomarkers associated with NSCLC prognosis have been discovered to help screen NSCLC patients with high risk.^12,13 However, to the best of our knowledge, none of the studies have developed molecular subtypes based on hypoxia for NSCLC. Compared with hypoxia-related biomarkers, molecular subtypes can directly evaluate whether a NSCLC patient could benefit from receiving hypoxia-target therapies or other personalized therapies such as immunotherapy.

Therefore, in this study, we aimed to classify a molecular subtyping system based on hypoxia-related genes for NSCLC patients. We analyzed functional pathways in different molecular subtypes and predicted the sensitivity of patients with different subtypes to immunotherapy and chemotherapy. In addition, we also developed a prognostic signature based on hypoxia-related genes, and the signature could serve as an independent risk factor for NSCLC patients.

Materials and Methods

Data Source

TCGA-LUAD and TCGA-LUSC datasets containing RNA-seq and clinical data were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). Expression profiles of gene chips, including GSE31210, GSE30219, GSE50081, GSE19188, GSE37745, and GSE29013 datasets, were acquired from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). These gene chips were all sequenced by GPL570 platform. Two hundred hypoxia-related genes were obtained from h.all.v7.3.symbols.gmt file in the Molecular Signatures Database v7.3 (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp).

Data Preprocessing

In TCGA datasets, genes with low expression (transcript per million (TPM) < 1 in over half the samples) were excluded. Ensembl ID was converted to gene symbol. Median expression value was selected when one gene symbol corresponded to multiple IDs. Expression was transformed to log2 (expression). Using “RemoveBatchEffect” function in limma R package,¹⁴ the batch effects of TCGA-LUAD and TCGA-LUSC datasets were removed and two datasets were combined to one dataset, named TCGA-NSCLC. Principal component analysis (PCA) described sample distribution before and after batch effect removal (Supplementary Figure S1A and B).

In GEO datasets, rma function in affy (V1.66.0) R package¹⁵ was applied to process and normalize raw data. Probes were converted to gene symbols, according to the GPL570 annotation file. Median expression was selected when multiple probes matched one gene symbol, whereas probes matching multiple gene symbols were eliminated. Samples without survival time or survival status were excluded. Using the “RemoveBatchEffect” function, the batch effects of the six datasets were eliminated and combined to one cohort. PCA showed that samples were clustered after removing batch effects (Supplementary Figure S1C and D).

Identifying Molecular Subtypes Based on Hypoxia-Related Genes

Firstly, univariate Cox regression analysis was conducted on the TCGA-NSCLC dataset and GSE cohorts to screen hypoxia-related genes significantly associated with prognosis (P<0.05). Intersected genes in the two datasets were retained for further analysis. ConsensusClusterPlus R package (v1.48.0)¹⁶ was applied to perform consensus clustering under the parameters of reps =100, pItem =0.8, pFeature =1, clusterAlg = “km”, distance = “euclidean”. Cluster number k =2–10 was analyzed to confirm the optimal cluster number.

Clinical Features and Known Immune Subtypes Analysis in Molecular Subtypes

The distribution of different clinical features, including survival status, age, gender, smoking status, and stages, in C1 and C2 were analyzed using Chi-squared test. Moreover, the distributions of six immune subtypes across 33 cancer types, namely C1 (wound healing), C2 (IFN-γ dominant), C3 (inflammatory), C4 (lymphocyte depleted), C5 (immunologically quiet), and C6 (TGF-β dominant), in C1 and C2 were analyzed using Chi-squared test.

T-Distributed Stochastic Neighbor Embedding (t-SNE)

T-SNE with two tree-based algorithms is commonly used for processing high-dimensional data and data visualization in scatter plots.¹⁷ The expression data of hypoxia-related genes in the TCGA dataset and GSE cohorts were used as an input in t-SNE visualization.

SubMap Analysis

Gene Pattern SubMap analysis¹⁸ was applied to compare the similarity of expression profiles between the two groups. This was performed here for evaluating the robustness of molecular subtypes in different datasets. The subtyping system was considered as robust if similar expression data of the same molecular subtype were shown in different datasets. In addition, we performed SubMap to predict the sensitivity of different molecular subtypes to immunotherapy (details in “Predicting Sensitivity to Immunotherapy and Chemotherapy” section).

Functional Analysis

Limma R package was employed to identify differential expressed genes between different molecular subtypes. False discovery rate (FDR) <0.05 and |fold change (FC)| >1.5 were set to screen significantly differentially expressed genes. R package of clusterProfiler (v3.16.0)¹⁹ was applied to analyze gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. Significantly enriched terms and pathways were screened under FDR <0.05. Gene set enrichment analysis (GSEA) in clusterProfiler was employed to enrich pathways (h.all.v7.3.symbols.gmt) obtained from MsigDB.

Assessment of Tumor Microenvironment

Microenvironment Cell Populations-counter (MCP-counter)²⁰ could characterize the immune cells in TME and calculate the abundance of 10 immune cells across mixtures. CIBERSORT,²¹ a computational method, is widely used to assess the enrichment of 22 immune cells through microarray or RNA-seq data. Estimation of STromal and Immune cells in MAlignant Tumours using Expression data (ESTIMATE)²² evaluates TME through calculating immune score and stromal score based on gene signatures. ESTIMATE score is the combined score of immune score and stromal score. Based on gene expression data, EPIC measurement²³ estimates the fraction of immune cells, cancer cells, and other non-malignant cells.

Predicting Sensitivity to Immunotherapy and Chemotherapy

Tumor Immune Dysfunction and Exclusion (TIDE) analysis²⁴ was performed to predict the immune response to immunotherapy through analyzing T-cell dysfunction and T-cell exclusion in LUAD. The proportion of positive and negative responders with different molecular subtypes to immune checkpoint blockade (ICB) was calculated. To more correctly evaluate the sensitivity of molecular subtypes to immunotherapy, we introduced a dataset treated by anti-programmed cell death protein 1 (PD-1) (GSE78220) and conducted SubMap analysis to compare expression profiles between TCGA-NSCLC and GSE78220 datasets. To analyze patients’ sensitivity to chemotherapy, pRRophetic R package²⁵ was applied to determine the estimated biochemical half maximal inhibitory concentration (IC50) of seven chemotherapeutic drugs, namely, bexarotene, doxorubicin, embelin, etoposide, gemcitabine, mitomycin C, and vinorelbine.

Constructing a Prognostic Model

Univariate Cox regression analysis was performed to filter hypoxia-related genes significantly associated with prognosis using the TCGA-NSCLC dataset and GSE cohorts. Intersection of the screened genes in the two datasets served as an input in LASSO analysis to reduce gene number.²⁶ StepAIC in MASS R package²⁷ and multiple Cox regression were applied to further optimize the prognostic model. Finally, the prognostic model was defined as: risk score = , where coef(i) represents the coefficient of gene(i). The risk score of each sample was calculated and converted to z-score. Surv_cutpoint function in survminer R package (http://www.sthda.com/english/rpkgs/survminer/) was conducted to determine the optimal cut-off of z-score to classify samples into high-risk and low-risk groups. Univariate and multivariate Cox regression analysis were implemented to evaluate the independence of risk score to predict LUAD prognosis. The reliability and practicability of risk score and nomogram was evaluated by decision curve analysis (DCA).

Statistical Analysis

Kaplan-Meier (KM) curves were used to show patients’ survival status. All statistical analyses were performed in R (v3.4.2) platform. Statistical methods were described in the corresponding sections and legends. P<0.05 was considered as significant. All parameters were default if not specified.

GitHub Page

https://github.com/jianwem/pro2022/blob/main/Hypoxia_LUAD_immune_model.R

Results

Construction of Molecular Subtypes for NSCLC Based on Hypoxia-Related Genes

The expression profiles of hypoxia-related genes were extracted from TCGA-NSCLC dataset and GSE cohorts. Genes related to hypoxia were obtained from the hallmark hypoxia pathway. Univariate Cox regression analysis was applied to identify the hypoxia-related genes associated with NSCLC prognosis, specifically, there were 45 and 56 genes from the TCGA-NSCLC dataset and GSE cohorts, respectively (P<0.05). Subsequently, NSCLC samples were categorized into two subtypes (C1 and C2) based on 101 prognosis-associated genes in the TCGA-NSCLC dataset by consensus clustering (Supplementary Figure S2). The two subgroups showed a significantly distinct prognosis, with more favorable overall survival (OS) being found in the C1 subtype (P=0.00011, Figure 1A). The two subtypes were separately visualized by t-SNE, the results of which proved that the classification system was valid (Figure 1B). Enrichment analysis of the hypoxia pathway revealed that, compared with the C1 subtype, the C2 subtype was significantly higher enriched to the hypoxia pathway (P<0.0001, Figure 1C and D). In GSE cohorts, we also observed similar results that 801 samples were clearly divided into the two subgroups, and that the hypoxia pathway was significantly enriched in C2 subtype (P<0.0001, Figure 1E–H). Furthermore, we performed SubMap analysis to compare the similarity between the two groups, which have been classified into C1 and C2 subtypes. The results showed significantly similar expression data between TGCA-NSCLC and GSE cohorts in the same subtype (Supplementary Figure S3). The above results demonstrated the effective classification system based on hypoxia-related genes.

Figure 1 Construction and validation of two molecular subtypes. (A) Kaplan-Meier survival curve between C1 and C2 subtypes in the TCGA-NSCLC dataset. Log rank test was conducted. (B) T-SNE analysis for describing the distribution of samples with two subtypes in the TCGA-NSCLC dataset. (C) Comparison of enrichment of hypoxia pathway between C1 and C2 subtypes in the TCGA-NSCLC dataset. Wilcoxon test was performed. (D) The enrichment score of hypoxia pathway ranking for samples in the TCGA-NSCLC dataset. (E) Kaplan-Meier survival curve between C1 and C2 subtypes in GSE cohorts. Log rank test was conducted. (F) T-SNE analysis for describing the distribution of samples with two subtypes in GSE cohorts. (G) Comparison of enrichment of hypoxia pathway between C1 and C2 subtypes in GSE cohorts. Wilcoxon test was performed. (H) The enrichment score of hypoxia pathway ranking for samples in GSE cohorts. Abbreviations: ES, enrichment score; NES, normalized enrichment score; FDR, false discovery rate.

Molecular Subtypes are Associated with Clinical Features

As for the relation between the two molecular subtypes and clinical features, we analyzed the distribution of C1 and C2 in different clinical features, including survival status, age, gender, smoking status, and stages. It has been found that the C1 subtype had a higher proportion of alive patients than C2 (P=0.0391, Figure 2A). There was no difference in age, gender, or smoking status (P>0.05, Figure 2B–D). In the staging status, it was observed that the C1 subtype consisted of a higher proportion of patients with mild stages, such as T0, N0, M0, and stage I (Figure 2E–H). The above data validated that C1 and C2 varied in clinical features and the subtyping was effective.

Figure 2 The distribution of different clinical features in C1 and C2 subtypes, including survival status (A), age (B), gender (C), smoking (D), T stage (E), N stage (F), M stage (G), and stage (H). Chi-square test was conducted.

Oncogenic Pathways are Highly Enriched in C2 Subtype

We next used limma R package to analyze the expression profiles of the TCGA-NSCLC dataset and GSE cohorts. FDR <0.05 and |FC| >1.5 were selected as the criteria to screen DEGs. In the TCGA-NSCLC dataset, the C2 subtype contained 610 DEGs incorporating 496 up-regulated and 114 down-regulated genes. GO function analysis revealed that extracellular matrix-related terms were the most annotated in up-regulated genes (Supplementary Figure S4A–C). Enrichment analysis of KEGG pathways demonstrated that up-regulated genes were highly enriched in tumor-related pathways, such as the PI3K-Akt signaling pathway, ECM-receptor interaction, and proteo glycans in the cancer pathway (Supplementary Figure S4D). In GSE cohorts, a total of 903 DEGs (552 up-regulated and 381 down-regulated) in the C2 subtype were screened. Analysis on the up-regulated genes showed that extracellular matrix-related terms were also enriched in up-regulated genes, and that tumor-related pathways including cell cycle, IL-17 signaling pathway, p53 signaling pathway, and DNA replication were highly enriched (Supplementary Figure S5). Moreover, GSEA was applied to evaluate pathways from hallmark genes (h.all.v7.3.symbols.gmt) in C1 and C2 subtypes. Oncogenic pathways, such as epithelial mesenchymal transition (EMT), hypoxia, NOTCH signaling, p53 pathway, TGF-β signaling, and WNT-β catenin signaling pathways, were significantly enriched in C2 subtype both in the TCGA-NSCLC dataset and GSE cohorts (Figure 3).

Figure 3 Assessment of enriched pathways from h.all.v7.3.symbols.gmt through GESA in TCGA-NSCLC dataset (A) and GSE cohorts (B). Horizontal axis represents the distribution of each sample and vertical axis represents the enrichment score of samples. Different pathways are labeled in different colors.

Comparison with Other Molecular Subtypes

In a pan-cancer study, Thorsson et al²⁸ proposed six immune subtypes across 33 cancer types, namely C1 (wound healing), C2 (IFN-γ dominant), C3 (inflammatory), C4 (lymphocyte depleted), C5 (immunologically quiet), and C6 (TGF-β dominant). Accordingly, we used the classification system to divide the samples in the TCGA-NSCLC dataset. The results showed that the C5 subtype was not present in these samples, and that C1, C2, and C3 formed the majority of them (Figure 4A). The distribution of our subtypes C1 and C2 was compared with these five subtypes, and we also analyzed the OS of the five subtypes. Here, C2 and C3 subtypes were found to have a more favorable prognosis than the rest, and they had a lower proportion in our subtype C2 than C1 (Figure 4B and C). The C3 subtype consisted of a markedly higher proportion (48%) in our subtype C2 than in C1 (29%), which was related to an unfavorable prognosis of our subtype C2. The above results indirectly demonstrated that our subtyping was effective.

Figure 4 Comparison with immune subtypes from a pan-cancer study. (A) Donut plot of comparing our molecular subtypes (outside circle) with immune subtypes (inside circle). (B) Histogram of the distribution of five immune subtypes in our molecular subtypes. Chi-square test was performed. (C) Kaplan-Meier survival plot of five immune subtypes. Log rank test was performed.

Distinct Immune Features Between C1 and C2 Subtypes

Next we compared the immune features of C1 and C2 subtypes in the TCGA-NSCLC dataset using CIBERSORT, MCP-counter, EPIC, and ESTIMATE. CIBERSORT analysis revealed a clear difference of 14 out of 22 immune cells between C1 and C2, especially M0 macrophages and CD8 T-cells (P<0.0001, Figure 5A). C2 had a higher proportion of M0 macrophages than C1 but a lower percent of CD8 T-cells. MCP-counter and EPIC methods showed similar results that the proportion of CD8 T-cells were all lower in C2, and a higher proportion of cancer-associated fibroblasts (CAFs) was shown in C2 (P<0.0001, Figure 5B and C). From the results of ESTIMATE analysis, it could be found that the C2 subtype had a significantly higher enrichment score of immune score and stromal score than C1, indicating that C2 involved more immune infiltrated than C1 (Figure 5D). In addition, the two subtypes exhibited a significant difference in immune checkpoint expression because 33 out of 47 immune checkpoints showed differential expression between C1 and C2 (P<0.05, Figure 5E). Furthermore, GSE cohorts demonstrated similar results to the above immune features, suggesting a distinct tumor immune microenvironment of two subtypes (Supplementary Figure S6).

Figure 5 Characterization of tumor microenvironment between two subtypes in TCGA-NSCLC dataset. (A) Enrichment of 22 immune cells evaluated by CIBERSORT. (B) Enrichment of 10 immune cells analyzed by MCP-Counter. (C) EPIC measurement for describing eight immune cells. (D) ESTIMATE method for calculating stromal score and immune score. (E) Expression of 47 immune checkpoints presented as log2 (TPM). Student t-test was performed. *P<0.05, **P<0.01, ***P<0.001, ****P<0.0001. Abbreviations: ns, no significance; TPM, transcript per million.

C1 Subtype is Predicted to Have a Favorable Prognosis in Immunotherapy and Chemotherapy

As C1 and C2 subtypes presented differential tumor microenvironments, whether they would show a differential response to targeted therapy was examined using TIDE measurement to predict patients’ response to immunotherapy. High TIDE score indicated a high immune escape to immunotherapy. In the TCGA-NSCLC dataset, the C1 subtype showed a lower TIDE score but a higher proportion of positive response to immunotherapy than C2, suggesting an insensitivity of C2 to immunotherapy (P<0.0001, Figure 6A–C). In GSE cohorts, similar results were obtained that C1 was more sensitive than C2 subtype, with a 57% positive response in C1 in comparison with 26% in C2 (Figure 6D–F).

Figure 6 Sensitivity of two subtypes to immunotherapy predicted by TIDE. (A) Kaplan-Meier survival plot of samples grouped by positive (true) and negative (false) responses to immunotherapy in the TCGA-NSCLC dataset. Log rank test was performed. (B) TIDE score of C1 and C2 subtypes in the TCGA-NSCLC dataset. Wilcoxon test was performed. (C) The distribution of positive and negative responses in C1 and C2 subtypes in the TCGA-NSCLC dataset. Chi-square test was performed. (D) Kaplan-Meier survival plot of samples grouped by positive and negative responses to immunotherapy in GSE cohorts. Log rank test was performed. (E) TIDE score of C1 and C2 subtypes in GSE cohorts. Wilcoxon test was performed. (F) The distribution of positive and negative responses in C1 and C2 subtypes in GSE cohorts. Chi-square test was performed.

Furthermore, we used an openly available dataset of melanoma patients (GSE78220) treated by anti-PD1 therapy to predict the response of NSCLC patients to immunotherapy. Submap analysis was used to compare the similarity of expression profiles between melanoma patients and NSCLC patients. In the TCGA-NSCLC dataset, C1 subtype showed similar expression compared with patients treated by anti-PD1 (P=0.011, Figure 7A), indicating that C1 was predicted to be actively responsive to anti-PD1 immunotherapy. In GSE cohorts, we also observed similar results, although there was no significant difference (Figure 7B). To know the response of two subtypes to chemotherapy, pRRophetic R package was applied to predict their sensitivity to seven chemotherapeutic drugs (bexarotene, doxorubicin, embelin, etoposide, gemcitabine, mitomycin C, and vinorelbine). A significant difference of response to these drugs was observed between C1 and C2 subtypes (Figure 7C and D). The C1 subtype had a higher estimated IC50 than the C2 subtype, which indicated that C1 was more sensitive to these chemotherapeutic drugs.

Figure 7 Predicting response of C1 and C2 subtypes to immunotherapy and chemotherapy. (A and B) Submap analysis between the anti-PD-1 treatment group and TCGA-NSCLC dataset, the anti-PD-1 treatment group and GSE cohorts. Anti-PD-1-NR indicates no response to anti-PD-1 and anti-PD-1-R indicates a positive response. P-value was corrected by Bonferroni correction. (C and D) Estimated IC50 of chemotherapeutic drugs including bexarotene, doxorubicin, embelin, etoposide, gemcitabine, mitomycin C, and vinorelbine in the TCGA-NSCLC dataset and GSE cohorts. **P<0.01, ****P<0.0001. Abbreviation: ns, no significance.

Construction of a Prognostic Model Based on Hypoxia-Related Genes

In the previous sections, we demonstrated that C1 and C2 subtypes identified based on hypoxia-related genes varied in OS, tumor microenvironment and sensitivity to immunotherapy and chemotherapy. Then, we screened prognostic genes related to hypoxia, and constructed a prognostic model for NSCLC patients. Here, 14 hypoxia-related genes significantly associated with prognosis were filtered using univariate Cox regression analysis. Then LASSO regression analysis was performed to reduce the number of genes with increased lambda value (Supplementary Figure S7A). Five-fold cross validation was applied to construct the model and analyze the confidence interval in each lambda value (Supplementary Figure S7B). When lambda =0.00978, the optimal model was generated, and here 11 genes remained. Subsequently, we performed stepAIC to further reduce the number of genes, and finally six genes were screened, but four genes showed a significant association with NSCLC prognosis (P<0.05, Figure 8). The prognostic model was defined as: risk score =0.173*ALDOA–0.139*EFNA1–0.071*GPC4+0.059*HOXB9+0.14*PGM2+0.157*PLAUR.

Figure 8 Multivariate Cox regression analysis of six prognostic genes. Log rank test was performed. *P<0.05, **P<0.01.

The risk score was calculated for each sample in TCGA-NSCLC dataset, and we used survminer R package to confirm the optimal cut-off in sample classification into high-risk and low-risk groups. Survival analysis revealed that the two groups had a significant differential prognosis (P<0.0001, Figure 9A). We also validated the 6-gene prognostic model in GSE cohorts, and observed that the samples were neatly classified into high-risk and low-risk groups with distinct OS (P<0.0001, Figure 9B). The above results indicated that the six prognostic genes could serve as predictors for predicting NSCLC prognosis.

Figure 9 Kaplan-Meier survival plots of high-risk and low-risk groups in TCGA-NSCLC dataset (A) and GSE cohorts (B). Log rank test was performed.

Risk Score is Associated with Clinical Features

To understand if there was a relation between risk score and clinical features, the distribution of risk score in different clinical features was analyzed (Figure 10). The results showed that risk score significantly varied in T1–T4 stages (P=1.9e-05, Figure 10A), N0–N2 stages (P=2.6e-05, Figure 10B), stage I–IV (P=2.4e-05, Figure 10D), and smoking status (P=0.0041, Figure 10F). Particularly, the risk score increased with the stages from I to IV. C2 subtype had a higher risk score than C1, which was consistent with the unfavorable prognosis in C2 (P<2.2e-16, Figure 10E). However, no significant difference was shown in M stage, genders, or ages (Figure 10C, G, and H). We then assessed the independence of the prognostic model in clinical use by univariate and multiple Cox regression analysis. Comparison with other clinical information demonstrated that risk score was an independent factor for clinical prediction (P<0.0001, Figure 11).

Figure 10 The relation between risk score and different features including T stage (A), N stage (B), M stage (C), stage (D), molecular subtypes (E), smoking (F), gender (G) and age (H). Kruskal-Wallis test was performed over two groups. Wilcoxon test was performed between two groups.

Figure 11 Univariate (A) and multivariate (B) Cox regression analysis of clinical features and risk score. Log rank test was performed. Abbreviation: HR, hazard ratio.

A Nomogram Based on Risk Score for Predicting Prognosis

To promote the clinical applicability of the risk score, we introduced a nomogram to serve the function of the prognostic model. T stage and risk score were used to construct the nomogram, and 1-, 3-, and 5-year survival were confirmed by the total points from T stage and risk score (Figure 12A). The predicted OS was corrected by the observed OS (Figure 12B). Then DCA was used to verify the robustness of the prognostic model and nomogram (Figure 12C). The result showed that nomogram was the most effective in predicting prognosis and was helpful in assisting clinical decision.

Figure 12 Construction of a nomogram based on risk score. (A) A nomogram based on risk score and T-stage for predicting 1-year, 3-year, and 5-year survival. (B) Correction for the nomogram by observed OS. (C) DCA for evaluating the efficiency of risk score, T-stage, and nomogram in predicting prognosis. Orange line indicates all patients receive no treatment and purple line indicates all patients receive treatment. Vertical axis represents standardized net benefit of positive treatment minus negative treatment.

Discussion

Hypoxia, which is one of the hallmarks in solid cancers, induces a tumor environment for cancer progression and metastasis. Hypoxia has been considered as a therapeutic target in lung cancer and achieved a promising outcome in clinical trials.⁸ To improve the efficiency and individuality of hypoxia-targeted therapy, we characterized the molecular features of NSCLC patients at a molecular level for helping design more personalized therapies. Therefore, in this study, we focused on hypoxia-related genes and constructed a molecular subtyping system for NSCLC.

In the two molecular subtypes, functional analysis exhibited that oncogenic pathways were highly enriched in the C2 subtype, indicating a higher activity of cancer progression in C2 compared with C1 subtype. Annotation on KEGG pathways showed that PI3K-Akt signaling, focal adhesion, proteoglycans in cancer, and p53 signaling pathways were significantly enriched in up-regulated genes in C2 subtype. The PI3K-Akt pathway can enhance the expression of hypoxia-inducible factor (HIF)-1 through growth factors such as transforming growth factor (TGF).²⁹ A number of studies suggested that hypoxia induces the EMT process, thereby promoting migration and metastasis in many cancer types.³⁰ NOTCH signaling and TGF-β signaling pathways play a modulation role in mediating EMT and hypoxia-induced cancer cell progression.^31–33 C1 and C2 subtypes manifested similar results of enriched pathways in different datasets, indicating that the subtyping based on hypoxia was reliable.

A pan-cancer research proposed that LUAD could be divided into five major immune subtypes, namely C1 (wound healing), C2 (IFN-γ dominant), C3 (inflammatory), C4 (lymphocyte depleted), and C6 (TGF-β dominant).²⁸ Using the TCGA dataset, we compared the distribution of the five immune subtypes with our subtypes. C1, C2, and C3 immune subtypes consisted of the majority, with C3 subtype showing the most favorable prognosis. Consistently, our subtype C1 with a better prognosis had a higher proportion of inflammatory (C3) subtype and a lower proportion of wound healing (C1) and IFN-γ (C2) subtypes, which indirectly proved the effectiveness of our molecular subtypes. In the above results, inflammatory features seemed to be a protective factor for NSCLC, and studies have demonstrated that high immune infiltration is correlated with a positive immune response to cancer cells in TME to some extent,³⁴ but chronic immune infiltration can inhibit effective immune response to cancer cells.³⁵

Hypoxia is a critical factor in TME regulation through modulating the expression of immune checkpoints such as PD-L1 and CD137.³⁶ Hypoxia has been found to play dual roles, such as stimulating anti-tumor effector T-cell functions and immune-suppressive functions of tumor-associated macrophages (TAMs).³⁶ We found that CD8 T-cells were higher-enriched in C1 subtype, while macrophages were more accumulated in C2 subtype. High infiltration of TAMs is considered to be associated with unfavorable prognosis in most cancers.³⁷ TAMs are widely studied in the relation to hypoxia, and evidence has revealed that hypoxia commonly co-exists with TAMs together with an increasing component of migratory factors such as CCL2, CCL5, and colony-stimulating factor 1 (CSF1).³⁸ Hypoxic TAMs release survival cytokines such as interleukin (IL)-6 and tumor necrosis factor (TNF)-α for cancer cells, thereby hindering the normal function of chemotherapeutic drugs. Furthermore, hypoxic TAMs could generate a forced tumor-promoted event in invasion, angiogenesis, and metastasis.³⁸

The above findings suggested a critical role of hypoxia in shaping TME and facilitating cancer progression. The status of TME could largely determine patients’ response to immunotherapy and chemotherapy. It could be concluded that hypoxia does play an important role in modulating TME, thereby leading to a differential response to immunotherapy and chemotherapy.

Based on hypoxia-related genes, we further identified six prognostic genes (ALDOA, EFNA1, GPC4, HOXB9, PGM2, and PLAUR) and constructed a prognostic signature. The signature manifested a strong performance in classifying NSCLC patients into high-risk and low-risk groups with a significantly distinct prognosis. Fructose-bisphosphate aldolase A (ALDOA) was reported to be an oncogene in pancreatic cancer,³⁹ and its high expression was associated with poor prognosis in colorectal cancer.⁴⁰ In lung cancer, another study found that ALDOA could induce lung cancer progression through activating the HIF-1α/MMP9 axis.⁴¹ Ephrin A1 (EFNA1) was found to be overexpressed in gastric cancer and colorectal cancer,^42,43 but it was confirmed as a protective factor to NSCLC in our study, which suggested its different roles in different cancer types. Glypicans are a family of heparan sulfate proteoglycans that interact with variable ligands such as chemokines, growth factors, and extracellular matrix.⁴⁴ Some glypicans were detected to be aberrantly expressed in cancer. Studies indicated that Glypican 4 (GPC4) was involved in tumorigenesis and migration in breast cancer and colorectal cancer.^45,46 Homeobox B9 (HOXB9) was considered as a biomarker in many cancer types, and it was commonly overexpressed in lung cancer to promote its metastasis.^47,48 Urokinase plasminogen activator receptor (PLAUR) was also identified as a biomarker in various cancers. Zhou et al⁴⁹ observed that PLAUR could induce gefitinib-resistance in NSCLC patients, indicating that it may be a potential target for chemotherapy. PGM2 was less reported in lung cancer, but it could be a new possibility to help explain the tumorigenesis mechanism in NSCLC.

In addition, a nomogram based on the prognostic signature showed promising results for its application in clinical practice. These identified prognostic genes related to hypoxia could serve as a signature for predicting NSCLC prognosis, and they were possibly involved in the cancer development and hypoxia in NSCLC.

Conclusion

In conclusion, this study proposed two novel immune subtypes (C1 and C2) and a 6-gene prognostic signature based on hypoxia-related genes. Analyses on multiple aspects demonstrated their effectiveness and validity for clinical management of NSCLC patients. Importantly, the study emphasized the important role of hypoxia in the tumorigenesis of NSCLC and the interactions of hypoxia-related pathways, TME, and targeted therapies. The immune subtypes could provide guidance to select suitable patients for immunotherapy, chemotherapy, or combined therapies. The six prognostic genes could serve as new targets for personalized therapies and to help further explain mechanisms in NSCLC development.