A novel prognostic model for colon adenocarcinoma based on cofactor and vitamin metabolism-related genes

Introduction

Colon adenocarcinoma (COAD) is a malignant neoplasm that predominantly arises in the proximal colon. Epidemiological data indicate that colorectal cancer (CRC) ranks third in global cancer incidence and second in cancer-related mortality, with COAD constituting a substantial proportion of cases (1). As a prevalent malignancy worldwide, COAD remains a significant clinical challenge due to the absence of effective early diagnostic strategies and reliable biomarkers, resulting in most patients being diagnosed at an advanced stage, which critically impacts prognosis (2,3). Current prognostic assessment of COAD primarily relies on clinical parameters such as the tumor-node-metastasis (TNM) staging system, pathological characteristics, and molecular subtyping. However, these conventional methods exhibit limitations in accurately predicting therapeutic response and long-term survival, necessitating the identification of novel biomarkers or predictive models for refined prognostic stratification (4). Metabolic reprogramming is a hallmark of cancer, and aberrant metabolic pathway alterations contribute to tumor progression, immune evasion, and therapeutic resistance (5,6). A deeper understanding of COAD-associated metabolic dysregulation may provide new insights into potential prognostic indicators and therapeutic targets.

Cofactors and vitamins are indispensable components of enzymatic reactions, playing crucial roles in cellular metabolism, redox homeostasis, and epigenetic regulation, thereby ensuring cellular stability (7,8). The availability and regeneration efficiency of cofactors directly influence the metabolic state of cancer cells. For instance, tumor cells regenerate nicotinamide adenine dinucleotide (NAD⁺) through lactate fermentation, a process known as the Warburg effect, to maintain a high NAD⁺/NADH ratio, which is essential for sustaining enhanced glycolytic metabolism(9). In addition, cancer cells adapt their metabolic pathways or upregulate the uptake of specific cofactors and vitamins through oncogenic signaling. For example, the oncogene MYC facilitates the acquisition of NAD⁺ and pantothenic acid, a precursor of coenzyme A (CoA), to promote metabolic reprogramming and support tumor progression (10,11). Beyond their metabolic functions, cofactors and vitamin-derived metabolites exert profound regulatory effects on epigenetic modifications. Vitamin C and flavin adenine dinucleotide (FAD), a coenzyme derived from riboflavin, serve as critical cofactors or modulators for DNA and histone demethylases, including ten-eleven translocation (TET) enzymes and lysine-specific demethylase 1 (LSD1), thereby influencing chromatin remodeling and gene expression. Furthermore, certain metabolic byproducts that accumulate during tumor progression, such as 2-hydroxyglutarate, succinate, and fumarate, act as competitive inhibitors of epigenetic enzymes, exacerbating aberrant epigenetic modifications and disrupting gene expression (12,13). Moreover, vitamin C-mediated post-translational modifications have been implicated in the regulation of tumor immune signaling pathways, thereby enhancing anti-tumor immune responses (14).

The regulation of cofactor and vitamin metabolism is not only fundamental for maintaining cellular homeostasis but also plays a pivotal role in tumor initiation, progression, and therapeutic response. A deeper understanding of these metabolic processes may provide novel insights into potential molecular targets and therapeutic strategies for cancer prevention and treatment. Existing studies have demonstrated that cofactor and vitamin metabolism play a critical role in the initiation, progression, and treatment response of colon cancer (15,16). Although substantial progress has been made in elucidating the underlying mechanisms, research on the prognostic implications of genes involved in these metabolic pathways remains limited. Given that conventional clinicopathological indicators often fail to fully capture individual patient prognosis and therapeutic response, developing a risk scoring model based on cofactor and vitamin metabolism-related genes (CVMRGs) may provide novel strategies and valuable insights for prognostic assessment and personalized treatment in colon cancer.

This study aims to construct a novel risk prediction model for prognostic evaluation in colon cancer patients based on genes associated with cofactor and vitamin metabolism and to validate the model using an independent cohort. The study systematically evaluates the predictive performance of the model in colon cancer prognosis and explores its potential association with the tumor immune microenvironment. Furthermore, the study preliminarily investigates the clinical relevance of the model in guiding therapeutic decision-making for colon cancer patients, with the goal of providing a scientific basis for accurate prognostic stratification and individualized treatment strategies. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1521/rc).

Methods

Data collection and processing

The study size was based on the availability of relevant public datasets. Transcriptomic profiles and clinical data were acquired from two public repositories: The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov/) provided RNA-sequencing (RNA-seq) data [raw counts and transcripts per kilobase million (TPM) values] with clinical annotations for 502 colon specimens (461 primary tumor samples and 41 normal controls) (17), while the Gene Expression Omnibus (GEO; sample series: GSE39582) contributed microarray data from 566 colon cancer patients (18). The transcriptomic and clinical data used in this study are derived from public databases. Due to the nature of such data, certain background information relevant to the study has not been detailed in the publicly available datasets. Predictors and outcomes were independently obtained from public datasets. No manual assessment was involved, and outcome labels of the test set were concealed during model development to avoid information leakage. After excluding 7 TCGA cases and 4 GEO samples with incomplete follow-up, the final cohorts comprised 454 tumor samples (TCGA discovery set) and 562 patients (GEO validation set). Raw sequencing counts were normalized via DESeq2 variance-stabilizing transformation, and microarray data underwent robust multi-array average (RMA) normalization with ComBat-based batch correction (19). This secondary analysis of pre-existing anonymized data from TCGA-COAD and GEO (GSE39582) databases complies with international ethics guidelines, exempting the requirement for institutional review board approval. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Identification of differentially expressed CVMRGs in TCGA

CVMRGs were initially identified through KEGG pathway annotations, yielding 244 candidate genes (20). Low-expressed genes (transcripts per million <1 in >50% of samples) were excluded, retaining 214 CVMRGs for cross-dataset analysis. Differential expression between 461 COAD tumors and 41 matched normal tissues was assessed using DESeq2 (R v4.4.3) with a lenient threshold (|log₂(fold change)| >0, false discovery rate (FDR)-adjusted P<0.05) to maximize candidate discovery.

Construction of prognostic metabolic signature

TCGA cohort was utilized as the discovery cohort for signature construction. Prognostic metabolic genes were initially screened through univariate Cox proportional hazards regression analysis using R packages “survival” and “survminer”, with overall survival (OS) as the primary endpoint. Genes showing statistically significant associations with OS (P<0.05) were selected for subsequent modeling. To establish a robust prognostic model, least absolute shrinkage and selection operator (LASSO) Cox regression with 10-fold cross-validation was implemented through 1,000 bootstrap resamples using the “glmnet” package. This approach ensured feature selection stability while controlling for overfitting through regularization parameter (λ) optimization.

Validation of prognostic signature performance

The metabolic risk model was applied in the TCGA cohort through stratified risk assessment. Patients were stratified into high- and low-risk groups using the median risk score as the predefined cutoff value. Kaplan-Meier survival curves with log-rank testing (α=0.05) were generated through the survival package, while risk score distribution and expression pattern visualization were implemented using the “pheatmap” package. Univariate and multivariate Cox proportional hazards regression analyses were systematically performed to evaluate independent prognostic factors, incorporating clinical variables (age, sex, TNM stage). Time-dependent receiver operating characteristic (ROC) analysis at 1-, 3-, and 5-year intervals was conducted using the “timeROC” and “pROC” packages, with area under the curve (AUC) values quantifying model discrimination capacity (21). A clinical nomogram integrating the risk score with established prognosticators (age, gender, TNM stage) was developed via the “rms” package. Survival probability calibration was assessed through bootstrap resampling (n=1,000 iterations). Model calibration was assessed using calibration plots comparing predicted probabilities with observed outcomes. No model updating or recalibration was performed after external validation; the model was evaluated in its original form as developed in the training dataset. The external validation dataset, sourced from the GEO database, was consistent with the development dataset in terms of study setting, eligibility criteria, outcome definitions, and predictor measurements, with no major differences identified.

Functional enrichment analysis in TCGA cohort

Differential gene expression analysis between high-risk and low-risk subgroups was performed using DESeq2 with default parameters. Significantly differentially expressed genes (DEGs) were defined by absolute log₂(fold change) >1.0 and Benjamini-Hochberg adjusted P<0.05. Functional annotation of DEGs was conducted through Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis using “clusterProfiler” (22,23). Gene set enrichment analysis (GSEA) implemented in “clusterProfiler” identified the top three significantly enriched pathways demonstrating positive/negative association with high-risk status (24). Pathway activity quantification was performed using gene set variation analysis (GSVA) with 20 HALLMARK gene sets from molecular signatures database (MSigDB), followed by differential pathway activity assessment through linear modeling with “limma” R package (25,26).

Construction of protein-protein interaction (PPI) network

PPI network was constructed by Search Tool for the Retrieval of Interacting Genes (STRING, https://string-db.org/) and the CytoHubba plug-in the Cytoscape software (version 3.10.1) (27).

Tumor immune microenvironment

The tumor immune microenvironment was comprehensively profiled through multi-modal computational approaches. Tumor immune evasion potential was quantified using the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu/), which evaluates immune checkpoint inhibitor (ICI) response likelihood through integrated modeling of T cell dysfunction and exclusion signatures. Immune cell composition was deconvoluted using Cell-type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT)with LM22 signature matrix, enabling quantification of 22 immune cell subtypes based on lineage-specific gene expression patterns. Tumor purity and stromal-immune infiltration were simultaneously assessed through Estimation of Stromal and Immune Cells in Malignant Tumors using Expression data (ESTIMATE) scoring, which calculates immune/stromal/tumor scores from transcriptomic profiles (28). Both analyses were implemented through the IOBR package with default parameters, ensuring methodological consistency across microenvironmental features (29).

Anticancer drug sensitivity analysis

Drug sensitivity prediction was conducted using “oncoPredict”, an R package implementing machine learning models trained on the Cancer Therapeutics Response Portal (CTRP) and Genomics of Drug Sensitivity in Cancer (GDSC) datasets (30). The algorithm estimates half-maximal inhibitory concentration (IC50) values by integrating tumor transcriptomic profiles with pharmacogenomic data from >1,000 cancer cell lines, utilizing ridge regression for robust prediction of drug response phenotypes.

Statistical analysis

The analysis based on the R package in this study was completed using R version 4.4.3. The log-rank test was used to compare OS rates between the TCGA and GEO cohorts. Two groups of data may be compared using the Wilcoxon test. P<0.05 is considered statistically significant.

Results

Flow chart of overall design

Transcriptomic expression data from COAD patients were obtained from TCGA, comprising a cohort of 495 patients, including tumor samples (n=454) and normal samples (n=41) (Table 1). Genes associated with cofactor and vitamin metabolism pathways (CVMRGs, n=214) were identified from the KEGG pathway database. Differential expression analysis of CVMRGs between tumor and normal tissues in the TCGA dataset revealed 165 differentially expressed CVMRGs (69 upregulated and 96 downregulated). Univariate Cox regression analysis identified 10 genes significantly correlated with OS in TCGA-COAD patients (P<0.05). Subsequent LASSO-Cox regression analysis with 1,000 iterations further refined these to 6 key genes, which were used to construct a prognostic signature. Finally, pathway enrichment analysis, immune microenvironment profiling, survival validation (GSE39582 dataset) (Table 2), and anticancer drug sensitivity analysis were performed to investigate the prognostic value of the signature in colon cancer and its associations with immune response and therapeutic efficacy (Figure 1).

Figure 1 Flow chart of overall design. COAD, colon adenocarcinoma; CVMRG, cofactor and vitamin metabolism-related gene; LASSO, least absolute shrinkage and selection operator; OS, overall survival; RNA-seq, RNA-sequencing; TCGA, The Cancer Genome Atlas.

Identification of differentially expressed CVMRGs in TCGA dataset

CVMRGs were systematically retrieved from the KEGG database. Following the hierarchical organization of KEGG, the screening process entailed accessing the KEGG PATHWAY module, navigating to the metabolism category, and further selecting the “Cofactor and Vitamin Metabolism” subcategory. Through this standardized procedure, a total of 244 CVMRGs were ultimately identified. Following intersection of CVMRGs across TCGA and GEO databases, 214 conserved metabolic genes were retained. Differential expression analysis between COAD and normal colon tissues in the TCGA cohort identified 165 significantly dysregulated CVMRGs (thresholds: |log₂ fold-change| >0 with Benjamini-Hochberg adjusted P<0.05, to retain more potential prognosis-related genes), comprising 69 upregulated and 96 downregulated genes (Figure 2A,2B).

Figure 2 Identification of prognosis-related CVMRGs by DEG analysis, univariate Cox analysis. (A) Volcano plot of DEGs in TCGA-COAD (P<0.05). Red: upregulated genes; blue: downregulated genes; gray: no-significant difference genes. (B) Venn diagram of CVMRGs and DEGs. (C) Forest plot displays univariate Cox-derived hazard ratios (95% CI) for CVMRGs (P<0.05). (D) Heatmap of 10 differentially expressed CVMRGs. Turquoise: normal; red: tumor. (E) PPI network of 10 CVMRGs constructed via STRING database. AK7, adenylate kinase 7; ALAD, aminolevulinate dehydratase; ALDH2, aldehyde dehydrogenase 2; ALPP, alkaline phosphatase; CI, confidence interval; COAD, colon adenocarcinoma; CSAD, cysteine sulfinic acid decarboxylase; CVMRG, cofactor and vitamin metabolism-related gene; CYP26A1, cytochrome P450 family 26 subfamily A member 1; DEG, differential expression gene; DLAT, dihydrolipoamide acetyltransferase; HR, hazard ratio; LRAT, lecithin-retinol acyltransferase; PPI, protein-protein interaction; STRING, Search Tool for the Retrieval of Interacting Genes; TCGA, The Cancer Genome Atlas; TH, tyrosine hydroxylase; TPK1, thiamine pyrophosphokinase 1.

Construction of prognostic risk signature based on 6 CVMRGs

Univariate Cox regression analysis in the TCGA-COAD cohort identified 10 CVMRGs significantly associated with OS (Figure 2C). Boxplots and heatmaps illustrated differential expression profiles of these 10 prognosis-linked CVMRGs between COAD tumors and normal tissues, including: Favorable prognosis genes: adenylate kinase 7 (AK7), aldehyde dehydrogenase 2 (ALDH2), dihydrolipoamide acetyltransferase (DLAT), thiamine pyrophosphokinase 1 (TPK1). Unfavorable prognosis genes: aminolevulinate dehydratase (ALAD), alkaline phosphatase, placental (ALPP), cysteine sulfinic acid decarboxylase (CSAD), cytochrome P450 family 26 subfamily A member 1 (CYP26A1), lecithin-retinol acyltransferase (LRAT), tyrosine hydroxylase (TH) (Figure 2D). PPI network analysis revealed DLAT and ALDH2 as pivotal hub genes, exhibiting extensive connectivity with other CVMRGs, suggesting central roles in metabolic regulation and neurotransmitter synthesis (Figure 2E). A prognostic risk model was subsequently developed through LASSO-penalized Cox regression, with 1,000 bootstrap resamples and 10-fold cross-validation, retaining 6 CVMRGs with non-zero coefficients. The risk score was calculated as: Exp ALDH2 × (−0.06739587) + Exp DLAT × (−0.02579058) + Exp TPK1 × (−0.01541147) + Exp ALAD × 0.07035241 + Exp ALPP × 0.04101327 + Exp TH × 0.04261918 (Figure 3A). Higher risk scores indicate a greater predicted risk of the outcome. Patients can be stratified into high- and low-risk groups based on the median risk score. Kaplan-Meier analysis of the TCGA cohort demonstrated that patients in the low-risk group had significantly better OS compared to those in the high-risk group (P<0.001) (Figure 3B).

Figure 3 Prognostic values of the six CVMRGs signature model in TCGA-COAD cohort. (A) The distribution of risk score, survival outcomes, and the expression pattern of 6 CVMRGs between high- and low-risk score group in the TCGA-COAD cohort. (B) Kaplan-Meier survival curve compares low- vs. high-risk groups (P=0.0005). (C) Violin-box plot depicts risk score distribution across TCGA-COAD CMS subtypes. (D) Heatmap illustrates associations between 6 CVMRGs expression profiles and clinical features. ALAD, aminolevulinate dehydratase; ALDH2, aldehyde dehydrogenase 2; ALPP, alkaline phosphatase; CMS, consensus molecular subtypes; COAD, colon adenocarcinoma; CVMRG, cofactor and vitamin metabolism-related gene; DLAT, dihydrolipoamide acetyltransferase; M, metastasis; N, node; TCGA, The Cancer Genome Atlas; T, tumor; TH, tyrosine hydroxylase; TPK1, thiamine pyrophosphokinase 1.

Risk signature correlates with consensus molecular subtypes (CMS) of TCGA-COAD

The molecular subtyping of CRC is based on the molecular characteristics of the tumor, mainly including chromosome instability, microsatellite instability, CpG island methylation, and KRAS gene mutations. Based on the synthesis of six independent classification systems, the International Colorectal Cancer Classification Consortium (CRCSC) proposed the CMS consensus subtyping, which includes CMS1 (MSI immune type), CMS2 (canonical type), CMS3 (metabolic type), and CMS4 (mesenchymal type) subtypes (31). To explore the biological relevance of the risk signature, we assigned TCGA-COAD cohort to CMS using RNA-seq data and “CMScaller” package. Notably, the risk score was significantly higher in CMS4 (the worst prognosis) than in CMS1, CMS2, and CMS3 (P=0.0001, 0.0003, and 4.8e−08, respectively) (Figure 3C). This indicates our model captures features linked to aggressive subtypes, supporting its biological relevance.

Relationship between clinicopathological characteristics and risk score

Using clinical data from the TCGA-COAD cohort, we investigated the association between risk score and prognostic factors (Figure 3D). The results revealed significant correlations between higher risk scores and advanced T stage (P=6.2e−07), N stage (P=1.2e−09), M stage (P=4.3e−06), tumor stage (P=8.4e−10), and postoperative tumor status (P=0.006). In contrast, no significant associations were observed with other clinical features, including gender (P=0.74), age (P=0.15), or histological type (P=0.64) (Figure 4A-4H).

Figure 4 Boxplots of risk scores across clinical variables. The clinical variables included T stage (A), N stage (B), M stage (C), overall tumor stage (D), gender (E), age (F), postoperative tumor status (G), and histological type (H). M, metastasis; N, node; T, tumor.

Independent prognostic factors and nomogram model construction

Based on the median risk score, we classified TCGA-COAD samples into high-risk and low-risk subgroups. Furthermore, univariate and multivariate analyses identified age, N stage, M stage, and risk score as independent prognostic factors (P<0.001) (Figure 5A,5B). A nomogram was constructed incorporating age, gender, T stage, N stage, M stage, and risk score (Figure 5C). The AUC for 1-, 3-, and 5-year OS predictions were 0.776, 0.771, and 0.759, respectively (Figure 5D) (32). The calibration curves further validated the superior predictive performance of the nomogram model (Figure 5E).

Figure 5 The results of prognostic factor analysis and the construction of a predictive nomogram model. (A,B) Univariate and multivariate Cox regression analyses of clinical variables (gender, age, T/N/M stage, risk score), presenting HRs (95% CIs) and P values. (C) Nomogram incorporating clinical factors (gender, T/N/M stage, age, risk score) to predict 1-, 3-, and 5-year survival probabilities. 0, female; 1, male. (D) ROC curves for 1-, 3-, and 5-year OS, with AUC values indicating the model’s predictive accuracy. (E) Calibration curves for survival probability assessment. AUC, area under the curve; CI, confidence interval; HR, hazard ratio; M, metastasis; N, node; OS, overall survival; ROC, receiver operating characteristic; T, tumor.

Functional enrichment analysis of the CVMRG-related signature in the TCGA cohort

A differential gene expression analysis was conducted on high-risk and low-risk subgroups of TCGA-COAD samples, identifying a total of 1,252 DEGs, of which 870 were upregulated and 382 were downregulated [threshold: |log₂(fold change)| >1, P<0.05]. GO annotation revealed enrichment in 329 GO terms, comprising 221 biological processes (BPs), 61 cellular components (CCs), and 47 molecular functions (MFs). The GO enrichment analysis highlighted key pathways, including modulation of chemical synaptic transmission, synapse organization, channel activity, and postsynaptic membrane function (Figure 6A). Additionally, KEGG pathway analysis identified enrichment in 227 pathways, with significant involvement in Cytoskeleton in muscle cells, Calcium signaling pathway, and Neuroactive ligand-receptor interaction (Figure 6B). Moreover, GSVA indicated that the high-risk group was significantly enriched in KRAS signaling DN, Hedgehog signaling, epithelial-mesenchymal transition (EMT), and Wnt beta catenin signaling, whereas the low-risk group exhibited enrichment in PI3K-AKT-mTOR signaling, G2M checkpoint, reactive oxygen species pathway, and oxidative phosphorylation (Figure 6C). GSEA further demonstrated significant enrichment of these gene sets across multiple biological pathways (Figure 7A-7F).

Figure 6 Functional and pathway enrichment analysis of DEGs. (A) GO enrichment analysis highlighting key biological processes and molecular functions. (B) KEGG pathway analysis showing significantly enriched pathways. (C) The heatmap of hallmark pathway activation states between low- and high-risk groups using GSVA enrichment analysis. BP, biological process; CC, cellular component; DEG, differentially expressed gene; GSVA, gene set variation analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function.

Figure 7 GSEA of key biological processes. The analysis highlights pathways related to modulation of chemical synaptic transmission (NES =1.99, adjusted P<0.001), synaptic membrane (NES =1.92, adjusted P<0.001), channel activity (NES =1.75, adjusted P<0.001), mitochondrial translation (NES =−2.41, adjusted P<0.001), DNA packaging complex (NES =−2.42, adjusted P<0.001), and catalytic activity, acting on a trim (NES =−1.91, adjusted P<0.001), all significantly enriched. GSEA, gene set enrichment analysis; NES, normalized enrichment score.

Comprehensive immune landscape analysis of prognostic characteristics of CVMRG risk score

TIDE is a computational tool used to evaluate tumor immune evasion mechanisms and predict responses to ICIs. By integrating multiple biomarkers, TIDE assesses tumor-immune interactions to estimate the potential responsiveness of tumors to immunotherapy. A higher TIDE score is associated with poorer response to immune checkpoint blockade therapy. Significant differences were observed between the high- and low-risk groups in TIDE, microsatellite instability expression signature (MSI Expr Sig), dysfunction, and cancer-associated fibroblasts (CAFs) (Figure 8A). The high-risk group exhibited elevated immune evasion or functional impairment, suggesting that tumors in this group may escape treatment by suppressing immune responses. In contrast, the low-risk group showed a higher proportion of myeloid-derived suppressor cells (MDSCs), indicating a stronger immunosuppressive characteristic in this group. CAFs, MDSCs, and tumor-associated macrophages (M2-type TAMs) act synergistically within the tumor microenvironment through multiple mechanisms, forming a potent immunosuppressive barrier that restricts T-cell infiltration and function.

Figure 8 Immune landscape analysis of high- and low-risk groups. (A) Comparison of immune-related biomarkers between high- and low-risk groups, including TIDE, IFNG, MSI Expr Sig, Merck18, CD274, CD8, dysfunction, exclusion, MDSC, CAF, and TAM M2. (B-D) ESTIMATE analysis of total, immune, and stromal scores between high- and low-risk groups. (E) Immune cell infiltration differences between high- and low-risk groups. *, P<0.05; **, P<0.01; ***, P<0.001; ns, P>0.05. CAF, cancer-associated fibroblasts; IFNG, interferon gamma; MDSC, myeloid-derived suppressor cell; MSI Expr Sig, microsatellite instability expression signature; TAM M2, tumor-associated macrophages M2; TIDE, Tumor Immune Dysfunction and Exclusion.

ESTIMATE is a computational method used to evaluate the composition of the tumor microenvironment by estimating the abundance of stromal and immune cells based on gene expression data. The analysis revealed that the stromal score (P=0.0002) and ESTIMATE score (P=0.03) were significantly higher in the high-risk group, indicating a more pronounced stromal component in these tumors (Figure 8B-8D). This suggests that the tumor microenvironment in the high-risk group may be characterized by a denser extracellular matrix, increased fibroblast activity, or enhanced interactions between stromal and malignant cells, potentially influencing tumor progression and therapeutic response.

CIBERSORT is a computational method that quantifies immune cell infiltration based on gene expression data by estimating the relative proportions of different immune cell types within tumor samples. In the TCGA-COAD cohort, significant differences in immune cell composition were observed between the high- and low-risk groups. Comparative analysis revealed that CD8⁺ T cells, CD4⁺ memory resting T cells, CD4⁺ memory activated T cells, gamma delta T cells, monocytes, M0 macrophages, activated dendritic cells, resting mast cells, activated mast cells, and neutrophils were differentially enriched across risk groups (Figure 8E). Specifically, the high-risk group exhibited greater enrichment of CD8⁺ T cells, monocytes, M0 macrophages, and resting mast cells, whereas the low-risk group had a higher abundance of CD4⁺ memory resting T cells, CD4⁺ memory activated T cells, gamma delta T cells, activated dendritic cells, activated mast cells, and neutrophils. These findings highlight distinct immune cell infiltration patterns between the two groups, suggesting potential differences in immune responses and tumor-immune interactions that may influence disease progression and therapeutic outcomes.

Drug sensitivity prediction analysis

Based on TCGA-COAD data, the “OncoPredict” package was used to estimate the IC50 values of multiple chemotherapeutic agents. The box plots illustrate significant differences in drug response between the high- and low-risk groups. For paclitaxel, fluorouracil, and gemcitabine, the low-risk group exhibited significantly lower IC50 values (P=1.7e−05, P=9.8e−10, and P=8.6e−07, respectively), suggesting higher sensitivity to these drugs (Figure 9A-9C). Conversely, for regorafenib, the high-risk group displayed a significantly lower IC50 value (P=0.007; Figure 9E), indicating greater sensitivity to this agent. In contrast, no significant difference in IC50 values was observed for oxaliplatin (P=0.51) and dabrafenib (P=0.55) between the two groups, suggesting a comparable drug response (Figure 9D,9F). These findings provide insights into potential personalized treatment strategies, indicating that different risk groups may exhibit distinct sensitivities to specific chemotherapeutic agents.

Figure 9 Anticancer drug sensitivity analysis between two groups. The IC50 of 6 anticancer drugs: paclitaxel (A), fluorouracil (B), gemcitabine (C), oxaliplatin (D), regorafenib (E), and dabrafenib (F). IC50, half-maximal inhibitory concentration.

Validation of the prognostic risk model

To evaluate the robustness of the constructed prognostic risk model, we applied it to an independent validation cohort (GSE39582). Kaplan-Meier survival analysis demonstrated that patients in the low-risk group had significantly better OS compared to those in the high-risk group (P=0.02), further supporting the predictive value of the risk score (Figure 10A). Time-dependent ROC curve analysis at 1, 3, and 5 years yielded AUC values of 0.672, 0.659, and 0.675, respectively, indicating a favorable predictive performance of the model (Figure 10B). Univariate Cox regression analysis in the validation cohort revealed that age, N stage, M stage, and risk score were significantly associated with prognosis (P<0.05). Furthermore, multivariate Cox regression analysis confirmed that the risk score remained an independent prognostic factor, even after adjusting for other clinical variables (Figure 10C,10D). Collectively, these results validate the reliability and clinical applicability of our prognostic risk model, reinforcing its potential utility for risk stratification and personalized prognosis prediction in colon cancer patients.

Figure 10 Validation of the prognostic model in the GSE39582 cohort. (A) Kaplan-Meier survival curve indicating poorer survival in the high-risk group compared to the low-risk group. (B) ROC curves assessing the predictive performance of the risk score at different time points. (C) Heatmap of gene expression across risk groups with clinical characteristics. Red and blue colors represent high and low expression. (D) Univariate and multivariate Cox regression analyses of clinical factors and risk score. 0, female; 1, male. ALAD, aminolevulinate dehydratase; ALDH2, aldehyde dehydrogenase 2; ALPP, alkaline phosphatase; AUC, area under the curve; CI, confidence interval; DLAT, dihydrolipoamide acetyltransferase; HR, hazard ratio; M, metastasis; N, node; ROC, receiver operating characteristic; T, tumor; TH, tyrosine hydroxylase; TPK1, thiamine pyrophosphokinase.

Discussion

COAD remains a clinically challenging malignancy due to its heterogeneous prognosis and limited biomarkers for personalized treatment (33). In this study, we developed a novel prognostic model based on CVMRGs, which demonstrated robust predictive accuracy for survival and therapeutic response. To our knowledge, this is the first prognostic signature integrating metabolic pathways involving enzymatic cofactors and vitamins, offering unique insights into the metabolic-immune interplay in COAD progression.

The 6-gene risk model (DLAT, TH, AK7, ALDH2, ALAD, CYP26A1) exhibited strong prognostic value, with high-risk patients showing significantly shorter OS and advanced TNM stages. Notably, DLAT and ALDH2 emerged as hub genes in the protein interaction network, aligning with their established roles in metabolic regulation. DLAT, a key component of the pyruvate dehydrogenase complex, facilitates acetyl-CoA production, which fuels lipid synthesis and epigenetic modifications in tumors (34,35). Similarly, TH regulates catecholamine synthesis, potentially influencing tumor angiogenesis and immune evasion through neurotransmitter-mediated signaling (36). These findings underscore the biological plausibility of CVMRGs as prognostic markers, linking metabolic dysregulation to aggressive tumor behavior.

Functional enrichment analysis revealed distinct pathway activation patterns between risk groups. High-risk tumors were enriched in EMT and Wnt/β-catenin signaling, both hallmarks of metastatic progression (37,38). This aligns with their immunosuppressive microenvironment characterized by elevated CAFs and TIDE scores, suggesting metabolic reprogramming may drive immune exclusion. Conversely, low-risk tumors showed activation of oxidative phosphorylation and reactive oxygen species pathways, indicative of a less aggressive metabolic phenotype. These observations resonate with recent studies highlighting mitochondrial metabolism as a determinant of tumor immunogenicity (39), though our model uniquely extends this paradigm to cofactor and vitamin dependencies.

The differential drug sensitivity between risk groups further supports the clinical utility of this model. High-risk patients exhibited greater sensitivity to fluorouracil and gemcitabine, possibly due to heightened dependence on nucleotide synthesis pathways targeted by these agents (40). In contrast, low-risk patients responded better to regorafenib, a multi-kinase inhibitor, potentially reflecting their reliance on kinase-driven metabolic adaptations (41). This heterogeneity underscores the need for risk-stratified chemotherapy regimens, a strategy increasingly emphasized in precision oncology (42).

Despite these advances, several limitations warrant consideration. First, retrospective design and reliance on public datasets may introduce selection bias. Second, while the model was validated in an independent cohort, prospective clinical trials are needed to confirm its utility in guiding treatment decisions. Third, the exact mechanisms linking specific CVMRGs such as CYP26A1 (involved in retinoic acid metabolism) to immune modulation remain speculative and require experimental validation (43). Future studies integrating multi-omics data, including metabolomics and single-cell sequencing, could elucidate how cofactor availability shapes tumor-immune crosstalk. Additionally, in vitro and in vivo models targeting these genes may reveal therapeutic vulnerabilities for high-risk COAD.

Conclusions

Our CVMRG-based model provides a metabolic lens for prognostic stratification in COAD, bridging the gap between tumor metabolism and clinical outcomes. By highlighting the interplay of metabolic enzymes, immune evasion, and drug response, this work paves the way for tailored therapies combining metabolic inhibitors with immunotherapy, ultimately improving outcomes for patients with this lethal malignancy.

Acknowledgments

We thank all database of the studies for providing the free data and the free R software which was used for analysis.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1521/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1521/prf

Funding: This work was supported by the grant from Hubei Provincial Department of Science and Technology (No. 2023AFB503).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1521/coif). All authors report the grant from Hubei Provincial Department of Science and Technology (No. 2023AFB503). The authors have no other conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Bray F, Laversanne M, Sung H, et al. Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2024;74:229-63. [Crossref] [PubMed]
2. Morgan E, Arnold M, Gini A, et al. Global burden of colorectal cancer in 2020 and 2040: incidence and mortality estimates from GLOBOCAN. Gut 2023;72:338-44. [Crossref] [PubMed]
3. Martelli V, Pastorino A, Sobrero AF. Prognostic and predictive molecular biomarkers in advanced colorectal cancer. Pharmacol Ther 2022;236:108239. [Crossref] [PubMed]
4. Viscaino M, Torres Bustos J, Muñoz P, et al. Artificial intelligence for the early detection of colorectal cancer: A comprehensive review of its advantages and misconceptions. World J Gastroenterol 2021;27:6399-414. [Crossref] [PubMed]
5. Martínez-Reyes I, Chandel NS. Cancer metabolism: looking forward. Nat Rev Cancer 2021;21:669-80. [Crossref] [PubMed]
6. Faubert B, Solmonson A, DeBerardinis RJ. Metabolic reprogramming and cancer progression. Science 2020;368:eaaw5473. [Crossref] [PubMed]
7. Dai Z, Ramesh V, Locasale JW. The evolving metabolic landscape of chromatin biology and epigenetics. Nat Rev Genet 2020;21:737-53. [Crossref] [PubMed]
8. Chornyi S, IJlst L, van Roermund CWT, et al. Peroxisomal Metabolite and Cofactor Transport in Humans. Front Cell Dev Biol 2020;8:613892. [Crossref] [PubMed]
9. Liu W, Wang P. Cofactor regeneration for sustainable enzymatic biosynthesis. Biotechnol Adv 2007;25:369-84. [Crossref] [PubMed]
10. Kreuzaler P, Inglese P, Ghanate A, et al. Vitamin B(5) supports MYC oncogenic metabolism and tumor progression in breast cancer. Nat Metab 2023;5:1870-86. [Crossref] [PubMed]
11. Sun L, Zhang H, Gao P. Metabolic reprogramming and epigenetic modifications on the path to cancer. Protein Cell 2022;13:877-919. [Crossref] [PubMed]
12. Li X, Egervari G, Wang Y, et al. Regulation of chromatin and gene expression by metabolic enzymes and metabolites. Nat Rev Mol Cell Biol 2018;19:563-78. [Crossref] [PubMed]
13. Camarena V, Wang G. The epigenetic role of vitamin C in health and disease. Cell Mol Life Sci 2016;73:1645-58. [Crossref] [PubMed]
14. He X, Wang Q, Cheng X, et al. Lysine vitcylation is a vitamin C-derived protein modification that enhances STAT1-mediated immune response. Cell 2025;188:1858-1877.e21. [Crossref] [PubMed]
15. Lyon P, Strippoli V, Fang B, et al. B Vitamins and One-Carbon Metabolism: Implications in Human Health and Disease. Nutrients 2020;12:2867. [Crossref] [PubMed]
16. Gao R, Wu C, Zhu Y, et al. Integrated Analysis of Colorectal Cancer Reveals Cross-Cohort Gut Microbial Signatures and Associated Serum Metabolites. Gastroenterology 2022;163:1024-1037.e9. [Crossref] [PubMed]
17. Colaprico A, Silva TC, Olsen C, et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res 2016;44:e71. [Crossref] [PubMed]
18. Marisa L, de Reyniès A, Duval A, et al. Gene expression classification of colon cancer into molecular subtypes: characterization, validation, and prognostic value. PLoS Med 2013;10:e1001453. [Crossref] [PubMed]
19. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 2014;15:550. [Crossref] [PubMed]
20. Kanehisa M, Furumichi M, Tanabe M, et al. KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res 2017;45:D353-61. [Crossref] [PubMed]
21. Robin X, Turck N, Hainard A, et al. pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinformatics 2011;12:77. [Crossref] [PubMed]
22. Yu G, Wang LG, Han Y, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 2012;16:284-7. [Crossref] [PubMed]
23. Ashburner M, Ball CA, Blake JA, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet 2000;25:25-9. [Crossref] [PubMed]
24. Subramanian A, Tamayo P, Mootha VK, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A 2005;102:15545-50. [Crossref] [PubMed]
25. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43:e47. [Crossref] [PubMed]
26. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013;14:7. [Crossref] [PubMed]
27. Szklarczyk D, Gable AL, Lyon D, et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res 2019;47:D607-13. [Crossref] [PubMed]
28. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612. [Crossref] [PubMed]
29. Zeng D, Ye Z, Shen R, et al. IOBR: Multi-Omics Immuno-Oncology Biological Research to Decode Tumor Microenvironment and Signatures. Front Immunol 2021;12:687975. [Crossref] [PubMed]
30. Yang W, Soares J, Greninger P, et al. Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res 2013;41:D955-61. [Crossref] [PubMed]
31. Guinney J, Dienstmann R, Wang X, et al. The consensus molecular subtypes of colorectal cancer. Nat Med 2015;21:1350-6. [Crossref] [PubMed]
32. Iasonos A, Schrag D, Raj GV, et al. How to build and interpret a nomogram for cancer prognosis. J Clin Oncol 2008;26:1364-70. [Crossref] [PubMed]
33. Miller KD, Nogueira L, Devasia T, et al. Cancer treatment and survivorship statistics, 2022. CA Cancer J Clin 2022;72:409-36. [Crossref] [PubMed]
34. Sutendra G, Michelakis ED. Pyruvate dehydrogenase kinase as a novel therapeutic target in oncology. Front Oncol 2013;3:38. [Crossref] [PubMed]
35. Wang N, Lu S, Cao Z, et al. Pyruvate metabolism enzyme DLAT promotes tumorigenesis by suppressing leucine catabolism. Cell Metab 2025;37:1381-1399.e9. [Crossref] [PubMed]
36. Entschladen F, Drell TL 4th, Lang K, et al. Tumour-cell migration, invasion, and metastasis: navigation by neurotransmitters. Lancet Oncol 2004;5:254-8. [Crossref] [PubMed]
37. Pastushenko I, Blanpain C. EMT Transition States during Tumor Progression and Metastasis. Trends Cell Biol 2019;29:212-26. [Crossref] [PubMed]
38. Nusse R, Clevers H. Wnt/β-Catenin Signaling, Disease, and Emerging Therapeutic Modalities. Cell 2017;169:985-99. [Crossref] [PubMed]
39. Vasan K, Werner M, Chandel NS. Mitochondrial Metabolism as a Target for Cancer Therapy. Cell Metab 2020;32:341-52. [Crossref] [PubMed]
40. Longley DB, Harkin DP, Johnston PG. 5-fluorouracil: mechanisms of action and clinical strategies. Nat Rev Cancer 2003;3:330-8. [Crossref] [PubMed]
41. Sastre J, Argilés G, Benavides M, et al. Clinical management of regorafenib in the treatment of patients with advanced colorectal cancer. Clin Transl Oncol 2014;16:942-53. [Crossref] [PubMed]
42. Di Nicolantonio F, Vitiello PP, Marsoni S, et al. Precision oncology in metastatic colorectal cancer - from biology to medicine. Nat Rev Clin Oncol 2021;18:506-25. [Crossref] [PubMed]
43. Yoo HS, Cockrum MA, Napoli JL. Cyp26a1 supports postnatal retinoic acid homeostasis and glucoregulatory control. J Biol Chem 2023;299:104669. [Crossref] [PubMed]