Novel Prognostic Signatures of Hepatocellular Carcinoma Based on Metabolic Pathway Phenotypes

Hepatocellular carcinoma is a disastrous cancer with an aberrant metabolism. In this study, we aimed to assess the role of metabolism in the prognosis of hepatocellular carcinoma. Ten metabolism-related pathways were identified to classify the hepatocellular carcinoma into two clusters: Metabolism_H and Metabolism_L. Compared with Metabolism_L, patients in Metabolism_H had lower survival rates with more mutated TP53 genes and more immune infiltration. Moreover, risk scores for predicting overall survival based on eleven differentially expressed metabolic genes were developed by the least absolute shrinkage and selection operator (LASSO)-Cox regression model in The Cancer Genome Atlas (TCGA) dataset, which was validated in the International Cancer Genome Consortium (ICGC) dataset. The immunohistochemistry staining of liver cancer patient specimens also identified that the 11 genes were associated with the prognosis of liver cancer patients. Multivariate Cox regression analyses indicated that the differentially expressed metabolic gene-based risk score was also an independent prognostic factor for overall survival. Furthermore, the risk score (AUC = 0.767) outperformed other clinical variables in predicting overall survival. Therefore, the metabolism-related survival-predictor model may predict overall survival excellently for HCC patients.

Introduction

Hepatocellular carcinoma (HCC) is one of the most prevalent primary cancers worldwide and ranks third in all cancer-related mortality (1, 2). Numerous therapeutic strategies for treating liver cancer have been developed, including surgical resection, radiofrequency ablation, liver transplantation, and targeted therapy (3, 4). Despite the fact that great progress has been made in clinical treatment, the survival rates of liver cancer patients within 5 years are still as low as 18% because of the highly malignant tumors, the high recurrence rate, and drug resistance (5, 6). Several factors have been identified to affect and predict HCC prognosis, such as microRNAs and blood groups (7, 8). However, these factors still could not predict the prognosis accurately. Therefore, it is paramount to exploring how specific cellular tumor progression pathways contribute to HCC prognostic stratification for cancer treatment development.

It is proposed that cancer cells must modify their metabolic programs to obtain energy and macronutrient during rapid growth (9, 10). Metabolism regulated by oncogenes allow tumor cells to survive and proliferate in the tumor microenvironment (11). In fact, metabolic reprogramming is a well-established hallmark of cancer (12, 13). Many studies suggested that in order to adapt to the growth and proliferation of HCC cells, the aberrant metabolism of cells develops, which is related to the prognosis of patients (14–16). Strikingly, the functional importance of metabolic alterations often diverges on tumor subtypes, leading to visible therapeutic vulnerability discrepancies in cancer therapy (17, 18). However, metabolic heterogeneity within different HCC subtypes defined by distinct metabolic pathways, which may further result in differences in oncogenes and tumor immunity, has not been well implemented.

In this study, we classified HCC into two different clusters by metabolism-related pathways profiling: Metabolism_H and Metabolism_L. Then, we explored the relationship between the classification and mutation of oncogenes and tumor immunity. Differentially expressed metabolic genes (DEMGs) were identified according to metabolism status and the DEMG-based survival-predictor model was also developed for predicting survival rates of HCC patients, as shown in Figure 1. Moreover, immunohistochemistry staining showed that, compared with normal tissues, 8 of the 11 genes were differentially expressed in cancer tissues, while 3 genes revealed no significant differences. Most of these differentially expressed genes (DEGs) were consistent with our prognostic model, which further verified the reliability of the model. Therefore, the DEMG-based survival-predictor model might have the great potential to predict survival rates of HCC patients.

FIGURE 1

Figure 1 The workflow of this work.

Materials and Methods

Patient Datasets

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.

Clustering

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.

Multi-Omics Analysis

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.

Evaluation of Immune Cell Infiltration Level, Tumor Purity, and Stromal Content Between the Two Clusters

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.

Assessment of Tumor-Infiltrating Immune Cells

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.

Identification of Metabolism-Related Pathways–Immune-Related Pathways Networks

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/).

Gene Set Enrichment and Functional Enrichment Analyses

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).

DEMG-Based Classifiers for Overall Survival

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.

Predictive Performance of the DEMG-Based Classifiers

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.

Immunohistochemical Staining

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).

Data Analysis

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.

Results

Metabolism-Related Pathways Profiling Identified Two HCC Clusters

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).

FIGURE 2

Figure 2 Metabolism-related pathways profiling identified two HCC clusters. (A) In TCGA, univariate Cox analyses showed that there were 14 significant metabolic pathways contributing to OS in HCC. (B) In LIRI_JP, there were 19 significant pathways. (C) A Venn diagram showed 10 metabolism-related pathways shared by TCGA and LIRI_JP. (D, E) Based on the different gene expressions of the 10 pathways, HCC was clustered: Metabolism_H and Metabolism_L in TCGA. (G, H) HCC was also clustered: Metabolism_H and Metabolism_L in LIRI_JP; Kaplan–Meier survival analysis results indicated that the two clusters had significantly different survival rates in both (F) TCGA (log-rank test p-value = 6.02e−05) and (I) LIRI_JP (log-rank test p-value = 2.527e−04). HCC, hepatocellular carcinoma; OS, overall survival.

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).

TABLE 1

Table 1 Correlations between risk score of the metabolism-related pathways classifier with overall survival and clinicopathological characteristics in the TCGA-LIHC cohort.

The Mutant Oncogenes, Copy Number Variation, and DNA Methylation Analysis

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).

FIGURE 3

Figure 3 Different mutant oncogenes between HCC Metabolism_H and Metabolism_L. (A) Waterfall plots suggested that TP53, TTN, and MUC16 were main mutant oncogenes in HCC patients, but they were different between (B) Metabolism_H and (C) Metabolism_L. (D) Metabolism_H was more likely to have TP53 mutations (Student’s t-test, p < 0.001). (E) The mutated TP53 gene was less expressed than the wild type (Wilcoxon test, p = 0.037). HCC, hepatocellular carcinoma. *p < 0.05, **p < 0.01, ***p < 0.001.

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).

FIGURE 4

Figure 4 The DNA methylation level analysis. (A) The heatmap of 30 methylation driving genes associated with prognosis. (B) K–M survival analysis of HCC patients with hypermethylated or hypomethylated PDK4 and TMEM165.

Immunological Evaluation of the Metabolism-Related Pathway-Based HCC Subtypes

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).

FIGURE 5

Figure 5 Distribution of immune cells in Metabolism_H and Metabolism_L. (A) In TCGA, Metabolism_L contained more naive B cells, gamma delta T cells, and resting mast cells (Wilcoxon test), the same as in (B) LIRI−JP. (C) A heatmap showed that different clusters led to significantly different gene expressions in immune pathways. (D) A metabolism-related pathways–immune-related pathways network indicated that metabolic pathways were associated with tumor immunity.

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).

GSEA-Based KEGG Analysis and GO Analysis

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.

FIGURE 6

Figure 6 Gene set enrichment and functional enrichment analyses. (A) KEGG pathways enriched in Metabolism_H and Metabolism_L, respectively. (B) A Venn diagram showed 135 DEMGs were overlapped among TCGA_ DEGs, LIRI_JP_DEGs, and MRGs. (C) KEGG pathways in Immunity_H and Immunity_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).

Prognostic Value of DEMGs

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).

FIGURE 7

Figure 7 The survival-predictor model based on eleven DEMGs. Univariate Cox analyses showed that 60 DEMGs contributed to the OS in (A) TCGA and (B) LIRI_JP, respectively. (C, D) The LASSO regression model identified the 11 most accurate predictive DEMGs in TCGA. (E) HCC patients were divided into two groups by the median risk score cutoff point in TCGA: high risk and low risk. (F) According to the same cutoff point, HCC patients were also divided into two groups in LIRI_JP; Kaplan–Meier Survival analysis results indicated that the two groups had significantly different survival rates in both (G) TCGA (log-rank test p-value = 1.217e−06) and (H) LIRI_JP (log-rank test p-value = 9.355e−10). DEMGs, differentially expressed metabolic genes; HCC, hepatocellular carcinoma.

Differential Expression Levels of the Eleven DEMGs

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.

FIGURE 8

Figure 8 The differential expression of 8 genes in cancer and matched normal tissues. (A) Representative image of immunohistochemistry staining and (B–I) immunohistochemical staining scores. (× 200, scale bar = 100 µm). Data are presented as the means ± SD. *p < 0.05 and **p < 0.01.

The DEMGs-Based Risk Score Outperforms Other HCC Prognostic Factors

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).

FIGURE 9

Figure 9 Comparison of prediction accuracy among the classifiers and other factors. The (A) univariate Cox analysis and (B) multivariate Cox analysis of DEMG-based risk score and other factors in TCGA. (C) The DEMG-based risk score was more accurate with an AUC of 0.767 in TCGA. The (D) univariate Cox analysis and (E) multivariate Cox analysis of risk score and other factors in LIRI_JP. (F) The DEMG-based risk score was more accurate with an AUC of 0.785 in LIRI_JP. DEMGs, differentially expressed metabolic genes; AUC, area under the curve.

Discussion

As the center of human metabolism, the liver is engaged in the metabolic interchange of water-soluble and lipid metabolites all the time; no other organ can match the metabolic rate of the liver (30). Moreover, the biological processes of toxic substance decomposition and metabolism processes, P450 pathway, p53 pathway, and alcohol dehydrogenase activity have been reported to be related to HCC in previous studies (31–33). Thus, metabolism pathways in the liver may provide a way to predict the prognosis of HCC (34).

In this study, we divided HCC patients into two groups by analyzing the activation of metabolic pathways (Figures 2E, F). The two groups mainly showed significant differences in expressions in 10 metabolic pathways, such as purine metabolism, pyrimidine metabolism, and fructose and mannose metabolism. Tumor cells generate energy by glycolysis despite the presence of sufficient oxygen to support the proliferation and differentiation of cancer cells, which is called the Warburg effect (35). Glutamines, which provide the nitrogen required for the biosynthesis of purine and pyrimidine nucleotides, are also needed in the survival and growth of tumor cells (35). These studies are consistent with our conclusions, suggesting that the Metabolism_L subtype with a lower expression of purine metabolic, pyrimidine metabolic, fructose and mannose metabolic, and amino sugar and nucleotide sugar metabolic pathways will likely have a better survival prognosis (Figures 2F, I). The two subtypes also showed significant differences on taurine and hypotaurine metabolism, retinol metabolism, and fatty acid metabolism. Free fatty acid was confirmed as an independent risk factor for cancer (36), and statin can improve survival outcomes and increases overall survival (37–39). Furthermore, the retinol and retinal are also associated with the diagnosis and prognosis of HCC (16). In brief, the metabolic pathways selected in the current study were all related to survival and proliferation of tumors.

Then, we analyzed mutated genes and tumor immunity in two subtypes, respectively. Figure 3 shows that the mutation rate of the TP53 gene was significantly different in Metabolism_H and Metabolism_L, and Metabolism_H expressed more TP53. The wild-type TP53 protein plays an important role in apoptosis after DNA damage and in cell cycle regulation (40). However, the mutant TP53 protein loses its wild-type function and accumulates in the nucleus, which is considered to be a highly specific marker of malignant tumors (41). Similarly, TP53 also plays an important role in HCC (42, 43). Mutant TP53 proteins simultaneously lose their tumor-suppressive functions and obtain new capacities to advance tumorigenesis (44). Ling et al. (45) indicated that HCC patients with non-functional mutant genes of TP53 tended to have a worse survival prognosis. CNV analysis suggested that the two HCC subtypes were different in chromosomal structural variation. As we have mentioned, several genes promoting the development of liver cancer, such as SCAMP3 and CCT3 (46), amplified more frequently in Metabolism_H. However, some gene expressions were negatively correlated with methylation. For instance, as shown in Figure 4, the levels of methylation of the MRG PDK4 were higher in Metabolism_H, and low expression of PDK4 promoted proliferation and metastasis of HCC (47). In conclusion, these findings explain why survival rates in Metabolism_H are lower. Interestingly, the two subtypes also showed significant differences in tumor immunity, and the immune-related pathways interacted with metabolic pathways. Some studies have proved that some immune cells associate with the prognosis of HCC, such as CD8+ T cells, regulatory T cells (Tregs), and B cells (48). HCC produced lactic and carbonic acids excessively by exacerbating glycolysis to change the tumor immunity microenvironment (49). Figure S1 shows that Metabolism_L contained more tumor cells and less immune cells compared with Metabolism_H. Recent studies reported that immune cell infiltration could affect the prognosis of HCC and intratumoral infiltration by dendritic cells and neutrophils may result in poor prognosis in HCC patients (50, 51). Neutrophils not only were involved in the activation and regulation of immune cells, but also promoted the progression of HCC by releasing cytokines (52). Moreover, macrophages, especially M2 macrophages, contributed to the poor prognosis of HCC (53), whose infiltration within the tumor microenvironment could facilitate tumor growth, angiogenesis, invasion, as well as metastasis (54). According to our results, tumor tissue of Metabolism_H contained more M2 macrophages and neutrophils. Conversely, tumor tissue of Metabolism_H contained less CD8 T cells, which were the primary cytotoxic lymphocytes exerting antitumor effects (55). Although p-values of M2 macrophages, neutrophils, and CD8 T cells were less than 0.05 in only one of the two databases, which may be due to insufficient data, the difference in the two immune cells was consistent between the two databases. Moreover, a lower proportion of CD8 T cells indicated an immune-suppressive state in Metabolism_H, and the Metabolism_H subtype had higher mRNAsi and higher expressions of immune checkpoint-related genes such as PD-L1 (56, 57). Therefore, it was no surprise that the Metabolism_H had a lower survival rate.

According to the LASSO regression method, we determined 11 DEMGs: RRM2, PYCR1, PTGES, LPCAT1, HK2, G6PD, CYP3A5, CYP2C9, BDH1, ADH4, and ADH1C. Reports have indicated that most of the 11 DEMGs were closely related to the OS of tumor (58–61). Among the 11 DEMGs, PYCR1 plays a vital part in the promotion of HCC cell proliferation by increasing proline biosynthesis effectively (62). Additionally, LPCAT1 participates in cell proliferation, migration, and invasion by modulating phospholipid composition, in HCC (63). In our study, we divided the HCC patients into two groups based on the 11-DEMG-based classification in TCGA. Then, we verified the correctness of this grouping method in TCGA. The two groups all showed significantly different survival rates in TCGA and LIRI-JP. Moreover, the differential expression of the 11 genes in cancer and matched normal tissues was observed by immunohistochemistry staining. Compared with normal tissues, 8 genes showed remarkable differential expression between cancer tissues and normal tissues, while 3 genes revealed no significant differences. The limitations may come from differential expression patterns of genes, population differences, or statistical noise, and require more experiments to be verified. In addition, when compared with the clinicopathological risk factors, the 11-DEMG-based risk score was better at predicting survival in both TCGA and LIRI-JP, which was the highlight of this study. There is absolutely no doubt that our 11-DEMG-based classifiers possessed their own unique prediction. When the classifiers are combined with clinicopathological risk factors, it would provide a more accurate prediction for OS at different times for HCC patients. Therefore, the DEMG-based survival-predictor model has shown a favorable effect on survival prediction, which will contribute to therapeutic decision-making.

However, there are several limitations in this study. Firstly, we found that tumor metabolism was associated with tumor immunity, but, regrettably, this study mainly focused on the association between the MRGs and the OS of HCC. It will be interesting to combine metabolic genes with immune genes to predict HCC OS in the future. Secondly, this study was a retrospective study utilizing the TCGA and LIRI-JP databases. Therefore, more prospective studies were still needed. Third, if we could discover tumor biomarker detection in a more accessible blood sample, it would be more clinically valuable. Finally, a study at the single-cell level would be better in entangling the heterogeneity among the cells (64–66), which will be the subject of a future work.

Conclusions

In summary, we identified two metabolism-based classifiers associated with OS in HCC and confirmed that the differences in survival rates in the two clusters may be related to mutated genes and tumor immunology. According to the LASSO regression method, we determined 11 DEMGs. Notably, the DEMG-based survival-predictor model could accurately predict the OS of HCC patients, and the results may contribute to the development of individual therapy.

Data Availability Statement

The datasets used during the current study are available from the corresponding author on reasonable request.

Ethics Statement

This study was reviewed and approved by the Medical Ethical Committee of the First Affiliated Hospital of Wenzhou Medical University (2021-R084).

Author Contributions

ZZ and YS were responsible for study conception and design. TY, LL, WH, and SW were responsible for data collection and analysis. TY, LC, and WH were responsible for drafting the manuscript. TY and ZZ were responsible for revision of the manuscript. All authors contributed to the article and approved the submitted version.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors would like to thank the investigators and staff who contributed to this study.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2022.863266/full#supplementary-material

Supplementary Figure 1 | Analysis of Copy number variation. Copy number variation in (A) Metabolism_H and (B) Metabolism_L; (C) Focal copy number alterations in several genes.

Supplementary Figure 2 | Different tumor immunoinfiltration between Metabolism_H and Metabolism_L. (A) the two clusters showed significantly different Stromal Score in LIRI−JP (Kruskal–Wallis test), while no significantly different Stromal Score in TCGA (Kruskal–Wallis test); (B) Immune Score was higher in Metabolism_H (Kruskal–Wallis test) both in TCGA and LIRI−JP; (C) Tumor Purity was higher in Metabolism_L (Kruskal–Wallis test) both in TCGA and LIRI−JP; (D) ESTIMATE Score was higher in Metabolism_H (Kruskal–Wallis test) both in TCGA and LIRI−JP. *P < 0.05, **P < 0.01, ***P < 0.001.

Supplementary Figure 3 | The immune checkpoints related gene expression and the cancer stemness analyses. (A) PD-L1, CTLA-4 and TIM-3 expressed more in Metabolism_H; (B) mRNA expression-based stemness index (mRNAsi) was higher in the Metabolism_H. ***P < 0.001.

Supplementary Figure 4 | The expression of 3 genes in cancer and matched normal tissues. (A) Representative image of immunohistochemistry staining and (B–D) immunohistochemical staining scores. (× 200, Scale bar = 100 µm). Data are presented as the means ± SD. ns, not significant.

Abbreviations

HCC, Hepatocellular carcinoma; TCGA, The cancer genome atlas; ICGC, International cancer genome consortium; DEMGs, Differentially expressed metabolic genes; LASSO, Least absolute shrinkage and selection operator; tdROC, Time-dependent receiver operating characteristics; OS, Overall Survival; CNV, Copy number variation; TACE, Transarterial chemotherapy with embolization; ssGSEA, Single-sample gene-set enrichment analysis; AUC, Area under the curve; IHC, Immunohistochemistry.

References

1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin (2021) 71:209–49. doi: 10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Liao Z, Wang X, Zeng Y, Zou Q. Identification of DEP Domain-Containing Proteins by a Machine Learning Method and Experimental Analysis of Their Expression in Human HCC Tissues. Sci Rep (2016) 6:39655. doi: 10.1038/srep39655

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Anwanwan D, Singh SK, Singh S, Saikam V, Singh R. Challenges in Liver Cancer and Possible Treatment Approaches. Biochim Biophys Acta Rev Cancer (2020) 1873:188314. doi: 10.1016/j.bbcan.2019.188314

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Huo Q, Ma Y, Yin Y, Qin G. Biomarker Identification for Liver Hepatocellular Carcinoma and Cholangiocarcinoma Based on Gene Regulatory Network Analysis. Curr Bioinf (2021) 16:31–43. doi: 10.2174/1574893615666200317115609

CrossRef Full Text | Google Scholar

5. Villanueva A. Hepatocellular Carcinoma. N Engl J Med (2019) 380:1450–62. doi: 10.1056/NEJMra1713263

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Liu H, Qiu C, Wang B, Bing P, Tian G, Zhang X, et al. Evaluating DNA Methylation, Gene Expression, Somatic Mutation, and Their Combinations in Inferring Tumor Tissue-Of-Origin. Front Cell Dev Biol (2021) 9:619330. doi: 10.3389/fcell.2021.619330

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Su XF, Li N, Meng FL, Chu YL, Li T, Gao XZ. MiR-16 Inhibits Hepatocellular Carcinoma Progression by Targeting FEAT Through NF-kappaB Signaling Pathway. Eur Rev Med Pharmacol Sci (2019) 23:10274–82. doi: 10.26355/eurrev_201912_19665

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Oral A, Sahin T. Prognostic Role of ABO Blood Group and Rhesus Factor in Cirrhotic Patients With Hepatocellular Carcinoma. Sci Rep (2019) 9:19087. doi: 10.1038/s41598-019-55685-8

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Llovet JM, Kelley RK, Villanueva A, Singal AG, Pikarsky E, Roayaie S, et al. Hepatocellular Carcinoma. Nat Rev Dis Primers (2021) 7:6. doi: 10.1038/s41572-020-00240-3

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Craig AJ, von Felden J, Garcia-Lezana T, Sarcognato S, Villanueva A. Tumour Evolution in Hepatocellular Carcinoma. Nat Rev Gastroenterol Hepatol (2020) 17:139–52. doi: 10.1038/s41575-019-0229-4

PubMed Abstract | CrossRef Full Text | Google Scholar

11. DeBerardinis RJ, Lum JJ, Hatzivassiliou G, Thompson CB. The Biology of Cancer: Metabolic Reprogramming Fuels Cell Growth and Proliferation. Cell Metab (2008) 7:11–20. doi: 10.1016/j.cmet.2007.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Faubert B, Solmonson A, DeBerardinis RJ. Metabolic Reprogramming and Cancer Progression. Sci (2020) 368:eaaw5473. doi: 10.1126/science.aaw5473

CrossRef Full Text | Google Scholar

13. Tang W, Wan S, Yang Z, Teschendorff AE, Zou Q. Tumor Origin Detection With Tissue-Specific miRNA and DNA Methylation Markers. Bioinformatics (2018) 34:398–406. doi: 10.1093/bioinformatics/btx622

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Vander Heiden MG, Cantley LC, Thompson CB. Understanding the Warburg Effect: The Metabolic Requirements of Cell Proliferation. Science (2009) 324:1029–33. doi: 10.1126/science.1160809

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Hu B, Lin JZ, Yang XB, Sang XT. Aberrant Lipid Metabolism in Hepatocellular Carcinoma Cells as Well as Immune Microenvironment: A Review. Cell Prolif (2020) 53:e12772. doi: 10.1111/cpr.12772

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Han J, Han ML, Xing H, Li ZL, Yuan DY, Wu H, et al. Tissue and Serum Metabolomic Phenotyping for Diagnosis and Prognosis of Hepatocellular Carcinoma. Int J Cancer (2020) 146:1741–53. doi: 10.1002/ijc.32599

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Kim J, DeBerardinis RJ. Mechanisms and Implications of Metabolic Heterogeneity in Cancer. Cell Metab (2019) 30:434–46. doi: 10.1016/j.cmet.2019.08.013

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Peng X, Chen Z, Farshidfar F, Xu X, Lorenzi PL, Wang Y, et al. Molecular Characterization and Clinical Relevance of Metabolic Expression Subtypes in Human Cancers. Cell Rep (2018) 23:255–69.e4. doi: 10.1016/j.celrep.2018.03.077

CrossRef Full Text | Google Scholar

19. Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring Tumour Purity and Stromal and Immune Cell Admixture From Expression Data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust Enumeration of Cell Subsets From Tissue Expression Profiles. Nat Methods (2015) 12:453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-Wide Expression Profiles. Proc Natl Acad Sci USA (2005) 102:15545–50. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Yu G, Wang LG, Han Y, He QY. Clusterprofiler: An R Package for Comparing Biological Themes Among Gene Clusters. OMICS (2012) 16:284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Ito K, Murphy D. Application of Ggplot2 to Pharmacometric Graphics. CPT Pharmacometrics Syst Pharmacol (2013) 2:e79. doi: 10.1038/psp.2013.56

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Walter W, Sanchez-Cabo F, Ricote M. GOplot: An R Package for Visually Combining Expression Data With Functional Analysis. Bioinformatics (2015) 31:2912–4. doi: 10.1093/bioinformatics/btv300

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Zhang J, Shi K, Huang W, Weng W, Zhang Z, Guo Y, et al. The DNA Methylation Profile of non-Coding RNAs Improves Prognosis Prediction for Pancreatic Adenocarcinoma. Cancer Cell Int (2019) 19:107. doi: 10.1186/s12935-019-0828-8

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Chen C, Sun X, Xie W, Chen S, Hu Y, Xing D, et al. Opposing Biological Functions of the Cytoplasm and Nucleus DAXX Modified by SUMO-2/3 in Gastric Cancer. Cell Death Dis (2020) 11:514. doi: 10.1038/s41419-020-2718-3

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Zheng L, Li S, Zheng X, Guo R, Qu W. AHNAK2 is a Novel Prognostic Marker and Correlates With Immune Infiltration in Papillary Thyroid Cancer: Evidence From Integrated Analysis. Int Immunopharmacol (2021) 90:107185. doi: 10.1016/j.intimp.2020.107185

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Zhou M, Zhang Z, Bao S, Hou P, Yan C, Su J, et al. Computational Recognition of lncRNA Signature of Tumor-Infiltrating B Lymphocytes With Potential Implications in Prognosis and Immunotherapy of Bladder Cancer. Brief Bioinform 22 (2021) 22:bbaa047. doi: 10.1093/bib/bbaa047

CrossRef Full Text | Google Scholar

29. Malta TM, Sokolov A, Gentles AJ, Burzykowski T, Poisson L, Weinstein JN, et al. Machine Learning Identifies Stemness Features Associated With Oncogenic Dedifferentiation. Cell (2018) 173:338–54.e15. doi: 10.1016/j.cell.2018.03.034

CrossRef Full Text | Google Scholar

30. Holmes E, Wijeyesekera A, Taylor-Robinson SD, Nicholson JK. The Promise of Metabolic Phenotyping in Gastroenterology and Hepatology. Nat Rev Gastroenterol Hepatol (2015) 12:458–71. doi: 10.1038/nrgastro.2015.114

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Elizondo-Alvarez MA, Uribe-Salas A, Bello-Teodoro S. Chemical Stability of Xanthates, Dithiophosphinates and Hydroxamic Acids in Aqueous Solutions and Their Environmental Implications. Ecotoxicol Environ Saf (2021) 207:111509. doi: 10.1016/j.ecoenv.2020.111509

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Malla S, Kadimisetty K, Jiang D, Choudhary D, Rusling JF. Pathways of Metabolite-Related Damage to a Synthetic P53 Gene Exon 7 Oligonucleotide Using Magnetic Enzyme Bioreactor Beads and LC-MS/MS Sequencing. Biochemistry (2018) 57:3883–93. doi: 10.1021/acs.biochem.8b00271

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Vaughn MB, Zhang J, Spiro TG, Dyer RB, Klinman JP. Activity-Related Microsecond Dynamics Revealed by Temperature-Jump Forster Resonance Energy Transfer Measurements on Thermophilic Alcohol Dehydrogenase. J Am Chem Soc (2018) 140:900–3. doi: 10.1021/jacs.7b12369

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Zhu GQ, Zhou YJ, Qiu LX, Wang B, Yang Y, Liao WT, et al. Prognostic Alternative mRNA Splicing Signature in Hepatocellular Carcinoma: A Study Based on Large-Scale Sequencing Data. Carcinogenesis (2019) 40:1077–85. doi: 10.1093/carcin/bgz073

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Pavlova NN, Thompson CB. The Emerging Hallmarks of Cancer Metabolism. Cell Metab (2016) 23:27–47. doi: 10.1016/j.cmet.2015.12.006

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Zhang L, Han L, He J, Lv J, Pan R, Lv T. A High Serum-Free Fatty Acid Level is Associated With Cancer. J Cancer Res Clin Oncol (2019) 146:705–10. doi: 10.1007/s00432-019-03095-8

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Yokomichi H, Nagai A, Hirata M, Tamakoshi A, Kiyohara Y, Kamatani Y, et al. Statin Use and All-Cause and Cancer Mortality: BioBank Japan Cohort. J Epidemiol (2017) 27:S84–91. doi: 10.1016/j.je.2016.12.011

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Nayan M, Punjani N, Juurlink DN, Finelli A, Austin PC, Kulkarni GS, et al. Statin Use and Kidney Cancer Survival Outcomes: A Systematic Review and Meta-Analysis. Cancer Treat Rev (2017) 52:105–16. doi: 10.1016/j.ctrv.2016.11.009

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Peled Y, Klempfner R, Kassif Y, Kogan A, Maor E, Sternik L, et al. Preoperative Statin Therapy and Heart Transplantation Outcomes. Ann Thorac Surg (2020) 110:1280–5. doi: 10.1016/j.athoracsur.2020.02.005

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Geisinger JM, Stearns T. CRISPR/Cas9 Treatment Causes Extended TP53-Dependent Cell Cycle Arrest in Human Cells. Nucleic Acids Res (2020) 48:9067–81. doi: 10.1093/nar/gkaa603

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Dowell SP, Wilson PO, Derias NW, Lane DP, Hall PA. Clinical Utility of the Immunocytochemical Detection of P53 Protein in Cytological Specimens. Cancer Res (1994) 54:2914–8. doi: 10.1111/j.1365-2303.1994.tb00530.x

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Shimada S, Mogushi K, Akiyama Y, Furuyama T, Watanabe S, Ogura T, et al. Comprehensive Molecular and Immunological Characterization of Hepatocellular Carcinoma. EBioMedicine (2019) 40:457–70. doi: 10.1016/j.ebiom.2018.12.058

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Rao CV, Asch AS, Yamada HY. Frequently Mutated Genes/Pathways and Genomic Instability as Prevention Targets in Liver Cancer. Carcinogenesis (2017) 38:2–11. doi: 10.1093/carcin/bgw118

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Brosh R, Rotter V. When Mutants Gain New Powers: News From the Mutant P53 Field. Nat Rev Cancer (2009) 9:701–13. doi: 10.1038/nrc2693

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Ling S, Shan Q, Zhan Q, Ye Q, Liu P, Xu S, et al. USP22 Promotes Hypoxia-Induced Hepatocellular Carcinoma Stemness by a HIF1alpha/USP22 Positive Feedback Loop Upon TP53 Inactivation. Gut (2019) 69:1322–34. doi: 10.1136/gutjnl-2019-319616

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Skawran B, Steinemann D, Weigmann A, Flemming P, Becker T, Flik J, et al. Gene Expression Profiling in Hepatocellular Carcinoma: Upregulation of Genes in Amplified Chromosome Regions. Mod Pathol (2008) 21:505–16. doi: 10.1038/modpathol.3800998

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Qin YJ, Lin TY, Lin XL, Liu Y, Zhao WT, Li XY, et al. Loss of PDK4 Expression Promotes Proliferation, Tumorigenicity, Motility and Invasion of Hepatocellular Carcinoma Cells. J Cancer (2020) 11:4397–405. doi: 10.7150/jca.43459

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Kurebayashi Y, Ojima H, Tsujikawa H, Kubota N, Maehara J, Abe Y, et al. Landscape of Immune Microenvironment in Hepatocellular Carcinoma and its Additional Impact on Histological and Molecular Classification. Hepatology (2018) 68:1025–41. doi: 10.1002/hep.29904

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Kuchuk O, Tuccitto A, Citterio D, Huber V, Camisaschi C, Milione M, et al. pH Regulators to Target the Tumor Immune Microenvironment in Human Hepatocellular Carcinoma. Oncoimmunology (2018) 7:e1445452. doi: 10.1080/2162402X.2018.1445452

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Zhou ZJ, Xin HY, Li J, Hu ZQ, Luo CB, Zhou SL. Intratumoral Plasmacytoid Dendritic Cells as a Poor Prognostic Factor for Hepatocellular Carcinoma Following Curative Resection. Cancer Immunol Immunother CII (2019) 68:1223–33. doi: 10.1007/s00262-019-02355-3

CrossRef Full Text | Google Scholar

51. Kuang DM, Zhao Q, Wu Y, Peng C, Wang J, Xu Z, et al. Peritumoral Neutrophils Link Inflammatory Response to Disease Progression by Fostering Angiogenesis in Hepatocellular Carcinoma. J Hepatol (2011) 54:948–55. doi: 10.1016/j.jhep.2010.08.041

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Lu C, Rong D, Zhang B, Zheng W, Wang X, Chen Z, et al. Current Perspectives on the Immunosuppressive Tumor Microenvironment in Hepatocellular Carcinoma: Challenges and Opportunities. Mol Cancer (2019) 18:130. doi: 10.1186/s12943-019-1047-6

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Tian Z, Hou X, Liu W, Han Z, Wei L. Macrophages and Hepatocellular Carcinoma. Cell Bioscience (2019) 9:79. doi: 10.1186/s13578-019-0342-7

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Anderson NM, Simon MC. The Tumor Microenvironment. Curr Biol (2020) 30:R921–5. doi: 10.1016/j.cub.2020.06.081

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Ye L, Li Y, Tang H, Liu W, Chen Y, Dai T, et al. Cd8+Cxcr5+T Cells Infiltrating Hepatocellular Carcinomas are Activated and Predictive of a Better Prognosis. Aging (Albany NY) (2019) 11:8879–91. doi: 10.18632/aging.102308

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Ford K, Hanley CJ, Mellone M, Szyndralewiez C, Heitz F, Wiesel P, et al. NOX4 Inhibition Potentiates Immunotherapy by Overcoming Cancer-Associated Fibroblast-Mediated CD8 T-Cell Exclusion From Tumors. Cancer Res (2020) 80:1846–60. doi: 10.1158/0008-5472.CAN-19-3158

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Romano S, Tufano M, D'Arrigo P, Vigorito V, Russo S, Romano MF. Cell Stemness, Epithelial-To-Mesenchymal Transition, and Immunoevasion: Intertwined Aspects in Cancer Metastasis. Semin Cancer Biol (2020) 60:181–90. doi: 10.1016/j.semcancer.2019.08.015

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Wei RR, Zhang MY, Rao HL, Pu HY, Zhang HZ, Wang HY. Identification of ADH4 as a Novel and Potential Prognostic Marker in Hepatocellular Carcinoma. Med Oncol (2012) 29:2737–43. doi: 10.1007/s12032-011-0126-3

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Xu J, Zhang J, Shan F, Wen J, Wang Y. SSTR5AS1 Functions as a ceRNA to Regulate CA2 by Sponging Mir15b5p for the Development and Prognosis of HBVrelated Hepatocellular Carcinoma. Mol Med Rep (2019) 20:5021–31. doi: 10.3892/mmr.2019.10736

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Yu T, Wang X, Zhu G, Han C, Su H, Liao X, et al. The Prognostic Value of Differentially Expressed CYP3A Subfamily Members for Hepatocellular Carcinoma. Cancer Manag Res (2018) 10:1713–26. doi: 10.2147/CMAR.S159425

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Tsunedomi R, Iizuka N, Hamamoto Y, Uchimura S, Miyamoto T, Tamesa T, et al. Patterns of Expression of Cytochrome P450 Genes in Progression of Hepatitis C Virus-Associated Hepatocellular Carcinoma. Int J Oncol (2005) 27:661–7.

PubMed Abstract | Google Scholar

62. Ding Z, Ericksen RE, Escande-Beillard N, Lee QY, Loh A, Denil S, et al. Metabolic Pathway Analyses Identify Proline Biosynthesis Pathway as a Promoter of Liver Tumorigenesis. J Hepatol (2019) 72:725–35. doi: 10.1158/1538-7445.AM2018-3477

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Morita Y, Sakaguchi T, Ikegami K, Goto-Inoue N, Hayasaka T, Hang VT, et al. Lysophosphatidylcholine Acyltransferase 1 Altered Phospholipid Composition and Regulated Hepatoma Progression. J Hepatol (2013) 59:292–9. doi: 10.1016/j.jhep.2013.02.030

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Peng L, Tian X, Tian G, Xu J, Huang X, Weng Y, et al. Single-Cell RNA-Seq Clustering: Datasets, Models, and Algorithms. RNA Biol (2020) 17:765–83. doi: 10.1080/15476286.2020.1728961

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Xu J, Cai L, Liao B, Zhu W, Yang J. CMF-Impute: An Accurate Imputation Tool for Single-Cell RNA-Seq Data. Bioinformatics (2020) 36:3139–47. doi: 10.1093/bioinformatics/btaa109

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Yang J, Liao B, Zhang T, Xu Y. Editorial: Bioinformatics Analysis of Single Cell Sequencing Data and Applications in Precision Medicine. Front Genet (2019) 10:1358. doi: 10.3389/fgene.2019.01358

PubMed Abstract | CrossRef Full Text | Google Scholar