Novel hypoxia-related gene signature for predicting prognoses that correlate with the tumor immune microenvironment in NSCLC

Abstract

Background: Intratumoral hypoxia is widely associated with the development of malignancy, treatment resistance, and worse prognoses. The global influence of hypoxia-related genes (HRGs) on prognostic significance, tumor microenvironment characteristics, and therapeutic response is unclear in patients with non-small cell lung cancer (NSCLC). Method: RNA-seq and clinical data for NSCLC patients were derived from The Cancer Genome Atlas (TCGA) database, and a group of HRGs was obtained from the MSigDB. The differentially expressed HRGs were determined using the limma package; prognostic HRGs were identified via univariate Cox regression. Using the least absolute shrinkage and selection operator (LASSO) and multivariate Cox regression, an optimized prognostic model consisting of nine HRGs was constructed. The prognostic model's capacity was evaluated by Kaplan‒Meier survival curve analysis and receiver operating characteristic (ROC) curve analysis in the TCGA (training set) and GEO (validation set) cohorts. Moreover, a potential biological pathway and immune infiltration differences were explained. Results: A prognostic model containing nine HRGs (STC2, ALDOA, MIF, LDHA, EXT1, PGM2, ENO3, INHA, and RORA) was developed. NSCLC patients were separated into two risk categories according to the risk score generated by the hypoxia model. The model-based risk score had better predictive power than the clinicopathological method. Patients in the high-risk category had poor recurrence-free survival in the TCGA (HR: 1.426; 95% CI: 0.997-2.042; p = 0.046) and GEO (HR: 2.4; 95% CI: 1.7-3.2; p < 0.0001) cohorts. The overall survival of the high-risk category was also inferior to that of the low-risk category in the TCGA (HR: 1.8; 95% CI: 1.5-2.2; p < 0.0001) and GEO (HR: 1.8; 95% CI: 1.4-2.3; p < 0.0001) cohorts. Additionally, we discovered a notable distinction in the enrichment of immune-related pathways, immune cell abundance, and immune checkpoint gene expression between the two subcategories. Conclusion: The proposed 9-HRG signature is a promising indicator for predicting NSCLC patient prognosis and may be potentially applicable in checkpoint therapy efficiency prediction.

Keywords: NSCLC; hypoxia; immune microenvironment; immunotherapy; prognosis.

FIGURE 1

Diagrammatic representation of the analytical process used in this study.

FIGURE 2

Building of the HRG model in the TCGA-NSCLC cohort. (A) Volcano graphic showing how HRG expression differs in tumor and normal samples based on both fold change and statistical significance. (B) Results of univariate Cox regression analysis presented in a forest plot, with gene name in the first column and box plot depicting the hazard ratio (HR) and its corresponding 95% confidence interval in the second column. (C,D) LASSO analysis for estimating the number of contributing components. (E) Correlation analysis of HRG expression in the risk model using Spearman’s correlation. *p < 0.05; **p < 0.01; ***p < 0.001. (F) Results of protein interaction analysis of the HRG model.

FIGURE 3

Prediction ability of the HRG model for the recurrence and stage of NSCLC patients. (A–C) TCGA cohort Fisher test outcomes for T and N stages and pathological stages comparing two risk categories. (D) Analyzing the differences between two risk categories on three known hypoxia scores. *p < 0.05; **p < 0.01; ***p < 0.001. (E,F) Recurrence survival curves for two risk categories. E: TCGA cohort, F: GEO cohort.

FIGURE 4

Prognostic data analysis from the training and validation cohorts as well as verification of the model’s precision. (A,B) The risk score curve is displayed in the top part. The distribution of the risk score, survival duration, and patient status are presented in the middle part. A heatmap of HRGs in the classifier is shown in the bottom part. (A) TCGA cohort, (B) GEO cohort. (C,D) NSCLC patient Kaplan‒Meier survival curve comparing two risk categories. (C) TCGA cohort, (D) GEO cohort. (E,F) ROC curve in the TCGA cohort (E) and GEO cohort (F).

FIGURE 5

Independent prognostic factor determination and predictive accuracy comparison. (A,B) Results of TCGA cohort (A) and GEO cohort’s multivariate Cox regression analysis (B). (C,D) Concordance index curve of 3 clinical parameters and risk scores for OS time from 1 to 5 years, (C) TCGA cohort, (D) GEO cohort. (E,F) Multi-index ROC curve of the risk score and other clinical parameters for 3-year OS time, TCGA cohort (E) and GEO cohort (F).

FIGURE 6

Clinical prognostic nomogram was created and validated. (A) Nomogram, taking into account risk score, tumor stage, age, and sex, predicted the likelihood of 1-, 3-, and 5-year OS. (B–D) 1-, 3-, and 5-year OS calibration curves; predicted survival probability graphed along the x-axis, whereas actual survival probability is represented along the y-axis.

FIGURE 7

Difference in biological pathways and cancer-related gene sets between the two risk categories. (A) Heatmap of GSVA enrichment of the KEGG pathway. Red represents high enrichment scores, while blue represents low enrichment scores. (B) Bar graph of hallmark gene set enrichment score. The color indicates the significance of the difference, and the x-axis represents the enrichment score fold change in this hallmark gene set between high- and low-risk categories.

FIGURE 8

Association analysis of HRGs with TIME in the TCGA cohort. (A) Correlation analysis of HRGs in risk signature and immune cell abundance. (B) Comparison of immune score in the two subcategories. (C) Quantitative analysis of immune cell abundance between patients in the two risk categories. (D) HRGs and immunomodulator-related gene correlation analysis. Red represents a positive correlation, while blue represents a negative correlation in the heatmap. *p < 0.05; **p < 0.01; ***p < 0.001.

FIGURE 9

Gene expression comparison relating to immunological checkpoints in the two categories. Middle line of the box represents the median of the data, while the upper and lower limits of the box represent the upper and lower quartiles of the data, the line extending from the box represents 1.5 times the interquartile range (IQR) from the upper and lower quartiles. *p < 0.05; **p < 0.01; ***p < 0.001.