A Novel Metabolism-Related Signature as a Candidate Prognostic Biomarker for Hepatocellular Carcinoma

Introduction

According to the 2018 Global Cancer Statistics, liver cancer is the sixth most common human malignancy worldwide (841,080 cases/year) and the fourth leading cause of cancer-related deaths (781,631 cases/year).¹ Hepatocellular carcinoma (HCC) originates from hepatocytes and accounts for the vast majority (approximately 75–85%) of primary liver cancers. Despite recent advances in early detection and management, mortality from HCC still remains high due to the low rate of early diagnosis.² Therefore, it is necessary to urgently identify effective prognostic biomarkers and therapeutic targets to complement and improve current screening strategies for HCC diagnosis and prognosis.

Recently, metabolic reprogramming has been recognized as a novel and essential hallmark of cancer cells.³ Cancer cells are able to regulate their metabolism to promote the production of energy, the synthesis of macromolecules and the maintenance of redox balance so as to support their rapid cell growth and proliferation.^4,5 Increased glycolysis under normoxic condition (Warburg effect) and glutamine metabolism are the best characterized metabolic changes seen in tumor cells.^5,6 Accumulating evidence has indicated that abnormal metabolism is associated with poor prognosis in many tumor types, including HCC.^7,8 Luo et al⁹ screened the combination of phenylalanyl-tryptophan and glycine as liver cancer metabolic markers, and confirmed that their diagnostic sensitivity was superior to AFP. A series of combined markers, such as lactic acid, glycolamine, phenylalanine, aconitic acid and ribose, were screened for identification of patients with and without recurrence.¹⁰ A recent study revealed thousands of metabolism-relevant genes encoding all known human metabolic enzymes and transporters,¹¹ however, the definitive role and underlying mechanisms of them in the development and prognosis of HCC remains poorly understood. More importantly, there are few models based on metabolism-related genes to predict the prognosis of HCC patients. Thus, it is of great clinical significance to establish a novel metabolism-related prognostic signature that can reliably predict the prognosis of HCC.

In this study, we identified the differentially expressed MRGs in HCC patients by associating the gene expression profiles with MRGs in two database systems, namely, TCGA and the ICGC. The functional enrichment analysis and prognostic value of the differentially expressed MRGs were also determined. The prognostic model was finally established based on Cox regression and Lasso regression analyses. To verify the accuracy of the model, the Kaplan-Meier (KM) estimator and the receiver operating characteristic (ROC) curve were applied. Moreover, the prognostic value of our metabolism-related prognostic model was further validated in TCGA database.

Materials and Methods

Data Collection

The transcriptomic and the corresponding clinical data of patients with hepatocellular carcinoma (HCC) were downloaded from the International Cancer Genome Consortium (ICGC; https://icgc.org/) and The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov/) databases. The RNA-seq data, including 243 HCC and 202 adjacent non-tumor cases from ICGC database, and 374 HCC and 50 adjacent non-tumor cases from TCGA database were examined. The current study did not require ethical approval, because all data had been collected from publicly available data in the ICGC and TCGA databases.

The Differentially Expressed MRGs Screening

In total, 2752 metabolism-related genes (MRGs), which have been reported to encode all known human metabolic enzymes and transporters, were obtained from Possemato’s research.¹¹ The differentially expressed MRGs in HCC and normal tissues were detected using the Wilcoxon test method.¹² |logFC|>1 and adjusted p<0.05 were considered as significant. MRGs with common differential expression were identified from the ICGC and the TCGA databases by using FunRich software. These intersection MRGs were selected for further analyses.

Functional Enrichment Analysis of MRGs

Gene Ontology (GO)¹³ and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis¹⁴ were performed by “clusterprofiler” R package and “enrichplot” R package in R software.¹⁵ Functional categories with false discovery rate (FDR) of less than 0.05 were considered statistically significant.

Survival-Associated MRGs

We collected clinical information of 232 HCC patients from the ICGC database. A total of 3 patients who were followed for less than 30 days or more than 2000 days were excluded to avoid the interference of irrelevant factors. Survival analysis was performed on 229 patients. To screen out the MRGs associated with prognosis in HCC patients, univariate Cox analysis was implemented by the R “survival” package, and data were visualized using forest plots. Only differentially expressed MRGs with a p value <0.001 were considered for subsequent analyses. To screen candidate small molecules that might be useful for HCC patients, the Broad Institute’s Connectivity Map (CMap) was used in predicting which compounds on the base of the prognostic MRGs.¹⁶

Construction of Metabolism-Related Signature for HCC

HCC patients in ICGC were used for constructing the COX prognostic signature, and patients in TCGA dataset were used for validating the signature. Lasso and multivariate Cox regression analyses were performed to construct a prognostic model.¹⁷ To avoid overfitting of the prognostic signature and consequent removal of highly related survival-related prognostic MRGs, Lasso Cox regression was carried out using R “survival” and “glmnet” package. MRGs detected via Lasso algorithm were evaluated by stepwise multivariate Cox regression analysis. By weighting the estimated Cox regression coefficients, the model of tumor-related metabolism genes risk was constructed.¹⁸ The prognostic metabolism-related gene signatures were shown as risk score = Ʃ (βi × Expi), where βi referred to the coefficients and represented the weight of the respective signature, whereas Expi represented the expression value. Based on the risk score formula, patients were assigned into the low-risk group and high-risk group with the median risk score as the cut-off point. The Kaplan-Meier (K-M) survival curve was used with the Log rank test to evaluate the differences in survival rate between the two groups. Furthermore, the receiver operating characteristic (ROC) curve was implemented by R “survival ROC” package¹⁹ and the corresponding Area Under the ROC Curve (AUC) was measured to assess the sensitivity and specificity of the metabolism-related signature.

Validation of Metabolism-Related Signature

To verify the prognostic value of metabolism-related signature, we used the TCGA database as a validation cohort. The same formula was used to calculate the risk scores for each patient. Survival and ROC curve analyses were implemented as described above. In addition, univariate and multivariate analyses were used to estimate the effect of risk scores on OS and the clinicopathologic features. We also explored the correlation between the expression of these MRGs and several clinical features. For further validation of our analyses, the Human Protein Atlas (HPA) online database (http://www.proteinatlas.org/) was applied to identify the expression of these MRGs at a translational level.²⁰

Cell Culture

Human normal hepatocyte cell line (LO2) and HCC cell lines (HepG2, Hep3B, HLF and PLC/PRF/5) which were purchased from National Collection of Authenticated Cell Cultures (https://www.cellbank.org.cn) were maintained in the DMEM medium (Gibco, Wuhan, China) with 10% heat-inactivated fetal bovine serum (FBS, GibcoBRL) and antibiotics of 1% streptomycin and penicillin at 37° C in an atmosphere of 5% CO2.

Quantitative Real-Time PCR

Total RNA from collected cells was extracted using TRIzol reagent. Quantitative real-time PCR was performed as described previously.²¹ The primer pairs used were as follows: human SLC5A11 (forward 5′-CCAAGGACATTCAGCACTGG-3′ and reverse 5′-GAAGCCATAGAAACAGCGGG-3′), human FLVCR1 (forward 5′-AGCTCTTCAAGACAGTCCCC-3′ and reverse 5′-TTAGCCCAATCCTTCCAGCA-3′), human RRM2 (forward 5′- CCAAGGACATTCAGCACTGG-3′ and reverse 5′-GAAGCCATAGAAACAGCGGG-3′), and human GAPDH (forward 5′-TCAAGAAGGTGGTGAAGCAGG-3′ and reverse 5′-TCAAAGGTGGAGGAGTGGGT-3′). The results were normalized to GAPDH.

Western Blotting Analysis

The protein extraction and Western blot analysis were performed as described previously.²²

Statistical Analysis

All statistical analyses were performed by version 3.6.1 of R software (https://www.r-project.org/) and version 3.1.3 of FunRich software (http://www.funrich.org/), along with version 3.7.2 of Cytoscape software (http://www.cytoscape.org/). If not otherwise stated, data were considered to be statistically significant with p value <0.05.

Results

Identification of Differentially Expressed MRGs

The flowchart of our study was shown in Figure S1. A total of 2752 MRGs that encoded all known human metabolic enzymes and transporters were obtained from Possemato’s study¹¹ (Table S1), and we matched these genes with the sequence data of HCC related mRNA in the ICGC and TCGA databases. The differentially expressed MRGs from HCC and normal samples were identified by Wilcoxon rank-sum test. Considering the cutoff criteria (adjusted p value <0.05 and |log FC| > 1.0), 475 differentially expressed MRGs (consisting of 94 downregulated and 381 upregulated genes) were extracted from the ICGC database (Figure 1A and C, Table S2), and 251 differentially expressed MRGs (consisting of 36 downregulated and 215 upregulated genes) were extracted from the TCGA database (Figure 1B and D, Table S3). Finally, a total of 178 differentially expressed MRGs that were common for both databases (consisting of 28 downregulated and 150 upregulated genes) (Figure 1E, Table S4) were selected for subsequent analysis.

Figure 1 Differentially expressed metabolism-related genes (MRGs) in hepatocellular carcinoma (HCC). Heatmap of MRGs between HCC and nontumor tissues in ICGC database (A) and TCGA database (B). The color from blue to red represents the progression from low expression to high expression; Volcano plot of MRGs in ICGC database (C) and TCGA database (D). The red dots in the plot represents upregulated genes and blue dots represents downregulated genes with statistical significance. Black dots represent no differentially expressed genes. (E) Venn diagram showing the gene numbers of the MRGs of HCC in ICGC and TCGA database.

GO, KEGG and PPI Analysis of Metabolism-Related DEGs

To evaluate the potential molecular mechanisms of MRGs in HCC, the 178 differentially expressed MRGs were further analyzed by GO functional annotation and KEGG pathway enrichment. The results showed the top 10 biological processes including GO terms, cellular component GO terms, molecular function GO terms (Figure 2A), and the top 10 KEGG pathway terms (Figure 2C). The correlation between the intersection genes and the top 5 biological processes (organic anion transport, organic acid transport, carboxylic acid transport, lipid catabolic process, and monovalent inorganic cation transport) was shown in Figure 2B. The KEGG analysis revealed that the intersection genes were associated with lipid metabolism and amino acid metabolism, especially arachidonic acid metabolism.

Figure 2 The GO and KEGG analysis of differentially expressed MRGs. (A) the top 10 of biological processes GO terms, cellular component GO terms, molecular function GO terms. (B) The correlation between intersection genes and top 5 biological processes GO terms. (C) The KEGG pathway analysis of differentially expressed MRGs.

Prognostic Values and Candidate Drugs of Survival-Related MRGs

To better define the characteristics of survival-related MRGs, we explored the differentially expressed MRGs associated with survival in HCC samples. In the univariate Cox proportional hazards regression analysis, 17 genes with significant effects on prognosis were identified. The resulting forest plot shown in Figure 3A demonstrated that 15 genes have the characteristics of HR>1, whereas 2 genes have the characteristics of HR<1. This suggested that most survival-related MRGs were risk factors for HCC prognosis. According to CMap database analysis, Cmap mode-of-action (MoA) analysis of 21 compounds indicated 11 mechanisms of action shared by the above drugs (Figure 3B), five drugs (BRD-K34437622, carmofur, raltitrexed, tegafur, capecitabine) shared the MoA of thymidylate synthase inhibitor, four drugs (CO-102862, flecainide, mexiletine, zonisamide) shared the MoA of sodium channel blocker, three drugs (cladribine, clofarabine, gemcitabine) shared the MoA of ribonucleoside reductase inhibitor, and two drugs (floxuridine, fludarabine) shared the MoA of DNA synthesis inhibitor. Our study identified candidate drugs targeting the survival-related MRGs, and might provide therapeutic targets for further analysis (Figure 3B). Furthermore, scatter plots were used to display the expression patterns of 17 differentially expressed MRGs from HCC and normal liver tissues in ICGC database (Figure 3C) and TCGA database (Figure 3D). Scatter plot showed expression patterns of 1 down-regulated genes (MOGAT2) and 16 up-regulated genes (ASRGL1, COX7B2, MTHFD1L, B3GALNT1, ITPKA, SCN4A, PLCB1, NQO1, RRM2, FLVCR1, TMCO3, CAD, TK1, TYMS, SLC5A11 and SPNS1) (Figure 3C and D).

Figure 3 Identification of survival-related differentially expressed MRGs and candidate drugs. (A) Forest plot of hazard ratios showing survival-related MRGs. p values <0.001 are considered to be statistically significant. (B) CMap database analysis identifies novel candidate drugs targeting the prognostic genes. The expression of 17 metabolism-related prognostic genes between HCC and normal tissues in ICGC database (C) and TCGA database (D), the red box shows the gene expression in HCC and the blue box shows the gene expression in normal liver tissue, ***p < 0.001.

Construction of Metabolism-Related Prognostic Signature for HCC

The clinical data of 229 HCC patients in ICGC database and 193 HCC patients in TCGA database are shown in Table 1. To establish the metabolism-associated prognostic signature, the Lasso regression and multivariate Cox proportional hazards regression analyses were conducted. The Lasso regression analysis was applied to exclude genes that might be highly correlated with other genes (Figure 4A and B, Table S5). Then, a prognostic signature model was established based on six MRGs selected from further multivariate Cox regression analysis (Figure 4C). Finally, a prognostic model was constructed to evaluate the prognosis of each patient as follows: Risk score = (0.031755×expression value of COX7B2) + (−2.27007×expression value of SCN4A) + (−0.15919 ×expression value of MOGAT2) + (0.164689×expression value of FLVCR1) + (0.02915×expression value of RRM2) + (0.068413×expression value of SLC5A11).

Table 1 Clinical Characteristics of Hepatocellular Carcinoma Patients in the ICGC and TCGA Database

Figure 4 Establishment of metabolism-related prognostic signature. (A) Screening of optimal parameter (lambda) at which the vertical lines were drawn. (B) Lasso coefficient profiles of the seventeen MRGs with non-zero coefficients determined by the optimal lambda. (C) Multivariate analyses assessing relationship between expression levels of MRGs and OS in patients with HCC.

Then, the risk score of each patient was calculated according to this prognostic model. Based on the median risk score, 229 HCC patients were classified into a high-risk group (n = 114) and a low-risk group (n = 115). The risk score, survival status and gene expression heatmap of these prognostic MRGs are presented in Figure 5A–C. Kaplan-meier Log rank test indicated that patients in the high-risk group showed markedly poorer OS than those in the low-risk group (Figure 5D). Areas under the curve with the value of the signature predicting 1-, 3- and 4-year OS rates were 0.805, 0.803 and 0.94, respectively, indicating that this prognostic model exhibited a good sensitivity and specificity (Figure 5E). Since only one patient had a follow-up period of 5 years, we did not plot the ROC curve for 5 years.

Figure 5 Construction of the metabolism-based prognostic risk signature in the ICGC cohorts. (A) The risk score distribution of HCC patients. (B) Survival status and duration of patients; (C) Heatmap of the metabolism‐related genes expression; (D) Survival curves for the low-risk and high-risk groups. (E) Time-independent receiver operating characteristic (ROC) analysis of risk scores for prediction the overall survival in the ICGC set.

Validation of the Metabolism-Related Prognostic Signature by TCGA Database

The TCGA database including 193 HCC samples was used for the validation of metabolism-related signature. According to the median risk score, we divided patients into high-risk (n = 80) and low-risk groups (n = 113). Consistent with the results derived from the ICGC dataset, the Kaplan–Meier curve demonstrated that patients in the high-risk group exhibited markedly poorer OS than those in the low-risk group (p< 0.001; Figure 6D). The risk score, survival status and gene expression heatmap of these prognostic MRGs are shown in Figure 6A–C. The AUCs for 1-, 3- and 5-year OS rates were 0.721, 0.693 and 0.737, respectively (Figure 6E). Univariate and multivariate Cox regression analyses were conducted to evaluate the independent prediction ability of metabolism-related prognostic signature between the signature and other common prognostic factors, such as age, gender, histological grade, pathological stage and TNM stage. Although univariate Cox analysis indicated that pathologic stage, T stage and our model were markedly associated with OS (Figure 6F, p<0.001), after the multivariate analysis, only metabolism-related prognostic signature (p<0.001) could be used as an independent prognostic factor (Figure 6G). These results demonstrated that this prognostic model exhibited great applicability and stability in predicting the prognosis of HCC patients.

Figure 6 Validation of the metabolism-based prognostic risk signature in the TCGA cohorts. (A) The risk score distribution of HCC patients. (B) Survival status and duration of patients. (C) Heatmap of the metabolism-related genes expression. (D) Survival curves for the low-risk and high-risk groups. (E) Time-independent receiver operating characteristic (ROC) analysis of risk scores for prediction the overall survival in the TCGA set. (F) Univariate Cox regression analysis of discrete clinical factors. (G) Multivariate Cox regression analysis of discrete clinical factors.

To validate the MRGs in this model, the protein expression levels were analyzed using the HPA database. The results showed that FLVCR1, SLC5A11, MOGAT2, and RRM2 protein levels matched their mRNA expression levels (Figure 7). However, representative images of the SCN4A and COX7B2 protein levels were not available in the HPA database.

Figure 7 Verification of hub MRGs expression in HCC and normal liver tissue using the HPA database. (A) FLVCR1, (B) SLC5A11, (C) RRM2, (D) MOGAT2.

Clinical Value of Prognostic Signature

To further evaluate the clinical value of MRGs, the relationship between MRGs’ prognostic indicators and clinical features were investigated, and the results indicated that FLVCR1, MOGAT2, RRM2, SCN4A, and COX7B2 were differentially expressed in patients with various clinical features (Figure 8). To validate the clinical value of the metabolism-related prognostic signatures, the association between the risk score and clinical characteristics was subsequently assessed, and the results demonstrated that high-risk scores were positively associated with histological grade, and survival status in patients with HCC (Figure 8).

Figure 8 Relationships between MRGs expression and clinicopathological factors in HCC (p < 0.05).

Validation of Hub Genes by WB and qRT-PCR

In order to further validate the hub genes expression in the lab, FLVCR1, SLC5A11 and RRM2 were tested by WB and qRT-PCR in normal hepatocyte cell line (LO2) and HCC cell lines (HepG2, Hep3B, HLF and PLC/PRF/5). As illustrated in Figure 9, the protein and mRNA expression levels of FLVCR1, SLC5A11 and RRM2 were significantly increased in human hepatocellular carcinoma cell lines compared with LO2, which were consistent with our bioinformatics analysis results.

Figure 9 Validation of hub genes by WB and qRT-PCR. The protein and mRNA expression of FLVCR1, SLC5A11 and RRM2 were further validated in HCC cell lines by WB (A) and qRT-PCR (B). *p < 0.05, ***p < 0.001.

Discussion

The metabolic changes and regulatory mechanisms involved in hepatocellular carcinoma are very complex and thus need further elucidation. HCC is different from normal liver tissue in terms of glucose metabolism, lipid metabolism and protein metabolism, and also different from other tumors. These disparate metabolic phenotypes provide a theoretical basis for therapeutic strategies targeting tumor metabolism in patients with HCC. To date, there are only few reports on MRGs-based prognostic signatures for HCC patients. Therefore, there is an urgent need to determine reliable metabolic biomarkers and predictive models to predict the prognosis of patients with HCC.

In the current study, based on the analysis of ICGC dataset and TCGA dataset, 178 differentially expressed MRGs were screened. In order to better understand how differentially expressed MRGs were involved in biological processes and signal transduction processes, GO and KEGG enrichment analyses were conducted. The results of the GO enrichment revealed that the above genes were primarily related to metabolic response. The KEGG pathways were mainly focused on metabolism-related pathways. In univariate regression analysis, performed for the differentially expressed metabolism-related genes, 17 genes were detected to display significant association with OS. Then, the metabolism-related prognostic model based on six prognostic MRGs (FLVCR1, MOGAT2, RRM2, SCN4A, SLC5A11, and COX7B2) was constructed by the Lasso and multivariate Cox regression analyses. Using this model, every HCC patient was assigned a particular risk score. The differences in survival between patients with low and high scores were significant in both the ICGC cohorts and the TCGA cohorts. The ROC curves and AUCs indicated that both models of the ICGC and TCGA cohorts performed well. Furthermore, we confirmed the important roles of MRGs in HCC through immunohistochemistry. Finally, FLVCR1, SLC5A11 and RRM2 were successfully validated with their high expression in HCC cell lines, which were in agreement with the results of bioinformatics analysis, indicating that bioinformatics analysis method was reliable.

Six metabolism-related genes that constituted the prognosis model were identified as potential biomarkers of HCC. SLC5A11, an inositol-specific sodium-dependent glucose cotransporter responsible for inositol uptake, has been reportedly related to anaplastic thyroid carcinoma, but its relationship with HCC has not been described previously.^23,24 In this study, our results indicated SLC5A11 might act as a potential oncogene in HCC. Several previous reports indicated that Ribonucleotide Reductase M2 subunit (RRM2) played an important role in the proliferation, invasion and metastasis of tumor cells, and thereby participating in many types of malignancies tumor including HCC.^25–27 A study done by Liu et al revealed that HBV induced upregulation of RRM2 and promoted the development of HCC, while suppression of RRM2 inhibited the growth of HCC.²⁸ In line with these studies, we demonstrated a high expression of RRM2 that might serve as a poor prognostic factor for HCC. Monoacylglycerol O-acyltransferase 2 (MOGAT2), a member of MOGAT gene family, played an important role in the monoacylglycerol pathway for triacylglycerol synthesis and was highly expressed in the small intestine and liver of humans.^29,30 Recent studies have shown that MOGAT2 was an attractive target for the treatment of obesity, type 2 diabetes mellitus, and nonalcoholic steatohepatitis.^31,32 Lian et al reported that MOGAT2 could have applications as a novel biomarker affecting the metastasis of colorectal cancer treatment.³³ So far, little is known about the existence of molecular links between MOGAT2 and HCC. In our study, however, MOGAT2 was found to be highly decreased which might indicate its protective role as a prognostic marker in HCC. Sodium voltage-gated channel alpha subunit 4 (SCN4A), encoding a skeletal muscle voltage-gated sodium channel, Nav1.4,³⁴ has been reportedly associated with neuromuscular disorders like hypo- and hyperkalaemic periodic paralyses, paramyotonia congenita, sodium channel myotonias and congenital myasthenic syndrome.³⁵ However, its underlying molecular mechanisms in cancers have not been explored. Liang et al reported that the single nucleotide polymorphism of the cytochromooxidase VIIb 2 (COX7B2) gene was associated with familial nasopharyngeal carcinoma.³⁶ On the other hand, COX7B2 was found to be a susceptibility gene for type 2 diabetes in an extended Arab family.³⁷ Nevertheless, little is known about the roles of COX7B2 in human cancers and other diseases. Herein, we demonstrated that COX7B2 was overexpressed in patients with HCC. Feline leukemia virus subgroup C receptor 1 (FLVCR1) encoding a ubiquitous heme exporter was demonstrated to promote the proliferation and tumorigenicity through inhibiting apoptosis and autophagy.^38,39 Recently, Wei et al reported that FLVCR1 showed the intracellular overproduction in HCC cell lines and that FLVCR1 expression was significantly associated with survival in HCC patients, which was consistent with our study.⁴⁰

In the current study, we not only identified six powerful MRGs prognostic signature for HCC, but also provided potential drugs for targeted therapy by targeting these six hub genes. As predicted, part of the 21 targeted drugs were already used in clinical for HCC treatment. Hepatic arterial infusion chemotherapy with raltitrexed and oxaliplatin post-TACE was a safe and efficacious therapy for patients with HCC.⁴¹ Capecitabine could be considered an active and safe option in advanced HCC, especially for patients unfit for other treatments.⁴² On the other hand, we developed six hub genes expression in cholangiocarcinoma (another common liver cancer), only FLVCRL1 and RRM2 were significantly overexpressed in cholangiocarcinoma, but all six hub genes expression could not be prognostic signatures for CholangiocarcinomaD patients (Figure S2). They were also not overexpressed in liver tissue from pathologies, such as cholestasis, NASH, cirrhosis (data were not shown). Herein, six powerful MRGs prognostic signatures were specific to HCC.

Although our research provides new insights into the treatment of HCC, there are some limitations. Firstly, the diagnostic efficiency and prognostic value of the key genes were analyzed and verified only in TCGA dataset. Therefore, predictive performance needs to be validated in other independent databases. Secondly, further experimental studies are needed to validate the predictive value of the six MRGs based on tumor samples and clinical data. In addition, more in vivo and in vitro experiments should be conducted in order to explore the functional role of the critical genes in the setting of HCC.

Conclusions

In this study, we assessed the expression profiles of metabolism-related genes based on both ICGC and TCGA databases. Moreover, a novel metabolism-related prognostic model was constructed, which could serve as an independent prognostic predictor for patients with HCC. The prognostic value of the metabolism-related prognostic model was further validated in TCGA database. Candidate drugs that aimed at hub ERGs were identified. Finally, the protein and mRNA expression were further validated in human hepatocellular carcinoma cell lines. Taken together, our findings may provide a novel prognostic model that would be effective in predicting the prognosis of HCC patients and providing targeted candidates drug for treatment.