Development and verification of a hypoxia- and immune-associated prognosis signature for esophageal squamous cell carcinoma

Abstract

Background: Esophageal cancer is one of the most common gastrointestinal malignancies worldwide, with high morbidity and mortality in China. The clinical importance of the interaction between hypoxia and immune status in the tumor microenvironment has been established in esophageal squamous cell carcinoma (ESCC). This study aims to develop a new hypoxia- and immune-based gene signature to predict the survival of ESCC patients.

Methods: The RNA-sequencing and clinical data of 173 cases of ESCC and 271 normal tissues were obtained from The Cancer Genome Atlas (TCGA) data portal and the Genotype-Tissue Expression (GTEx) database. Hypoxia-related genes (HRGs) and immune-related genes (IRGs) were retrieved from publicly shared data. Differentially expressed gene (DEG) analyses were carried out by the DESeq2 method using the edgeR package in R. Based on the intersection of the DEGs and HRGs/IRGs, differentially expressed HRGs (DEHRGs) and differentially expressed IRGs (DEIRGs) were obtained. DEHRGs and DEIRGs associated with prognosis were evaluated using univariate Cox proportional hazards analysis. A prognostic risk score model was constructed according to the genes acquired through Cox regression. Univariate analysis and Cox proportional hazards analysis were used to determine the independent prognostic factors related to prognosis. A nomogram was developed to predict the 1-, 2-, and 3-year overall survival (OS) probability.

Results: A total of 73 intersecting genes were obtained as DEHRGs and a total of 548 intersecting genes were obtained as DEIRGs. The risk score was established using 8 genes (FABP7, TLR1, SYTL1, APLN, OSM, EGFR, IL17RD, MYH9) acquired from univariate Cox analysis. Based on this 8-gene-based risk score, a risk prognosis classifier was constructed to classify the samples into high- and low-risk groups according to the median risk score. The nomogram model was constructed to predict the OS of ESCC patients.

Conclusions: The hypoxia- and immune-based gene signature might serve as a prognostic classifier for clinical decision-making regarding individualized management, follow-up plans, and treatment strategies for ESCC patients.

Keywords: Esophageal cancer; hypoxia; immune; microenvironment; prognosis.

Figure 1

Differentially expressed hypoxia-related genes (DEHRGs) and differentially expressed immune-related genes (DEIRGs) in patients with esophageal squamous cell carcinoma (ESCC). Heat map (A) and volcano plot (B) of differentially expressed genes (DEGs) in ESCC and normal esophageal tissues. Heat map of DEHRGs (C) and DEIRGs (D) in ESCC and normal esophageal tissues.

Figure 2

Functional enrichment analysis of differentially expressed hypoxia-related genes (DEHRGs) and differentially expressed immunerelated genes (DEIRGs). (A-C) Functional enrichment of Gene Ontology (GO) terms of DEHRGs; (D) functional enrichment of Kyoto Encyclopedia of Genes and Genomes (KEGG) terms of DEHRGs; (E-G) functional enrichment of GO terms of DEIRGs; (H) functional enrichment of KEGG terms of DEIRGs.

Figure 3

Hypoxia- and immune-associated prognosis signature. (A) Forest plot of hazard ratios for 10 hypoxia- and immune-associated prognostic differentially expressed genes (DEGs). (B) Least absolute shrinkage and selection operator (LASSO) coefficient profiles of the hypoxia- and immune-associated prognostic DEGs. Each curve in the figure represents the changing trajectory of each independent variable coefficient. The Y axis indicates the value of the coefficient. The lower X axis shows log (λ) and the upper X axis shows the number of non-zero coefficients in the model. (C) Three-fold cross-validation of the tuning parameter selection in the LASSO model. The lower X axis indicates log (λ), and the upper X axis indicates the average number of genes associated with prognosis. Partial likelihood deviance values are shown, with error bars representing standard error (SE). The vertical black dotted lines are drawn at the optimal values by minimum criteria and 1 − SE criteria, which provides the best fit; (D) Forest plot of hazard ratios for 8 hypoxia- and immune-associated prognostic DEGs; (E) Distributions of risk score and survival status. The predicted value of event 1 (death) is significantly higher than event 0 (survival).

Figure 4

Validation of the effectiveness of the hypoxia- and immune-associated prognosis prediction model. (A) The high- and low-risk groups according to the median risk score; (B) Kaplan-Meier curves of overall survival (OS) in 172 patients of the training cohort based on risk score; (C) distributions of risk score and survival status; (D) survival-dependent receiver operating characteristic (ROC) curve validation of the model for prognosis; (E) heat map of the expression distribution of 8 genes in the prognosis model. Red indicates high expression and blue indicates low expression.

Figure 5

Validation of the stability of the hypoxia- and immune-associated prognosis prediction model. Kaplan-Meier plot of overall survival by risk groups according to smoking status (A); more than 1 cigarette per day (B); body mass index (BMI) ≤24 kg/m² (C); patients who received radiotherapy (D); N0 staging (E); N1 staging (F); M0 staging (G); and male gender (H).

Figure 6

Relationship between the risk score and immune cell infiltration. (A,B) Resting dendritic cells; (C,D) naive CD4⁺ T cells; (E,F) activated mast cells; (G,H) regulatory T cells; (I,J) neutrophils. *P<0.05, **P<0.01, and ****P<0.0001. NES, normalized enrichment score.

Figure 7

Development and evaluation of the nomogram for predicting overall survival. (A) Nomogram based on the risk score and clinicopathological characteristics to predict the 1-, 2-, and 3-year overall survival probability. Calibration of the nomogram according to the consistency between the predicted and the actual results. (B) The nomogram depicts curves relative to the black line, suggesting perfect prediction.