Metabolic-cell-death gene trio predicts survival and cuproptosis sensitivity in colorectal cancer

Introduction

Colorectal cancer (CRC) represents a formidable global oncologic burden, accounting for approximately 10% of cancer incidence and 9% of cancer-related mortality worldwide (1,2). In 2020 alone, this disease contributed to over 1.9 million new cases and 935,000 deaths, constituting nearly one-tenth of the global cancer burden (1). Although the overall incidence has plateaued or declined in developed countries, an alarming epidemiological shift has emerged: a steady rise in both incidence and mortality among individuals younger than 50 years (3). This early-onset CRC frequently presents at advanced clinical stages, characterized by distant organ involvement and high recurrence rates. Compounding this challenge, the prognosis for advanced-stage patients remains dismal, with 5- and 10-year survival rates of merely 65% and 58%, respectively (4,5).

Current diagnostic paradigms predominantly rely on colonoscopy and biopsy; however, conventional biomarkers such as carcinoembryonic antigen (CEA) and carbohydrate antigen 19-9 (CA19-9) exhibit inadequate sensitivity and specificity for effective diagnosis and recurrence surveillance (6). This biomarker deficiency precludes early detection for most patients, resulting in a substantial proportion presenting with advanced disease at initial diagnosis. While therapeutic modalities—including surgery, chemotherapy, radiotherapy, and immunotherapy—have evolved considerably, clinical benefits for CRC patients remain limited (7). These unmet clinical needs underscore the urgent imperative to identify novel molecular targets and establish robust prognostic risk models that could enable personalized therapeutic strategies and improve patient outcomes.

Metabolic cell death (MCD) represents a distinct form of regulated cell death triggered by overload or depletion of specific nutrients (e.g., glucose and amino acids) or metals (e.g., iron and copper), leading to cellular metabolic imbalance. Although the mechanisms of MCD are complex, several recently discovered pathways have been characterized, including ferroptosis, cuproptosis, disulfidptosis, lysosomal zinc death, and alkaliptosis (8,9). Cancer cells frequently undergo metabolic reprogramming, which exposes unique metabolic vulnerabilities that can be therapeutically exploited by inducing specific MCD pathways (10). In CRC, accumulating evidence demonstrates that ferroptosis induction—achieved by elevating intracellular Fe²⁺ and reactive oxygen species levels, depleting glutathione (GSH), or inactivating glutathione peroxidase 4 (GPX4)—holds therapeutic promise, while ferroptosis inhibition promotes tumor progression and chemoresistance (11,12). Additionally, reduced expression of cuproptosis-related gene FDX1 correlates with diminished disease-free and overall survival (OS) in CRC patients and increased metastatic risk (13). Furthermore, loss of disulfidptosis-related MSH3 protein has been identified as a prevalent feature in MLH1-deficient CRC, where it associates with adverse prognosis (14). Collectively, these findings establish cuproptosis and disulfidptosis as potential prognostic determinants in CRC (15). However, the comprehensive prognostic value and biological functions of MCD in CRC remain inadequately investigated.

To systematically delineate the key targets of MCD in CRC, we integrated The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) transcriptomic datasets and constructed an MCD-related gene (MRG)-based prognostic risk signature via comprehensive bioinformatics. Based on this model, we characterized the immune microenvironment, mutational landscape, and drug-sensitivity profiles across distinct risk strata. Furthermore, we validated the expression of signature genes in clinical CRC tissues and functionally demonstrated that silencing the representative prognostic gene CDKN2A exerts tumor-suppressive effects by enhancing colon cancer cell sensitivity to cuproptosis, as evidenced by reverse transcription quantitative real-time polymerase chain reaction (RT-qPCR), Western blot (WB), 5-ethynyl-2’-deoxyuridine (EdU), and cell viability assay (CCK-8) assays. Our findings offer novel insights and therapeutic implications for MCD-guided prognostication and immunotherapy in CRC. We present this article in accordance with the MDAR and TRIPOD reporting checklists (available at https://jgo.amegroups.com/article/view/10.21037/jgo-2026-0419/rc).

Methods

Data collection

The data of CRC patients [colon cancer (COAD) and rectal cancer (READ)] with their RNA-seq, survival information, and clinical characteristics were downloaded from the TCGA database (https://portal.gdc.cancer.gov/). The TCGA-COAD dataset contained 41 normal tissue samples and 455 COAD tissue samples, and the TCGA-READ dataset contained 10 normal tissue samples and 163 READ tissue samples (visit time: January 2, 2025). The TCGA-COAD and TCGA-READ datasets were subsequently combined and designated as the TCGA-CRC dataset, which included a total of 618 CRC tissue samples and 51 normal tissue samples. This dataset was utilized as the training set. Among the 618 CRC samples, survival information was available for 585 patients. Additionally, in the GEO database (https://www.ncbi.nlm.nih.gov/geo/), 556 CRC tissue samples were picked from the GSE39582 dataset on the GPL570 platform. The 556 CRC samples, which included survival details and gene expression data, were available for use as a validation set. MCD encompassed various forms, including ferroptosis, cuproptosis, disulfidptosis, lysozincrosis, and alkaliptosis (8). Genes associated with the 5 categories of MCD were collected from the FerrDb database (http://www.zhounan.org/ferrdb/current/) and relevant literature (8,16), respectively. The five categories of genes were combined, resulting in a total of 628 unique genes after the removal of duplicates.

Weighted gene co‑expression network analysis (WGCNA)

To obtain important modules genes related to CRC from the training set, the “WGCNA” (v 1.18.0) was used. The hierarchical clustering analysis was conducted on all samples using the Euclidean distance metric to identify and exclude outliers in the sample expression profiles. In order to determine an appropriate soft threshold, the pick Soft Threshold function was utilized to calculate the average number of connections and the R² (scale-free fitting index) for various soft thresholds. Subsequently, the genes were classified into distinct modules in accordance with the specifications of the hybrid dynamic tree cutting algorithm, with a minimum of 100 genes per module and a module merge parameter of 0.4. The cor function was essential for conducting a Spearman correlation analysis, which aimed to determine between modules and CRC (P<0.05, |correlation (r)| >0.3). The modules exhibiting the highest positive and negative correlations with CRC were identified as key modules based on the evaluation of correlation coefficients and significance levels. Key module genes were combined and documented as module genes of CRC for subsequent analysis.

Differential expression analysis

To identify differentially expressed genes (DEGs) in the TCGA-CRC set, the “DESeq2” (v 1.38.0) (17) was used to perform differential expression analysis on CRC and normal samples (|log₂ fold change (FC)| >1, adj. P<0.05). Moreover, the “ggplot2” (v 3.4.1) (18) was employed to construct a volcano plot of DEGs, and the top 10 genes (sorted by |log₂FC| from high to low) with significant up- and down-regulation differences were labeled, and then the heatmap of these genes was drawn using the “pheatmap” (v 1.0.12) (19).

Identification, functional enrichment, and protein-protein interaction (PPI) network of candidate genes

To identify the genes associated with MCD-RGs that were differentially expressed in CRC, the “ggvenn” (v 0.1.9) (20) was employed to determine the intersection of DEGs, MCD-RGs, and module genes, and these genes were recorded as candidate genes. Subsequently, an analysis of these genes was conducted to investigate their biological pathways and functional enrichment, utilizing Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) through the application of the “clusterProfiler” (v 4.2.2) (21) (P<0.05). Candidate genes were submitted to the STRING database (https://string-db.org/) for the purpose of constructing a PPI network, utilizing an interaction score threshold greater than 0.4. This approach aimed to elucidate the interactions among genes at the protein level. Subsequently, the identified interactions were visualized employing Cytoscape software (v 3.9.1) (22).

Identification of prognostic genes

To identify candidate prognostic genes that exhibit a significant correlation with OS in CRC, a univariate Cox regression analysis was performed on the candidate genes using the “survival” (v 3.5.3) (23) in the 585 CRC samples with survival information [hazard ratio (HR) ≠1, P<0.1]. Candidate prognostic genes were further subjected to the proportional hazards (PH) assumption test. Genes with a P<0.05 were filtered out, while those with a P>0.05 were retained. Additionally, the “glmnet” (v 4.1.4) (24) was employed to conduct least absolute shrinkage and selection operator (LASSO) regression analysis on the remaining genes. The optimal lambda value, which corresponded to the minimal error, was identified through cross-validation in order to select prognostic genes.

Construction and validation of a risk model

Based on prognostic genes, the “random Forest SRC” (v 3.2.2) (25) was utilized to create random survival forest (RSF) model from 585 CRC samples with survival information. The model parameters, ntree and mtry, were adjusted to calculate the risk score for each patient when the RSF model exhibited the lowest error rate. Through the computation of the best cutoff value derived from the risk score, the cohort of 585 CRC patients was classified into 2 distinct categories: a high-risk group (HRG) and a low-risk group (LRG). The implementation of the “pheatmap” (v 1.0.12) enabled the visualization of prognostic gene expression. The Kaplan-Meier (K-M) survival curve was generated utilizing the “survminer” (v 0.4.9) (23) to investigate the differences in OS rates between HRG and LRG (log-rank test, P<0.05). Additionally, “survival ROC” (v 1.0.3.1) (26) was utilized to create receiver operating characteristic (ROC) curves [with an area under the curve (AUC) >0.6] corresponding to the 1-, 2-, and 3-year time intervals. Finally, an additional assessment was conducted to evaluate the predictive accuracy and robustness of the risk model, with validation conducted in the validation set.

Independent prognostic analysis and construction of the nomogram

The distributions of risk score were compared among 585 CRC patients with clinical characteristics from the training set. Specifically, Wilcoxon test was employed to investigate the variations in risk score across the age and gender subgroups independently (P<0.05). Additionally, Kruskal-Wallis test was utilized to analyze the differences in risk score among the subgroups categorized by T stage, N stage, and overall stage (P<0.05).

Furthermore, in order to identify independent risk factors, clinical characteristics and risk score were analyzed through univariate/multivariate Cox regression analysis (HR ≠1, P<0.05), as well as PH assumption test (P>0.05). Following this, a nomogram was constructed utilizing the independent risk factors, employing “rms” (v 6.5.0) (https://CRAN.R-project.org/package=rms), to estimate the probability of survival for CRC at 1-, 2-, and 3-year. To assess the reliability of this nomogram, calibration curves were generated with the utilization of “rms” (v 6.5.0). Additionally, ROC curves at time spans of 1-, 2-, and 3-year were plotted by the “survivalROC” (v 1.0.3.1) to evaluate the nomogram (AUC >0.7).

Gene set enrichment analysis (GSEA) and gene set variation analysis (GSVA)

GSEA was carried out with the objective of elucidating the biological functions in HRG and LRG of CRC patients. The reference gene set (c2.cp.kegg.v7.5.1.symbols.gmt) were selected from the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb). In the analysis of the 585 CRC samples, the “DESeq2” (v 1.40.2) was employed to assess the differential expression between HRG and LRG. Subsequently, the genes were ranked in descending order based on their log₂ FC values. GSEA was performed using “clusterProfiler” (v 4.2.2), with a threshold of |NES| >1 and P. adjust <0.05.

To further investigate the biological functional differences between HRG and LRG, GSVA was performed. Specifically, the reference gene set (c2.cp.kegg.v2023.1.Hs.symbols.gmt) were downloaded from MSigDB database. Of the 585 CRC samples, GSVA score for each sample was calculated using the single-sample gene set enrichment analysis (ssGSEA) algorithm from the “GSVA” (v 1.46.0) (27). The “limma” (v 3.54.1) (17) was then applied to compare the differences in GSVA scores between HRG and LRG, with the thresholds established at |t|>2 and P adjust <0.05. The top 10 activated pathways and top 10 inhibited pathways (sorted by |t| from high to low) were identified as significantly enriched pathways for presentation. Spearman correlation analyses were conducted using the “psych” (v 2.2.9) (28) to investigate the relationships among risk score, prognostic genes, and significantly enriched pathways (|r|>0.3, P<0.05).

Analysis of metabolic cell death pathway activity and enrichment

The activity scores of MCD pathways, including ferroptosis, cuproptosis, disulfidptosis, and alkaliptosis, were quantified for each sample using ssGSEA implemented in the “GSVA” (v 1.46.0) package. Differences in pathway scores were compared between the high-risk and low-risk groups using the Wilcoxon rank-sum test, and P values were adjusted using the Benjamini-Hochberg (BH) method. An adjusted P value <0.05 was considered statistically significant. To investigate the association between the prognostic genes and MCD pathway scores, Spearman correlation analysis was performed. Furthermore, GSEA was conducted to determine the enrichment of MCD-related pathways between the two risk groups. Adjusted P values < 0.05 were considered statistically significant.

Immune microenvironment analysis

To evaluate the tumor microenvironment in CRC patients, we calculated tumor microenvironment scores—including the immune score, ESTIMATE score, and stromal score—for each sample using the ESTIMATE algorithm. This analysis was conducted on 585 CRC samples. Wilcoxon test was then utilized to compare the differences of these scores between 2 groups (P<0.05). Spearman correlation analyses were conducted using the “psych” (v 2.2.9) to investigate the relationship between risk score and tumor microenvironment scores (|r|>0.3, P<0.05).

According to the ssGSEA algorithm, the “GSVA” (v 1.46.0) was utilized to determine the infiltration scores of 28 kinds of immune cells (29) between the HRG and LRG in the 585 CRC samples with survival information. Next, the Wilcoxon test was applied to make a comparison of the differences in immune-infiltrating cells between the 2 groups (P<0.05). The “psych” (v 2.2.9) was further utilized to calculate correlations and significance among risk score, prognostic genes, and different immune cells (|r|>0.3, P<0.05).

To further understand the relationship between risk score, prognostic genes, immune indicators, and immune checkpoint genes, cytotoxic activity (CYT) score and tumor inflammatory signature (TIS) score for each patient of 585 CRC samples were calculated. Additionally, the correlations of risk score and prognostic genes with 2 immune indicators (CYT and TIS) and immune checkpoint genes (30) were analyzed using the “psych” (v 2.2.9) (|r|>0.3, P<0.05).

Cancer immunity cycle analysis

A compilation of genes associated with the cancer immunity cycle was sourced from the Immunophenotype (TIP) database (http://biocc.hrbmu.edu.cn/TIP). Based on the 585 CRC samples, each sample was given an enrichment score based on the cancer immune cycle using the ssGSEA algorithm. We then employed Wilcoxon test to compare the differences in enrichment scores of the cancer immune cycle between 2 groups (P<0.05).

Chemotherapeutic drug sensitivity analysis

To investigate the therapeutic effects of chemotherapeutic drugs in CRC patients, a total of 198 chemotherapeutic drugs were obtained from the Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/, v Release 8.5). Drug response prediction was performed using the “oncoPredict” package (v 1.2) (31), which estimates the half-maximal inhibitory concentration (IC₅₀) of each agent based on gene expression profiles. Lower IC₅₀ values represent higher sensitivity to the corresponding drug. Differences in predicted IC₅₀ values between the high-risk and low-risk groups were evaluated using the Wilcoxon rank-sum test, and drugs with P<0.05 were considered significantly differentially sensitive agents. The top 10 chemotherapeutic drugs with the most significant P values were displayed. Additionally, the correlations between the risk score and the top 10 chemotherapeutic drugs were analyzed using Spearman correlation analyses from the “psych” (v 2.2.9) (|r|>0.3, P<0.05).

Analysis of genetic variation in prognostic genes

To assess the genetic variation of prognostic genes during the development of CRC, the genetic variation of these genes in the TCGA dataset was analyzed using the cBioPortal database (https://www.cbioportal.org/). This included data on mutation types, mutation frequencies, and copy number alteration (CNA) for the gastrointestinal dataset, as well as an analysis of OS (P<0.05).

Protein and mRNA levels analysis of the prognostic genes

To explore the differential protein expression of prognostic genes between CRC and normal tissues, immunohistochemistry (IHC)-based protein expression profiles of these prognostic genes were retrieved from the Human Protein Atlas (HPA) database (http://www.proteinatlas.org/). To evaluate the expression levels of prognostic genes in CRC and normal samples, the relative expression levels of these genes were analyzed using UALCAN (http://ualcan.path.uab.edu) (from highest to lowest) in the different organizations, including comparisons between tumours of different stages and subtypes. The percentage of patients exhibiting high and moderate expression of prognostic gene proteins in each type of cancer was further illustrated. The Wilcoxon test (P<0.05) was employed to analyze the mRNA expression of these genes within both the CRC and normal samples derived from the TCGA-CRC dataset.

Human CRC and adjacent normal tissues

Tissue samples, including cancer tissues and their corresponding adjacent non-cancerous tissues, were collected from patients with CRC who were treated at the Department of Gastroenterology, The Second People’s Hospital of Guiyang (Jinyang Hospital), between January 2024 and January 2025. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Ethics Committee of The Second People’s Hospital of Guiyang (Jinyang Hospital) (No. JYYY-2025-XM-31). Written informed consent was obtained from all individual participants included in the study before sample collection. We enrolled five male patients with advanced colon cancer, aged between 40 and 50 years. All patients were classified as stage II according to the American Joint Committee on Cancer staging system and had no other underlying diseases.

Cell culture and siRNA transfection

The human CRC cell lines SW480 (RRID: CVCL_0546), SW620 (RRID: CVCL_1763), HCT116 (RRID: CVCL_0291), and the normal colonic epithelial cell line NCM460 (RRID: CVCL_9267) were obtained from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China) and maintained in DMEM (Gibco) supplemented with 10% FBS and 1% penicillin-streptomycin at 37 ℃ in 5% CO₂. For transient knock-down, 1×10⁵ cells were seeded in 6-well plates and transfected with 50 nM siRNA (GenePharma) using Lipofectamine 3000 (Invitrogen) according to the manufacturer’s instructions.

RT-qPCR

Total RNA was extracted with TRIzol reagent (Invitrogen) and reverse-transcribed using the PrimeScript RT Reagent Kit (Takara). qPCR was performed on a CFX96 system (Bio-Rad) with TB Green Premix Ex Taq II (Takara). Relative expression was calculated by the 2^-ΔΔCT method with GAPDH as endogenous control. Primers are provided in Table S1.

Western blotting (WB)

Cells or snap-frozen tissues were lysed in RIPA buffer containing protease and phosphatase inhibitors. Protein (30 µg) was resolved on 10–12% ulfatepolyacrylamide gel electrophoresis (SDS-PAGE), transferred to polyvinylidene difluoride (PVDF) membranes (Millipore), blocked with 5% non-fat milk and probed overnight at 4 ℃ with primary antibodies against CDKN2A [Proteintech, catalogue number (Cat No.) 10883-1-AP], Lip-DLAT (Abcam, Cat No. ab58724), Lip-DLST (Abcam, Cat No. ab187699), CTR1 (Solarbio, Cat No.: K900084P) and GAPDH (Proteintech, Cat No. 60004-1-Ig). After incubation with HRP-conjugated secondary antibodies, bands were visualised using enhanced chemiluminescence (ECL) substrate (Bio-Rad) and quantified with ImageJ software.

Cell viability assay (CCK-8)

Cells (3×10³ per well) were seeded in 96-well plates, transfected the next day, and incubated for 0–5 days. 10 µL CCK-8 reagent (Dojindo) was added for 2 h and absorbance at 450 nm was recorded on a microplate reader (BioTek). Each time-point contained six replicates.

EdU proliferation assay

Forty-eight hours post-transfection, cells were exposed to 50 µM EdU (RiboBio) for 2 h, fixed with 4% paraformaldehyde, permeabilised with 0.3% Triton X-100 and stained with Azide 488 and Hoechst following the kit protocol. EdU-positive nuclei were counted in five random fields per well using a fluorescence microscope (Olympus).

Transwell invasion assay

Matrigel (BD Biosciences) was diluted 1:8 and coated on 8-µm pore inserts (Corning). 2×10⁴ serum-starved cells in 200 µL DMEM were added to the upper chamber; 600 µL complete medium served as chemoattractant. After 48 h, invaded cells were fixed with methanol, stained with 0.1% crystal violet and counted in five fields per insert.

Wound-healing (scratch) migration assay

Confluent monolayers were scratched with a 200-µL pipette tip, washed twice to remove debris, and cultured in serum-reduced (1% FBS) medium. Images were captured at 0, 2, 24 and 48 h using an inverted microscope. Wound area was measured with ImageJ and percentage closure calculated.

Copper treatment and cuproptosis induction

Where indicated, 5 µM CuCl₂ (Sigma) was added 12 h after siRNA transfection and cells were harvested 24 h later for Western blot analysis of cuproptosis-related proteins.

Statistical analysis

All bioinformatics analyses were carried out with the R software (v 4.2.2). Data are presented as mean ± standard deviation from at least three independent experiments. Comparisons were performed using two-tailed Student’s t-test or one-way ANOVA followed by Tukey’s post-hoc test. P<0.05 was considered statistically significant.

Results

The 3,274 module genes were associated with CRC

In the WGCNA analysis, the results indicated that no obvious outliers were present in the samples (Figure S1A). In the process of building the scale-free network, the ideal soft threshold was identified as 12, corresponding to a scale-free R² of 0.85, and the mean connectivity approached but did not reach 0 (Figure S1B). Based on hierarchical cluster analysis, 21 gene modules (removing the grey module) were identified in total (Figure S1C). Then, MEgreen module (r=0.42, P<0.0001) and MEpurple module (r=−0.85, P<0.0001) were selected as the key modules associated with CRC based on the results of the correlation analysis between the modules and the CRC (Figure S1D). Among them, there were 2,401 genes in the MEgreen module and 873 genes in the MEpurple module, resulting in a total of 3,274 module genes related to CRC.

Differential expression analysis showed that there were 9,375 DEGs between the CRC group and the normal group. Among them, 5,611 genes in the CRC group were identified as upregulated genes, and 3,764 genes were identified as downregulated genes (Figure 1A,1B). Moreover, a total of 39 shared genes between DEGs, module genes, and MCD-RGs were identified and selected as candidate genes (Figure 1C). The enrichment analysis of the 39 candidate genes revealed associations with 15 GO terms, including 391 BPs, 108 MFs, and 25 CCs. These terms encompassed a wide range of functions such as phospholipid metabolic process, lipid homeostasis, melanosome, and kinase inhibitor activity among others (Figure 1D). Moreover, the KEGG enrichment analysis of the candidate genes revealed that a total of 12 KEGG pathways were enriched, for instance MicroRNAs in cancer, ferroptosis, p53 signaling pathway, sulfur metabolism, and central carbon metabolism in cancer (Figure 1E). Next, a PPI network consisting of 28 interaction relationships corresponding to 39 candidate genes was constructed (Figure 1F), among which 14 genes formed isolated targets. In this network, SLC3A2, CDKN2A, FGFR4, and RPL8 had frequent protein-level interactions with other genes.

Figure 1 Integrated multi-omics and functional enrichment dissect the core signatures of differentially expressed genes in colorectal cancer. (A,B) Differentially expressed genes between colorectal cancer and adjacent normal tissues. (A) Volcano plot displaying log2 fold change versus −log10 P value. (B) Heatmap showing expression patterns of top DEGs. (C) Venn diagram illustrates the intersection of DEGs, WGCNA module genes, and MRGs. (D) GO enrichment analysis of candidate genes. (E) Top significant terms in BP, CC, and MF categories. (F) KEGG pathway enrichment analysis of candidate genes. BP, biological process; CC, cellular component; CRC, colorectal cancer; DEGs, differentially expressed genes; FC, fold change; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; MRGs, mitochondrial-related genes.

A risk model based on 3 prognostic genes demonstrated superior predictive power in CRC

Next, a univariate Cox regression analysis was performed. The results showed that 3 candidate prognostic genes were significantly associated with OS of CRC (HR ≠1, P<0.1) (Figure 2A). Among the genes examined, MPC1 and AHCY were found to be associated with a more favourable prognosis (HR <1), indicating that they might exert an inhibitory effect on CRC progression. In contrast, CDKN2A was linked to a worse prognosis (HR >1), indicating that they might contribute to CRC progression. Furthermore, all 3 genes successfully met the PH assumption test criteria (P>0.05 (Figure S2A-S2C). Following LASSO regression analysis, 3 prognostic genes were identified at optimal lambda 0.003676226, including CDKN2A, MPC1, and AHCY (Figure 2B,2C). Subsequently, the RSF model exhibited the lowest error rate with ntree set to 12 and mtry set to 1, resulting in the calculation of risk score (Figure S2D). The 585 CRC samples were classified into HRG (n=176) and LRG (n=409) according to the best cutoff value for the risk score of 23.44623. It was noteworthy that the mortality rate of CRC patients increased in a progressive manner with increasing risk score (Figure 2D,2E). In the HRG, CDKN2A was highly expressed, while the expression of MPC1 and AHCY was reduced (Figure 2F). A substantial survival disparity was observed between HRG and LRG via K-M curve, with the HRG demonstrating a marked reduction in survival probability (P<0.0001) (Figure 2G). The 1-, 2-, and 3-year AUCs were 0.83, 0.83, and 0.85, respectively (Figure 2H). The AUC values were all found to be greater than 0.6, and the risk model demonstrated a more robust predictive power with CRC prognostic value. Similarly, in the GSE39582 set, best cutoff value for the risk score of 23.44094 was used to classify CRC patients into HRG (n=193) and LRG (n=363).

Figure 2 Three metabolic-related gene (CDKN2A/MPC1/AHCY) features predict colorectal cancer prognosis. (A) Identification of prognosis-related genes by univariate Cox regression analysis. (B,C) LASSO regression model construction for survival-associated genes using the glmnet package. (B) Coefficient profiles of selected genes across regularization parameters. (C) Ten-fold cross-validation error curve; vertical lines denote λ.min and λ.1se. (D,E) Risk stratification of colorectal cancer patients. (D) Distribution of risk scores with optimal cutoff value indicated. (E) Survival status and time in high-risk and low-risk groups. (F) Expression heatmap of prognostic genes in high-risk versus low-risk groups. (G) Time-dependent ROC curve analysis evaluating the predictive performance of the three-gene signature. (H) Kaplan-Meier survival curves comparing high- and low-risk groups (log-rank test). AUC, area under the curve; CI, confidence interval; CRC, colorectal cancer; FP, false positive; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic; TP, true positive.

The pattern of risk score distribution, the corresponding survival status and prognostic genes expression (Figure S3A-S3C), the K-M curves (P=0.00029) (Figure S3D), and ROC curve analysis (Figure S3E) confirmed the robustness of the constructed risk model. These findings indicated that the risk model served as a valuable tool for personalized prognostic assessment in the clinical management of CRC.

A nomogram predicting CRC survival outcome was created using risk score, age, and stage as independent prognostic factors

With respect to clinical characteristics, the data revealed significant variations in risk score among the T stage, N stage, and stage groups (P<0.05) (Figure 3A-3C). The univariate Cox regression analysis and PH assumption test identified risk score, age, T stage, N stage, and stage as risk factors for OS in patients with CRC (Figure S4). Subsequently, risk score, age, and stage remained independent prognostic factors in the multivariate Cox analysis for OS of patients with CRC (Figure 3D). Then a nomogram incorporating these factors was devised to forecast 1-, 2-, and 3-year OS in CRC (Figure 3E). Furthermore, the calibration curve demonstrated that the actual 1-, 2-, and 3-year survival exhibited high degree of concordance with the predicted values (Figure 3F). The ROC curves verified the desirable efficacy of nomogram in predicting CRC patients’ 1-, 2-, and 3-year survival outcomes (Figure 3G). The AUC values were all found to be greater than 0.7, and the risk model demonstrated a more robust predictive power with CRC prognostic value.

Figure 3 Prognostic model development and validation based on risk scoring system. (A-C) Differential analysis of risk scores among groups stratified by different clinical characteristics (Stage, T stage, and N stage). The results demonstrate significant variations in risk scores across different staging classifications, highlighting the prognostic relevance of our risk scoring system. (D) Multivariate Cox regression analysis of factors meeting the PH assumption. Only variables that satisfied the PH assumption were included in the final multivariate model to ensure statistical validity. (E) Nomogram construction based on prognostic genes for predicting 1-, 2-, and 3-year OS. The nomogram integrates multiple prognostic factors to provide individualized survival probability estimates. (F) Calibration curves for 1-, 2-, and 3-year survival predictions using the nomogram constructed with R package “rms”. The calibration plots demonstrate good agreement between predicted and observed survival probabilities across all time points. (G) ROC curves generated using R package “survivalROC” for evaluating the predictive performance of the prognostic model. AUC, area under the curve; CI, confidence interval; N, node; OS, overall survival; PH, proportional hazards; ROC, receiver operating characteristic; T, tumor.

Exploring the functional and immune infiltration differences between HRG and LRG

GSEA and GSVA are advantageous method for analyzing biological differences between different risk groups. Then, the signaling pathways between the HRG and LRG were explored by GSEA and GSVA. GSEA results showed that a total of 82 signalling pathways were significantly different between the two risk groups, such as ribosome (Figure 4A). Among these, glycosaminoglycan biosynthesis chondroitin sulfate and dilated cardiomyopathy were found to be activated in HRG, while Selen amino acid metabolism and pyruvate metabolism were observed to be inhibited in HRG (Figure 4B). Correlation results indicated that the risk score associated with the peroxisome and Selen amino acid metabolism exhibited a significant negative correlation. In contrast, the prognostic genes MPC1 and AHCY were significantly and positively associated with terpenoid backbone biosynthesis and Selen amino acid metabolism, respectively (Figure 4C).

Figure 4 Multi-level characterization of high- vs. low-risk groups. (A) GSEA visualization of hallmark pathways enriched between high- and low-risk groups using the clusterProfiler R package. (B) Top 10 activated (blue) and suppressed (green) Gene Ontology biological-process terms ranked by NES in the high-risk versus low-risk contrast. (C) Spearman correlation matrix (psych R package) among risk score, prognostic signature genes and the differentially enriched pathways. (D) Box-plots comparing ESTIMATE-derived immune, stromal and composite scores between the two risk strata. (E) Correlation heatmap linking individual risk score and signature genes to the three ESTIMATE components. (F) Relative abundance of 28 immune cell subsets quantified by ssGSEA in CRC patients stratified by risk. (G) Differential infiltration of immune cell populations between high- and low-risk groups; coloured bars indicate adjusted P values. (H) Spearman correlations among risk score, signature genes and the immune-cell subsets that differed significantly between groups. CRC, colorectal cancer; GSEA, gene set enrichment analysis; NES, normalized enrichment score; ssGSEA, single-sample gene set enrichment analysis.

To further analyze differences in MCD-related pathway activity between the high- and low-risk groups. The results showed that ferroptosis and disulfidptosis scores were significantly higher in the high-risk group than in the low-risk group (P<0.05), whereas cuproptosis and alkaliptosis did not differ significantly between the two groups (Figure S5A). Further correlation analysis showed that MPC1 was significantly positively correlated with cuproptosis (cor=0.34, P<0.05), whereas AHCY was significantly negatively correlated with disulfidptosis (cor=−0.414, P<0.05) (Figure S5B). The GSEA results further showed that the ferroptosis pathway had a significantly negative NES, suggesting that, under the current ranking direction, ferroptosis-related genes were more likely to be enriched in the low-risk group (Figure S5C,S5D). Collectively, these findings suggested that the three-gene signature was closely associated with MCD-related pathways, especially ferroptosis and disulfidptosis, while ferroptosis-related regulation might be context-dependent in CRC.

Immune score, stromal score, and ESTIMATE score exhibited significant differences between HRG and LRG with scores being elevated in the HRG (Figure 4D). Correlation results indicated that risk score was significant but weakly correlated with all three scores (Figure 4E). For HRG and LRG samples, the proportions of 28 kinds of immune cell score were displayed in Figure 4F. A comparative analysis of immune cell infiltration score revealed significant disparities among 17 distinct immune cell types between two groups (Figure 4G). The 17 distinct immune cells, including activated CD8 T cells, exhibited elevated levels of abundance in the HRG. Figure 4H described the correlation among risk score, prognostic genes, and differential immune cells. Myeloid-derived suppressor cells (MDSCs) demonstrated a significant positive correlation with macrophages. AHCY demonstrated significant negative correlation with type 2 T helper cells. Moreover, risk score demonstrated a significant positive correlation with CDKN2A and a significant negative correlation with MPC and AHCY. The complexity and activity of the tumor microenvironment might make the HRG of tumors more aggressive.

Significant differences were identified between HRG and LRG in immunotherapy and chemotherapy drug sensitivity

In the correlation analyses of risk score and prognostic genes with immune indicators and immune checkpoint genes, the results revealed significant but weak correlations (Figure S6A). The anti-tumor immune cycle represents the steps and interactions that generate an effective immune response against tumors. The results demonstrated that the enrichment scores of genes for cancer immune cycle in the HRG and LRG exhibited significant disparities, with the exception of stage 4 (Figure S6B). Furthermore, these enrichment scores were consistently higher in the HRG. Drug sensitivity prediction revealed that the high-risk group exhibited significantly lower predicted IC50 values for AZ960, AZD8186, JAK inhibitors, JQ1, PCI-34051, PRT062607, WIKI, and XAV939, suggesting that these patients may be more sensitive to these agents (Figure S6C). Moreover, the prognostic gene AHCY was significantly and positively associated with the chemotherapy drug AZD8186 (Figure S6D). To summarize, the findings of the analysis indicated that the risk score exhibited a degree of predictive value with regard to the efficacy of chemotherapy in patients diagnosed with CRC.

Analysis of prognostic genetic variation and association with gastrointestinal cancer

The most prevalent type of genetic alteration observed in AHCY across nearly all datasets was amplification (Figure 5A). In contrast, the most common genetic alterations in CDKN2A were deep deletion and mutation (Figure 5B). For MPC1, the predominant types of genetic alterations were mutation and amplification (Figure 5C). Compared to proteins of AHCY and MPC1, the CDKN2A protein had a greater number of mutation sites (Figure 5D). In addition, diploid status with mutations in the prognostic genes AHCY, CDKN2A, and MPC1 was the most prevalent type of mutation found in gastrointestinal cancers (Figure 5E-5G). Finally, a significant difference was observed between mutations in the prognostic genes AHCY, CDKN2A, and MPC1 and the survival of patients with gastrointestinal cancers (Figure 5H-5J). Genetic variants in the prognostic genes AHCY and MPC1 were associated with an increased survival rate for patients with gastrointestinal cancers, while genetic variants in CDKN2A were linked to a decreased survival rate in these patients.

Figure 5 Genetic landscape and clinical relevance of the three prognostic genes across gastrointestinal cancers. (A-C) Mutation frequencies of AHCY, CDKN2A and MPC1 in colorectal, gastric and esophageal cohorts; histograms depict the percentage of samples harbouring somatic mutations, amplifications or deep deletions. (D) Lollipop diagram summarising protein domains and recurrent somatic mutation hotspots of AHCY, CDKN2A and MPC1; the height of each lollipop corresponds to mutation frequency. (E-G) Bar charts of mutation frequencies for three prognostic genes. (H,I) Kaplan-Meier survival curves comparing overall survival of patients with versus without alterations in AHCY, CDKN2A or MPC1. CNA, copy number alteration.

CDKN2A is transcriptionally and translationally regulated in CRC

We first quantified MPC1, CDKN2A and AHCY mRNA in 3 freshly resected CRC specimens and matched adjacent normal mucosa. Consistent with our in-silico findings, CDKN2A expression was significantly higher in tumors, whereas MPC1 and AHCY were down-regulated (Figure 6A-6C). To corroborate the clinical relevance of CDKN2A our top candidate linking metabolic cell death to colorectal carcinogenesis we extended the analysis to an independent cohort of five pairs. RT-qPCR and Western blot revealed that both CDKN2A transcript and protein levels were markedly elevated in tumor tissue (Figure 6D-6F). Finally, we compared CDKN2A abundance between normal colonic epithelial cells (NCM460) and three CRC cell lines (SW480, SW620, and HCT116). All cancer-derived lines exhibited significantly higher CDKN2A mRNA and protein than NCM460 (Figure 6G-6I), recapitulating the patient-derived data and underscoring CDKN2A as a consistently activated metabolic-oncogenic node in CRC.

Figure 6 Transcriptional and translational profiling of AHCY, CDKN2A and MPC1 in colorectal cancer. (A-C) Comparing mRNA levels of MPC1 (A), CDKN2A (B) and AHCY (C) in tumor (T) versus matched adjacent normal tissues (N) (n=3). (D,E) Representative Western blot images of CDKN2A protein in five randomly selected T/N pairs. I Quantification of CDKN2A protein abundance normalized to GAPDH; band intensities were analysed by ImageJ and expressed as fold-change (T/N). (F) RT-qPCR validation of CDKN2A mRNA in the same five T/N pairs. (G) Basal CDKN2A mRNA expression in normal colonic epithelial cell line NCM460 and three CRC cell lines (SW480, SW620 and HCT116) determined by RT-qPCR. (H) Western blot showing CDKN2A protein levels in the four cell lines. (I) Densitometric quantification of (H) normalized to GAPDH. *, P<0.05; **, P<0.01; ***, P<0.001. CRC, colorectal cancer; RT-qPCR, reverse transcription quantitative real-time polymerase chain reaction.

CDKN2A silencing suppresses CRC aggressiveness and sensitises cells to cuproptosis

To dissect the functional role of CDKN2A in CRC, we transfected SW480 and HCT116 cells with three independent siRNAs (si-CDKN2A#1–3). All duplexes efficiently reduced CDKN2A mRNA and protein, with si-CDKN2A#3 producing the strongest knock-down (Figure 7A-7F); this reagent was used for all subsequent assays. EdU and CCK-8 experiments revealed that si-CDKN2A#3 markedly attenuated proliferation in both lines (Figure 7G-7L). Transwell and wound-healing assays further showed the reduction in invasion and migration after CDKN2A depletion (Figure 8A-8H). We next asked whether CDKN2A influences the newly characterized cuproprotein pathway. CuCl₂ alone modestly elevated the lipoylated forms of DLAT/DLST and the copper importer CTR1. Combined treatment with CuCl₂ and either si-CDKN2A#2 or si-CDKN2A#3 amplified these changes, indicating enhanced cuproptotic priming (Figure 8I-8L). Collectively, the data indicate that CDKN2A loss not only curbs CRC cell fitness but also lowers the threshold for CuCl₂-induced cuproptosis, identifying CDKN2A as a metabolic-oncogenic dependency that can be exploited through copper-based therapeutics.

Figure 7 Oncogenic dependency on CDKN2A in colorectal cancer cells. (A,B) RT-qPCR validation of CDKN2A knock-down efficiency 24 h after transfection with three independent siRNAs (siCDKN2A-#1, -#2, -#3) or non-targeting control (siNC) in SW480 and HCT116 cells. (C-F) Representative Western blots (C,E) and densitometric quantification (D,F) showing CDKN2A protein levels under the same conditions; GAPDH served as loading control. (G-J) EdU incorporation assays performed 72 h post-transfection with the most effective siRNA (siCDKN2A-#3); representative images (G,H) and quantitative percentages of EdU-positive cells (I,J) are presented (n=3). (K-L) CCK-8 proliferation curves of SW480 and HCT116 cells transfected with siCDKN2A-#3 or siNC; absorbance at 450 nm was measured daily for 72 h and normalized to day 0. *, P<0.05; **, P<0.01; ***, P<0.001. CCK-8, Cell Counting Kit-8; EdU, 5-ethynyl-2′-deoxyuridine; OD, optical density; RT-qPCR, reverse transcription quantitative real-time polymerase chain reaction.

Figure 8 Silencing CDKN2A impairs cell invasiveness and wound-healing capacity while enhancing CuCl₂-induced death in colorectal cancer cells. (A-D) Transwell invasion assays 24 h after siCDKN2A-#3 transfection; representative crystal-violet images (A,D) and quantification of invaded cells (B,C) for SW480 and HCT116 (n=3). (E-H) Scratch-wound migration assays; photomicrographs at 0,12, 24 and 48 h (E,H) and percentage wound closure (F,G) following CDKN2A knock-down. (I,J) HCT116 cells pre-treated with 5 nM CuCl₂ for 12 h after transfection with siCDKN2A-#3 or -#2; Western blot (I) and densitometry (J) of lipoylated DLAT/DLST, CTR1 and CDKN2A (GAPDH loading control). (K-L) Parallel experiment in SW480 cells; Western blot (K) and quantitative analysis (L) performed as in (I,J).

Discussion

MCD is no longer regarded as a passive by-product of oncogenic stress; instead, it represents a critical and targetable driver of CRC progression. Leveraging >1,100 transcriptomes, we constructed the first MCD-centric prognostic signature for CRC and identified CDKN2A, MPC1, and AHCY as three master genes that concurrently govern tumor-cell survival, immune escape, and drug sensitivity. Furthermore, we validated the expression of signature genes in clinical CRC tissues and functionally demonstrated that silencing the representative prognostic gene CDKN2A exerts tumor-suppressive effects by enhancing colon cancer cell sensitivity to cuproptosis, as evidenced by RT-qPCR, Western blot, EdU, and CCK-8 assays. Our findings not only delineate the mechanistic link between MCD dysregulation and adverse prognosis but also deliver a readily translatable biomarker platform for precision therapy in CRC.

CDKN2A is a canonical cell-cycle brake that has recently been re-classified as an anti-cuproptotic gene (32). Together with MTF1 and GLS, it governs the copper-dependent cell death pathway whereby Cu²⁺ binds lipoylated tricarboxylic-acid-cycle enzymes, causing their aggregation and proteotoxic stress (33). By restraining this cascade, CDKN2A is thought to shelter tumour cells from copper-induced lethality. Pharmacological degradation of CDKN2A with the copper ionophore elesclomol curtails proliferation of acute myeloid leukaemia cells, hinting at similar therapeutic value in CRC (34). Consistently, in silico data sets—including ours—link high CDKN2A expression to poor CRC prognosis (35). Extending these correlations, we demonstrate that siRNA-mediated CDKN2A silencing sensitises two CRC cell lines to cuproptosis and markedly impairs their viability, providing functional evidence that targeting CDKN2A could be exploited for copper-based precision therapy in CRC.

MPC1 and MPC2 assemble into a hetero-oligomeric channel that gates cytosolic pyruvate into the mitochondrial matrix (36). Loss of MPC1 compromises this gate, locking cells into a glycolytic state—an archetypal Warburg shift implicated in cancer and other disorders (37). In CRC, downregulation of MPC1 expression promotes histone lactylation and tumor invasiveness (38). Furthermore, inhibiting MPC1 expression confers ferroptosis tolerance via the KDM5A-MPC1 axis, thereby promoting head and neck cancer cell proliferation (39). Our large-scale transcriptomic and functional data now extend this paradigm by demonstrating that MPC1 down-regulation similarly impairs ferroptotic vulnerability in CRC, positioning MPC1 re-activation or combination with ferroptosis inducers as a novel therapeutic leverage.

AHCY is one of the most evolutionarily conserved proteins in eukaryotes and the only known mammalian enzyme that reversibly hydrolyzes S-adenosyl-L-homocysteine (SAH) into adenosine (Ado) and homocysteine (Hcy) (40). Initially characterized as a tumour suppressor, AHCY mRNA is frequently downregulated in multiple solid malignancies. Functional analyses revealed that this suppressive activity is context-dependent: AHCY silencing impairs invasion and migration in breast cancer (41), whereas reduced AHCY activity induces cell-cycle arrest and attenuates proliferation in hepatocellular carcinoma (42). Recent work further links AHCY to CRC by showing that its depletion suppresses proliferation, clonogenicity, invasion, and tumor angiogenesis while simultaneously dampening inflammation and oxidative stress (43). In non-small-cell lung cancer, AHCY inhibition enhances ferroptosis sensitivity (44). Collectively, these findings implicate AHCY as a metabolic-immunologic node that may modulate CRC progression via ferroptosis, meriting further therapeutic exploration.

The dismal and heterogeneous prognosis of CRC demands robust stratification tools to deliver precision therapy. Here, we built a three-MCD-gene risk signature and, for the first time, combined GSEA and GSVA to uncover that the high-risk group is specifically enriched for extracellular matrix (ECM)-receptor interaction and related micro-environmental pathways, offering actionable targets beyond canonical oncogenes. Immune-infiltration profiling further revealed a myeloid-biased, CD8⁺ T-cell-excluded landscape in high-risk tumours, providing a biological rationale for combining immune-checkpoint blockade with myeloid-targeting agents. Drug-sensitivity prediction based on oncoPredict and the GDSC2 database suggested that AZ960, AZD8186, and JAK inhibitors may represent potential therapeutic candidates for high-risk patients. Collectively, our model not only predicts patient outcomes with high accuracy but also provides a basis for the identification of risk-specific therapeutic strategies in CRC. Nevertheless, these findings are based on computational predictions and require further experimental and clinical validation.

Although this study conducted comprehensive bioinformatics analyses and experimental validations, several limitations must be acknowledged. First, our study was primarily based on retrospective transcriptomic datasets and in vitro experimental validation. Due to time and resource constraints, in vivo animal experiments were not performed. Therefore, the functional role of CDKN2A in tumor progression and cuproptosis regulation warrants further validation using appropriate animal models in future studies. Second, although we observed a significant association between the risk signature and tumor immune microenvironment features (including MDSC infiltration), these findings were derived solely from bioinformatics correlation analyses rather than mechanistic confirmation. Future studies will employ immunohistochemistry/immunofluorescence staining on clinical samples and MDSC co-culture systems to further elucidate the underlying mechanisms by which the CDKN2A/MPC1/AHCY three-gene signature remodels the tumor immune microenvironment and regulates MDSCs.

Conclusions

In summary, we developed and validated a robust, three-gene MCD-related signature that accurately forecasts prognosis and guides therapeutic decision-making in CRC. Mechanistic dissection identified CDKN2A as a metabolic-oncogenic dependency whose inhibition augments cuproptosis and cripples tumour aggressiveness. The signature provides an immediately translatable framework for risk-adapted treatment, including copper-based and immune-metabolic combination strategies. Prospective multicentre trials are warranted to confirm its clinical utility.

Acknowledgments

We thank all the patients who provided tissue samples for this study.

Footnote

Reporting Checklist: The authors have completed the MDAR and TRIPOD reporting checklists. Available at https://jgo.amegroups.com/article/view/10.21037/jgo-2026-0419/rc

Data Sharing Statement: Available at https://jgo.amegroups.com/article/view/10.21037/jgo-2026-0419/dss

Peer Review File: Available at https://jgo.amegroups.com/article/view/10.21037/jgo-2026-0419/prf

Funding: This work was supported by Guizhou Provincial Key Laboratory for Digestive System Diseases [No. ZSYS (2025)021].

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://jgo.amegroups.com/article/view/10.21037/jgo-2026-0419/coif). The authors have no 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. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Ethics Committee of The Second People’s Hospital of Guiyang (Jinyang Hospital) (No. JYYY-2025-XM-31). Written informed consent was obtained from all individual participants included in the study before sample collection.

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. Huo Z, Liu G, Li J. Recent research progress and clinical status of immunotherapy for colorectal cancer. J Adv Res 2026;82:729-50. [Crossref] [PubMed]
2. Chen Y, Zhang J, Ding Y, et al. Colorectal cancer pathogenesis, oncogenic signaling networks and targeted therapeutic advances. Mol Biomed 2026;7:32. [Crossref] [PubMed]
3. Siegel RL, Wagle NS, Star J, et al. Colorectal cancer statistics, 2026. CA Cancer J Clin 2026;76:e70067. [Crossref] [PubMed]
4. Siegel RL, Miller KD, Wagle NS, et al. Cancer statistics, 2023. CA Cancer J Clin 2023;73:17-48. [Crossref] [PubMed]
5. Yang L, Fang C, Zhang R, et al. Prognostic value of oxidative stress-related genes in colorectal cancer and its correlation with tumor immunity. BMC Genomics 2024;25:8. [Crossref] [PubMed]
6. Yalçıner M, Örüncü MB, Kayaalp M, et al. CEA and CA-19-9 Dynamics Associate with Survival in Regorafenib-Treated Metastatic Colorectal Cancer: A Real-World Analysis. J Clin Med 2026;15:2599. [Crossref] [PubMed]
7. Fan A, Wang B, Wang X, et al. Immunotherapy in colorectal cancer: current achievements and future perspective. Int J Biol Sci 2021;17:3837-49. [Crossref] [PubMed]
8. Mao C, Wang M, Zhuang L, et al. Metabolic cell death in cancer: ferroptosis, cuproptosis, disulfidptosis, and beyond. Protein Cell 2024;15:642-60. [Crossref] [PubMed]
9. Hao Y, Shao J, Lian N, et al. Metabolic cell death in cancer: mechanisms and therapeutic potential. Apoptosis 2025;30:2588-611. [Crossref] [PubMed]
10. Xue ZR, Xin YY, Jin WL. Exploiting metabolic vulnerabilities in cancer: From mechanisms to therapeutic opportunities. Cancer Lett 2025;634:218067. [Crossref] [PubMed]
11. Guo C, Liu P, Deng G, et al. Honokiol induces ferroptosis in colon cancer cells by regulating GPX4 activity. Am J Cancer Res 2021;11:3039-54.
12. Chaudhary N, Choudhary BS, Shah SG, et al. Lipocalin 2 expression promotes tumor progression and therapy resistance by inhibiting ferroptosis in colorectal cancer. Int J Cancer 2021;149:1495-511. [Crossref] [PubMed]
13. Xie J, Yang Y, Gao Y, et al. Cuproptosis: mechanisms and links with cancers. Mol Cancer 2023;22:46. [Crossref] [PubMed]
14. Xu K, Zhang Y, Yan Z, et al. Identification of disulfidptosis related subtypes, characterization of tumor microenvironment infiltration, and development of DRG prognostic prediction model in RCC, in which MSH3 is a key gene during disulfidptosis. Front Immunol 2023;14:1205250. [Crossref] [PubMed]
15. Gong X, Wu Q, Tan Z, et al. Identification and validation of cuproptosis and disulfidptosis related genes in colorectal cancer. Cell Signal 2024;119:111185. [Crossref] [PubMed]
16. Li J, Song Z, Chen Z, et al. Association Between Diverse Cell Death Patterns Related Gene Signature and Prognosis, Drug Sensitivity, and Immune Microenvironment in Glioblastoma. J Mol Neurosci 2024;74:10. [Crossref] [PubMed]
17. 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]
18. Gustavsson EK, Zhang D, Reynolds RH, et al. ggtranscript: an R package for the visualization and interpretation of transcript isoforms using ggplot2. Bioinformatics 2022;38:3844-6. [Crossref] [PubMed]
19. Wu C, Li Y, Luo Y, et al. Analysis of glutathione Stransferase mu class 5 gene methylation as a prognostic indicator in low-grade gliomas. Technol Health Care 2024;32:3925-42. [Crossref] [PubMed]
20. Chen H, Boutros PC. VennDiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R. BMC Bioinformatics 2011;12:35. [Crossref] [PubMed]
21. Wu T, Hu E, Xu S, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) 2021;2:100141. [Crossref] [PubMed]
22. Shannon P, Markiel A, Ozier O, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res 2003;13:2498-504. [Crossref] [PubMed]
23. Lei J, Qu T, Cha L, et al. Clinicopathological characteristics of pheochromocytoma/paraganglioma and screening of prognostic markers. J Surg Oncol 2023;128:510-8. [Crossref] [PubMed]
24. Friedman J, Hastie T, Tibshirani R. Regularization Paths for Generalized Linear Models via Coordinate Descent. J Stat Softw 2010;33:1-22.
25. Taylor JM. Random Survival Forests. J Thorac Oncol 2011;6:1974-5. [Crossref] [PubMed]
26. Navarro G, Gómez-Autet M, Morales P, et al. Homodimerization of CB(2) cannabinoid receptor triggered by a bivalent ligand enhances cellular signaling. Pharmacol Res 2024;208:107363. [Crossref] [PubMed]
27. 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]
28. Robles-Jimenez LE, Aranda-Aguirre E, Castelan-Ortega OA, et al. Worldwide Traceability of Antibiotic Residues from Livestock in Wastewater and Soil: A Systematic Review. Animals (Basel) 2021;12:60. [Crossref] [PubMed]
29. Liu ZY, Huang RH. Integrating single-cell RNA-sequencing and bulk RNA-sequencing data to explore the role of mitophagy-related genes in prostate cancer. Heliyon 2024;10:e30766. [Crossref] [PubMed]
30. Ye Y, Xu G. Construction of a new prognostic model for colorectal cancer based on bulk RNA-seq combined with The Cancer Genome Atlas data. Transl Cancer Res 2024;13:2704-20. [Crossref] [PubMed]
31. Maeser D, Gruener RF, Huang RS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform 2021;22:bbab260. [Crossref] [PubMed]
32. Liu H. Pan-cancer profiles of the cuproptosis gene set. Am J Cancer Res 2022;12:4074-81.
33. Li L, Zhou H, Zhang C. Cuproptosis in cancer: biological implications and therapeutic opportunities. Cell Mol Biol Lett 2024;29:91. [Crossref] [PubMed]
34. Du A, Yang Q, Luo X. Cuproptosis-related lncRNAs as potential biomarkers of AML prognosis and the role of lncRNA HAGLR/miR-326/CDKN2A regulatory axis in AML. Am J Cancer Res 2023;13:3921-40.
35. Jiang T, Wang Z, Sun Z, et al. Single-cell and Multi-omics Analysis Confirmed the Signature and Potential Targets of Cuproptosis in Colorectal Cancer. J Cancer 2025;16:1264-80. [Crossref] [PubMed]
36. Bricker DK, Taylor EB, Schell JC, et al. A mitochondrial pyruvate carrier required for pyruvate uptake in yeast, Drosophila, and humans. Science 2012;337:96-100. [Crossref] [PubMed]
37. Herzig S, Raemy E, Montessuit S, et al. Identification and functional expression of the mitochondrial pyruvate carrier. Science 2012;337:93-6. [Crossref] [PubMed]
38. Zhang B, Xu A, Wang H, et al. MPC-mediated lactate production drives histone lactylation in dendritic cells to affect tumor progression and immunotherapy. Cell Mol Life Sci 2025;82:371. [Crossref] [PubMed]
39. You JH, Lee J, Roh JL. Mitochondrial pyruvate carrier 1 regulates ferroptosis in drug-tolerant persister head and neck cancer cells via epithelial-mesenchymal transition. Cancer Lett 2021;507:40-54. [Crossref] [PubMed]
40. Vizán P, Di Croce L, Aranda S. Functional and Pathological Roles of AHCY. Front Cell Dev Biol 2021;9:654344. [Crossref] [PubMed]
41. Park SJ, Kong HK, Kim YS, et al. Inhibition of S-adenosylhomocysteine hydrolase decreases cell mobility and cell proliferation through cell cycle arrest. Am J Cancer Res 2015;5:2127-38.
42. Belužić L, Grbeša I, Belužić R, et al. Knock-down of AHCY and depletion of adenosine induces DNA damage and cell cycle arrest. Sci Rep 2018;8:14012. [Crossref] [PubMed]
43. Vande Voorde J, Steven RT, Najumudeen AK, et al. Metabolic profiling stratifies colorectal cancer and reveals adenosylhomocysteinase as a therapeutic target. Nat Metab 2023;5:1303-18. [Crossref] [PubMed]
44. Zheng H, Chen H, Cai Y, et al. Hydrogen sulfide-mediated persulfidation regulates homocysteine metabolism and enhances ferroptosis in non-small cell lung cancer. Mol Cell 2024;84:4016-4030.e6. [Crossref] [PubMed]