Identification of a six-gene metabolic signature predicting overall survival for patients with lung adenocarcinoma

Data expression datasets

The transcriptome expression profiles and corresponding clinical information for LUAD were downloaded from The Cancer Genome Atlas (TCGA; http://portal.gdc.cancer.gov/) and Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) databases. From the TCGA, gene expression data were of the HTSeq-FPKM type, obtained from 497 LUAD and 54 non-tumor samples. From the GEO, the GSE68465 dataset included 443 LUAD and 19 non-tumor samples, using the GPL96 platform (Affymetrix Human Genome U133A Array). The metabolic genes in the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway were extracted from Gene Set Enrichment Analysis (GSEA) (https://www.gsea-msigdb.org/gsea/index.jsp), and the overlapping metabolism-related genes were identified from TCGA and GSE68465 (Possemato et al., 2011; Zhu et al., 2020).

Construction and validation of the prognostic metabolic gene signature

The clinical cases from the TCGA database were used to assess the prognostic associations of the metabolic genes with clinical outcomes. The differentially expressed genes (DEGs) between LUAD and paired non-tumor samples were obtained by the Wilcoxon rank sum test using the R package called “limma”, and the adjusted P-value <  0.05 and absolute log2 fold change (FC) >1 were considered as the selection criterion. Univariate Cox regression analysis was used to identify prognosis-related metabolic genes, and adjusted P-values <  0.001 were considered statistically significant. The lasso penalty for Cox proportional hazards model (1,000 iterations) was used to construct the prognostic gene-expression signature utilizing an R package called “glmnet.” The prognostic gene-expression signature was designed using a risk scoring method with the following formula: ${\displaystyle \mathit{Risk}\mathit{score}=\sum _{\mathit{i}}^{\mathit{n}}\left({\mathit{x}}_{\mathit{i}}\ast {\beta }_{\mathit{i}}\right)}$

where x_i indicates the expression of gene i and β_i indicates the coefficient of gene i generated from the Cox multivariate regression.

The R package “survminer” was used to explore the cutoff point of the risk score, which divided patients into high- and low-risk groups. The R package “survival” was used to draw the Kaplan–Meier survival curves to demonstrate the overall survival (OS) in the high- and low-risk groups. The R package “survival ROC” was used to evaluate the prognostic value of the gene-expression signature.

Independence of the prognostic gene signature from other clinical characteristics

To determine whether the predictive power of the prognostic gene-expression signature could be independent from other clinicopathological variables in patients with LUAD (including age, sex, TNM stage, T stage, N stage, and M stage), univariate and multivariate Cox regression analyses were performed. The hazard ratio (HR), 95% confidence intervals (Cis), and P-values were calculated.

Construction and validation of a predictive nomogram

The nomogram was constructed using all the independent prognostic factors of the Cox regression analyses using R package “rms.” Validation of the nomogram was assessed by discrimination and calibration using the concordance index (C-index) by Harrell, Lee & Mark (1996) (bootstraps with 1,000 resamples) and the calibration plot, respectively.

External validation of the prognostic metabolic gene signature

To verify the prognostic metabolic-gene molecular signature in the GEO dataset, the risk score of patients was calculated directly with the gene-expression signature constructed from the TCGA dataset for further analysis. The receiver operating characteristic (ROC) and Kaplan–Meier analyses were performed identically with the gene signature in the TCGA dataset. The mRNA expression levels of the signature genes were analyzed further using online databases (the Oncomine database (http://www.oncomine.org/) and TIMER database (http://cistrome.shinyapps.io/timer/)). The protein expression levels associated with the signature genes were validated using the Human Protein Atlas database (http://www.proteinatlas.org/). The known genetic alterations of the signature genes were investigated using cBioPortal for Cancer Genomics (http://www.cbioportal.org/).

Gene set enrichment analysis

Enrichment analysis of the KEGG pathways of the signature genes was performed using GSEA on the TCGA dataset. The nominal (NOM) P-value <  0.05 and the False Discovery Rate (FDR) q-value <0.25 indicated statistical significance.

Statistical analysis

All analyses were performed using R software v3.6.3 (R Foundation for Statistical Computing, Vienna, Austria). Two-tailed P-values <  0.05 were considered statistically significant.

Clinical characteristics

The TCGA dataset included 486 patients with LUAD (Table S1). The GEO dataset included 443 patients with LUAD (Table S1). Patients with a survival time of less than 30 days were omitted. For the study, 454 and 439 patients remained in the TCGA and GEO datasets, respectively. The detailed clinical characteristics of all patients are listed in Table 1.

Building and validation of the prognostic metabolic gene signature

To clarify our study design, a flow chart of the analysis procedure is presented in Fig. 1. A list of 994 genes in the KEGG pathway was identified from GSEA (Table S2), and 633 overlapping metabolism-related genes were abstracted from TCGA and GSE68465 (Table S3). The 96 DEGs (72 up-regulated genes and 24 down-regulated genes) between LUAD and paired non-tumor samples were identified from the further analysis (Fig. 2; Table S4). Seven significant genes associated with OS were identified using univariate analysis (Table S4). Furthermore, six genes were selected to build the prognostic model using a lasso-penalized Cox analysis (Table 2). The six genes were phosphofructokinase platelet (PFKP), pyruvate kinase muscle (PKM), triosephosphate isomerase 1 (TPI1), lactate dehydrogenase A (LDHA), prostaglandin E synthase (PTGES), and thymidylate synthase (TYMS). Risk score = (0. 00005 × PFKP mRNA level) + (0. 00173 ×PKM mRNA level) + (0. 00038 ×TPI1 mRNA level) + (0. 00379 ×LDHA mRNA level) + (0. 00292 ×PTGES mRNA level) + (0. 02490 × TYMS mRNA level).

The 445 patients with LUAD were divided into the high-risk or low-risk group based on the median risk score of 0.861 in the TCGA dataset. Patients in the high-risk group had significantly poorer OS than those in the low-risk group (P < 0.001; Fig. 3A). The distribution of the risk score and survival status of the patients is presented in Fig. 3C, which showed a higher mortality in the high-risk group than in the low-risk group. The expression of the six prognostic genes is shown in the heatmap. All the six genes had a significant positive correlation with the high-risk group (Fig. 3E). The area under the curve (AUC) of the time-dependent ROC curve was used to identify the prognostic ability of the six-gene molecular signature. The AUCs of the six-gene signature model were 0.693, 0.655, and 0.565 for the 1-, 3-, and 5-year OS, respectively, suggesting that the prediction model had a good performance in predicting the OS in patients with LUAD (Fig. 3G).

The prognostic model was validated in the GSE68465 dataset. The 439 patients with LUAD were divided into the high-risk or low-risk group based on the median risk score of 0.861. Patients in the high-risk group had a poor OS compared with those in the low-risk group (P < 0.001; Fig. 3B). The distribution of the risk score and survival status showed a higher mortality in the high-risk group than in the low-risk group (Fig. 3D). The expression heatmap of the six prognostic genes showed that all the six genes had a significant positive correlation with the high-risk group (Fig. 3F). The AUCs of the six-gene signature model were 0.728, 0.654, and 0.618 for the 1-, 3-, and 5-year OS, respectively (Fig. 3H). Taken together, these results suggested that the prognostic model had a high sensitivity and specificity in predicting the OS in patients with LUAD.

The prognostic gene signature was independent from other clinicopathological factors

Univariate and multivariate Cox regression analyses were conducted to assess the independent predictive value of the six-gene prognostic signature. In the TCGA dataset, univariate Cox regression analysis demonstrated that the prognostic model (HR: 2.845, P < 0.001), TNM stage (HR: 1.666, P < 0.001), T stage (HR: 1.605, P < 0.001), and N stage (HR: 1.806, P < 0.001) had a prognostic value for OS (Fig. 4A). Multivariate Cox regression analysis demonstrated that the only prognostic model (HR: 2.448, P <0.001) and TNM stage (HR: 1.950, P <0.01) were independent prognostic factors for OS (Fig. 4A). In the GSE68465 dataset, the prognostic model, T stage, N stage, and age had a prognostic value in the univariate and multivariate Cox regression analyses (Fig. 4B). Gender was the only independent prognostic factor for OS in the univariate Cox regression analysis (Fig. 4B).

In addition, the time-dependent ROC curve was used to identify the predictive ability of the prognostic model compared with the other clinicopathological characteristics. In the TCGA dataset, the AUCs of the prognostic model were 0.693, 0.655, and 0.565 for the 1-, 3-, and 5-year OS, respectively, which were higher than most of the other clinicopathological characteristics including age (0.498, 0.511, 0.485), gender (0.579, 0.485, 0.451), T stage (0.673, 0.613, 0.608), N stage (0.685, 0.666, 0.628), and M stage (0.508, 0.527, 0.530) (Fig. 5A). Furthermore, in the GSE68465 dataset, the AUCs of the prognostic model were 0.728, 0.654, and 0.618 for the 1-, 3-, and 5-year OS, respectively, which were higher than most of the other clinicopathological characteristics including age (0.593, 0.568, 0.581), gender (0.539, 0.549, 0.547), grade (0.580, 0.571, 0.548), T stage (0.647, 0.606, 0.606), and N stage (0.690, 0.680, 0.655) (Fig. 5B). The prognostic model had a larger AUC value compared with other clinicopathological characteristics. These results indicated that the model was an excellent prognostic model for LAUD patients, especially for the 1- and 3-year OS.

These results suggested that our prognostic model could be an independent predictor of prognosis in patients with LAUD.

Building and validating a predictive nomogram

A nomogram was built to predict the survival probability in patients with LAUD in the TCGA dataset. The nomogram was constructed using four prognostic factors (the TNM stage, T stage, N stage, and prognostic model; Fig. 6A). The C-index was calculated to evaluate the predictive ability of the nomogram for OS. The C-index for the nomogram was 0.754 (95% CI [0.561–0.947]). Calibration plots indicated that the nomogram had a good accuracy in predicting the 1- and 3-year OS (Fig. 6B).

To predict the survival probability more accurately, the combined prognostic model was built based on the nomogram. The combined prognostic model consisted of the TNM stage, T stage, N stage, and prognostic model. A time-dependent ROC curve was used to identify the predictive ability of the combined prognostic model. The AUCs of the combined prognostic models were 0.782, 0.717, and 0.688 for the 1-, 3-, and 5-year OS, respectively, which were higher than other clinical models including the TNM stage model (0.732, 0.687, 0.681), T stage model (0.671, 0.612, 0.613), N stage model (0.686, 0.661, 0.648), and the prognostic model (0.692, 0.634, 0.576). The combined model had the largest AUC value compared with other factors, which indicated that the combined model had a good predictive accuracy for survival. These results suggested that the predictive ability of the combined model built with the nomograms is better than other models, especially for predicting 1- and 3-year survival (Fig. 6C).

Gene set enrichment analysis

To recognize signaling pathways that are differentially activated in LUAD, a GSEA was used, and a total of 49 significantly enriched KEGG pathways were found in the high-risk group and low-risk group (Table S5) of the TCGA dataset (FDR q-val < 0.25, NOM p-val < 0.05). Among them, many enriched pathways were related to metabolism and some highly dysregulated pathways including cell cycle, p53 signaling pathway, and basal transcription factors were also contained in these results (Table S5). We chose the top five significantly enriched metabolism-signaling pathways depending on the normalized enrichment score from the high-risk group or low-risk group. We found that the top five most significantly enriched metabolism-related pathways of the high-risk group were the cysteine and methionine, fructose and mannose, glyoxylate and dicarboxylate, purine, and pyrimidine pathways (Fig. 7A). The top five most significantly enriched metabolism-related pathways of the low-risk group were the alpha linolenic acid, arachidonic acid, ether lipid, glycerophospholipid, and linoleic acid pathways (Fig. 7B). Most of the metabolism-related pathways in the high-risk group mainly focused on amino acid and glycolysis metabolism, while the pathways in the low-risk group mainly focused on lipid metabolism. The results of the ten representative enriched metabolism-related KEGG pathways are given in Table 3. Furthermore, all the six metabolic genes of the prognostic model enriched these metabolism pathways significantly. LDHA enriched the cysteine and methionine pathway (Table S6); PFKP and TPI1 enriched the fructose and mannose pathway (Table S6); PKM enriched the purine pathway (Table S6); TYMS enriched the pyrimidine pathway (Table S6); and PTGES enriched the arachidonic acid pathway (Table S6). The results further elucidated the role of metabolism in LUAD and the value of the six-gene signature in predicting the prognosis of LUAD.

External validation using online databases

To further identify the role of the six metabolic genes in LUAD, we compared the mRNA expression levels of the six metabolic genes (PFKP, PKM, TPI1, LDHA, PTGES, and TYMS) in the LAUD tissues with those in the normal lung tissues using data from the Oncomine database (Fig. 8). Obviously, all the six genes were overexpressed in lung cancer in all the datasets from the Oncomine database with the threshold of fold change = 2, P-value = 0.001 (Fig. 8A). Furthermore, the mRNA levels of all the six genes in LUAD were significantly upregulated than those in normal tissues in the combined LUAD datasets from the Oncomine database (Fig. 8B; Table 4). To further validate the overexpression of the six genes in LUAD, we analyzed the expression of the six genes using TIMER databases (Fig. 9). The results revealed that all the mRNA expression of the six genes in LUAD were significantly higher than in normal tissues. All the results from the Oncomine and TIMER databases were consistent with our results for the TCGA and GEO datasets. In addition, the mRNA expression of the six genes was also higher in esophageal carcinoma, head and neck squamous cell carcinoma, lung squamous cell carcinoma, and stomach adenocarcinoma from the TIMER databases (Fig. 9). The protein expressions of these six genes were analyzed using clinical specimens from the Human Protein Profiles (Figs. 10A and 10B; Table 5). The representative images of the six gene protein levels from the Human Protein Profiles are shown in Fig. 10A. Compared with the expression level in normal lung tissue, LDHA (100%, n = 7) and TYMS (80%, n = 5) showed a significantly higher percentage of high/medium expression levels in the LAUD tissue (Fig. 10B; Table 5). PKM (50%, n = 6), PFKP (33.33%, n = 6), and PTGES (16.67%, n = 6) showed a significant moderate percentage of high/medium expression levels in the LAUD tissue (Fig. 10B; Table 5). However, TPI1 showed no detected expression both in the LAUD and normal lung tissue (Fig. 10B; Table 5). The genetic alterations were explored in the cBioPortal database. Amplifications and mutations were the most common alterations in the six metabolic genes (Fig. 10C). The aberrant genetic alterations might elucidate the overexpression of these six genes in LUAD.

Altogether, the correlation of the aberrant expression of these six genes with LAUD cancer was further validated using multiple online databases.