The Cancer Genome Atlas (TCGA) Colon and Rectal Cancer (COADREAD) dataset and GDC TCGA Colon Cancer (COAD) dataset was downloaded from the UCSC Xena server. TCGA, a landmark cancer genomics program based on American population, molecularly characterized over 20,000 primary cancer and matched normal samples spanning 33 cancer types. TCGA was built by National Cancer Institute and National Human Genome Research Institute. COADREAD in UCSC Xena server was derived from TCGA Data coordinating Center in January 2016 and had 434 samples. COAD data in UCSC Xena server was derived from TCGA data coordinating Center in August 2019 and had 512 samples. Gene expression in COADREAD and COAD was converted to log2 (1 + counts). COADREAD data were used as the training set in this study. After removing the same samples from training set, the remaining samples in COAD were used as the validation set.

Principal component analysis (PCA) methods (Trozzi et al., 2021) were used to evaluate gene expression patterns in cancer samples and normal samples in training set.

In the training set, multivariate Cox regression with forward stepwise selection was used to identify CRGs related to OS. The cuproptosis-related prognostic risk score (CPRS) of each sample was calculated using the formula CPRS = Σ Exp (mRNAί) × Coefficient (mRNAί). All tumor samples were grouped by the CPRS median and named the low-risk group and the high-risk group. Kaplan–Meier survival analysis and time-dependent receiver operating characteristic (ROC) curves between the low-risk and high-risk groups were used to evaluate the prognostic predictive performance in training sets. Furthermore, survival analysis and time-dependent ROC analysis were also used in the validation set (R package: timeROC). The independent predictive variables identified by multivariate Cox regression analyses were used to construct the predictive nomogram and the corresponding calibration curves (R package: rms) in training set.

Differentially expressed genes (DEGs) between the low-risk and high-risk groups were identified with the limma (version 3.40.6) package in training set (Ritchie et al., 2015). Genes that met the criteria of p < 0.05 and log2Fold change >1.5 were regarded as DEGs between the two groups.

GO and KEGG enrichment analyses of DEGs in training set were performed by clusterProfiler (version 3.14.3) (Yu et al., 2012). The minimum gene set was set as 5, and the maximum gene set was set as 5000. A p value <0.05 and an FDR <0.01 were considered statistically significant.

ESTIMATE analysis was employed to estimate the proportion of immune cells and stromal cells in the tumor tissue in training set. MCP-Counter analysis was employed to further estimate human immune and mesenchymal cell subsets according to the gene expression data in training set. The TIDE score and predicted response to immunotherapy were obtained via the TIDE database (http://tide.dfci.harvard.edu/) in training set.

Immunohistochemical images of CRGs in CRC and normal colorectal tissues were downloaded from the Human Protein Atlas (https://www.proteinatlas.org/). The H-score was used to quantify the staining.

All statistical analyses and visualizations were carried out by R version 4.0.2 (http://www.r-project.org). For detecting CRG differential expression, the Wilcoxon signed-rank test was used to estimate the differences between normal and tumor groups, and the Kruskal–Wallis test was used to compare more than two groups in different clinical characteristics. Also, ESTIMATE score, MCP-Counter score and TIDE score were compared between low and high risk groups by Wilcoxon signed-rank test. For identifying relationship between CRG and CRC prognosis, Kaplan-Meier analyses and the log-rank test were used to assess the survival differences between low and high levels of CRGs. Prognosis between high and low risk groups were compared by Kaplan-Meier analyses and the log-rank test. Spearman analysis was used to compute the correlation coefficients between immune score and DLAT/CDKN2A. p < 0.05 was regarded as statistically significant.

We obtained all CRGs including FDX1, LIAS, LIPT1, DLD, DLAT, PHDA1, PDHB, MTF1, GLS and CDKN2A from previous literature (Tsvetkov et al., 2022). We performed PCA on the TCGA CRC dataset, and one sample in the tumor group was removed (Supplementary Figure S1A). Subsequently, we detected the expression of these CRGs in tumor and normal samples and found that 6 genes, FDX1, LIAS, DLD, DLAT, PDHB and MTF1, were significantly downregulated in tumors, while 2 genes, GLS and CDKN2A, were significantly upregulated (Figure 1A). Immunohistochemical staining showed that among these CRGs with p < 0.05 in TCGA, their protein expression trend was consistent with the RNA expression in TCGA (Figure 1B).

In tumors, we found that compared with adenocarcinoma, PDHA1 and GLS were expressed at low levels in mucinous adenocarcinoma, while CDKN2A was highly expressed (Figure 2A). In addition, between the recurrent and nonrecurrent tumor samples, only LIAS expression was different (Supplementary Figure S1B). Meanwhile, there was no difference in any CRGs between M0 and M1 (Supplementary Figure S1C), indicating that these genes may not be related to distant metastasis of the tumor. In terms of lymphatic metastasis, we found that 5 genes, including LIAS, DLD, DLAT, PDHB and CDKN2A, had significantly different expression between tumors with and without lymphatic metastasis (Figure 2B). LIAS, DLD, DLAT and PDHB were expressed at low levels in the group with lymphatic metastasis, while CDKN2A was highly expressed. Regarding T stage, we found that the expression of most genes showed a positive or negative trend with the T stage (Figure 2C). Among them, DLD, PDHA1 and CDKN2A were significantly differentially expressed among different T stages. DLD and PDHA1 expression decreased gradually with T stage, while CDKN2A expression increased gradually. In addition, we observed a similar phenomenon for stage. The expression of LIAS, DLD, DLAT and CDKN2A decreased with stage (Figure 2), while the expression of CDKN2A increased. Different tumor pathological types have different molecular expression characteristics. Therefore, we divided CRC into adenocarcinoma and mucinous adenocarcinoma and then observed the expression of CRGs among different clinical features. The results showed that the expression of CRGs (Supplementary Figures S2A–C) in adenocarcinomas was similar to that before in different clinical features including lymphatic metastasis, T stage and stage (Figures 2B–D). However, there were no differences in different clinical features in mucinous adenocarcinoma (Supplementary Figures S2D–F). The above studies suggest that there is a close relationship between CRG expression and the clinical characteristics of CRC.

After we found that CRG expression was closely related to the clinical characteristics of CRC, we sought to determine whether CRG was related to prognosis. We grouped CRC patients according to the median expression value of each gene. The results showed that not all CRGs were related to prognosis. Three genes, PDX1, DLD1 and DLAT, were associated with OS (all p < 0.05) (Figure 3). The lower expression of these genes was associated with worse OS. Three genes, LIPT1, FDX1 and PDHA1, were related to DSS (all p < 0.05) (Figure 3). Higher LIPT1 expression was associated with worse DSS, while lower PDX1 and PDHA1 expression was associated with worse DSS. Additionally, 4 genes, LIAS, GLS, CDKN2A and LIPT1, were related to PFI (Figure 3). The lower expression of LIAS was associated with higher PFI, and the lower expression of GLS, CDKN2A and LIPT1 was associated with worse PFI. The above results suggest that there is a certain relationship between CRG expression and prognosis in CRC.

To further evaluate the relationship between CRG and prognosis, we first performed univariate Cox regression analysis of CRGs for OS. DLAT, CDKN2A, DLD, and PDHB were highly correlated with OS (Figure 4A). Then, multivariate Cox regression analysis was performed to further evaluate the relationship. The results showed that only DLAT and CDKN2A had a significant relationship with OS (Figure 4B). We believe that these two genes are most closely related to the prognosis of CRC. Subsequently, we used these two genes to construct a prognostic risk score using their regression coefficients: score = −0.641310371 × DLAT +0.130766422 × CDKN2A. Kaplan–Meier survival analysis showed that higher risk scores were associated with poorer OS (Figure 4C). With the increase in risk score, the survival rate decreased significantly. As a protective factor, DLAT expression decreased gradually as the risk score increased, and the expression level of CDKN2A increased gradually as a risk factor (Supplementary Figure S1D). The prediction accuracy of the area under the curve (AUC) assessment for 1 year, 3 years and 5 years were all more than 70% (Figure 4D). To further evaluate the prognostic risk model, we calculated the risk scores of each sample in the validation set. After the samples were divided into high-risk and low-risk groups according to the median score, the results were consistent with the previous results in the training set. Patients with high risk scores had a worse prognosis (Figure 4E). The ROC curve showed that the prediction accuracies at 1 year and 3 years were 0.65 and 0.62, respectively (Figure 4F). The above results suggested that the risk score we constructed is robust for predicting the prognosis of CRC patients.

Next, we conducted Cox multivariate analysis with forward stepwise selection of clinical characteristics and the risk score for OS. Age (young represented age<65 and older represented age≥65), T stage, pathological stage and risk score were identified as risk factors for the prognosis of CRC (Supplementary Table S1). Then, we built a nomogram for visualization (Figure 4G). C-index was calculated to be 3.59494602157775e-24 for OS, suggesting a relatively excellent predictive performance of the nomogram. Additionally, calibration plots demonstrated favorable concordance between the predicted OS and the observed OS at 1 and 3 years (Figure 4H).

To identify the potential differences between the high-risk and low-risk groups, we conducted differential expression analysis. The DEGs between the high-risk and low-risk groups were identified (Figure 5A), of which there were more upregulated DEGs than downregulated DEGs. The heatmap shows the top 50 DEGs with fold change >2. The expression of these genes was significantly different between the low- and high-risk groups. KEGG pathways (Figure 5B) were enriched in several classical tumor signaling pathways, including the Wnt signaling pathway, PI3K Akt signaling pathway, mTOR signaling pathway, etc. In addition, we found that the immune-related signaling pathways included the B-cell and T-cell receptor signaling pathways, Th17 differentiation, and the PDL1 and PDK1 checkpoint signaling pathways. We speculated that immune infiltration may be closely related to tumor prognosis.

Based on the above results, we used the ESTIMATE algorithm to evaluate the proportion of immune, stromal and tumor cells in tumor tissues (Figure 5C). The results showed that there was no difference in the immune score between the high-risk and low-risk groups, but the stromal score and ESTIMATE score were significantly higher in the high-risk group, indicating that high-risk tumors had more stromal cells and tumor cells. Next, we used the MCP-Counter algorithm to further evaluate immune cells and stromal cells (Figure 5D). The results showed that stromal cells, including endothelial cells and fibroblasts, were significantly more abundant in the high-risk group than in the low-risk group. In addition, most immune cells did not differ between the two groups, except that cytotoxic lymphocytes and monocytes were higher in the high-risk group than in the low-risk group. We analyzed the correlations between DLAT and CDKN2A and the above scores and found that DLAT was negatively correlated with the above scores, while CDKN2A was positively correlated with the above scores (Figure 5E). Both DLAT and CDKN2A are related to endothelial cells and fibroblasts. Additionally, DLAT was mainly related to CD8 T cells, cytotoxic lymphocytes and B cells, while CDKN2A was mainly related to T cells, cytotoxic cells, monocytes and neutrophils. Next, we analyzed the TIDE score between the high-risk and low-risk groups to evaluate the immune escape ability (Figure 5E). The results showed that the TIDE score in the high-risk group was significantly higher than that in the low-risk group, and we predicted the response to immunotherapy of the two groups (Figure 5F). The response rate in the high-risk group was significantly lower than that in the low-risk group. This may be attributed to the difference in immune escape ability.