Identification and validation of a two-gene metabolic signature for survival prediction in patients with kidney renal clear cell carcinoma

Differentially expressed MRGs and functional enrichment analysis

The expression levels of 944 MRGs, obtained from The Molecular Signatures Database (MSigDB, were assessed in KIRC samples and normal kidney samples from the TCGA database using the Wilcoxon signed-rank test. A total of 123 differentially expressed MRGs, including 60 upregulated and 63 downregulated genes, were eventually identified based on the criteria of |log2FC| > 2 and FDR < 0.05 (Figure 1A, 1B).

Functional enrichment analysis was next performed to explore potential molecular mechanisms related to the differentially expressed MRGs. The most enriched GO terms in the biological process (BP) category were ‘small molecule catabolic process’, ‘cellular amino acid metabolic process’, and ‘organic acid catabolic process’. Significantly enriched GO terms related to the cellular component (CC) category included ‘microbody part’, ‘peroxisomal part’, and ‘peroxisomal matrix’. In the molecular function (MF) category, the differentially expressed MRGs were highly enriched in the terms ‘oxidoreductase activity’, ‘coenzyme binding’, and ‘iron ion binding’ (Figure 1C). In addition, KEGG pathway analysis revealed that these genes were notably associated with pathways in retinol metabolism, metabolism of xenobiotics by cytochrome P450, and drug metabolism − cytochrome P450 (Figure 1D).

Identification of a metabolism-related prognostic signature for KIRC

To explore the prognostic value of MRGs in renal cancer progression, we performed univariate Cox regression analysis to examine potential relationships between the expression levels of the 123 MRGs and patient OS. Results demonstrated that 15 MRGs were significantly associated with OS (P < 0.01) (Figure 2A). Among those, P4HA3, IL4I1, RRM2, ITPKA, PSAT1, TYMP, HK3, PLCB2, and AANAT were considered as risk genes (HR > 1), while AGMAT, GATM, HAO2, FBP1, ADH6, and ALDH6A1 were considered as protective genes (HR < 1). We then used least absolute shrinkage and selection operator (LASSO) Cox regression on the above-mentioned 15 MRGs to identify the most optimal risk score model for predicting survival in KIRC patients (Figure 2B, 2C). Eventually, RRM2 and ALDH6A1 were retained as target genes, and their respective coefficients were calculated to construct the metabolism-related prognostic signature.

We constructed the OS prognostic signature based on the expression of the 2 target genes and their prognostic coefficients using the following formula: Risk score= (0.0361 × expression level of RRM2) + (–0.0184 × expression level of ALDH6A1). According to the median risk score, 253 and 254 KIRC patients were sorted into a high-risk group and a low-risk group, respectively. The Kaplan-Meier curve displayed a significant difference in OS between the high- and the low-risk groups (5-year survival rate, 49.3% vs. 72.6%, respectively; P <0.001) (Figure 2D). We then applied the ROC curve to evaluate the predictive accuracy of the signature. The area under the ROC curve was 0.705, suggesting a moderate prognostic value (Figure 2E). In addition, the distributions of risk scores and survival status of patients were ranked according to the risk scores (Figure 2F, 2G).

Evaluation of the prognostic signature as an independent prognostic factor

We performed univariate and multivariate Cox regression analyses to further determine whether the prognostic signature could serve as an independent prognostic factor. Univariate analysis revealed that age, grade, AJCC stage, T stage, M stage, and risk score were significantly associated with OS (Figure 2H). Subsequent results showed that age (P < 0.001), grade (P = 0.024), AJCC stage (P = 0.043), and risk score (P < 0.001) were still significantly correlated with OS in multivariate analyses (Figure 2I). These data indicate that our metabolism-related prognostic signature is an independent prognostic factor for KIRC patients.

Clinical utility of the prognostic signature

We further explored the relationship between the metabolism-related prognostic signature and various clinical parameters. The expression levels of the two signature MRGs in the high- and low-risk groups are shown in a heatmap (Figure 3A). This analysis showed that RRM2 and ALDH6A1 were expressed at high and low levels, respectively, in the high-risk group. Additionally, we observed significant differences between low- and high-risk groups in grade (P = 4.841e−08), AJCC stage (P = 1.648e−08), T stage (P = 3.826e−07), and M stage (P = 5.918e−05) (Figure 3B–3E).

To better evaluate the survival outcomes and define the broad applicability of the prognostic signature, we next performed survival analyses stratified by age, gender, tumor grade, AJCC stage, T stage, and M stage. As shown in Figure 4, patients in the high-risk group had significantly shorter OS than those in the low-risk group for cases with age ≤ 60 (P = 1.236e−02), age > 60 (P = 2.404e−05), female gender (P = 1.139e−04), male gender (P =3.673e−04), G1-2 (P = 3.163e−02), G3-4 (P = 2.236e−02), AJCC stage I and II (P =1.166e−02), T1-2 (P = 2.623e−04), and M0 (P = 5.778e−03). In contrast, no significant differences were observed for OS between high- and low-risk groups for KIRC patients with AJCC stage III and IV (P = 7.808e−02), T3-4 (P = 1.673e−01), and M1 (P = 7.726e−02).

We then analyzed the relationship between the expression level of each target gene in the prognostic signature and clinicopathological features to assess the function of the two genes in disease progression. The results indicated that enhanced RRM2 expression was significantly associated with advanced tumor stage and high-grade tumor, suggesting in turn a positive correlation between RRM2 expression and poor prognosis in KIRC. Indeed, the highest RRM2 expression was detected in the most progressive clinicopathological stage, that is, G4 and stage IV, T4, and M1 (Figure 5A). In contrast, ALDH6A1 expression tended to gradually decrease with KIRC progression (Figure 5B), which suggested that ALDH6A1 is a protective factor for KIRC.

Construction and validation of a predictive nomogram

A predictive nomogram was constructed by incorporating the prognostic signature and several clinical parameters to generate individual numerical probabilities of OS (Figure 6A). The C-index of the developed nomogram was 0.774. DCA demonstrated that the nomogram provided a higher net benefit in predicting OS if the threshold probability was larger than 3% (Figure 6B). Additionally, the nomogram displayed an obviously higher net benefit than did tumor grade and AJCC stage. The calibration curves indicated that the nomogram performed well in predicting 1-, 3- and 5-year OS compared with the ideal model (Figure 6C–6E). These results are thus indicative of the reliability and predictability of the nomogram.

External validation of the expression and prognostic value of RRM2 and ALDH6A1

Consistent with our results, analysis of the Oncomine (Figure 7A), TIMER (Figure 7B), and GEPIA (Figure 7C) databases revealed that RRM2 expression was significantly upregulated, while ALDH6A1 expression was significantly downregulated, in KIRC samples compared with normal ones. Moreover, similar expression trends for these two genes were frequently detected in various cancer types (Figure 7B). In addition, the prognostic value of each gene was further confirmed by Kaplan Meier analysis in the GEPIA database, which indicated favorable prognosis for the RRM2 low-expression and the ALDH6A1 high-expression groups (Figure 7D).

To determine protein expression levels of RRM2 and ALDH6A1, we interrogated the HPA database. The information available indicated that ALDH6A1 expression in kidney cancer is significantly lower than in normal kidney tissue (Figure 7E). In contrast, no significant difference was found for RRM2 expression between kidney cancer and normal tissues. In parallel, gene expression levels of RRM2 and ALDH6A1 were further verified in two independent cohorts (GSE53757, Figure 7F, and GSE66270, Figure 7G) in the GEO database. Results showed that RRM2 was significantly overexpressed, while ALDH6A1 was significantly underexpressed, in KIRC samples compared with normal ones. Therefore, our assessment of multiple datasets was highly consistent, especially at the transcriptional level, with our expression data for the two prognostic signature genes.

GSEA analysis of signature genes

The top 50 genes with significant correlation with our two prognostic genes were retrieved using GSEA. Subsequently, we performed MSigDB Hallmark analysis for RRM2 and ALDH6A1. Results indicated that the most significant pathways related with RRM2 included ANTIGEN PROCESSING AND PRESENTATION, ALLOGRAFT REJECTION, SYSTEMIC LUPUS ERYTHEMATOSUS, LEISHMANIA INFECTION, and AUTOIMMUNE THYROID DISEASE (Figure 8A). A heatmap showing transcriptional expression profiles of the top 50 RRM2-correlated genes for each phenotype (Figure 8B). Correspondingly, GSEA showed that significantly enriched pathways for ALDH6A1 included VALINE LEUCINE AND ISOLEUCINE DEGRADATION, FATTY ACID METABOLISM, PROPANOATE METABOLISM, CITRATE CYCLE TCA CYCLE, and PYRUVATE METABOLISM (Figure 9A). Transcriptional expression profiles of the top 50 genes associated with ALDH6A1 for each phenotype are shown in a heatmap on Figure 9B.

Study cohort

RNA sequencing transcriptome data and clinical information were downloaded for 539 KIRC and 72 normal tissues from The Cancer Genome Atlas (https://cancergenome.nih.gov/). The extracted data were normalized and processed by log2 transformation. KIRC samples with incomplete data or corresponding to survival times shorter than 30 days were excluded. Metabolism-related gene (MRG) sets were obtained from the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways in The Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb/) [33]. A total of 944 genes were selected for further analysis after removing duplicates.

Identification of differentially expressed MRGs and enrichment analysis

Differentially expressed MRGs between cancer samples and normal control samples were identified by LIMMA package in R language [34] with the following thresholds: |log2 fold change (FC)| > 2.0 and a false discovery rate (FDR) corrected p-value < 0.05 (calculated by the Benjamini and Hochberg procedure). A series of functional enrichment analyses was performed to explore potential roles for the differentially expressed MRGs using gene ontology (GO) and KEGG. The ggplot2 and enrichplot R packages were used to supply the visual enrichment maps of annotation analysis results to help interpretation. P < 0.05 was regarded as statistically significant.

Construction and evaluation of the prognostic signature

The relationship between the expression levels of the differentially expressed MRGs and OS was first examined via univariate Cox regression analysis. Genes with significant prognostic value (P < 0.05) were included for subsequent validation. To obtain the most optimal MRGs and to control the complexity of the model, we carried out a LASSO Cox regression analysis with the glmnet package. LASSO Cox regression is a robust model building method that reduces the dimension of the model and prevents data overfitting [35]. After that, several target genes were obtained to develop the prognostic signature. Based on the signature, the risk score for each patient was calculated using the following formula:$$\text{Risk}\text{Score}={\sum }_{i=1}^{n}Coef(i)\times x(i),$$ where n represents the number of module genes, Coef(i) denotes the estimated regression coefficient by LASSO analysis, and x(i) indicates the relative expression level of each MRG.

The median risk score was taken as the cut-off point to separate all KIRC patients into high-risk and low-risk groups. The Kaplan–Meier method and the log-rank test were used to compare the difference in OS between high- and low-risk groups. Time-dependent receiver operating characteristic (ROC) curve analysis was conducted to assess the predictive accuracy of the prognostic signature; the value of the area under the ROC curve (AUC) ranged from 0.5 (no predictability) to 1 (perfect predictability).

Univariate and multivariate Cox regression analyses were applied to assess whether the metabolism-related signature could be an independent predictor of OS for the TCGA-KIRC cohort. Age, gender, tumor grade, American Joint Committee on Cancer (AJCC) stage, T stage, and M stage were used as covariates. N stage was not analyzed because of a large proportion of missing data.

To detect the prognostic value of the risk score in different subgroups, stratified survival analysis was carried out according to clinical characteristics related to the prognosis: age (≤ 65 and > 65 years), gender, tumor grade (G1-2 and G3-4), AJCC stage (I/II and III/IV), T stage (T1-2 and T3-4), and M (M0 and M1) stage. The relationship between the expression of each target gene and clinical parameters was also compared to further understand the impact of each individual target gene in our prognostic signature.

Development of a predictive nomogram

A nomogram incorporating several clinical variables (age, gender, grade, AJCC stage, T stage, and M stage) and the risk score calculated from the prognostic signature was developed to evaluate the probability of 1-, 3-, and 5-year OS for TCGA-KIRC patients. Subsequently, the discriminative and predictive abilities of the nomogram were evaluated by determining the concordance index (C-index) and by generating calibration plots. C-index values range from 0.5 to 1.0, indicating respectively no discriminating power and perfect discriminating ability. A reference line with a slope of one in the calibration plots represents perfect calibration. Decision curve analysis (DCA) was performed to evaluate the clinical utility of the signature-based nomogram model by quantifying the net benefits under different threshold values.

External verification of the prognostic genes

To further evaluate the reliability of the prognostic gene signature set, we investigated the expression of RRM2 and ALDH6A1 in both tumor and normal tissues using Oncomine (https://www.oncomine.org/resource/main.html) and Tumor Immune Estimation Resource (TIMER, https://cistrome.shinyapps.io/timer/) databases. In addition, the Gene Expression Profiling Interactive Analysis (GEPIA, http://gepia.cancer-pku.cn/) online tool was accessed to verify the prognostic value of the above target genes through survival analysis. The expression of the prognostic genes in the gene signature was further investigated at the protein level using the Human Protein Atlas (HPA, https://www.proteinatlas.org/) database. Then, differences in the expression of these target genes between tumor and normal tissue samples were verified in the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/).

Gene set enrichment analysis (GSEA)

GSEA (http://www.broadinstitute.org/gsea) was performed using 1000 permutations for each task to identify significantly enriched pathways. The minimum and maximum criteria for selection of gene sets from the collection were 15 and 500 genes, respectively. Significantly related genes were defined with a nominal P value < 0.05 and FDR < 0.25.

Statistical analyses

All statistical analyses were performed with R 3.6.2 (https://www.r-project.org/). All analyses were two-sided and statistical significance was defined as P < 0.05.