Identification of differentially expressed hypoxia genes associated with OS
After data preprocessing, a total of 3166 DEGs in the GEO set and 4808 DEGs in TCGA were obtained. These are illustrated in the volcano plots in Figure 1A and 1B and heatmaps in Supplementary Figure 1. We then intersected the DEGs with the list of 270 HRGs obtained from the collection of the CancerSEA database and MSigDB database to obtain differentially expressed hypoxia genes. As shown in the Venn diagrams in Figure 1C, 30 unique elements (15 upregulated and 15 downregulated) were identified and consistently differentially expressed in the GEO training set and TCGA validation set (Supplementary Table 1).
A subsequent PPI network was constructed by uploading the aforementioned 30 target genes to the STRING online database and visualized using Cytoscape to elucidate their interactions (Figure 1D). In the functional analysis of DEGs, the top 6 enriched GO annotations and KEGG pathways in each category were visualized intuitively as Circos plots (Figure 1E). The expression data of 5 real hub genes extracted from univariate analysis (P < 0.05) were used to assess the impact of these variables on the survival outcome and were visualized with forest plots. As shown in Figure 1F, SLC2A1, PGM2, SULT2B1, and CA9 all had a lower risk of death except BGN, suggesting that BGN may be expressed at higher levels in tumor tissues than in normal tissues.
Consensus clustering to identify distinct subgroups and intercluster prognosis analysis
A total of 179 tumor samples from the GSE53625 dataset were consistently clustered based on expression similarity, and the correlation between hypoxic conditions and clinical phenotypes was investigated. After executing ConsensusClusterPlus, we obtained the cluster consensus and item-consensus results. Considering the minimum descent slope of the CDF and Delta area plots (Figure 2A-B), although K = 2 is not the last inflection point that minimizes the area under the CDF curve, the consensus matrix helps us to determine the clearest division when divided into two clusters, named C1 and C2. Compared to those of the other categorical numbers (Supplementary Figure 2), the consensus matrix graph corresponding to K = 2 showed that the distribution of 2 blue blocks on the diagonal along the white background was well defined (Figure 2C).
The results of Kaplan-Meier survival analysis revealed significant differences in prognosis among cluster 1 and cluster 2 (P = 0.014). As shown in Figure 2D, the samples in cluster 2 had better performance for overall survival than those in cluster 1. Sorting the samples by cluster produced the gene expression heatmap shown in Supplementary Figure 3, which indicated the composition and quantity of clustering. Gene expression patterns differed among the subgroups, suggesting the credibility of the 2 structural clusters. A heatmap corresponding to the dendrogram in Figure 2E, annotated by grade, T stage, N stage, stage, sex, and 5 real hub genes as mentioned previously, shows the distribution of these features based on cluster 1 or 2 intuitively. This result demonstrated the heterogeneity between the two clusters.
Identifying the correlations between the obtained clusters and immune infiltration
We calculated the levels of 22 immune cell types in each sample and compared their differences in the C1 and C2 subgroups (Figure 2F). The results showed that 11 types of immune cells, including CD8 T cells, resting memory CD4 T cells, follicular helper T cells, regulatory T cells (Tregs), activated NK cells, monocytes, M0 macrophages, activated dendritic cells and activated mast cells, had significantly different infiltration levels in different subgroups. Then, we performed GSEA to explore the underlying regulatory mechanism that contributed to the differences in the tumor immune microenvironment between the 2 clusters, and the results showed that they were mainly related to the activation of P53 signaling pathway molecules (Figure 2G). These results indicated that each cluster had unique gene expression and pathway characteristics, and hypoxic cells may potentially regulate the tumor immune microenvironment between the two clusters through the P53 signaling pathway.
Identification of modules associated with hypoxia by WGCNA
After sample clustering to detect outliers, the WGCNA was then restricted to 178 patients from GSE53625. Different power values (1–20) were analyzed, and the best power value of β = 6 (scale-free R2= 0.95) was chosen for soft thresholding for subsequent coexpression scale-free network construction (Supplementary Figure 4A-C). Eight gene modules were obtained and visualized for the next analysis (Figure 3A). Module-feature relationship analysis revealed that the green module was significantly associated with consensus subgroups (C1, r = -0.46, P = 1e-10; C2, r=0.46, P = 1e-10, Supplementary Figure 4D-E), suggesting that the module is suitable for identifying the hub genes associated with C1/C2. Therefore, the green module was chosen for further analysis.
Constructing and evaluating a hypoxia-related prognosis signature
The results of differential expression analysis among the obtained clusters are shown in Figure 3B. After further screening for the most significant hub genes, 259 differentially expressed HRGs (Figure 3C) resulted from the intersection of 443 green module genes and 731 DEGs from consensus clustering analysis. The hypoxia gene signature with the best prognostic performance was constructed using the stepAIC algorithm based on the expression level of the above genes, identifying the risk score by using the 6 most relevant genes, including PNPLA1, CARD18, IL-18, SLC37A2, ADAMTS18 and FAM83C (Figure 4). The risk score was calculated as follows: risk score = (0.066265*PNPLA1 expression) + (-0.149270*CARD18 expression) + (-0.183367*IL18 expression ) + (-0.037724079*SLC37A2 expression) + (0.119388782*ADAMTS18 expression) + (-0.031834954*FAM83C expression). The samples were assigned a risk score and ordered to determine whether the expression level varied systematically with the risk score (Figure 4A-B). The expression levels of the 6 genes distinctly decreased as the risk score increased, except for ADAMTS18 (Figure 4C). After the risk curves were plotted based on the risk score versus patient survival status and divided into the high- and low-risk groups based on the median risk score, a higher percentage of patient deaths was associated with high-risk patients, and a higher percentage of long-term survival was associated with low-risk patients (P = 0.0063). Furthermore, the prognoses differed significantly between the two groups, as shown in Figure 4D. The results of ROC curve analysis using the risk score calculated for each sample are shown in Figure 4E. The AUCs for 1-, 3-, and 5-year prognostic prediction were 0.71, 0.68 and 0.71, respectively, indicating the relatively excellent predictive efficacy of the model.
Then, the established prognostic signature was further validated in the test group, including the internal validation set and all GSE53625 and TCGA datasets. Expression level profiles for the 6 selected genes were obtained from 3 testing group samples, and the risk scores were calculated using the abovementioned method for each patient. Consistent with the above results, sorting the samples by risk score produced the heatmap shown in Figure 5A-C and Supplementary Figures 5A-C and 6A-C, which indicated that the risk score increased as the expression levels decreased, except for ADAMTS18. The Kaplan-Meier survival curves (Figure 5D) showed that this risk model could effectively distinguish the survival of the high- and low-risk groups in the validation set (P = 0.0014 for the internal validation set, P < 0.001 for all datasets, and P = 0.043 for the TCGA dataset). In the external validation group, the AUCs of the 1-, 3-, and 5-year overall survival predictions for the risk scores were 0.64, 0.78, and 0.79, respectively (Figure 5E). The efficiency results of the risk score classification for prognosis prediction at 1, 3, and 5 years are presented. These results were consistent with those obtained from the training dataset, demonstrating the relatively excellent predictive accuracy and stability of our risk score model.
Evaluating the independent role of the prognostic signature and building a predictive nomogram for OS prediction
To verify the use of these candidate prognostic genes as independent biomarkers, univariate and multivariate Cox regression analyses were used to evaluate whether the predictive value of the model-incorporated prognostic signature was affected by other clinical factors (Figure 6A-B). The results indicated that age (P = 0.010, HR = 1.011), stage, location and risk score (P < 0.001, HR = 1.869) were independent prognostic factors for OS. Taking these results of the univariate and multivariate Cox analyses into consideration, the risk score integrated with age and stage was finally chosen together to construct a nomogram model, as presented in Figure 6C. The survival probability for the individuals at 1, 3, and 5 years was obtained through the function conversion relationship of the total scores. The calibration plot of the nomogram (Figure 6D) showed better consistency between the predicted OS outcomes and actual observations, indicating a good predictive performance of the hypoxia-related prognosis model.
Correlation analysis of the risk score with clinicopathological features and immune infiltration
The characteristics of the clinicopathological features in the high- and low-risk groups are shown by heat maps (Figure 7A), which were closely correlated with the clinical phenotypes. We then analyzed the relationship between the risk score and immune cell infiltration as well as the tumor microenvironment in ESCC. The results of the correlation analysis between the immune-related score and the risk score (Figure 7B, Supplementary Table 2) showed that the stromal and immune scores were comparable between the high- and low-risk groups. However, the estimated score was higher in the high-risk group than in the low-risk group, and the difference was statistically significant (P = 0.03). Differences in the infiltration of 22 immune cell subtypes in the high- and low-risk groups are shown in Figure 7C. GSEA between the high- and low-risk groups showed that the most significantly enriched pathway was the P53 signaling pathway (Figure 7D). Notably, we found consistent results for the Wnt pathway and hypoxia pathway (Figure 7E and 7F). The results reconfirmed that esophageal cancer cells in a hypoxic state could affect the tumor immune microenvironment through an underlying potential regulatory mechanism.
Validation of the expression of selected HRGs
To verify the accuracy of the abovementioned HRGs obtained from previous analyses, we further detected the protein expression levels of BGN and IL-18 according to previous publications and antibody availability. The clinical details of the 232 patients involved are presented in Table 1. The results of IHC-based staining of TMA showed that BGN was highly expressed in these ESCC samples and mainly localized to the cytoplasm of cancer cells (Figure 8A-D). IL-18 is normally expressed in cancerous tissues, with a significantly higher percentage of expression in normal tissues than in tumor tissues (Figure 8E-H). Moreover, the survival analysis suggested that patients with higher BGN expression were predicted to have poorer survival (Figure 8I). In contrast, patients with high expression of IL-18 showed better overall survival (Figure 8J). The chi–squared test indicated that high BGN expression was notably associated with tumor size, tumor invasion, and lymph node metastasis (P < 0.05 for all) of patients with ESCC (Table 1). For IL-18, no significant correlation between the IL-18 level and clinicopathological factor except lymph node metastasis (P < 0.001), was observed.