Comprehensive analysis reveals a metabolic ten-gene signature in hepatocellular carcinoma

Identification and acquisition of TCGC, GEO and ICGC data

Gene expression profiles and clinical data associating with HCC were identified and acquired from TCGA project (https://cancergenome.nih.gov/), GEO database (https://www.ncbi.nlm.nih.gov/gds/) and ICGC project (https://icgc.org/). A total of 415 candidates (365 patients and 50 normal candidates) and the clinical information (age, gender, grade, stage, myometrial invasion, lymph node status, distant metastasis status) were identified and acquired from TCGA dataset. A total of 220 patients and the clinical information (age, gender, ALT, AFP, stage) were identified and acquired from GEO database. A total of 243 patients and the clinical information (gender, age, stage, prior Malignancy) were identified and acquired from ICGC dataset.

Patients characteristics and grouping

A total of 365 patients from TCGA dataset, used as the training cohort. Using the R package Caret, with the ratio of 1:1 in a random manner, 184 patients from TCGA dataset were selected as the internal testing cohort. A total of 243 patients from ICGC dataset and 220 patients from GSE14520 dataset, used as the external testing cohorts. Finally, we integrated all the 828 patients from the TCGA, GEO and ICGC datasets to use as the entire testing cohort. The clinical information of the four cohorts is summarized in Table 1.

Identification of metabolism related genes

Metabolism-related genes in the KEGG pathway associated with metabolism were screened from the GSEA (http://software.broadinstitute.org/gsea/index.jsp), and the overlapping metabolism-related genes were identified from TCGA, GSE14520, and ICGC gene expression profile.

Prognostic genes were identified and gene signature were constructed by utilizing training cohort

Using the Limma version 3.36.2 R package, DEGs were calculated between tumors and paired non-tumor samples of 50 patients from TCGA dataset, the adjusted P-value < 0.05 and absolute log2 fold change (FC) > 1.5 were considered as the selection criterion. By using the R package survival and Coxph function, Univariate Cox proportional hazard regression analysis was performed to discover the prognostic genes in the training cohort, with the adjusted P-value < 0.05 as the significance cutoff. We further narrowed the gene range to construct a gene signature by performing LASSO analysis, R package glmnet was used to perform LASSO analysis.

The performance of gene signature

With the risk-formula: Risk score = expression of gene1 × β1gene1 + expression of gene2 × β2gene2 +…… expression of genen × βngenen, patients in training cohort were divided into low or high risk group basing on the median risk score. The expression level was compared between low and high risk group. Kaplan–Meier analysis was performed to compare the survival rate between the two groups. ROC curve analysis for OS was performed to assess the clinically predictive ability of the gene signature. Next, Univariate and Multivariate Cox proportional hazards analysis were performed to investigate whether the gene signature could be independent of other clinical parameters, including age, gender, grade, stage, myometrial invasion, lymph node status, distant metastasis status. Furthermore, stratification analysis were used to assess the prognostic value of the gene signature in different subgroups stratified clinical variables. Moreover, the diagnostic capability of the gene signature was investigated to further explore the clinical significance, including differentiating normal tissue and HCC, different stages and grades.

Validation of the gene signature

Internal testing cohort, GSE14520 testing cohort, ICGC testing cohort and entire testing cohort were used to validate the reliability and validity of the gene signature. According to the risk-formula and the median risk score, the patients from the four testing cohorts were divided into low or high risk group. And the same analyses were performed to validate the performance, including Kaplan–Meier analysis, the ROC curve analysis, Univariate and Multivariate Cox proportional hazards analysis, stratification analysis.

Go enrichment analysis and Gene set enrichment analysis

Metabolism related DEGs were identified between high- and low-risk groups, with corrected P-value < 0.05 and absolute log fold change (FC) > 1.5 being considered as the cutoff criterion. Next, gene ontology processes were considered as enriched by using the DAVID database (https://david.ncifcrf.gov/). Furthermore, we generated an ordered list of all genes according to their correlation with two subtypes and elucidated the significant survival difference between high- and low-risk groups by GSEA. A total of 1,000 times were performed for gene set permutations. The nominal P value was used to sort the pathways enriched in each phenotype.

Statistical analysis

P < 0.05 was considered statistically significant, The statistical analyses were conducted by employing the R (version 3.4.3) and GraphPad Prism 7.

Identification of prognostic genes

We carry on our study as described in the flow chart (Fig. 1). A total of 944 metabolism-related genes were identified from the online website GSEA (Table S1) and 918 overlapping metabolism-related genes were obtained between TCGA, GSE14520 and ICGC gene expression profile (Table S2). A total of 71 up-regulated genes and 25 down-regulated genes were calculated for subsequent analysis (Fig. 2; Table S3). By conducting Univariate Cox proportional hazards regression analysis, 11 significant genes associated with OS were obtained for further analysis (adjusted P < 0.05) (Table S4).

Establishment of gene signature from the training cohort

LASSO analysis was performed to narrow the gene range to construct a ten-gene signature from 11 significant genes (Risk-formula: Risk score = expression of G6PD × 0.0015558631423743 + expression of LPCAT1 × 0.00310356967612766 + expression of ME1 × 0.00632401235412166 + expression of PRIM1 × 0.0026285256475841 + expression of RRM2 × 0.00979851158278189 + expression of TXNRD1 × 0.00783495109084 + expression of UCK2 × 0.057281553558209 + expression of CAD × 0.06164749118803 + expression of DTYMK × 0.00849922964649245 + expression of ENTPD2 × 0.03691738437542). The risk score was computed for each patient in the training cohort, 365 patients from training cohort were divided into low risk group (183 patients) and high risk group (182 patients) according to the median risk score: 0.73. A total of 243 patients from ICGC testing cohort were divided into low risk group (100 patients) and high risk group (143 patients). A total of 220 patients in GSE14520 testing cohort were divided into low risk group (112 patients) and high group (108 patients). A total of 184 patients in the internal testing cohort were divided into low risk group (95 patients) and high risk group (89 patients), 828 patients in the entire testing cohort were divided into low risk group (356 patients) and high risk group (432 patients). The 10 prognostic genes expression level distribution between low and high risk group of the training cohort was showed in Fig. 3A, the expression of ten prognostic genes in high risk group was higher than in low risk group, the result was consisted with internal testing cohort in Fig. 3B, GSE14520 testing cohort in Fig. 3C and ICGC testing cohort in Fig. 3D.

The performance of gene signature

As showed in Fig. 4A, with the increasing risk score, patients in the training cohort have a worse OS, the expression ten prognostic genes increased. The ROC curve was presented to assess the clinically predictive ability of the gene signature, the AUC for 1-year (Fig. 4B), 3-year (Fig. 4C), and 5-year (Fig. 4D) OS were 0.805, 0.756, 0.716 for training cohort, which was higher than other clinical characteristics, including age (0.512, 0.536, 0.540), gender (0.504, 0.535, 0.565), grade (0.522, 0.505, 0.527), stage (0.659, 0.688, 0.669). Besides, we found the gene signature also performed more specific and sensitive than any single gene (Table 2). The OS in the training cohort was significantly different between low and high risk group, in the 1-year, 3-year and 5-year, the OS in the high risk group were 0.576, 0.181, 0.076, the OS in the low risk group were 0.814, 0.306, 0.142, the result indicated that patients with a high-risk score have more poor OS than the patients with low-risk score (P < 0.001), detail was presented in Fig. 4E. By conducting Univariate and Multivariate Cox regression analysis, we noted that gene signature have a significant correlation with worse OS, the HR of the gene signature was 3.919 (95% CI [2.766–5.555]) with P-value < 0.001 in Univariate Cox regression analysis (Fig. 4F), 3.490 (95% CI [2.401–5.073]) with P-value < 0.001 in Multivariate Cox regression analysis (Fig. 4G), thus the gene signature was an independent prognostic factor of other clinical variables.

Validation of the gene signature

To validate the predictive ability in different HCC populations, we applied the gene signature to ICGC testing cohort, the result was similar to the training cohort. Figure 5A showed the distribution of risk scores for each patients, patients in high risk group had a worse OS than patients with a low-risk group. In addition, the AUC for 1-year (Fig. 5B), 3-year (Fig. 5C), and 5-year (Fig. 5D) OS were 0.775, 0.754, 0.778, which was higher than other clinical characteristics, including age (0.542, 0.520, 0.595), gender (0.587, 0.574, 0.545), prior Maliganancy (0.526, 0.572, 0.508). Even though the AUC of gene signature was a little less than the TNM stage at 1-year OS (0.775 vs 0.809), the AUC of gene signature was much larger than the TNM stage at 3-year OS (0.754 vs 0.658), 5-year OS (0.778 vs 0.564). The AUC of ROC for gene signature was obviously greater than single gene (Table 2). Moreover, the OS for patients in the high risk group was 0.576 at 1-year, 0.181 at 3-year, 0.076 at 5-year, compared with 0.814, 0.306, 0.142 in the low risk group (P < 0.001, Fig. 5E). Further Univariate Cox regression analysis and Multivariate Cox regression analysis displayed gene signature was a powerful and independent factor in external testing cohort (Figs. 5F and 5G). For GSE14520 testing cohort (Fig. 6), internal testing cohort (Fig. 7) and entire testing cohort (Fig. 8), the gene signature had the similar predictive ability. The distribution of risk scores, gene expression were evaluated in GSE14520 testing cohort (Fig. 6A), internal testing cohort (Fig. 7A) and entire testing cohort (Fig. 8A). The ROC curve demonstrated that gene signature was more specific and sensitive than any clinical characteristics and any single gene in GSE14520 testing cohort (Figs. 6B–6D; Table 2), internal testing cohort (Figs. 7B–7D; Table 2) and entire testing cohort (Figs. 8B–8D; Table 2). Patients in a high-risk group have poorer OS than the patients in low-risk group for GSE14520 testing cohort (Fig. 6E), internal testing cohort (Fig. 7E), entire testing cohort (Fig. 8E) (P < 0.001). Univariate and Multivariate Cox regression analysis indicated the gene signature was an independent prognostic factor for GSE14520 testing cohort (Figs. 6F and 6G), internal testing cohort (Figs. 7F and 7G) and entire testing cohort (Figs. 8F and 8G).

Stratification analysis

To further demonstrate the clinical significance of the gene signature in HCC, we perform the survival analysis stratified by clinical variables (age, gender, grade, stage) in training cohort and internal testing cohort, by clinical covariates (age, gender, stage) in ICGC testing cohort, GSE14520 testing cohort and entire testing cohort. Patients of stage I–II, stage III–IV, grade 1–2, grade 3–4, age <60, age >=60, female, male were stratified into high risk group and low risk group. The log-rank test indicated that HCC patients in high risk group still had obviously worse OS than patients in low risk group for training cohort (Fig. 9A), internal testing cohort (Fig. 9B), ICGC testing cohort (Fig. 9C), GSE14520 testing cohort (Fig. 9D) and entire testing cohort (Fig. 9E), and the high-risk patients of subgroup subdivided by the signature had poorer survival than the low-risk patients in training cohort (Figs. S1A–S1H), internal testing cohort (Figs. S1I–S1P), GSE14520 testing cohort (Figs. S2A–S2C, S2E and S2F), ICGC testing cohort (Figs. S2G–S2L) and entire testing cohort (Figs. S3A–S3F). There was no different trend between high and low risk group for female patients in GSE14520 testing cohort, small number female patients may be an important reason. Interestingly, subgroup stage I, stage II, stage III in entire testing cohort were stratified into high risk group and low risk group, HCC patients in high risk group had poorer OS than patients in low risk group (Figs. S3G–S3I).

Bidkhori et al. (2018) published a article that was proposing 3 sub-types of HCC, including iHCC1, iHCC2 and iHCC3. They reported that the expression of iHCC3 tumors are markedly distinct from those of iHCC2 and iHCC1, and a larger number of genes are differentially expressed between iHCC3 and iHCC1/iHCC2 compared with iHCC1 vs iHCC2. Consistent with the result, we revealed that the expression level of majority of genes of the gene signature were higher in iHCC3 than iHCC1/iHCC2, and the expression level of iHCC1 was similar to iHCC2 (Fig. S4).

Diagnostic capability

The gene signature have been validated the predictive ability in different HCC populations. To further explore diagnostic capability of the gene signature, we compared risk score between normal liver and HCC in training cohort (Fig. S5A, P < 0.0001), ICGC testing cohort (Fig. S5B, P < 0.0001) and GSE14520 testing cohort (Fig. S5C, P < 0.0001), the risk score in HCC was obviously higher in than normal liver. The AUC of ROC curve was 0.98 in training cohort (Fig. S6A), 0.77 in ICGC testing cohort (Fig. S6B) and 0.97 in GSE14520 testing cohort (Fig. S6C), which revealed the strong diagnostic capability for HCC. In addition, we also investigated the distribution of high-risk and low-risk patients in different stages and grades, the proportion of high risk patients is higher in advanced tumor grade (grade 3–4 ) than early tumor grade (grade 1–2) in training cohort (Fig. S5D). The proportion of high risk patients is higher in late stage (TNM III–IV) than early stage (TNM I–II) in training cohort (Fig. S5E), ICGC testing cohort (Fig. S5F) and GSE14520 testing cohort (Fig. S5G). Furthermore, as the TNM stage and tumor grade increased, the risk score increased (Figs. S5H–S5K). For early and advanced tumor grade in training cohort (Figs. S6D and S6E), early and advanced TNM stage in training cohort (Figs. S6F and S6G), ICGC testing cohort (Figs. S6H and S6I), GSE14520 testing cohort (Figs. S6J and S6K), the AUC of ROC curve indicated the modest diagnostic capability.

Comparing the performance of the gene signature with other gene signatures and TMN stage

TMN stage is still significant to predict the survival of HCC patients. In the training cohort, we found the AUC of gene signature was larger than the TNM stage at 1-year, 3-year, 5-year (0.805 vs 0.659 at 1-year, 0.756 vs 0.688 at 3-year, 0.716 vs 0.669 at 5-year ) (Figs. 4B–4D), the result was consistent with GSE14520 testing cohort (Figs. 6B–6D), internal testing cohort (Figs. 7B–7D) and entire testing cohort (Figs. 8B–8D). For ICGC testing cohort (Figs. 5B–5D), even though the AUC of gene signature was a little less than the TNM stage at 1-year OS (0.775 vs 0.809), the AUC of gene signature was much larger than the TNM stage at 3-year OS (0.754 vs 0.658), 5-year OS (0.778 vs 0.564). Thus we believed the gene signature was more specific and sensitive than TNM stage. Next, we also compared the AUC of ROC between the gene signature and all single genes, the AUC of the gene signature was larger than any single gene (Table 2). Moreover, we further compared the gene signature with other gene signatures, the ROC analysis indicated our model had better performance, 0.805 at 1-year, 0.756 at 3-year, 0.716 at 5-year, the AUC of Long et al. (2018) model is 0. 7674, 0.7040, 0.6919 at 1, 3 and 5-year, the AUC of Qiao et al. (2019) model is 0.71, 0.69 at 3 and 5-year, the AUC of Li et al. (2017) model is 0.67, 0.67 at 3 and 5-year, the AUC of Xiao-Hong Xiang et al. (2019) model is 0.708, 0.699, 0.678 at 1, 3 and 5-year, the AUC of Li et al. (2017) model is 0.727, 0.709, 0.604 at 1, 3 and 5-year, the AUC of Yan et al. (2019) model is 0.712, 0.661 at 1, 3-year, indicating that the model had a high sensitivity and specificity for predict the survival.

Go enrichment analysis and Gene set enrichment analysis

We performed a differential expression analysis between low and high risk group to reveal the association between the metabolic genes and two subtypes. We identified 18 upregulated and 14 downregulated genes between patients in the low vs high risk group. Go enrichment analysis revealed that upregulated genes in high risk group were mainly involved in cell cycle regulation, such as DNA replication, chromatid/chromosome segregation and regulation, spindle organization and recombination. Upregulated genes in low risk group were mainly involved metabolic and energy regulation, including metabolism of amino acids, oxidative phosphorylation, fatty acid β oxidation and catabolism (Fig. S7A).

In order to further explore the mechanism of prognostic genes in patients with hepatocellular carcinoma, we conducted GSEA between low and high risk group to identify the significant pathways (FDR < 0.05, NOM P-value < 0.05). The most meaningful pathway were identified according to the FDR standard, we uncovered the most meaning pathways which were active in the high-risk group, including KEGG_CELL_CYCLE, KEGG_NUCLEOTIDE_EXCISION_REPAIR, KEGG_OOCYTE_MEIOSIS, KEGG_PURINE_METABOLISM, KEGG_PYRIMIDINE_METABOLISM. And the most meaning pathways which were active in the low-risk group, including KEGG_DRUG_METABOLISM_CYTOCHROME_P450, KEGG_FATTY_ACID_METABOLISM, KEGG_GLYCINE_SERINE_AND_THREONINE_METABOLISM, KEGG_PRIMARY_BILE_ACID_BIOSYNTHESIS, KEGG_VALINE_LEUCINE_AND_ISOLEUCINE_DEGRADATION. The result was consistent with go enrichment analysis, significant pathways which were active in high risk group were mainly related with cell cycle regulation. In contrast, meaningful pathways which were active in low risk group were mainly associated with metabolic and energy regulation (Fig. S7B).