Transcriptome analysis reveals a reprogramming energy metabolism-related signature to improve prognosis in colon cancer

Colon cancer datasets preparation

Three independent and publicly available gene expression microarray datasets (GSE39582 (Marisa et al., 2013), GSE17538 (Freeman et al., 2012; Smith et al., 2010; Williams et al., 2015a), and GSE41258 (Martin et al., 2018; Sheffer et al., 2009)), which also contained patient outcome and clinicopathological features, were retrieved from the Gene Expression Omnibus (GEO). In total, 879 colon cancer patients were included. For the dataset GSE39582, we excluded 86 samples which received treatment but not without corresponding drug information or did not contain the exact treatment information. The 316 patients who didn’t receive adjuvant chemotherapy and the 19 normal subjects from GSE39582 were treated as training set, and the other two datasets (GSE17538 and GSE41258) were served as independent validation cohorts. In the training cohort, 3 tumor samples that had no available clinical survival information and the 19 normal samples were excluded in survival analysis. For the two validation cohorts, we used the similar strategies to select samples. Besides, 164 colon cancer patients who received fluorouracil-based ACT in GSE39582 were further used for stratification analysis.

Identification of reprogramming energy metabolism-related genes associated with prognosis in colon cancer

Differential expression analyses to detect transcriptome differences between tumor and normal samples were performed by the “limma” R package (3.42.2) in GSE39582. We assigned the false discovery rate (FDR) as 0.1% and fold change (FC) threshold as 2 to screen differentially expressed genes (DEGs). Among DEGs, we focused on genes involved in reprogramming energy metabolism as well as genes related to malignant tissue development. The reprogramming energy metabolism-related genes were identified by hallmark-associated GO terms in our previous research (Deng et al., 2018). Then, we performed the univariable Cox regression analysis to select genes associated with overall survival in the training cohort. Finally, a risk score model was established: Risk score = (0.4911 × expression level of AARS) + (0.2992 × expression level of COL4A2) + (0.2865 × expression level of COL4A1) + (−0.3543 × expression level of PCNA) + (−0.2377 × expression level of MMP12). The formula was a linear combination of the five genes weighted by its corresponding regression coefficient in the univariable Cox regression analysis. Based on this model, each sample could be assigned an assessment score and then assigned to the high-risk or low-risk group according to the median score in the corresponding cohort.

Statistical analysis

The Kaplan–Meier method was used to assess the survival time in each cohort (Dudley, Wickham & Coombs, 2016). The difference of overall survival between the high-risk group and the low-risk group in each cohort was evaluated by the log-rank test. We applied the univariate Cox regression model to assess the prognostic value of the signature and performed the multivariate Cox regression models to estimate whether it was an independent predictor of overall survival in colon cancer (Zhang et al., 2016b). The Cox proportional hazards regression model was used to calculated hazard ratio (HR) and 95% confidence intervals (CI). The statistically significant result was defined as the two-tailed P-value being less than 0.05. We further performed ROC analysis to assess the sensitivity and specificity of the signature in overall survival prediction. The prediction models of the five-gene risk score and the TNM stage were also assessed, respectively. And then, the two variables were regarded as predictors and combined in a logistic regression model. The pROC package was used to plot ROC curves of the predictions (Robin et al., 2011), and the areas under the curve (AUC) of these models were computed to compare their abilities. All the statistical analyses were conducted under the R program 3.5.2 (http://www.r-project.org).

Transcriptome analysis reveals a prognostic signature related to reprogramming energy metabolism

Transcriptional profiles of colon cancer patients and corresponding clinical information were obtained from Marisa et al. (2013) After removing patients who received fluorouracil-based ACT or without available clinical survival information, a total of 316 tumor samples and 19 normal samples were recruited in the training cohort (Table 1). Based on the transcriptional profile, we identified 1,489 differentially expressed genes in colon cancer patients (FDR < 0.001 and FC > 2 by “limma” method), including 752 up-regulated and 737 down-regulated genes. Then, by subjecting the differentially expressed genes to reprogramming energy metabolism as well as malignant tissue development, we obtained 43 metabolism-related genes for subsequent analysis (Table S1). To extract a core prognostic gene set related to reprogramming energy metabolism, we applied univariate Cox proportional hazards regression analysis. We obtained a five-gene signature that was significantly related to overall survival (Table 2, P < 0.05). Among the five genes, high expression levels of two genes (PCNA and MMP12) were strongly associated with better overall survival, and the other three genes (COL4A1, COL4A2, and AARS) showed opposite correlations. Based on the risk score formula of the five genes (see Materials and Methods), each sample assigned a risk score (Fig. S1). Then, samples could be divided into high-risk (n = 156) and low-risk (n = 157) groups according to the median score 3.71 in the training cohort in which 3 tumor samples without available clinical survival information were excluded. For the high-risk group, COL4A1, COL4A2, and AARS showed relatively higher expression and PCNA and MMP12 showed relatively lower expression than the low-risk group. Moreover, the high-risk samples presented a trend of short survival time (Fig. S2).

Validation of the prognostic value of the reprogramming energy metabolism-related signature

To evaluate the prognostic value of the signature, we performed Kaplan–Meier analysis to evaluate the survival outcomes of patients with colon cancer. Our result showed that the overall survival time of the patients in the high-risk group was significantly shorter than that in the low-risk group (Fig. 1A, log-rank test, P = 0.0011). By using univariate Cox regression analysis, we also found that patients with higher five-gene risk scores had shorter overall survival time (hazard ratio (HR) = 1.92, 95% CI [1.29–2.87]).

To further validate the prognostic ability of the five-gene signature, we collected another two independent validation cohorts, i.e., GSE17538 (n = 232) and GSE41258 (n = 167). Consistent with the results in the training cohort, patients classified as “high-risk” suffered significantly poorer overall survival than those defined as “low-risk” in both two independent validation cohorts (Figs. 1B–1C, P = 5.2e−04 and 3.8e−04, respectively). Through the univariate Cox regression model analysis, the correlation between the signature and the overall survival was consistently strong in the two cohorts (HR = 2.09, 95% CI [1.37–3.2] for GSE17538 and HR = 2.08, 95% CI [1.37–3.16] for GSE41258). In addition, we collected 373 colon cancer samples containing RNA-seq expression and clinical information from The Cancer Genome Atlas (TCGA) (https://cancergenome.nih.gov/). Our signature was still able to significantly predict overall survival of colon cancer patients in TCGA (Fig. S3A, P = 0.038, HR = 1.57, 95% CI [1.02–2.4]). We further explored our risk model in other gastro-intestinal tumors from TCGA but not of colonic origin. We found that the risk model did not have the ability to predict patient survival in other gastro-intestinal tumors (Figs. S3B–S3C, P = 0.8, HR = 0.85, 95% CI [0.23–3.17] for pancreatic cancer and P = 0.096, HR = 0.69, 95% CI [0.45–1.07] for liver cancer). Our signature might be specific for colon tissue.

The five-gene signature predicts patient survival independently

To further test whether the prognostic value of the five-gene signature was independent, we performed the multivariate Cox regression model analysis according to clinicopathological parameters of the colon cancer patients, including age, sex, and tumor stage. Our results proved that the signature remained to have prognostic power of overall survival when considering those factors in the training cohort (Table 3, P = 0.009, HR = 1.72, 95% CI [1.15–2.59]). The multivariate Cox regression model analysis based on the two validation cohorts also showed that the risk score was significant among all clinicopathological parameters and significantly associated with overall survival (Table 3, P = 0.006, HR = 1.84, 95% CI [1.19–2.82] for GSE17538 and P = 0.01, HR = 1.76, 95% CI [1.14–2.71] for GSE41258, respectively). While for other factors, only the TNM stage was independent in all cohorts (Table 3). In addition, we further explored the connection between the signature and the TNM stages. The risk scores had a positive Spearman’s correlation with the TNM stages (Fig. S4, Spearman’s correlation, 0.26 for the training cohort; 0.24 for the validation cohort GSE17538; 0.271 for the validation cohort GSE41258).

Prognostic prediction of the signature within clinical stages

The disease-specific survival of colon cancer is impacted by several factors, including the tumor–node–metastasis (TNM) staging. The stage of disease at diagnosis has been proved to be a major predictor to help evaluate the prognosis (Crooke et al., 2018). In the multivariate survival analysis, our results showed that the prognostic value of the signature was significant in all cohorts and independent of TNM stage. To examine the significance of our reprogramming energy metabolism-related signature in patients with the same stage, a stratified analysis was performed. Our results showed that the signature could further divide the subgroups of patients with the same stage. Those patients with high-grade (Stage III–IV) tumors were divided into two subgroups with significantly different survival times (Fig. 2B, log-rank test, P = 0.042). For patients with low-stage (Stage I–II) tumors, our signature could also divide them into two different prognostic groups (Fig. 2A, log-rank test, P = 0.1). Similar results were also found in the other two cohorts (Figs. 2C–2D, Fig. S5).

The signature is associated with adjuvant chemotherapy and can improve the predictive value of TNM stage

Adjuvant chemotherapy is frequently considered for patients with colon cancer. We found that the signature could significantly predict overall survival and relapse-free survival for all samples with and without fluorouracil-based adjuvant chemotherapy (ACT) in GSE39582 (Fig. S6A, Fig. 3A, log-rank test, P = 0.0013 for overall survival; P = 0.0013 for relapse-free survival). The high-risk samples who did not receive ACT, had a significantly worse overall survival and were more likely to relapse (Figs. 1A, 3B, log-rank test, P = 0.0011 for overall survival; P = 0.00031 for relapse-free survival). For the samples who received ACT, there was a similar phenomenon, but without significance (Fig. S6B, Fig. 3C, log-rank test, P = 0.14 for overall survival; P = 0.29 for relapse-free survival). The results implied that the samples could benefit from adjuvant chemotherapy without a significant difference in risk scores. To assess the combinative prognostic value of the signature with clinical factors, we further analyzed the benefit of the signature by evaluating the risk score together with the TNM stage. As shown in Table 3, the risk score and the TNM stage were both independent prognostic factors. So, it is reasonable to generate a new score combining the risk score and the TNM stage. When considering the reprogramming energy metabolism-related signature together with TNM stage, the prognostic power of TNM stage could be significantly enhanced: the AUC was 0.60 (95% CI [0.55–0.65]) for TNM stage alone, and the AUC became 0.67 (95% CI [0.61–0.74]) when adding the risk score (Fig. 4A, P = 0.0093). The other two independent cohorts also showed significant improvement after adding the risk score (P = 0.0072 for GSE17538 and P = 0.0025 for GSE41258, Figs. 4B and 4C).