Identification and validation of a glycolysis-related gene signature for depicting clinical characteristics and its relationship with tumor immunity in patients with colon cancer

Clinical characteristics of CC patients

Overall, 566 CC and 19 normal control data stored on the Affymetrix Human Genome U133 Plus 2.0 Array platform were downloaded from a GEO dataset. The clinical characteristics of the samples, including sex, age, survival status, survival time, clinical TNM stage, clinical T stage, clinical N stage, and clinical M stage, are shown in Table 1.

Differentially expressed glycolysis genes between the CC and normal control samples

Six glycolysis-related gene sets were obtained from the Molecular Signatures Database: REACTOME_GLYCOLYSIS, WP_COMPUTATIONAL_MODEL_OF_AEROBIC_GLYCOLYSIS, WP_GLYCOLYSIS_IN_SENESCENCE, REACTOME_REGULATION_OF_GLYCOLYSIS_BY_FRUCTOSE_2_6_BISPHOSPHATE_METABOLISM, HALLMARK_GLYCOLYSIS, and GO_GLYCOLYTIC_PROCESS. A total of 301 glycolysis-related genes were selected from the six gene sets for additional research, and 217 glycolysis-related genes were found to be significantly and differentially expressed in the CC samples compared with the controls (Figure 1).

Figure 1. Significantly differentially expressed glycolysis-related genes between CC and normal controls.

Construction of the glycolysis-related gene prognostic prediction model for CC patients

UniCox and multiCox regression analyses were performed on the differentially expressed glycolysis-related genes to construct a novel biomarker model for predicting the prognosis of CC patients. Thirteen risk genes (NUP107, SEC13, ALDH7A1, ALG1, CHPF, FAM162A, FBP2, GALK1, IDH1, TGFA, VLDLR, XYLT2, and OGDHL) were identified (P < 0.05, Table 2). Then, the prognostic gene model based on the 13 glycolysis-related genes was used to divide patients into low- and high-risk groups as follows: Risk score = (-0.29×NUP107 Expression + [-0.35×SEC13 Expression] + 0.30×ALDH7A1 Expression + [-0.29×ALG1 Expression] + 0.45×CHPF Expression + [-0.46×FAM162A Expression] + 0.91×FBP2 Expression + [-0.65×GALK1 Expression] + [-0.34×IDH1 Expression] + 0.36×TGFA Expression + 0.18×VLDLR Expression + [-1.28×XYLT2 Expression] + [-0.41×OGDHL Expression]). The alterations in the 13 genes in the CC tissues were analyzed, and the results showed that the alteration rates in NUP107, SEC13, ALDH7A1, ALG1, CHPF, FAM162A, FBP2, GALK1, IDH1, TGFA, VLDLR, XYLT2, and OGDHL were 3%, 0.8%, 1.9%, 1%, 1.9%, 0.2%, 2.3%, 1.5%, 1.7%, 0.4%, 4%, 4%, and 4%, respectively (Figure 2A). Many mutations occurred in the gene domains, and the specific mutation sites are presented in Figure 2B. Finally, the expression of the 13 genes in the CC patient and normal tissues was further analyzed. NUP107, SEC13, ALDH7A1, ALG1, CHPF, GALK1, XYLT2, and OGDHL were highly expressed in CC tissues, while FAM162A, FBP2, IDH1, TGFA, and VLDLR were downregulated (P < 0.05, Figure 2C).

Figure 2. Mutations and expression levels of model genes in CC. (A) Mutation frequency of model genes. (B) Model gene-specific mutation domains. (C) Model gene expression between CC and normal control samples.

Efficacy of the risk score in CC patients

Based on the gene model, 280 CC patients and 281 CC patients were classified into high- and low-risk groups by the median risk score (Figure 3A). KM analysis was carried out between the low- and high-risk groups, and the high-risk group had a significantly poorer prognosis (P < 0.05, Figure 3B). ROC analysis was carried out, and the AUC was calculated to be 0.716, which showed a good prediction effect (Figure 3C). The expression of the 13 genes was calculated; TGFA FBP2, CHPF, and VLDLR were significantly higher in the high-risk group, and NUP107, OGDHL, SEC13, IDH1, GALK1, FAM162A, ALG1, and XYLT2 ALDH7A1 were downregulated (P < 0.05, Figure 3D). The risk plot indicates that the high-risk group was closely related to poor prognosis (Figure 3E). PCA showed that high- and low-risk groups could be significantly distinguished using our model, which reduced the dimension of the expression of multiple genes (Figure 3F). All the analysis results indicated that our gene model had good efficacy based on the risk score.

Figure 3. Text efficacy of the gene model in CC patients. (A) CC patient distribution according to the median risk score. (B) Kaplan–Meier analysis between the low- and high-risk groups. (C) ROC analysis of the risk score. (D) Model gene expression in the low- and high-risk groups. (E) The relationship between patient survival status and risk score. (F) PCA between the low- and high-risk groups.

Independent risk factor for risk score in CC patients and its relationship with clinical characteristics of the patients

Univariate and multivariate independent prognostic analyses were carried out to identify independent prognostic factors, including age, sex, TNM stage, and risk score. The results showed that age, TNM stage, and risk score were independent prognostic factors in CC patients, and they were all positively correlated with poor survival prognosis (P < 0.05, Figure 4A, 4B). For a more accurate analysis of the relationship between risk factors and patient OS, a series of Kaplan–Meier curve analyses was carried out. The results revealed that an age of > 65 years, TNM III–IV, T3-4, N1–3 and M1 were positively related to poor OS (P < 0.05, Figure 4C). The relationship between patients’ clinical characteristics and OS in the low- and high-risk groups was determined, and high-risk patients were positively associated with a poor OS when age ≤ 65 years, age>65 years, female sex, male sex, TNM I–II, TNM III–IV, T1–2, T3–4, M0, M1, N0, and N1-3 subgroups (P < 0.05, Figure 4D). These results demonstrate that our glycolysis-related gene model could distinguish clinical features and show the reliability of our glycolysis-related gene model in CC.

Figure 4. Analysis of risk factors. Univariate (A) and multivariate Cox regression analyses (B) of the relationship between risk core and other clinical features. (C) Kaplan–Meier analysis of clinical features, including age, sex, TNM stage, T stage, N stage, and M stage. (D) Kaplan–Meier analysis of clinical features in low- and high-risk subgroups.

Role of EMT in a high glycolysis tumor environment

Studies have demonstrated that glycolysis can promote EMT in cancers. To explore the role of EMT in CC, we analyzed the expression of EMT biomarkers (SNAI1, SNAI2, TWIST1, TWIST2, and ZEB2) in our study. SNAI1 and TWIST1 showed significantly higher expression in CC tissues (P < 0.05, Figure 5A). Moreover, SNAI1, SNAI2, TWIST1, TWIST2, and ZEB2 expression was upregulated in the high-risk group (P < 0.05, Figure 5B). This suggests that EMT plays a role in the high-glycolytic tumor environment of CC patients. The relationship between glycolysis and EMT in CC requires further research.

Figure 5. Expression of EMT biomarkers. (A) SNAI1 and TWIST1 expression in CC patients and normal controls. (B) SNAI1, SNAI2, TWIST1, TWIST2, and ZEB2 expression in the low- and high-risk groups.

Tumor immunity microenvironment between the low- and high-risk groups

We investigated the fractions of 22 types of infiltrated immune cells in every single cancer among all 561 CC samples (Figure 6A). The results showed that every sample had its own immune cell infiltration characteristics. The heatmap and vioplot of infiltrated immune cells were performed between the low- and high-risk groups according to the glycolysis-related gene model (Figure 6B, 6C). In low-risk group tissues, the top five most abundant fractions of infiltrating immune cells were M2 macrophages, CD8 T cells, plasma cells, follicular helper T cells, and resting memory CD4 T cells, and this was also true for the high-risk group (Figure 6B, 6C). However, no significant fractions of memory B cells, CD8 T cells, naïve CD T cells, resting CD4 memory T cells, activated CD4 memory T cells, follicular helper T cells, activated NK cells, monocytes, M0 macrophages, resting dendritic cells, activated mast cells, or eosinophils were observed in either the low- or high-risk group. (P > 0.05, Figure 6C). Naïve B cells, M1 macrophages, M2 macrophages, resting mast cells, and neutrophils were significantly higher in the high-risk group than in the low-risk group (p < 0.05, Figure 6C). Plasma cells, Tregs, gamma delta T cells, resting NK cells, and activated dendritic cells were higher in the low-risk group (p < 0.05, Figure 6C).

Figure 6. Immune infiltration between the low- and high-risk groups. (A) Fractions of immune cells in every single CC sample. (B, C) Heatmap and vioplot of immune cells between the low- and high-risk groups. Correlation between immune cells in low-risk samples (D), high-risk samples (E), and all colon cancer samples (F).

Finally, the correlation of immune cells in CC tissues in both the low- and high-risk groups was analyzed. In the low-risk group, a moderate to high negative correlation was observed between the T cells CD4 memory resting and T cells CD8, M1 macrophages and activated mast cells, T cells follicular helper and Tregs, activated mast cells and resting mast cells, and resting dendritic cells and M0 macrophages; however, a positive correlation was observed between the T cells follicular helper and M1 macrophages, activated T cells CD4 memory and M1 macrophages, T cells CD8 and activated NK cells, activated NK cells and resting mast cells, resting NK cells and activated mast cells, and M0 macrophages and neutrophils (Figure 6D). In the high-risk group, a moderate to high negative correlation was observed between the resting memory CD4 T cells and CD8 T cells, M1 macrophages and activated mast cells, follicular helper T cells and Tregs, resting mast cells and activated mast cells, and neutrophils and plasm cells. However, a naïve positive correlation was observed between neutrophils and M0 macrophages, resting NK cells and activated mast cells, follicular helper T cells and M1 macrophages, activated M1 macrophages and CD4 memory T cells, and activated mast cells and CD4 T cells (Figure 6E). A negative correlation was observed between resting memory CD4 T cells and CD8 T cells, M1 macrophages and activated mast cells, follicular helper T cells and Tregs, resting mast cells and activated mast cells, and resting dendritic cells and M0 macrophages. However, a positive relationship was observed in all CC patients between M1 macrophages and follicular helper T cells, M1 macrophages and activated memory CD4 T cells, activated mast cells and resting NK cells, CD8 T cells and activated NK cells, and activated NK cells and resting mast cells (Figure 6F). In the low-risk group, high-risk group, and all cancer samples, a negative correlation was observed between resting memory CD4 T cells and CD8 T cells, M1 macrophages and activated mast cells, follicular helper T cells and Tregs, and activated mast cells and resting mast cells; however, a positive correlation was observed between follicular helper T cells and M1 macrophages and activated memory CD4 T cells and M1 macrophages.

Validation of glycolysis-related gene prognostic prediction model in TCGA dataset

RNA-seq and corresponding clinical data of 437 colon cancer patients (Tumor=398, Normal=39) were downloaded from TCGA GDC (https://portal.gdc.cancer.gov/). Matching with clinical information, 189 patients and 190 patients were classified into high- and low-risk groups according to the above glycolysis-related gene prognostic prediction model, respectively (Figure 7A). The high-risk group had a significantly poorer prognosis than the low-risk group by KM analysis (P < 0.05, Figure 7B). Then, ROC analysis showed that the AUC of our model was 0.660 (Figure 7C). The risk plot showed that a higher risk score was related to poor prognosis (Figure 7D). Our model could significantly distinguish high- and low-risk group samples through PCA (Figure 7E).

Figure 7. Validation of the efficacy of the gene model in TCGA samples. (A) CC patient distribution according to the median risk score. (B) Kaplan–Meier analysis between the low- and high-risk groups. (C) ROC analysis of the risk score. (D) The relationship between patient survival status and risk score. (E) PCA between the low- and high-risk groups.

These results indicated that our gene model also had good efficacy in the TCGA validation dataset. Finally, the relationship between clinical characteristics and patient OS was examined. The results revealed that an age of > 65 years, M1, N1-3 and Stage III-IV were positively related to poor OS in the validation dataset (P < 0.05, Figure 8A). Then, the relationship between patients’ clinical features and OS in the low- and high-risk groups was analyzed, and the results indicated that high-risk patients were positively associated with poor OS when age ≤ 65 years, male sex, TNM stage I-II, T1-2 and N0 subgroups (P < 0.05, Figure 8B). These results showed that our genetic model could distinguish some clinical features and showed good efficacy. Although not all clinical features could be distinguished, this may be because TCGA colon cancer data were not as rich as GSE39582. It is also a vital reason that we used the GSE39582 dataset as a training cohort to construct the glycolysis-related gene prognostic prediction model.

Figure 8. Clinical features and risk factors in the validation cohort. (A) Kaplan–Meier analysis of clinical features, including age, sex, TNM stage, T stage, N stage, and M stage. (B) Kaplan–Meier analysis of clinical features in low- and high-risk subgroups.

Tumor immune sensitive drug prediction

These studies have confirmed that glycolysis is closely related to colon cancer immunity. To further study colon cancer-sensitive tumor immune drugs in the glycolytic microenvironment, the top 400 genes of the absolute value of significantly differential expression multiples between the high- and low-risk groups were screened (P < 0.05, Figure 9A). GO function analysis of the differentially expressed genes suggested that the functions of the differentially expressed genes were enriched in DNA replication, rRNA metabolic process, mitochondrial matrix, mitochondrial inner membrane, catalytic activity and acting on RNA and tubulin binding, which are related to cell metabolism (P < 0.05, Figure 9B). The KEGG results demonstrated that the differentially expressed genes were enriched in metabolism-related pathways such as carbon metabolism, the cell cycle, apoptosis, DNA replication and the proteome (P < 0.05, Figure 9C). We screened tumor immune-related drugs through Dreimt [7], and the top 400 genes were input into the software.

Figure 9. Tumor immune sensitive drug prediction. (A) Top 200 up- and downregulated significantly expressed genes between the high- and low-risk groups. (B, C) GO and KEGG enrichment analyses of differentially expressed genes. (D) Pie chart of tumor immune-sensitive drug types.

The results showed that the top 3 drug types were topoisomerase inhibitors, dopamine receptor antagonists and regenerative receptor antagonists (Figure 9D). The drugs were further screened (tau absolute value > 90, top 15 of drug specificity score > 0.7), and pirarubicin/mitoxantrone (topoisomerase inhibitor), daunorubicin (topoisomerase inhibitor, RNA synthesis inhibitor), BMS-387032/CGP-60474/JNJ-7706621/AT-7519/alvocidib/purvalanol-a/PHA-793887/AZD-5438 (CDK inhibitor), bisindolylmaleimide-ix (PKC inhibitor), dactinomycin RNA (polymerase inhibitor), fludarabine (ribonucleotide reductase inhibitor) and PF-562271 (focal adhesion kinase inhibitor) were filtered out. Our data provide a theoretical basis for tumor immune targeted therapy related to glycolysis in patients with colorectal cancer.

Validation of model gene expression in CC patient samples

To lay the foundation for the next step of our basic research, we validated the expression of all the model genes through western blotting. Our results showed that NUP107, SEC13, ALDH7A1, ALG1, CHPF, GALK1, XYLT2, and OGDHL were highly expressed, whereas FAM162A, FBP2, IDH1, TGFA, and VLDLR were expressed at low levels in CC patients (P < 0.05, Figure 10). These gene expression levels were consistent with the gene expression in GSE39582.

Figure 10. Model gene expression confirmed using western blotting.

Patient and control data

The Affymetrix Human Genome U133 Plus 2.0 array gene expression profile and corresponding clinical information data of 566 CC patients and 19 normal control tissues were downloaded from a GEO dataset (GSE39582). The clinical information included survival time, survival status, sex, age, tumor-node-metastasis (TNM) stage, clinical T stage, clinical N stage, and clinical M stage.

Identification of differentially expressed glycolysis-related genes

A total of 301 mRNAs of six glycolysis-related gene sets were selected from the GSEA database (WP_COMPUTATIONAL_MODEL_OF_AEROBIC_GLYCOLYSIS, WP_GLYCOLYSIS_IN_SENESCENCE, REACTOME_REGULATION_OF_GLYCOLYSIS_BY_FRUCTOSE_2_6_BISPHOSPHATE_METABOLISM, REACTOME_GLYCOLYSIS, HALLMARK_GLYCOLYSIS, and GO_GLYCOLYTIC_PROCESS). The Wilcoxon test was used to determine the differentially expressed genes. Then, a prognostic gene model was constructed based on the 217 differentially expressed glycolysis-related genes. The cutoffs were set as Log2(fold change) > 1 and p value< 0.05. Heatmaps were generated using the Pheatmap package.

Construction of a glycolysis-related gene risk model

Based on the 217 differentially expressed glycolysis-related genes in GSE39582, uniCox and multiCox analyses were conducted to identify the model genes. The risk score was obtained from the results of multiplication of the Cox coefficient of multiCox analysis and gene expression values. The following formula was used for the calculation of the risk score: Risk score = (Cox coefficient of model gene A* A gene expression + Cox coefficient of model gene B* B gene expression+……). Then, all patients were divided into low- and high-risk subgroups based on the median risk score. Cbioportal (http://www.cbioportal.org/) was used to screen the model gene mutation status and gene mutation sites. Survival and survminer packages were used for survival analysis. Receiver operating characteristic (ROC) analysis was performed using the survival ROC package. An area under the curve value of > 0.65 was regarded as an acceptable cutoff value. The risk score, principal component analysis (PCA), and survival status plots were generated to establish a prediction model between the low- and high-risk groups.

Independent prognostic value of risk score

To assess the prognostic value of the model, we used univariate and multivariate Cox regression analyses to determine the prognostic factors. The factors (age, sex, TNM stage, and risk score) with p values less than 0.05 in univariate and multivariate Cox regression analyses were considered as independent factors.

Clinical application of the gene model

To examine the efficiency of our model in CC, we analyzed the correlation between the gene model and clinical features (age, sex, clinical TNM stage, clinical T stage, clinical N stage, and clinical M stage). All patients were classified into two subgroups based on age (> 65 and ≤ 65 years old), sex (female and male), TNM stage (TNM stage I-II and TNM stage III-IV), clinical T stage (T1-2 and T3-4), clinical M stage (M0 and M1), and clinical N stage (N0 and N1-3). Survival analyses were performed between the two groups.

Immune cell infiltration between low- and high-risk patients

Subsequently, CIBERSORT, including 22 immune cells (seven kinds of T cells, naïve B cells, memory B cells, plasma cells, resting NK cells, activated NK cells, monocytes, M0–M2 macrophages, resting dendritic cells, activated dendritic cells, resting mast cells, activated mast cells, eosinophils, and neutrophils), was used for immune cell infiltration evaluation.

Functional enrichment analyses between low- and high-risk groups

Gene ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were conducted using the org.Hs.e.g..db R package and clusterProfiler R package. A P value of < 0.05 was identified as statistically significant.

Model validation using western blot analysis

CC patient and control tissues were collected from the Affiliated Kunshan Hospital of Jiangsu University. The patient tissues were preserved in liquid nitrogen and lysed using RIPA lysis buffer (Strong) (MedChemExpress, China). A BCA Protein Quantification Kit (Vazyme, China) was used to measure the protein concentration. Proteins were transferred to nitrocellulose membranes (Millipore, Billerica, MA, USA) after separation through 10% sodium dodecyl sulfate–polyacrylamide gel electrophoresis. Then, the membranes were incubated with primary antibodies against the following: NUP107/SEC13/ALDH7A1/ALG1/CHPF/FAM162A/FBP2/GALK1/IDH1/TGFA/VLDLR/XYLT2/OGDHL (Abcam, USA, 1:1000) and GAPDH (Abcam, USA, 1:3000). The membranes were washed and incubated with horseradish peroxidase-labeled goat anti-rabbit IgG (1:5000; Sigma). The signals were visualized using BIO-RADXR.

Statistical analysis

Statistical analysis was performed using SPSS 21.0 software (SPSS, Chicago, IL, USA). Data were presented as the means ± standard deviation (SD). The difference between the groups was determined using ANOVA with repeated measures. To compare the difference between the two groups, the independent sample t test was used. P<0.05 was considered as significant.