Identification of Hypoxia-Related Subtypes, Establishment of Prognostic Models, and Characteristics of Tumor Microenvironment Infiltration in Colon Cancer

Abstract

Background: Immunotherapy is a treatment that can significantly improve the prognosis of patients with colon cancer, but the response to immunotherapy is different in patients with colon cancer because of the heterogeneity of colon carcinoma and the complex nature of the tumor microenvironment (TME). In the precision therapy mode, finding predictive biomarkers that can accurately identify immunotherapy-sensitive types of colon cancer is essential. Hypoxia plays an important role in tumor proliferation, apoptosis, angiogenesis, invasion and metastasis, energy metabolism, and chemotherapy and immunotherapy resistance. Thus, understanding the mechanism of hypoxia-related genes (HRGs) in colon cancer progression and constructing hypoxia-related signatures will help enrich our treatment strategies and improve patient prognosis. Methods: We obtained the gene expression data and corresponding clinical information of 1,025 colon carcinoma patients from The Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO) databases, respectively. We identified two distinct hypoxia subtypes (subtype A and subtype B) according to unsupervised clustering analysis and assessed the clinical parameters, prognosis, and TME cell-infiltrating characteristics of patients in the two subtypes. We identified 1,132 differentially expressed genes (DEGs) between the two hypoxia subtypes, and all patients were randomly divided into the training group (n = 513) and testing groups (n = 512). Following univariate Cox regression with DEGs, we construct the prognostic model (HRG-score) including six genes (S1PR3, ETV5, CD36, FOXC1, CXCL10, and MMP12) through the LASSO-multivariate cox method in the training group. We comprehensively evaluated the sensitivity and applicability of the HRG-score model from the training group and the testing group, respectively. We explored the correlation between HRG-score and clinical parameters, tumor microenvironment, cancer stem cells (CSCs), and MMR status. In order to evaluate the value of the risk model in clinical application, we further analyzed the sensitivity of chemotherapeutics and immunotherapy between the low-risk group and high-risk group and constructed a nomogram for improving the clinical application of the HRG-score. Result: Subtype A was significantly enriched in metabolism-related pathways, and subtype B was significantly enriched in immune activation and several tumor-associated pathways. The level of immune cell infiltration and immune checkpoint-related genes, stromal score, estimate score, and immune dysfunction and exclusion (TIDE) prediction score was significantly different in subtype A and subtype B. The level of immune checkpoint-related genes and TIDE score was significantly lower in subtype A than that in subtype B, indicating that subtype A might benefit from immune checkpoint inhibitors. Finally, an HRG-score signature for predicting prognosis was constructed through the training group, and the predictive capability was validated through the testing group. The survival analysis and correlation analysis of clinical parameters revealed that the prognosis of patients in the high-risk group was significantly worse than that in the low-risk group. There were also significant differences in immune status, mismatch repair status (MMR), and cancer stem cell index (CSC), between the two risk groups. The correlation analysis of risk scores with IC₅₀ and IPS showed that patients in the low-risk group had a higher benefit from chemotherapy and immunotherapy than those in the high-risk group, and the external validation IMvigor210 demonstrated that patients with low risk were more sensitive to immunotherapy. Conclusion: We identified two novel molecular subgroups based on HRGs and constructed an HRG-score model consisting of six genes, which can help us to better understand the mechanisms of hypoxia-related genes in the progression of colon cancer and identify patients susceptible to chemotherapy or immunotherapy, so as to achieve precision therapy for colon cancer.

Keywords: HRG-score; colon cancer; hypoxia-related genes; immune checkpoint blockade; immunotherapy; molecular subtype; tumor microenvironment.

FIGURE 1

HRG subtypes and clinical parameters and biological characteristics of two distinct subtypes of samples divided by consistent clustering. (A) Consensus matrix heatmap defining two subtypes (k = 2) and their correlation area. (B) PCA showing a remarkable difference in transcriptomes between the distinct HRG-subtypes, and each dot represents a single sample. (C) KM survival curve analysis showed that the overall survival time of the distinct HRG-subtypes was different (log-rank tests, p < 0.001). (D) Differences in clinical parameters and HRG expression levels between the two distinct HRG-subtypes.

FIGURE 2

Correlations of tumor immune cell microenvironments and two HRG-subtypes. (A) GSVA of biological pathways between two distinct subtypes, in which red and blue represent activated pathways and blue represents inhibited pathways, respectively. (B) Relative abundance of 23 infiltrating immune cell types in the two HRG-subtypes. (C)Correlations between the two CRC subtypes and TME score. (D–I) Expression levels of PD-L1, PD-L2, PDCD1, LAG3, TIGIT, and CTLA4 in two distinct HRG-subtypes. (*p < 0.05; **p < 0.01; ***p < 0.001).

FIGURE 3

Construction of the HRG-score in the training set based on the differentially expressed genes of two distinct HRG-subtypes. (A) Volcano plot of differentially expressed genes between the two distinct HRG-subtype. Gray dots represent not significant genes, green dots represent upregulated genes in HRG-subtype A, and red dots represent upregulated genes in HRG-subtype B (B–C) GO and KEGG enrichment analyses of DEGs among two distinct HRG-subtypes. (D–E) LASSO regression analysis and partial likelihood deviance on the prognostic genes. (F) Forest plot of multivariate cox regression analysis for prognostic genes. (G–H) Ranked dot and scatter plots showing the HRG-score distribution and patient survival status. (I) Heatmap of the expression of six genes involved in the HRG-score in low- and high-risk groups. (J) Survival analysis of the patients in low- and high-risk groups. (K) PCA based on the prognostic signature. (L) ROC curves to predict the sensitivity and specificity of 1-, 3-, and 5-year survival according to the HRG-score.

FIGURE 4

Validation of the HRG-score signature in the testing set. (A,B) Ranked dot plot indicates the PRG-score distribution, and the scatter plot presents the patients’ survival status. (C) Heatmap of the expression of six genes involved in the HRG-score in low- and high-risk groups. (D) KM analysis of the OS between the low- and high-risk groups. (E) PCA demonstrated that the patients in the different risk groups were distributed in two directions. (F) ROC curves to predict the sensitivity and specificity of 1-, 3-, and 5-year survival according to the PRG-score.

FIGURE 5

Correlation and independent prognosis analysis of HRG-score and clinical parameters in the training set. (A,B) Univariate and multivariate analyses of the prognostic value of the HRG-score. (C) Relationships between clinical parameters and the low- and high-risk groups. (D) Clinical application value of HRG-score in predicting T stage, N stage, M stage, and TNM stage, respectively (*p < 0.05; **p < 0.01; ***p < 0.001).

FIGURE 6

Evaluation of the TME and checkpoints between the two risk groups. (A) Correlations between HRG-score and immune cell types. (B) Correlations between HRG-score and TME score. (C) Expression of immune checkpoint-related genes in the low- and high-risk groups. (D) Correlations between the relative abundance of immune cells and six genes involved in the HRG-score. (*p < 0.05; **p < 0.01; ***p < 0.001)

FIGURE 7

Comprehensive analysis of the HRG-score in COAD. (A) Relationships between the HRG-score and MMR status. (B) Relationships between the HRG-score and CSC index. (C) Relationships between HRG-score and sensitivity of five chemotherapeutics. (D) Prediction of the response of different risk samples to the combination of anti-CTLA4 and anti-PD1 based on IPS. (E) Boxplot for assessing HRG-score in predicting anti-PD-L1 response through the IMvigor210 cohort.