Prognostic implications of metabolism-associated gene signatures in colorectal cancer

Data sources

Both the RNA-sequencing datasets and clinical data of CRC were obtained from the TCGA database (Data Release 22.0—January 16, 2020, https://tcga- data.nci.nih.gov/tcga/). The two main filter criteria for our data were as follows: (1) The keywords of cases are “colon, rectosigmoid junction, and rectum [Primary Site]”, “TCGA[Program]”, “TCGA-COAD, TCGA-READ [Project]”, “Adenomas and Adenocarcinomas [Disease Type]”.(2) The keywords of files are “Transcriptome Profiling [Data Category]”, “Gene Expression Quantification [Data Category]”, “RNA-Seq [Experimental Strategy]”, “HTSeq—FPKM [Workflow Type]”, “Clinical [Data Category]”, “BCR XML [Data Format]”. The matrix files of RNA-sequencing for different samples were collated and annotated onto the genome. The expression of mRNA was extracted from the matrix file obtained from the RNA-sequencing data. The data access used GSE39582 and GSE87211 as the validation data obtained from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/). Raw data of gene chips are normalized using the RMA algorithm provided by “limma” (Ritchie et al., 2015). Perl and R-package “sva” were used to merge microarray data and reduce heterogeneity between the two studies. Metabolism-related genes were obtained from KEGG symbols of the Gene Set Enrichment Analysis (GSEA) website (https://www.gsea-msigdb.org). R software (version 3.6.2) was used for data annotation and the extraction of metabolic gene expression for the TCGA and GEO data.

Identification of prognosis-associated metabolism-related genes

Data extraction and integration were performed using Perl software. We screened differentially expressed metabolism-related genes using the Wilcox Test with R-package “edgeR” (Robinson, McCarthy & Smyth, 2010), and “limma”. |Log Fold Change| > 1.5 and False Discovery Rate (FDR) < 0.05 were set as the cutoff. Bidirectional hierarchical clustering was analyzed and we drew a heatmap using R-package “pheatmap” (https://cran.r-project.org/web/packages/pheatmap/). R-package “ggplot2” was used to draw a volcano map. Prognosis-associated metabolism-related genes were identified by univariate Cox regression analysis.

Construction prognosis model of metabolism-related genes and survival analysis

We screened out the metabolism-related genes that had a significant correlation with overall survival (OS) of colorectal cancer cohorts using univariate Cox regression analysis (p < 0.01). The least absolute shrinkage and selection operator (Lasso) Cox regression analysis constructed the optimal model of metabolism-related genes using the R-package “glmnet” (Friedman, Hastie & Tibshirani, 2010). Genes expression and survival analysis were evaluated using the Kaplan–Meier method and the log-rank test (p < 0.05). The risk score for each patient calculated as: Risk score $={\sum }_{j=1}^{n}Coefj\ast Xj$, Coefj denotes the coefficient and Xj denotes the expression levels of each metabolism-related gene (Liu et al., 2019). Survival data were selected from each sample obtained from the clinical data downloaded from the TCGA and combined with the previously acquired expression profiling data. The median risk score was selected as a cutoff value to create colorectal cancer cohorts. The survival curve was drawn according to the high and low-risk value by R-package “survival” “survminer”. The R-package “survival ROC” drew the Receiver Operating Characteristic (ROC) curve, which was used to investigate the sensitivity and specificity of the survival prediction by the gene marker risk score (Le, Yapp & Yeh, 2019). Area Under Curve (AUC) served as an indicator of prognostic veracity (Le, 2019; Sachs, 2017). The risk curve, survival state diagram, and heatmap were drawn based on the different risk scores of the patients. Independent prognostic metabolic genes were recognized using univariate and multivariate Cox proportional risk regression analysis. We conducted clinical correlation an analysis using the R-package “beeswarm”. The prognostic metabolism-related gene was verified in the independent colorectal cancer cohort of GEO (GSE39582, GSE87211) (Hu et al., 2018; Marisa et al., 2013).

Functional enrichment analysis of metabolism-related genes

Functional enrichment analysis of metabolic genes was conducted based on the DAVID database (https://david.ncifcrf.gov/home.jsp) (Huang, Sherman & Lempicki, 2009), which identified Gene Ontology (GO) categories in the biological processes (BP), cellular components (CC), or molecular functions (MF). The DAVID database was also used to determine the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. FDR < 0.05 was the cutoff value. R-package “GOplot” and “ggplot2” were used to integrate expression data and functional and pathway enrichment analysis. Gene Set Enrichment Analysis (GSEA) was also used to uncover the signaling pathways and biological processes of differentially expressed genes between the high and low-risk subgroups. The number of permutations was set to 1,000, and FDR < 0.25 was considered to be statistically significant.

Exploitation of the nomogram

Age, gender, stage, TNM stage, and risk score were used to draw a nomogram using the R-packages “Hmisc”, “lattice”, “Formula”, “ggplot2”, “foreign” and “rms”. Calibration traces were used to assess the consistency between the actual and predicted survival rates. The Consistency index (C-index), ranging from 0.5 to 1.0, was calculated to evaluate the model’s capability for predicting an accurate prognosis. Measurements of 0.5 and 1.0 from the model represent a random probability or an excellent performance for predicting survival, respectively.

Statistical analysis was performed using the R software (Version 3.6.2; https://www.R-project.org). P-values for all of the analyses were: less than 0.05 (statistically significant), 0.01 (more statistically significant), and 0.001 (most statistically significant).

Characteristics of patients

The TCGA CRC cohort consisted of 624 patients. Patients who have not survival time and TNM stage were excluded, resulting in a total of 596 patients. Patients were further screened, and those whose survival time was less than 90 days were also excluded, resulting in a total of 545 patients. A total of 692 mRNA expression profiles of CRC were downloaded from TCGA. Among them, 51 (7.4%) derive from healthy samples, and 641 (92.6%) come from tumor samples.

The Marisa cohort consisted of 585 patients (GSE39582), including 566 with stage I-IV colon adenocarcinoma and 19 colon mucosa samples (Marisa et al., 2013). The Hu cohort included 363 patients (GSE87211), consisting of 203 rectal tumors and 160 rectal mucosa samples (Hu et al., 2018). A total of 948 patients were obtained after merging the two sets of data. Patients with incomplete information were excluded, resulting in a total of 720 tumor patients. Patients were further screened and those whose survival time was less than 90 days were also excluded, resulting in a total of 701 tumor patients. The demographic and clinical characteristics of patients are listed in Table S1.

Construction prognosis model for TCGA colorectal cancer cohort

A total of 102 differentially expressed metabolism-associated genes were screened containing 40 upregulated and 62 down-regulated genes. Gene expression profiles of 102 genes were displayed separately for colorectal cancer cohorts in normal and tumor tissues in Fig. S1. Univariate Cox regression analysis revealed seven metabolism-related genes that were significantly related to the OS of CRC. NAT2 and XDH were low-risk PRMGs (Hazard Ratio (HR) < 1, p < 0.05); GPX3, AOC3, AKR1C4, SPHK1, and ADCY5 were high-risk PRMGs (HR > 1, p < 0.05) (Fig. 2A).

Six PRMGs were filtered out by LASSO COX regression analysis and the coefficient was calculated. The model provided the best figure when six genes were included (Figs. 2B and 2C). The coefficient of each CRC gene is shown in Fig. 2D. The prognosis model was constructed based on these six genes (NAT2, XDH, GPX3, AKR1C4, SPHK1, ADCY5). The full names, pathways, and coefficients of these genes are shown in Table 1.

We investigated genetic alterations in the role of CRC risk-related metabolic genes in CRC using Gene Set Cancer Analysis (GSCA) Lite (http://bioinfo.life.hust.edu.cn/web/GSCALite/). The genes of interest in CRC were changed in 38 of 38 queried samples (100%) (Fig. 2E). The frequent genetic changes indicated the crucial roles of these genes in the development of CRC.

Each patient’s risk score was calculated based on the mRNA expression level and risk coefficient of each gene. The calculation of the risk score was: risk score = 0.0283 × expression of AKR1C4 + 0.0075 × expression of SPHK1 + 0.0110 × expression of GPX3 + (−0.0631) × expression of NAT2 + (−0.0293) × expression of XDH + 0.1465 × expression of ADCY5. The risk score was used to predict prognosis. Patients were divided into a high-risk group and a low-risk group with the median risk score as the cut-off value. A heatmap was developed to show the gene expression profiles of the high- and low-risk TCGA-CRC groups (Fig. 3A) and the GEO-CRC group (Fig. 3B). The risk score had a significant correlation with age, T, N, M, and clinical stage in TCGA-CRC, and age, gender, T, N, and clinical stage in GEO-CRC. The high-risk genes (GPX3, AKR1C4, SPHK1, ADCY5) were more likely to be expressed in the high-risk group patients, while low-risk genes (NAT2 and XDH) were expressed in the low-risk group patients of TCGA-CRC (Fig. 3A) and GEO-CRC (Fig. 3B). The association between survival time and risk score is shown in Fig. 3C (TCGA-CRC), 3D (GEO-CRC). The survival time decreased as risk score increased and the higher the risk score, the more deaths in TCGA-CRC (Fig. 3E) and GEO-CRC (Fig. 3F). Risk scores both were significantly related to OS in the univariate independent prognostic analysis of TCGA-CRC (HR = 3.324, 95% CI [1.956–5.650], P < 0.001) (Fig. 3G) and GEO-CRC (HR = 2.096, 95% CI [1.383–3.177], P < 0.001) (Fig. 3H).

The multivariate independent prognostic analysis showed that the risk score was an independent prognostic predictor in the TCGA-CRC group (HR = 2.639, 95% CI [1.413–4.928], P = 0.002) (Fig. 4A) and the GEO-CRC group (HR = 1.658, 95% CI [1.059–2.598], P = 0.027) (Fig. 4B). Kaplan–Meier cumulative curves indicated that the OS rate between the high-risk group and the low-risk group was significantly different. The survival time of patients with high-risk score was significantly shorter than that of patients with low-risk score in the TCGA-CRC group (P < 0.001) (Fig. 4C) and GEO-CRC group (P < 0.001) (Fig. 4D). The 1-, 3-, and 5-year AUC values of the risk score were 0.672, 0.608, and 0.648, respectively, and the prognostic accuracy of the stage was higher than other clinical characteristics in TCGA-CRC (Figs. 4E, 4G and 4I). Similarly, the 1-, 3-, and 5-year AUC values of the risk score were 0.612, 0.566, and 0.573, respectively, and the prognostic accuracy of the M-stage was higher than other clinical characteristics in GEO-CRC (Figs. 4F, 4H and 4J). These results confirm that six metabolic gene signatures can also predict survival in the independently validated GEO colorectal cohort.

The risk score and metabolic genes related to the clinicopathological features of CRC

The expressions of NAT2, ADCY5, SPHK1, GPX3, and the risk score were significantly correlated with the clinicopathological features of TCGA-CRC (Table S2, Figs. 5A–5K). NAT2 was highly expressed in stages I, II, N₀, and M₀, but had low expression in stages III, IV, N_1–2, and M₁ (Figs. 5A–5C) (p < 0.05). However, the expression level of ADCY5 was higher in stages III, IV, and N_1–2 than in stages I, II, and N₀ (Figs. 5D and 5E). SPHK1 and GPX3 were also highly expressed in T_3–4 but had low expressions in T_1–2 (Figs. 5F and 5G). The risk score was significantly correlated with the clinical stage, which was highly expressed in stages III, IV, T_3–4, N_1–2, and M₁, but had low expression in stages I, II, T_1–2, N₀, and M₀, (Figs. 5H–5K) (p < 0.05). These results showed that NAT2 was a low-risk metabolic gene, and ADCY5, SPHK1, GPX3, were high-risk metabolic genes, which is consistent with our previous results. More importantly, the risk score closely correlated with the malignant clinicopathological characteristics of CRC and is an independent prognostic factor.

Similar results were seen in the validation group (Table S3, Figs. 5L–5AA). NAT2 and XDH were highly expressed in the <65-year age group, M₀, and had low expression in the >65-year age group, M₁ (Figs. 5L–5O) (p < 0.05). The expression level of ADCY5 was the opposite of NAT2 and XDH (Figs. 5P and 5Q). SPHK1 and GPX3 were highly expressed in stage III, IV, T_3–4, and N_1–3, and had low expressions in stages I, II, T_1–2, and N₀ (Figs. 5R–5W). The expression of AKR1C4 was higher in stages III, IV, and N_1–3, than in stages I, II, and N₀ (Figs. 5X–5Y). The risk score was higher in the group of >65 years old and M₁, lower in the group of <65 years old and M₀ (Figs. 5Z and 5AA) (p < 0.05).

Identification involved signaling pathways of PRMGs

GO enrichment analysis of the seven metabolic genes showed that there were two GO terms in BP, one GO term in MF, and one GO term in CC, which was significant (FDR < 0.05). GO enrichment analysis showed that these seven genes could be classified into several essential processes, including oxidation–reduction, the xenobiotic metabolic process, electron carrier activity, and cytosol (Fig. 6A). KEGG pathway enrichment shown that PRMGs were mainly enriched in caffeine metabolism, metabolic pathways, and drug metabolism by other enzymes (Fig. 6B) (p < 0.05). GSEA analysis showed that changed genes were observably enriched in several common pathways. 98 of 178 gene sets were upregulated in the high-risk phenotype group, and 95 gene sets were remarkable at FDR < 25%. The top-five gene sets of the high-risk group were significantly related to basal cell carcinoma (NES = 2.31, P = 0.000), dilated cardiomyopathy (NES = 2.28, P = 0.000), vascular smooth muscle contraction (NES = 2.25, P = 0.000), glycosaminoglycan biosynthesis chondroitin sulfate (NES = 2.24, P = 0.000), and axon guidance (NES = 2.23, P = 0.000). 80 of 178 gene sets were upregulated in the low-risk phenotype group. 74 gene sets were markedly enriched at FDR < 25%, and the five most common gene sets of the low-risk group were negatively associated with peroxisome (NES = −2.37, P = 0.000), fatty acid metabolism (NES = −2.36, P = 0.000), butanoate metabolism (NES = −2.31, P = 0.000), valine, leucine, and isoleucine degradation (NES = −2.26, P = 0.000), and propanoate metabolism (NES = −2.26, P = 0.000) (Fig. 6C). The integration of the 5 most prevalent phenotypes of the high-risk and the low-risk groups was visualized in Fig. 6D.

A characterized prognostic prediction model

A nomogram is a powerful tool quantifying an individual’s risk in a clinical setting by integrating multiple risk factors (Liang et al., 2015; Won et al., 2015). We used the nomogram to predict the probabilities of 1-, 3- and 5-year OS by incorporating the age, gender, TNM stage, and risk score of TCGA-CRC (Fig. 7A). The results were verified in the GEO-CRC. Each factor was assigned a score in proportion to its contribution to the risk of survival. The calibration curve showed that the actual survival time is in agreement with the predicted survival time and the C-index is 0.8. (Figs. 7B, 7D and 7F). A 75-year old (60 points) female patients (6 points) would acquire a total of 224 points if she had stage III (41points), T3 (32 points), N1 (0 points), and M0 (0 points), with a risk score (85 points). Her 1-, 3-, and 5-year survival rate was approximately 67%, 45%, and 28%, respectively. The nomogram in the GEO colorectal cancer cohorts and the 1-, 3- and 5- year calibration curves are shown in Figs. 7C, 7E and 7G, respectively.