mRNA expression profiles and corresponding clinical information from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/) COAD dataset were used as the discovery cohort in this study. Eligible samples met the following criteria: (1) histologically confirmed COAD; (2) mRNA expression profile data and clinical information were both available; and (3) only samples with a patient survival time $>$30 days were included in the survival analysis. Adjacent normal tissues (n = 40) and COAD samples (n = 404) were used to identify DEGs. The FPKM expression of each gene was applied with Log (1+FPKM) and normalized with Z-score (mean-centered). For external validation, gene expression matrices for the GSE91061, GSE78220, GSE74602, GSE17536, and GSE17537 series were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). GSE91061 and GSE78220 datasets have high-throughput RNA sequencing information. GSE91061 contains 109 RNA-seq samples (58 on-treatment and 51 pre-treatment) from 65 patients, while GSE74602 contains 30 matched pairs of normal and colorectal tumor samples. GSE78220 contains the RNA-seq transcriptome of responsive (n = 15) and non-responsive (n = 13) pretreated tumors with high-quality RNA (27 pretreatment and 1 early treatment). GSE17536 and GSE17537 contain microarray data from 177 and 55 colorectal cancer patient samples, respectively. The R package “sva” was used to eliminate batch effects that can occur because of the different data sources. To explore the metabolic landscape in colon cancer more comprehensively, genes from all metabolic KEGG pathways were extracted from the Molecular Signature Database (MsigDb) v7.1 (https://www.gsea-msigdb.org/gsea/msigdb/human/genesets.jsp?collection=CP:KEGG). A total of 38 metabolic pathways were analysed, including alanine aspartate and glutamate metabolism, alpha linolenic acid metabolism, amino sugar and nucleotide sugar metabolism, arachidonic acid metabolism, arginine and proline metabolism, ascorbate and aldarate metabolism, beta alanine metabolism, butanoate metabolism, cysteine and methionine metabolism. The 38 metabolic pathways and the genes within them are listed in Supplementary Table 1.

This study was approved by the Institutional Ethics Review Board of the Third People’s Hospital of Chengdu and conducted following the Chinese ethical guidelines for human genome/gene research.

The R package “limma” was used to perform differentiation analysis of metabolism-related gene expression by comparing tumor and normal samples. The thresholds for DEGs were $|$logFC$|$ $\ge$1 and adjusted p 0.05. Heatmaps were drawn using the “pheatmap” package and used to visualize differential gene expression.

Univariate analysis and least absolute shrinkage and selector operation (LASSO) Cox regression analysis were performed using the R packages “glmnet” and “survival”. These were used to identify metabolism-related prognostic gene signatures and included only those genes with p 0.01. The risk score for each patient in the training and validation cohorts was calculated separately according to the following formula:

Risk score = $\sum coefficient(mRN{A}_i)\times expressionvalue(mRN{A}_i)$

The median value of the risk scores was set as the cutoff point, and subsequently used to divide patients into high- and low-risk groups. Survival analysis was performed using Kaplan-Meier (K-M) survival curves. Discrimination between the outcomes for observations and predictions were performed using the area under the curve (AUC) of ROC in the package “survivalROC”. The AUC value ranged 0.5 to 1.0, with 0.5 indicating a random probability and 1.0 indicating perfect predictive ability.

KEGG analysis (p-value $\le$ 0.01) was performed using the KEGG automatic annotation web service (https://www.genome.jp/tools/kaas/). GSVA was performed using the R package “GSVA” [14]. All genes in the purine metabolism pathway were extracted from the KEGG database using the R package “KEGGREST” (version 1.35.0). Consensus clustering was performed using the R package “ConsensusClusterPlus”. Patients were divided into two clusters according to the expression of DEGs in the purine metabolism pathway. CIBERSORT was performed using the webserver tool (https://cibersort.stanford.edu/). Gene Set Enrichment Analysis (GSEA) was performed using GSEA software (http://www.broadinstitute.org/gsea). The Gene Expression Profiling Interactive Analysis (GEPIA) database (http://gepia.cancer-pku.cn/detail.php) was used to assess the expression level of genes highly related to immune checkpoints. Images of immunohistochemical staining for the selected core genes in normal tissue and COAD tissue were obtained from the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/) [15].

Tissues were fixed in 4% PFA at room temperature, embedded in paraffin and sectioned at 3 µm, then incubated with H${}_2$O${}_2$ to block endogenous peroxidases. Primary antibodies used for IHC analysis included rabbit anti-PDE1A (Proteintech, 12442-2-AP, Wuhan, China), rabbit anti-PDE5A (Proteintech, 22624-1-AP, Wuhan, China), rabbit anti-GUCY1A1 (Proteintech, 12605-1-AP, Wuhan, China), and rabbit anti-GUCY1B1 (Proteintech, 19011-1-AP, Wuhan, China). The IHC images were captured using an OLYMPUS v200 ASW microscope. All human participant experiments were conducted in compliance with the ethical standards of the Declaration of Helsinki. This study was approved by the Institutional Ethics Review Board of the Third People’s Hospital of Chengdu and conducted following the Chinese ethical guidelines for human genome/gene research.

The software package R 3.6.0 (R Foundation for Statistical Computing, Vienna, Austria) was used for statistical analyses. Groupwise comparisons were performed by Wilcox-ranked sum testing. Pearson correlation analysis was performed to identify metabolism-related prognostic genes (MRPGs) and potential immune-metabolic checkpoints. The overall survival (OS) rate of patients was estimated using K-M analysis, and the log-rank test was used to evaluate differences between groups. p 0.05 indicated statistical significance.

Data for 945 metabolism-related genes was extracted from the KEGG dataset for the differential expression analysis. Using the criteria of adjusted p 0.05 and $|$log2FC$|$ $\ge$1, a total of 129 MRDEGs were identified (Fig. 1A). Univariate Cox regression analysis was performed to identify MRDEGs with prognostic value, which identified 7 genes with p 0.01 (Fig. 1B). The Lasso regression algorithm with non-zero regression coefficients was used to select the most valuable prognostic genes. Using lasso regression analysis (lambda = 6), a prognostic model comprised of 6 MRDEGs was constructed (EPHX2, GPX3, PTGDS, NAT2, ACOX1, and CPT2) (Fig. 1C,D). To further validate the prognostic value of this 6-gene signature in COAD, the risk score for each patient in the TCGA and GEO databases was calculated using the formula: (exprGPX3 $\times$ 0.0083) + (exprPTGDS $\times$ 0.0269) – (exprEPHX2 $\times$ 0.0314) – (exprNAT2 $\times$ 0.0249) – (exprACOX1 $\times$ 0.0584) – (exprCPT2 $\times$ 0.067). All patients were then divided into high- or low-risk groups according to the median risk score (Fig. 1E). The log-rank test found that patients in the high-risk group of the TCGA cohort had significantly shorter overall survival (OS) (p = 0.0012; Fig. 1F). The AUCs for 3- and 5-year OS were 0.693 and 0.822 respectively (Fig. 1G,H), and the predictive ability of the signature was better than that of TNM stage.

Fig. 1.

Construction of a metabolism-related prognostic signature for COAD. (A) Heatmap showing the expression of MRDEGs. (B) Forest plots of univariate analysis for OS in patients from the TCGA database. (C) LASSO coefficient profile for the common genes. (D) Cross-validation for turning parameter screening in the LASSO regression model. (E) Kaplan-Meier survival curves for patients from the TCGA dataset stratified according to their risk score. (F) Distribution of the risk scores and survival times for TCGA patients. (G,H) Receiver operating characteristic (ROC) curve analyses for the prediction of 3- and 5-year OS according to the risk score and other clinical features. COAD, colon adenocarcinoma; MRDEGs, metabolism-related differentially expressed genes; TCGA, The Cancer Genome Atlas; LASSO, least absolute shrinkage and selector operation; OS, overall survival.

To determine whether the 6-gene signature was an independent prognostic factor for COAD patients, the expression of these genes was determined using publicly available transcriptome data (TCGA). The expression of EPHX2, GPX3, NAT2, PTGDS, ACOX1, and CPT2 were lower in COAD tumor tissue than in normal tissues (Supplementary Fig. 1A–F). K-M survival analysis showed that three genes were positively correlated with the OS of COAD patients (Supplementary Fig. 1G–I). Patients with high CPT2, EPHX2 and NAT2 expression levels showed better OS. Additional independent cancer datasets on immunotherapy for melanoma (GSE91061, GSE78220) were used to evaluate the predictive value of these genes for sensitivity to immunotherapy. Patients with higher expression of PTGDS were more sensitive to immunotherapy (Supplementary Fig. 1J). Patients with high CPT2 and NAT2 expression also showed significantly better response to immunotherapy than those with low expression (Supplementary Fig. 1J). The above results suggest that patients with high expression levels of PTGDS, CPT2 and NAT2 are more likely to benefit from immunotherapy.

We next performed univariate and multivariate Cox regression analysis. TNM stage and the risk score showed significant prognostic value for OS in univariate analysis (Fig. 2A). Significant factors were then entered into a multivariate Cox regression model. The 6-gene signature risk score and patient age were found to be independent predictors of OS (Fig. 2B).

Fig. 2.

Validation of the prognostic value of the MRDEG signature. (A) Univariate analysis for the OS of patients from the TCGA database. (B) Multivariate analysis for the OS of patients from the TCGA database. (C,D) Distribution of risk scores and survival times for GEO patients stratified according to their risk score. (E) K-M survival curves for patients in the GEO dataset stratified according to their risk score. GEO, Gene Expression Omnibus; K-M survival curves, Kaplan-Meier survival curves.

To validate the prognostic value of the signature, patients in the GSE17536 and GSE17537 series were separated into high- and low-risk groups according to the 6-gene signature risk score (Fig. 2C,D). Patients in the high-risk group had worse OS than those in the low-risk group (Fig. 2E), thus confirming the prognostic value of the 6-gene signature based MRDEGs.

We next investigated the metabolic pathways involved in determining the prognosis of COAD. KEGG enrichment analysis was performed on the 129 MRDEGs to identify reprogrammed metabolic pathways in COAD. Purine metabolism was found to be the most significantly enriched pathway (Fig. 3A). Other significantly altered pathways were drug metabolism, retinol metabolism, tyrosine metabolism, and pentose and glucuronate interconversion. Pearson correlation analysis was used to identify genes that were co-expressed with the 6-gene signature, using p 0.01 and $|$R$|$ $>$0.5 as the cut-off criterion. A total of 615 genes were retained and these were defined as metabolism-related prognostic genes (MRPGs). KEGG enrichment analysis was performed on the MRPGs to identify metabolic pathways that could play a role in determining the prognosis of COAD (Fig. 3B). The results showed that purine metabolism, biosynthesis of antibiotics, pyrimidine metabolism, carbon metabolism, fatty acid degradation, and pyruvate metabolism were associated with the prognosis of COAD patients. Interestingly, and consistent with the results shown in Fig. 3A, the purine metabolism pathway was the most enriched for metabolism-related differentially expressed genes s (p = 3.20 $\times$ 10${}^{-32}$) and also contained the most metabolism-related prognostic genes (p = 3.14 $\times$ 10${}^{-74}$). Therefore, its significance in COAD cannot be ignored.

Fig. 3.

Purine metabolism is the most significantly altered pathway related to the prognosis of COAD, and plays an important role in the development of COAD. (A) The top 12 enriched pathways identified from KEGG analysis of the MRDEGs. The bar color changes gradually from blue to red in ascending order according to the adjusted p-values. (B) The top 15 enriched pathways identified from KEGG analysis of the MRPGs. The bar color changes gradually from blue to red in ascending order according to the adjusted p-values. The size of the node represents the number of counts. (C) Heatmap for the expression of purine metabolism genes. (D–G) Distribution of GSVA enrichment scores for purine metabolism stratified according to TNM stage. (H–K) Associations between purine metabolism and the P53, PI3K/AKT/mTOR, cell-cycle, and MYC signaling pathways, as assessed by GSVA enrichment scores. These pathways were activated together with purine metabolism. KEGG, The Kyoto Encyclopedia of Genes and Genomes; GSVA, Gene Set Variation Analysis; TNM, tumor-node-metastasis.

GSVA was used to confirm that purine metabolism played an important role in COAD and to further examine the biological role of purine metabolism in COAD. The GSVA enrichment score was calculated for each patient, who were then divided into two groups according to the median score. The heatmap for the two groups showed a significant difference in the expression of genes for purine metabolism (Fig. 3C). The Wilcox test was performed to compare the GSVA enrichment score between patients with different clinical characteristics (Fig. 3D–G). A lower purine metabolism enrichment score was associated with higher tumor-node-metastasis (TNM) stage (p 0.05). Pearson correlation analysis was performed to explore the relationships between purine metabolism and other pathways in the KEGG and hallmark gene sets. The P53 signaling pathway, cell cycle, PI3K/AKT/mTOR signaling, and MYC targets were found to be positively correlated with purine metabolism, indicating these pathways were activated (Fig. 3H–K). The above results indicate that purine metabolism plays an essential role in COAD development.

In view of the importance of purine metabolism in COAD, we next performed consensus clustering. The k = 2 was identified by optimal clustering stability from k = 2 to 10 based on the similarity displayed by the expression levels of the MRDEGs and the proportion of ambiguous clustering measure. Patients were divided into cluster-1 (n = 177) and cluster-2 (n = 139) (Fig. 4A–C). The OS of cluster-1 patients was significantly better (p 0.01) than that of cluster-2 patients (Fig. 4D). The cluster subtypes defined by differentially expressed purine metabolism genes were closely related to the survival of COAD patients, thus further confirming our hypothesis on the important role of this pathway in COAD.

Fig. 4.

The expression of purine metabolism genes is closely related to the survival of COAD patients. (A) Consensus clustering cumulative distribution function for k = 2 to k = 10. (B) Relative change in the area under the CDF curve for k = 2 to k = 10. (C) Consensus clustering matrix for k = 2. (D) K-M survival curves for patients in the two clusters.

GSEA analysis was performed to compare the biological features between the two subtypes. Samples in cluster-2 were more enriched in immune-related pathways, including the immune network for IgA production, natural killer cell-mediated cytotoxicity, antigen processing and presentation, graft versus host disease, and cytokine-receptor interaction (Fig. 5A). These results suggest that purine metabolism may have some interaction with tumor immunity.

Fig. 5.

Consensus clustering of purine metabolism genes with anti-tumor immunity. (A) GSEA showed that samples in cluster-2 were more likely to be enriched in immune-related pathways. (B–K) Expression levels for PD1, PD-L1, PD-L2, LGALS9, HAVCR2, CD86, CD80, CTLA4, TNFRSF9 and TNFSF9 in cluster-1 and cluster-2 subtypes from the TCGA cohort. (L) Infiltrating levels of 22 immune cell types in cluster-1 and cluster-2 subtypes from the TCGA cohort. * p 0.05, ** p 0.01 and *** p 0.001. GSEA, Gene Set Enrichment Analysis.

The expression of immune checkpoints are closely related to the efficacy of immunotherapy. For example, the overexpression of PD-L1 is widely used as a predictive biomarker for the response to immunotherapy with checkpoint inhibitors. The expression of immune checkpoints was therefore investigated in the two subtypes. Patients in cluster-2 expressed higher levels of PD1, PD-L1, PD-L2, LGALS9, HAVCR2, CD86, CD80, CTLA4, TNFRSF9 and TNFSF9 than patients in cluster-1 (Fig. 5B–K), indicating they may respond better to immunotherapy.

The distribution of 22 immune cell types in the two subgroups was compared using the CIBERSORT algorithm. Cluster-1 showed higher levels of infiltration with monocytes and CD4 memory-resting T cells, whereas cluster-2 showed more CD8 T-cells, M1 macrophages, and follicular helper T cells (Fig. 5L).

Purine metabolism status had closely related to the expression of immune checkpoints [16, 17]. Therefore, targeting of the purine metabolism pathway may be a potential strategy for improving the efficacy of immune blockade therapy. Pearson correlation analysis was used to compare the expression of genes in the purine metabolism pathway with immune checkpoints. With a cut-off point of $|$R$|$ $>$0.5, we found 8 genes that were highly correlated with the expression of immune checkpoints. The expression of four of these genes (GUCY1A1, GUCY1B1, PDE1A and PDE5A) were significantly lower in COAD than normal tissues (Fig. 6A–D). Moreover, a similar expression pattern for these genes was seen in an independent patient cohort (GSE74602) (Supplementary Fig. 2A). IHC staining also revealed the expression of GUCY1A1, GUCY1B1, PDE1A and PDE5A were downregulated in COAD tissues (Fig. 6E, Supplementary Fig. 2B). As shown in Fig. 6F–K, PD-L2 expression was positively correlated with the expression of GUCY1A1 (r = 0.55, p = 8.52 $\times$ 10${}^{-34}$), GUCY1B1 (r = 0.55, p = 2.41 $\times$ 10${}^{-33}$) and PDE1A (r = 0.52, p = 6.05 $\times$ 10${}^{-30}$). The expression of CD80 was positively correlated to that of GUCY1A1 (r = 0.55, p = 1.44 $\times$ 10${}^{-33}$) and PDE5A (r = 0.56, p = 3.63 $\times$ 10${}^{-35}$), while the expression of CD86 was positively correlated with GUCY1B1 expression (r = 0.51, p = 5.24 $\times$ 10${}^{-28}$). These findings suggest that targeting of GUCY1A1, GUCY1B1, PDE1A and PDE5A may improve the effectiveness of immunotherapy, especially in patients with higher PDE1A (Supplementary Fig. 2C).

Fig. 6.

Identification of potential immuno-metabolic checkpoints in COAD. (A–D) Expression levels for GUCY1A1, GUCY1B1, PDE1A and PDE5A in tumor and normal tissues from the GEPIA database. (E) IHC staining of GUCY1A1, GUCY1B1, PDE1A and PDE5A was significantly lower in COAD tissues compared to normal tissue. Scale bars represent 100 µm. (F,G) Correlation between the expression of immune checkpoints and GUCY1A1 in the TCGA database. (H,I) Correlation between the expression of immune checkpoints and GUCY1B1 in the TCGA database. (J) Correlation between the expression of PD-L2 and PDE1A in the TCGA database. (K) Correlation between the expression of CD80 and PDE5A in the TCGA database.

The heatmap shown in Supplementary Fig. 2D provides an overview of the correlation coefficients between the signature genes and the immune checkpoints. We found no significant difference in the expression of PDE6C, NT5C1B, ENTPD1 and AMPD3 between COAD and normal tissues (Supplementary Fig. 3A–G). However, the expression levels of PDE6C, NT5C1B, ENTPD1 and AMPD3 were also positively correlated to the expression of immune checkpoints (Supplementary Fig. 3H–M).