From the UCSC Xena database, NSCLC clinical data and transcriptome expression profiles for LUAD and LUSC patients in The Cancer Genome Atlas (TCGA) were downloaded. A total of 1,011 tumor and 108 normal samples were employed for HRG expression differential analysis and as a training set for the construction of the risk score model. All scores indicating the degree of hypoxia in the tumor, such as the Winter, Buffa, and Ragnum hypoxia scores (Winter et al., 2007; Buffa et al., 2010; Ragnum et al., 2015), were obtained from TCGA. The higher these scores, the more hypoxic is the tumor.

For external validation, four microarray datasets together with detailed clinical characteristics were retrieved from the Gene Expression Omnibus (GEO) (GSE30219, GSE31210, GSE37745, and GSE50081) using the GEOquery R package. A total of 628 tumor samples were acquired, and the clinical characteristics are listed in Supplementary Table S1. Batch effects were eliminated using the sva package (Leek et al., 2012), and the scale function was used to standardize the expression value before model construction and validation.

A total of 190 of 200 HRGs and available data were retrieved from the MSigDB Hallmark database (Supplementary Table S2), and the expression of these genes was extracted from integrated datasets (Liberzon et al., 2015). Differential analysis between tumor and normal samples of 190 HRGs was performed using the limma package in the TCGA database. Genes with p values less than 0.05 were defined as differentially expressed HRGs (DEHRGs). Volcano mapping was then performed using the ggplot2 package.

We conducted univariate Cox regression analysis based on overall survival to determine prognostic DEHRGs (p < 0.05). The DEHRGs related to OS were identified by least absolute shrinkage and selection operator (LASSO) regression analysis using the glmnet package to shrink the gene set associated with prognosis. To achieve the optimal model, we combined the remaining genes step by step and tested performance. Ultimately, nine genes were selected, and a multivariate Cox regression analysis was adopted to build a prognostic risk model. The algorithm used to determine the patient’s risk level was as follows: R i s k   s c o r e = 0.10515 × S T C 2 exp + 0.0234 × A L D O A exp + 0.02207 × M I F exp + 0.13426 × L D H A exp + 0.09928 × E X T 1 exp − 0.01151 × P G M 2 exp − 0.1475 × E N O 3 exp − 0.16294 × I N H A exp − 0.14694 × R O R A exp

A risk score was then obtained for all patients. Patients were separated into two risk categories by using the median risk score of patients.

To determine the predictive ability of this model for tumor recurrence, recurrence-free survival was compared between two risk categories using stage I and II cases in the TCGA and GEO databases.

Overall survival (OS) was also compared between the two risk categories using all patients. The survivalROC R package was employed to plot the receiver operating characteristic (ROC) curves of 1-, 3-, and 5-year OS in the training cohort, and the same method was used in the validation cohort.

Multivariate Cox regression analysis was applied to construct an independent prognostic model. A nomogram was used to predict patient survival rates, and the performance was verified by drawing calibration curves. ROC curves and the concordance index were introduced to detect the reliability of the HRG model. Clinical features were incorporated into multivariate Cox regression to determine whether the risk model was an independent indicator of prognostic prediction.

The gene sets of the Kyoto Encyclopedia of Genes and Genomes and the Hallmark gene sets were downloaded from MSigDB. The enrichment scores of individual patients were computed with the GSVA R package. Using the limma package, the enrichment score was then compared between two TCGA risk categories, and the difference in HRG expression between the two categories was also analyzed. A cluster heatmap was used to show the top 30 genes with the most significant p values.

Immunomodulators, a group of immunoregulatory genes containing antigen-presenting factors, ligands, and receptors, play crucial roles in cancer immunotherapy, and the expression of HRGs and immunomodulators was compared with the Spearman correlation method.

In addition, the differences in their infiltration abundance in the tumor microenvironment affected the patients’ outcomes and the efficiency of immunotherapy. Therefore, the ssGSEA algorithm was employed to compute immune cell abundance in individual patients. Using the Spearman correlation method, the correlation between the hypoxia-related model and immune cells was analyzed, and, finally, the difference in immune cell abundance between the two risk categories was assessed.

Immunomodulator-related genes and immune cell-type marker genes were obtained from the TCGA immune response working group (Charoentong et al., 2017; Thorsson et al., 2018). The immunological score was computed with the estimate package, and we compared the difference between the two risk categories.

Kaplan‒Meier curves and the log-rank test were used to compare the survival rate between various subcategories. The independent prognostic significance of the clinical features in OS was detected by univariate and multivariate Cox proportional hazard regression analyses. The prognostic prediction significance of the risk models for 1/3/5-year OS was evaluated by ROC curves and area under the curve (AUC) values. Differences in TMN stage and clinicopathologic stage composition between the two risk categories of the TCGA cohort were analyzed by Fisher’s test. The correlation analysis of HRGs with immune cell abundance and with immune regulator-related genes was performed using Spearman’s correlation analysis. The Wilcoxon test was applied to compare immune cell infiltration. p values less than 0.05 were defined as significant. R 3.6.1 was used to conduct all the analyses.

Genes that are expressed normally in both tumor and normal tissues do not typically play a significant role in the development and progression of tumors. Therefore, it is necessary to exclude such genes from the analysis beforehand.

Differential analysis between 1,011 tumor and 108 normal samples on 190 HRGs was performed, and 160 genes were significantly differentially expressed (Supplementary Table S3) in the TCGA cohort. Of these, 72 genes were downregulated in tumor tissues, while 88 were upregulated. The differentially expressed volcano map is shown in Figure 2A.

The analysis revealed that a majority of HRGs showed significant differences in expression, highlighting the importance of the hypoxia pathway as a mechanism underlying tumor growth and metastasis.

Although there may be aberrant gene expression, it may not necessarily affect a patient’s survival status or time. Thus, the identification of hypoxia-related genes that are linked to prognosis is crucial.

In total, 990 NSCLC cases in the TCGA cohort with detailed OS information were identified. To further identify prognosis-related DEHRGs, 160 DEHRGs were analyzed using univariate Cox regression analysis, and 41 genes were significantly correlated with prognosis (Supplementary Table S4). A forest map was used to display the prognosis-related HRGs’ hazard ratios and the matching confidence intervals (CIs), which demonstrated that the majority of candidates were at risk (Figure 2B). To ensure the stability and feasibility of the model, LASSO regression analysis was performed with the optimum λ value to reduce the scope of candidates, and 25 DEHRGs were identified (Figures 2C, D). Finally, by performing multivariate Cox regression, nine genes (STC2, ALDOA, MIF, LDHA, EXT1, PGM2, ENO3, INHA, and RORA) were selected to construct a risk prognostic model. In the TCGA training cohort, patients were split into two risk categories based on the median risk score. Notably, compared to patients in the low-risk category, those in the high-risk category were more likely to experience relapse (HR: 1.426; 95% CI: 0.997–2.042; p = 0.046; Figure 3E).

A linear regression analysis was conducted to determine the interaction of these nine genes; the results indicated that 25 candidates exhibited high correlation, such as ALDOA and LDHA, ALDOA and STC2, ALDOA and MIF, and ENO3 and EXT1 (Figure 2E). In addition, protein interaction analysis was performed on these nine genes using the STRING online database and visualized by Cytoscape. We observed that ALDOA, LDHA, ENO3, and PGR2 strongly interacted with multiple genes—in the end, these four candidates were regarded as the core genes (Figure 2F).

Fisher’s test showed that the percentage of stage I samples, T1 samples, and N0 samples in the low-risk category (stage: 59.7%; T: 36.1%; N: 72.5%) was greater than the percentage in the high-risk category (stage: 43.5%; T: 20.9%; N: 57.9%) (Figures 3A–C). The difference in M0 staging in the low-risk category was not significant in comparison with the low-risk category (Supplementary Figure S1). A higher risk score indicates higher expression levels of most HRGs, which are associated with a more advanced disease stage, a larger tumor size, and increased lymph node involvement. One possible explanation for the higher proportion of stage I, T1, and N0 samples in the low-risk category compared with the high-risk category is that HRGs are expressed more strongly in later-stage and more invasive tumors, which may lead to higher risk scores and worse prognoses. In addition, the lack of a significant difference in M0 staging between the low- and high-risk categories may be due to distant metastasis being less affected by hypoxia than local tumor growth and invasion. The Buffa, Winter, and Ragnum hypoxia scores in the high-risk category were greater than those in the low-risk category (Figure 3D), implying that the established model could effectively evaluate the degree of hypoxia.

Moreover, the scatter plot of OS status showed that, in comparison with low-risk category patients, more deaths were observed in the high-risk category (Figure 4A). More importantly, high-risk category patients (n = 495) had a worse prognosis than low-risk patients (n = 495) (p < 0.0001; HR: 1.778; 95% CI: 1.457–2.169; Figure 4C). To accurately evaluate the prognostic performance of the nine-gene signature, ROC analysis was conducted using the 1-, 3-, and 5-year cutoff values of OS—the AUC values were 0.679, 0.654, and 0.588, respectively, indicating that our prognostic model had good predictive performance in patients with NSCLC (Figure 4E).

Early-stage cancer patients are more highly represented in the low-risk category based on HRGs than in the high-risk category. Thus, when we specifically only analyzed early-stage patients, we discovered that significant differences in overall survival still existed between high- and low-risk categories in early-stage cancer prognostic indicators (Supplementary Figure S2). This finding suggests that risk stratification based on HRGs can still be used to effectively predict the overall survival of NSCLC patients, even when analyzing only the subset of patients in the early stages of the disease.

To summarize, a model based on the expression of HRGs that are prognostically relevant can effectively assess the extent of hypoxia, exhibit a close correlation with clinical staging, and reveal significant differences in recurrence, survival, and overall survival between high- and low-risk categories.

External validation plays a crucial role in improving the reliability and reproducibility of scientific research. By confirming the generalizability and applicability of research findings, external validation can help to identify and correct potential problems and limitations in study design and methodology.

To evaluate the prognostic capability of this novel HRG-based risk model, external cohorts were obtained from the GEO. The risk score was computed via the identical formula described before. The patients were separated into two risk categories using the median score. In recurrence-free survival analysis, patients in the high-risk category had a higher tendency to relapse than low-risk patients (HR: 2.374; 95% CI: 1.743–3.233; p < 0.0001; Figure 3F).

Scatter plots for death events are displayed in Figure 4B. High-risk category patients showed higher mortality rates than patients in the low-risk category. Additionally, we confirmed the performance of our HRG risk model for OS in the GEO cohort. As expected, a worse OS was indicated by a higher risk score (p < 0.0001, HR: 1.81, 95% CI: 1.428–2.295; Figure 4D). Figure 4F shows the 1-, 3-, and 5-year cutoff OS values in the validation dataset; the AUC values were 0.668, 0.643, and 0.623, respectively. ROC analysis further verified the effectiveness and reliability of the HRG risk model for prognostic prediction in NSCLC patients.

In our external dataset validation analysis, we found that our model demonstrated strong predictive performance for both overall survival and recurrence-free survival, indicating that it is both stable and generalizable.

Patients with complete clinical information were identified to construct risk prognostic models; Table 1 displays the clinical information. Along with risk scores, clinical stage, age, and sex were considered in a multivariate Cox regression analysis to determine independent prognostic markers. According to the outcome, the risk score was an independent prognostic factor in the TCGA (HR: 0.58, 95% CI: 0.47–0.71; Figure 5A) and GEO cohorts (HR: 0.65, 95% CI: 0.50–0.84; Figure 5B).

We also calculated the conformance index (C-index) of these three clinical factors and risk score within 1–5 years of OS time. For the TCGA cohort, the results suggested that the C-index of the risk score was greater than that of other clinical factors (Figure 5C), and that the C-index had a similar numerical approximation between the risk score and clinicopathological stage in the GEO cohort (Figure 5D).

The ROC curve of the risk score was detected with these three clinical factors for the 3-year OS. The AUC value of the risk score was greater than that of the other three clinical characteristics in both cohorts (Figures 5E, F).

Taken together, even when other factors that may affect patient prognoses were taken into account, the risk score had a significant and distinct impact on patient outcomes. It performed better than other clinical parameters in predicting overall survival.

To assist clinical practitioners with more accurate predictions of patient survival rates during treatment and management, we constructed a clinical nomogram that serves as an important auxiliary decision-making tool.

Using the “rms” and “survival” R packages, we generated a nomogram containing the risk score and valuable clinical parameters predicting 1-, 3-, and 5-year overall survival of NSCLC patients in the TCGA dataset and a calibration plot that analyzed estimated and actual overall survival to assess the effectiveness of the prognostic nomogram (Figure 6A). The calibration plot showed that the nomogram’s estimated survival and actual survival rates were highly correlated (Figures 6B–D). It is possible to evaluate the probability of a patient surviving at a certain time according to the patient’s clinical information and gene expression levels using this clinical prognosis nomogram. In conclusion, this nomogram system for predicting the overall survival of patients based on gene expression and clinical information is reliable.

Different biological mechanisms often lead to distinct prognostic outcomes in cancer patients. KEGG pathway and cancer-related hallmark gene set analyses were performed to clarify the biological mechanisms that were related to the risk model.

In the low-risk category, several metabolically related pathways, such as histidine, fatty acid, and ether lipid metabolism, were highly enriched (Figure 7A). The metabolism of these amino acids and lipids was reported to be closely associated with tumors (Currie et al., 2013; Li and Zhang, 2016; Bian et al., 2021). Interestingly, immune-related pathways were observed to cluster in the low-risk category, such as the T-cell and B-cell receptor signaling pathways. Conversely, for the high-risk category, cell proliferation-related pathways such as the G2/M checkpoint, E2F targets, and MYC targets were enriched. It is worth noting that hypoxia was also identified in the high-risk category, as anticipated (Figure 7B). The expression heatmap of the 30 most significant DEHRGs between the two risk categories is shown in Supplementary Figure S3. It is evident from the heatmap that there is a significant difference in the expression levels of these genes.

In summary, there are significant mechanistic differences between the low-risk and high-risk categories, with the immune features of the former and proliferation characteristics of the latter being consistent with their respective prognostic outcomes.

To further explore the risk score and immunity, we performed a correlation analysis of immune cell abundance and HRG expression in the risk model. The findings revealed an extensively significant negative correlation between immune cell abundance and HRGs, such as ENO3, INHA, and ALDOA. RORA was positively correlated with most immune cell types (Figure 8A).

In the TCGA cohort, immune cell abundance was calculated between the two risk categories. The abundance of activated CD8 T cells, activated B cells, effector memory CD8 T cells, and central memory CD4 T cells was higher in the low-risk category (Figure 8C). In addition, the immune score of the low-risk category was significantly greater than that of the high-risk category (Figure 8B).

To illustrate the relationship between the HRG model and the immune microenvironment, we generated the correlation pattern of the HRG model and immune-modulator-related genes. The results demonstrated that an extensively significant negative correlation existed, including STC2, ALDOA, MIF, LDHA, EXT1, PGM2, ENO3, and INHA. Only RORA was positively correlated with most immune modulators (Figure 8D).

The low-risk category was associated with high immune cell infiltration, immune-modulator-related gene expression, and immune scores, which are indicative of favorable clinical outcomes related to preexisting anticancer immunity.

In consideration of the therapeutic significance of treatment approaches related to immune checkpoint blockade in NSCLC, the correlation between the risk score and several immune checkpoints, such as PD1, PDL1, and CTLA4, were analyzed. In the low-risk category, we found a substantial increase in the expression of PD1 and CTLA4, while the expression of PDL1 was not significantly different (Figure 9).