The data of RNA-seq and clinical features in HCC patients were extracted from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/) and LIRI_JP from International Cancer Genome Consortium (ICGC) (https://icgc.org/). For further analysis, a total of three hundred and twelve patients with both gene expression and overall survival (OS) data were extracted. Moreover, a total of thirty pairs of paraffin-embedded specimens were collected for this study from the pathology department, including both HCC and matched para-carcinoma tissues. All the specimens were obtained from HCC patients who underwent tumor resection. Clinical and pathological information of these specimens are presented in Table S1 . This study was approved by the Ethics Department of the First Affiliated Hospital of Wenzhou Medical University.

We first downloaded KEGG pathways from GSEA (http://software.broadinstitute.org/gsea/index.jsp). Then, we conducted univariate Cox analyses by SPSS 19.0 in TCGA and LIRI_JP, respectively, to acquire significant metabolism-related pathways. We selected the 10 mutual metabolism-related pathways between TCGA and LIRI_JP. For each HCC dataset, we quantified the enrichment levels of the 10 mutual metabolism-related pathways in each HCC sample through the single-sample gene set enrichment analysis (ssGSEA) score. Based on the enrichment levels (ssGSEA scores) of the 10 mutual metabolism-related pathways, we performed hierarchical clustering analysis of HCC.

We obtained the mutation data of HCC patients from TCGA and LIRI_JP. The data containing somatic variants were stored in the form of Mutation Annotation Format (MAF). Mutation data were analyzed and summarized using “GenVisR” package in R software. Copy number variation (CNV) analysis was performed using data of TCGA through cBioPortal (www.cbioportal.org).

To screen out methylation-driven genes, we calculated correlation between gene methylation level and expression using the “MethylMix” package in R software with corFilter = −0.3 and adjusted p-value < 0.05 as the cutoff value. Then, based on the median of beta-values of methylation, we divided HCC patients into two groups and performed Kaplan–Meier (K-M) survival analysis. The genes were considered to be significantly associated with OS based on the threshold of p < 0.05.

We used ESTIMATE (19) to evaluate the immune cell infiltration level (immune score), tumor purity, and stromal content (stromal score) of each HCC sample in Metabolism_H and Metabolism_L cohorts.

After removing data with p-value >0.05 of the correlation between the samples and the immune cells, the RNA-sequencing data of TCGA and LIRI_JP were used to estimate the proportions of 22 types of infiltrating immune cells using the CIBERSORT algorithm following the procedure as previously reported (20). We calculated the expression level of 22 immune cells in each sample through a deconvolution algorithm to quantify the number of cells in each sample. The R packages of e1071, parallel, and preprocessCore were used. The estimate package could adopt the RNA-seq data to calculate the immune and mechanism scores of the specimen, and then evaluate the purity of the tumor. The principle was to evaluate the above content through the signature of the characteristic tumor’s RNA-seq and the input file needed RNA-seq. Common genes data were also required to calculate the inner matrix. We used the estimate package to calculate the score of each sample’s immunity and matrix, then obtained the purity of the tumor and used it for the next calculation.

We first quantified the enrichment levels of immune-related pathways in each HCC sample of Metabolism_H and Metabolism_L by the ssGSEA score. Then, metabolism-related pathways–immune-related pathways networks were drawn by cytoscape online software (http://www.cytoscape.org/).

We performed gene set enrichment analysis of the LIRI_JP and TCGA datasets by GSEA (R implementation) (21, 22). This analysis identified the KEGG (22) pathways that were upregulated or downregulated in Metabolism_H and Metabolism_L, respectively. Terms in KEGG with a false discovery rate (FDR) < 0.05 were considered significantly enriched and were visualized using R package “plyr”, “grid”, “gridExtra”, and “ggplot2” (23). Gene ontology (GO) analysis was performed using the R package “GOplot” (FDR < 0.05) (24).

The least absolute shrinkage and selection operator (LASSO)-Cox regression model (25) was used to identify the most accurate predictive DEMGs for OS. The risk score of each patient was determined based on the DEMG-based classifiers. The patients were categorized into two groups by median score. The survival estimation of patients was analyzed by the K-M method.

The univariate and multivariate Cox regression analyses were conducted to identify significant prognostic predictors associated with OS. The time-dependent receiver operating characteristics (tdROC) analysis by the “timeROC” package of R software was used to assess performance of clinical variables and classifiers. The area under the curve (AUC) of tdROC represented the predictive accuracy. In addition, p-values < 0.05 were considered statistically significant.

Paraffin-embedded liver tissue sections (4 µm) were deparaffinized in xylene and rehydrated in ethanol solutions. Then, the tissue sections were boiled in sodium citrate buffer using a microwave oven for 15 min to perform antigen retrieval, and 3% hydrogen peroxide was used for inhibiting the activity of endogenous peroxidase. Subsequently, to prevent nonspecific binding, the sections were blocked with 5% normal goat serum for 30 min at a temperature of 37°C. Primary antibodies for RRM2 (Proteintech, Wuhan, China, 11661-1-AP; 1:200), LPCAT1 (Proteintech, Wuhan, China, 16112-1-AP; 1:200), G6PD (Proteintech, Wuhan, China, 25413-1-AP; 1:500), CYP3A5 (Proteintech, Wuhan, China, 13737-1-AP; 1:200), CYP2C9 (Proteintech, Wuhan, China, 16546-1-AP; 1:500), BDH1 (Proteintech, Wuhan, China, 15417-1-AP; 1:200), ADH4 (Proteintech, Wuhan, China, 16474-1-AP; 1:200), PYCR1 (ABclonal, Wuhan, China, A13346; 1:100), PTGES (ABclonal, Wuhan, China, A18632; 1:100), HK2 (ABclonal, Wuhan, China, A0994; 1:100), and ADH1C (ABclonal, Wuhan, China, A8081; 1:100) diluted with antibody diluent buffer were added to corresponding tissues and then incubated overnight at 4°C. After washing, secondary antibodies (goat anti-rabbit antibody) (BioSharp Inc., China) were employed to incubate for 1 h. Finally, the sliders were stained with diaminobenzidine (DAB) (Beijing Zhongshan Golden Bridge Biotechnology, China) for color visualization and counterstained with hematoxylin. The staining results were evaluated by three blinded pathologists independently, and their median values were adopted as final score. Staining intensity was divided into four levels according to the following criteria: 0 point (negative), 1 point (weak), 2 points (moderate), and 3 points (strong). The proportion of positive staining area was scored as follows: 0 point (<5%), 1 point (5%–25%), 2 points (26%–50%), 3 points (51%–75%), and 4 points (>75%). The total scores consisted of multiplying the positive staining area percentage scores by staining intensity scores (26, 27).

The Student’s t-test, Wilcoxon test, and other data processing were completed by SPSS 19.0 and GraphPad Prism 7.0 software. The K-M log rank test was calculated by medcalc (Version 19.0). Continuous variables were expressed as the mean ± standard deviation (SD). When all the hypotheses have a p-value < 0.05, the difference is statistically significant.

Firstly, we acquired KEGG signaling pathways by GSEA and singled out 42 metabolism-related pathways ( Table S2 ). There were 14 significant pathways in TCGA ( Figure 2A ) and 19 significant pathways in LIRI_JP ( Figure 2B ), which were both associated with OS of HCC. Among the metabolism-related pathways, there were 10 pathways shared by the two datasets ( Figure 2C ).

Next, the 10 metabolism-related pathways were clustered in TCGA and LIRI_JP, respectively ( Figures 2D, G ). Interestingly, the two datasets showed similar clustering results, with two clusters being clearly separated ( Figures 2E, H ). The 10 metabolic pathways related to prognosis were clustered using the k-means method. As shown in Figure 2E , in the two classifications, 10 pathways have different expression trends. We defined the one with poor prognosis as Metabolism_H, and the other as Metabolism_L. The clusters significantly correlated with histologic grade, TMN stage, and AJCC pathological stage (p < 0.001) ( Table 1 ). Moreover, survival analyses showed that the different metabolic subtypes of HCC had distinct clinical outcomes. The Metabolism_L subtype likely had a better survival prognosis than the Metabolism_H subtype ( Figures 2F, I ).

To explore why the Metabolism_L subtype likely had a better survival prognosis than the Metabolism_H subtype, we analyzed three parts of the differences between the two clusters: the mutant oncogenes, copy number variation, and DNA methylation. Obviously, there were several oncogenes mutated in most HCC patients, especially TP53, TTN, and MUC16 mutated in more than 50 samples (20%) ( Figures 3A–C ). For most cancer types, the TP53, TTN, and MUC16 genes were found to mutate frequently (25). The waterfall plot illustrated that TP53 and TTN were mutated in different numbers of patients. Metabolism_H was more likely to have TP53 mutations (Student’s t-test, p < 0.001) ( Figure 3D ), while TTN showed the opposite. Figure 3E shows that the mutated TP53 gene was less expressed than the wild (Wilcoxon test, p < 0.05).

As demonstrated in Figure S1 , HCC patients showed significantly different CNV in the two subtypes. For instance, according to the blue representing deletions and the red representing amplifications, large chunks of DNA were removed in chromosome 4 of the Metabolism_H subtype. The detailed information of CNV of each gene is shown in Table S3 . Figure S1C shows ten genes with the most obvious differences in CNV between Metabolism_H and Metabolism_L: NUP210L, KCNN3, FAM189B, SCAMP3, CTSS, DPM3, EFNA1, GBA, GBAP1, and KRTCAP2, all of which were amplified on chromosome 1, and more frequent in Metabolism_H.

Through whole-genome DNA methylation analysis, we firstly screened out 240 methylation-driven genes whose genetic expressions were negatively correlated with methylation, and found that 30 of them were related to the HCC prognosis. Then, based on the 30 methylation-driven genes, the heatmap indicated that methylation levels of the 30 genes were significantly different between Metabolism_H and Metabolism_L ( Figure 4A ). For example, the methylation level of PDK4 was higher in Metabolism_H, while the methylation level of TMEM165 was higher in Metabolism_L. Moreover, the survival analyses also showed that the survival rates of HCC patients with hypermethylated PDK4 and hypomethylated TMEM165 were lower, which was consistent with differences in gene methylation levels and prognosis between Metabolism_H and Metabolism_L ( Figure 4B ).

Due to the fact that tumor infiltration lymphocytes were closely related to metabolism and prognosis of tumors (28), we explored the immune infiltration of the metabolism-related pathway-based HCC subtypes according to immune scores ( Figure S2 ). When comparing the tumor immunity microenvironment of the two HCC subtypes, we found that the two clusters were significantly different. Compared with Metabolism_L, the immune scores were significantly higher in Metabolism_H (Kruskal–Wallis test, p < 0.001) ( Figure S2B ). In addition, although stromal scores did not have significant difference between Metabolism_H and Metabolism_L in TCGA, we obtained opposite trends when comparing the tumor purities and stromal scores of the two HCC subtypes. The tumor purity was higher in Metabolism_L while the stromal score was higher in Metabolism_H (Kruskal–Wallis test, p < 0.05) ( Figures S2A, C ). The ESTIMATE Score was also higher in Metabolism_H (Kruskal–Wallis test, p < 0.01) ( Figure S2D ). In conclusion, these results indicated that Metabolism_H contained more immune cells and stromal cells, while Metabolism_L contained more tumor cells.

Therefore, we analyzed immune cell makeups of the two subtypes and found that they were obviously different. Metabolism_H contained more M2 macrophages in LIRI_JP and neutrophils in TCGA, while Metabolism_L contained more CD8 T cells in LIRI_JP and M1 macrophages in TCGA (Wilcoxon test, p < 0.05) ( Figures 5A, B ). Based on the two subtypes, we analyzed their immune pathways, which were also differently expressed ( Figure 5C ). Then, for identifying the relationship between metabolism and immune in HCC, a metabolism-related pathways–immune-related pathways network was conducted including 5 immune-related pathways and 9 metabolism-related signal pathways ( Figure 5D ).

Furthermore, we further analyzed the immune checkpoint-related gene expression and found that the expression levels of PD-L1, CTLA-4, TIM-3 were significantly higher in Metabolism_H ( Figure S3A ). Based on the previous study of the cancer stemness (29), mRNA expression-based stemness index (mRNAsi) was higher in Metabolism_H ( Figure S3B ).

A total of 110 prominent KEGG pathways including pathways expressed differently in Metabolism_H and Metabolism_L were selected ( Table S4 ). The 82 high-expression pathways in Metabolism_H, such as “Pathways in cancer” and “TOLL like receptor signaling pathways”, were related to tumor proliferation and metastasis, indicating worse survival prognosis of Metabolism_H. However, the 28 high-expression pathways in Metabolism_L were mainly concentrated on metabolic process, such as “Tryptophan metabolism”, “Primary bile acid biosynthesis”, and “Retinol metabolism”. Figure 6A shows GSEA enrichment plots of representative gene sets on several representative pathways of Metabolism_H and Metabolism_L.

There were 946 metabolism-related genes (MRGs) that expressed in 110 prominent KEGG pathways. There were 757 DEGs between Metabolism_H and Metabolism_L in TCGA and 2,468 DEGs in LIRI_JP. A total of 135 DEMGs overlapped among TCGA_DEGs, LIRI_JP_DEGs, and MRGs in the Venn diagram ( Figure 6B ). The 135 DEMGs were mainly associated with the following biological processes: small molecule catabolic process; carboxylic acid biosynthetic and catabolic processes; organic acid biosynthetic and catabolic process; cellular amino acid biosynthetic, metabolic, and catabolic process; and alpha-amino acid metabolic and catabolic process ( Figure 6C ). The results also indicated that the DEMGs were mainly associated with the following cellular contents: mitochondrial matrix, peroxisome, and microbody ( Figure 6C ). Moreover, the DEMGs were related to molecular functions, such as coenzyme binding, iron ion binding, and cofactor binding ( Figure 6C ).

In order to investigate the effect of DEMGs on HCC prognosis, we first conducted univariate Cox analyses in TCGA and LIRI_JP, respectively. Sixty DEMGs in TCGA and LIRI_JP related to mortality were identified ( Figures 7A, B ). Among the 60 DEMGs, there are 36 genes shared by the two datasets ( Figure 7C ). According to the results of the LASSO-Cox regression model, 11 prognostic DEMGs with non-zero regression coefficients were finally chosen as the potential prognostic biomarkers for the OS of HCC patients ( Figures 7C, D ). The detailed information of DEMGs for constructing the prognostic signature is summarized in Table S5 . The formula of the eleven-DEMG survival-predictor model was as follows: eleven-DEMG survival-predictor model score = (0.0177074543570851 * RRM2) + (0.000599168151290748 * PYCR1) + (0.000451238392995456 * PTGES) + (0.0000832126857397287 * LPCAT1) + (0.0154490143978134 * HK2) + (0.00824505291131197 * G6PD) − (0.00170896754573046* CYP3A5) − (0.000591128394736733* CYP2C9) − (0.00190395464082525 * BDH1) − (0.000171745083443705 * ADH4) − (0.000267396449305932 * ADH1C). Based on the survival-predictor model, we evenly divided HCC patients into two groups by the median risk score cutoff point, whose value is −0.0687, in TCGA: high risk and low risk ( Figure 7E ). The enrichment levels of the 11 genes in the two groups quantified by the ssGSEA was also significantly different. Then, K-M analysis showed that survival rates were significantly lower in the high-risk group (p < 0.001) ( Figure 7G ). Interestingly, we used the same eleven-DEMG survival-predictor model and cutoff point to cluster patients in LIRI_JP, in which the similar results were obtained ( Figure 7F ). The survival analysis also indicated that high risk had a worse OS (p < 0.001) ( Figure 7H ).

To further verify the bioinformatics analysis results, we collected both HCC and matched para-carcinoma tissues for IHC. Compared with normal tissue, the PYCR1, LPCAT1, and G6PD significantly expressed more in HCC tissue, while CYP3A5, CYP2C9, BDH1, ADH4, and ADH1C expressed less ( Figure 8A ), and the expression of RRM2, PTGES, and HK2 has no significant difference between two tissues ( Figure S4A ). The statistical analysis results are shown in Figures 8B–I and Figures S4B–D . Most of these results were consistent with our prognostic model, further indicating that the gene-based classifier had great value in predicting the mortality for HCC patients.

For identifying the clinical significance of the DEMG-based survival-predictor model, we conducted the univariate Cox analysis in TCGA. The results indicated that tumor stage, T classification, M classification, and risk score were correlated with the survival rates (p < 0.05) ( Figure 9A ). Moreover, in multivariate Cox analysis, the risk score was significant (p < 0.01) while other factors were not associated with OS (p > 0.05) ( Figure 9B ). More importantly, the time-dependent ROC curves suggested that the DEMG-based risk score with an AUC of 0.767 could predict mortality more accurately than other HCC prognostic factors: age (AUC = 0.527), gender (AUC = 0.501), grade (AUC = 0.501), stage (AUC = 0.661), T (AUC = 0.667), N (AUC = 0.494), and M (AUC = 0.506) ( Figure 9C ). In LIRI-JP, we acquired the same results ( Figures 9D–F ).