A model of twenty-three metabolic-related genes predicting overall survival for lung adenocarcinoma

Data collections

First, we got clinical information for patients with LUAD from TCGA (https://portal.gdc.cancer.gov/). It included 497 LUAD patients with mRNA expression profiles and clinical follow up information was available for our study. The number of obtainable clinical cases for the selected subjects was 445 after removing 52 patient samples from the study due to a lack of clinical information (such as survival time, T stage, N stage, and so on) or survival time less than 30 days (avoiding non-cancer-related death samples). A total of 445 LUAD patients and their information were utilized to build a training cohort for identifying prognostic metabolic-related genes and building a prognostic risk model. Next, we downloaded the LUAD gene expression data from GEO (https://www.ncbi.nlm.nih.gov/geo/) in two accessed datasets GSE30219 and GSE72094. Removing other cancer pathological types, such as lung squamous cell carcinoma, 83 LUAD samples, and 393 LUAD samples were utilized to build a testing cohort for validating the prognostic value of the TCGA-LUAD prognostic risk model (up to April 01, 2020).

Identification of metabolic-related genes in TCGA-LUAD

First, we obtained 944 hub metabolism-related genes from the intersection of the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb) and TCGA-LUAD. Then, a Wilcoxon signed-rank test was performed on normal and cancer tissues in the training cohort by “limma” R package (—log FC—>  0.5; FC: fold change; a false discovery rate (FDR) P <  0.05) (Diboun et al., 2006). The heatmap was plotted by the ”pheatmap” R package and we obtained 336 metabolic-related differentially expressed genes. Second, after univariate Cox regression analysis, 59 metabolic genes were retained (P <  0.05) by using the method that the correlation between expression values of metabolic genes and survival of samples in the training cohort. Last, we performed lasso-penalized Cox regression analysis to identify more important metabolic genes for OS prediction through the “glmnet, survival” R package (Zhang et al., 2019). We obtained twenty-three metabolic-related genes for risk model building. The three-step screening method was robust and performed via Perl (https://www.perl.org/) and R (version 3.6.1).

Building the prognostic metabolic gene signature

To construct the prognostic model, we utilized lasso-penalized Cox regression analysis to select the prognostic metabolic-related gene (Tibshirani, 1997). We obtained a risk score for each patient by their coefficient. Risk score= (Coef AKR1A 1 × expression of AKR1A1) + (Coef NT5E × expression of NT5E) + ……(Coef TYMS × expression of TYMS) (Liu et al., 2020). R software packages “survival” and “survminer” were used to calculate the optimal cut-off value for risk scores and plot Kaplan–Meier survival curves (Chan et al., 2018). Using the median as a point of differentiation, we differentiated patients into two groups: high-risk and low-risk. The R package “survivalROC” was used to plot time-dependent ROC curves for predicting the diagnostic value (Heagerty, Lumley & Pepe, 2000). The concordance index (C-index) was used to evaluate the predictive ability of the risk model.

Verification of the prognostic signature as an independent risk factor and correlation analysis between the clinical characteristics and risk scores

Patients with complete information on the corresponding clinical data were available for univariate and multivariate analysis. P <  0.05 symbolizes statistically significant (Liu et al., 2020). We performing the Student’s t-test to verify the correlation between clinical characteristics and risk scores.

Construction and verification of the predictive nomogram

The nomogram was built by the “rms” R package according to training cohort data (Iasonos et al., 2008). In our study, the tumour-node metastasis (TNM) model and the prognostic signature were integrated into the predictive nomogram. We performed the calibration plot and C-index to investigate the predictive ability of the nomogram. The calibration plot was used to assess whether the numerical value of the predicted value of the model and the probability of the occurrence of the ending event were consistent (Fenlon et al., 2018). We used C-index to assess the predictive ability of the nomogram. It estimates the probability that the predicted result is consistent with the actual observed result. We compared the TNM model, prognostic model, and the nomogram model through ROC analysis and C-index (Liu et al., 2020). Next, we verified the nomogram by C-index, ROC analysis, and calibration plot in the testing cohort data.

KEGG and GO pathways enrichment analyses

To study the biological role of mRNA markers in LUAD patients, we utilized the Kyoto Encyclopedia of Genes and Genomics (KEGG) and Gene Ontology (GO) pathway enrichment analysis to explore which pathways the differentially expressed genes were mainly enriched in (up to April 11, 2020). Gene Set Enrichment Analysis (GESA) (https://www.gsea-msigdb.org/gsea/index.jsp) was utilized to find enriched terms in the training cohort or testing cohort (Subramanian et al., 2005). We choose “c2. cp. kegg. v6.2. symbols. gmt. gene sets” as a reference gene set from the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb/). P <  0.05, FDR q-value <  0.25, and normalized enrichment score —NES— ≥1 suggested statistically significant. We plotted the results by “ggplot2, gridExtra, grid, plyr” R package. All operations are carried out in GSEA_4.0.3. GO pathway enrichment analysis of metabolic genes was performed by “clusterProfiler, org.Hs.eg.db, plot, ggplot2” in R p ackage (Pathan et al., 2015).

Statistical analyses

logFC: logarithmic value of FC; positive/negative logFC indicates the logarithmic foldness of upregulation/downregulation; —log FC—>0.5 indicates multiple differences in the gene expression greater than 0.5 between normal tissues and cancer tissues. A coefficient is a number that expresses a measurement of a particular quality of a substance or object under specified conditions (Bøvelstad et al., 2007). The risk score was calculated according to the formula: ${\displaystyle \mathrm{Risk score}\left(\mathrm{patients}\right)=\sum _n\left(\mathrm{coefficient}\left(\mathrm{mRNAn}\right)\ast \mathrm{expression (mRNAn)}\right).}$ C-index had a lower accuracy from 0.50–0.70, medium accuracy between 0.71–0.90, and higher accuracy when greater than 0.90 degrees (Kim, Schaubel & Mccullough, 2018). The Area Under Curve (AUC) value of the time-dependent ROC lay in the range of 0.5–0.9 viewed as statistically significant. All statistical analyses were conducted through R 3.6.1.

Construction of the prognostic signature from the training cohort

The workflow of the study was shown in Fig. 1. In our study, 445 LUAD patients (Table S1, clinical information in Table S1) from the TCGA data set were assigned to the training sample cohort; 83 LUAD patients (Table S3, clinical information in Table S4) and 393 LUAD patients (Table S5, clinical information in Table S6) from the GEO data set were assigned to the testing sample cohort for batch processing. After the Wilcoxon signed-rank test was applied to the training set, we obtained 336 meaningful metabolic genes (Fig. 2, Table S7). From the 336 meaningful metabolic genes, we got 59 mRNAs which were considered to be significantly associated with OS in LUAD patients, 42 high-risk genes, and 17 low-risk genes (Fig. 3, Table S8). Finally, lasso-penalized Cox analysis identified 23 genes (AKR1A1, NT5E, PTGIS, GMPS, MBOAT1, ADCY9, B4GALT1, MAOB, INPP4B, NEU1, ALDOA, ENTPD2, GNPNAT1, GSTA3, PKM, HK3, ALDH2, AK2, LDHA, CHPT1, SMS, CTPS2, and TYMS) to construct the prognostic model (Table 1).

The twenty-three metabolic genes signature and predictability assessment in the training cohort

Finding that there is a significant and independent correlation between the expression of twenty-three prognostic mRNAs and OS, we believed that combining the 23 metabolic genes to form a twenty-three gene signature could predict a patient’s prognosis. According to the optimal cut-off of 1.53, we classified training cohort samples into two groups: a high-risk group (n = 222) and a low-risk group (n = 223) (Table S9). The OS rate between the two risk groups was significantly different (P = 5.543e−10) (Fig. 4A).

The AUC was 0.798, 0.747, and 0.734 for the 1, 3, and 5-year OS, separately (Fig. 4B). We plotted patients’ risk curves in the training cohort and analysed their distribution in Figs. 4D–4E, and the heat map reveals the prognostic mRNA expression patterns between two distinct prognostic patient groups (Fig. 4C).

Validation of the twenty-three metabolic gene signatures

To test the robustness of the prognostic signature, according to the risk score cut-off of the training cohort, GSE30219 was structured into a high-risk group (n = 25) and a low-risk group (n = 58) (Table S10); and GSE72094 was structured into a high-risk group (n = 196) and a low-risk group (n = 197) (Table S11). In the GSE30219 cohort, the Kaplan–Meier survival curves of the prognostic signature have a statistically significant difference in the two predicted risk groups (Fig. 5, P = 1.176e−02). The AUC was 0.694, 0.645, and 0.637 for the 1, 3, and 5-year OS, separately (Fig. 5B). The risk curves and expression of twenty-three metabolic genes in the GSE30219 cohort were shown in Figs. 5C–5E Fig 5C-E. And in the GSE72094 cohort, Kaplan–Meier survival curves of the prognostic signature has a statistically significant difference in the two predicted risk groups (Fig. 6, P = 1.417e−10). The AUC was 0.695, 0.725, and 0.742 for the 1, 3, and 5-year OS, separately (Fig. 6B). The risk curves and expression of the twenty-three metabolic genes in the GSE30219 cohort were shown in Figs. 6C–6E. We also compared our results with the two published gene signature studies (Xie & Xie, 2019; Zhao, Li & Tian, 2018) and found that the C-index results were better than their signatures (Table 2).

Taken together, these results indicated a better predictive performance about our prognostic signature.

Correlation analysis between the prognostic signature and clinical characteristics

445 patients with their information from TCGA-LUAD cohort were utilized for the correlation analysis. Being male (P = 0.017) and TNM stage (P < 0.01) have a significant correlation with a higher risk score. Samples with higher T, N, and M grading were also significantly correlated with a higher risk score (Table 3).

83 patients with their information from GSE30219 cohort were utilized for the correlation analysis. TNM stage was significantly associated with a higher risk score (P <  0.001). Samples with higher T, N, and M grading were also significantly correlated with a higher risk score (P <  0.01) (Table 3). And 393 patients with their information from GSE72094 cohort were utilized for the correlation analysis. TNM stage was significantly correlated with a higher risk score (P <  0.001) (Table 3).

Validation of the independent prognostic factor

We analysed 445 patients, with a median age of 68, grouping them by gender (Table 4). The results identified that our prognostic signature was an independent OS prognostic indicator (Figs. 7A and 7B). We also verified that the prognostic signature can serve as an independent prognostic indicator in the GSE30219 cohort (Figs. 7C and 7D) and GSE72094 cohort (Figs. 7E and 7F). The multivariate Cox analysis indicating that the prognostic signature was significantly associated with OS in each cohort when adjusted for the TNM stage (Table 4). And stratification analysis showed that high-risk group was significantly correlated with a poorer OS (Figs. 8A–8F). However, in the GSE30219 testing cohort, patients in the high-risk group from TNM stages III and IV show no significantly correlated with OS (Fig. 8D).

Construction and verification of the predictive nomogram

We built the nomogram by including the independent prognostic roles (Fig. 9A). Calibration plots verified the performance of the nomogram (Fig. 9B). The C-index of the TNM model, prognostic signature, and nomogram model were 0.654, 0.730, and 0.793, separately (Table 5). The AUC was 0.838, 0.785, and 0.779 for the 1, 3, and 5-year OS, respectively. The nomogram model showed a better AUC predicting 1, 3, and 5-year OS than the prognostic model in the training cohort (Figs. 9C–9E).

Next, we verified the clinical application of the nomogram in the GSE30219 and GSE72094 cohorts. In the GSE30219 and GSE72094 testing cohorts, calibration plots verified the performance of the nomogram (Figs. 10A and 10E). In the GSE30219 testing cohort, the C-index of the TNM model, prognostic signature, and nomogram model were 0.612, 0.682, and 0.684 (Table 5). The AUC was 0.731, 0.686, and 0.722 for the 1, 3, and 5-year OS, separately (Figs. 10B–10D). In the GSE72094 testing cohort, the C-index of the TNM model, prognostic signature, and nomogram model were 0.579, 0.709, and 0.713, separately (Table 5). The AUC was 0.708, 0.794, and 0.771 for the 1, 3, and 5-year OS, separately (Figs. 10F–10H). The nomogram model identified the better AUC predicting 1, 3, and 5-year OS than the prognostic signature in both testing cohorts.

Taken together, the nomogram model increased the predicting ability of the prognostic signature. These results indicated the better predictive performance of the nomogram model.

KEGG and GO pathways enrichment analyses

KEGG enrichment analyses suggesting that a majority of the metabolism-related pathways such as the metabolism of fatty acid, arachidonic acid, glycerophospholipid, alpha-linolenic acid, and pyrimidine were associated with the low-risk group, while the cell cycle, mismatch repair, ubiquitin-mediated proteolysis, and p53 signalling pathways were associated with the high-risk group (Fig. 11A, Table 6, Table S12). Besides, we performed GO pathway enrichment analysis on 23 genes. These 23 genes were statistically significant in metabolic-related pathways, such as nucleotide biosynthetic processes, secretory granule lumen, and monosaccharide binding in Biological Processes (BP), cellular component (CC), and molecular function (MF) (Fig. 11B, Table S13).