Expression profile and survival data in the GSE14520 dataset were obtained from tissue samples of 247 HCC patients and 239 normal controls (NC) (20). Expression profile in the GSE76427 dataset were obtained from 115 primary tumors tissues and 52 adjacent non-tumor tissues of HCC patients (21). Expression profile in the GSE174570 dataset were obtained from 57 paired tumors tissues and adjacent non-tumor tissues of HCC patients (22). Expression profile and survival data in the GSE141198 dataset were obtained from tumors tissues of 148 HCC patients (23). In addition, expression, survival and somatic mutation data of 373 HCC patients and 50 normal controls were obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). Genes in four metabolic pathways (citrate cycle TCA cycle, fatty acid metabolism, glycerolipid metabolism, and glycolysis/gluconeogenesis) were collected from Molecular Signatures Database (MSigDB) (24). PXD006512 was collected from PRIDE database (25), including proteomic data with 124 paired tumors tissues and adjacent non-tumor tissues of HCC patients. Single-cell RNA sequencing (scRNA-seq) data were collected from GSE149614 for 10 HCC patients with tumor and non-tumor liver (26).

Differential analysis between HCC and controls were assessed using limma package in R (27) for GSE14520, GSE76427, GSE174570, and PXD006512, and using DEseq2 package in R (28) for TCGA, then differentially expressed genes (DEGs) and differentially expressed proteins (DEPs) were obtained with |logFC (log fold change)|> 1 and P < 0.05.

Gene Set Enrichment Analysis (GSEA) gene set (29) was used to evaluate the activation of four metabolic pathways in HCC. The scores of four metabolic pathways were calculated separately with gene set variation analysis (GSVA) (30), and the effect of the median score on patient OS was analyzed by Kaplan-Meier (K-M) curves. Subsequently, DEGs in metabolic pathways with prognostic significance were screened out as candidate genes.

Candidate genes were used to perform consensus cluster analysis using ConsensusClusterPlus package in R (31). HCC samples were therefore clustered as different clusters. The K-M curves were used to analyze and compare the OS of different clusters. DEGs between clusters in GSE14520, GSE141198, and TCGA were obtained with P < 0.05. Somatic mutation was calculated to evaluate the tumor mutation burden (TMB) in different clusters using Maftools (32). The stromal score, immune score, ESTIMATE score, tumor purity, and glycolysis/gluconeogenesis in TCGA were assessed between different clusters using GSVA.

Weighted gene co-expression network analysis (WGCNA) was performed to construct co-expression networks in TCGA for intersecting DEGs between clusters in GSE14520, GSE141198, and TCGA using WGCNA package in R (33). Correlations between pairs of genes were first calculated using gene expression profiles and transformed into adjacency matrices. Then the optimal β was set such that the connections between genes in the network obeyed a scale-free network distribution, and the adjacency matrix was transformed into a topological overlap matrix. Hierarchical clustering trees were subsequently constructed, with different branches (colors) representing different modules. Correlations between modules and clinical trait were calculated by Pearson correlation.

Enrichment analysis was used to identify the Gene Ontology (GO) analysis, and Kyoto Encyclopedia of Genes and Genomes (KEGG) in which module genes were involved using clusterProfiler package in R (34). Significant pathways were determined by P < 0.05.

CIBERSORT (35) was used to determine the immune cell infiltration based on the gene expression data in TCGA. The proportion of infiltrating immune cells was estimated by LM22 signal between HCC and control or between different clusters. Tumor Immune Dysfunction and Exclusion (TIDE; http://tide.dfci.harvard.edu/) was utilized to predict the responsiveness of samples in different clusters to immunotherapy.

Univariate Cox regression was used to analyze the prognostic role of candidate genes in GSE14520, GSE141198, and TCGA. Forest plot was plotted with hazard ratio (HR) and 95% confidence interval (CI). Genes that significantly affected patients’ overall survival were identified as prognostic genes. Nomogram was generated using rms package in R. The calibration curves were established to illustrate the agreement between the nomogram-predicted and the observed probabilities of HCC. Then a risk score was calculated and the HCC samples in TCGA were divided into high- and low-risk groups according to the median risk score. OS in different groups was predicted by K-M curve. The receiver operating characteristic (ROC) curves was generated by survivalROC package in R to evaluate time-dependent OS of HCC patients.

Multivariate Cox regression was used to analyze the prognostic role of candidate genes in protein level of PXD006512. Proteins that significantly affected patients’ OS were used to calculate a risk score. The HCC samples in PXD006512 were divided into high- and low-risk groups according to the median risk score. OS in different groups was predicted by K-M curve. The ROC curves were also generated by survival ROC package in R to evaluate time-dependent OS of HCC patients.

To filter out low-quality cells, cells with more than 8000 or fewer than 200 expressed genes were removed. Single cell sequencing data were normalized using a standardized data algorithm, and variable genes were filtered using the FindVariableFeatures function. Highly variable genes of top 2000 were selected to perform clustering analysis using unified manifold approximation and projection (UMAP). Marker genes for cell types were collected from known cell specific marker genes.

A total of 30 paired tumors tissues and adjacent normal tissues from HCC patients were collected in the First Affiliated Hospital of Xinjiang Medical University. Peripheral blood samples of 10 HCC patients and 10 healthy volunteers were also collected in the First Affiliated Hospital of Xinjiang Medical University. Samples in this study were obtained with approval by Ethics Committee of the First Affiliated Hospital of Xinjiang Medical University (NO. K202304-20), and consent were obtained for all participants.

Tumors tissues and normal tissues samples were used to extract total RNA using was TRIzol reagent (Invitrogen, CA, USA). Complementary DNA was obtained through reverse-transcribe using total RNA with PrimeScript™ RT reagent Kit (Takara, Dalian, China). Expression of target genes were detected using RT-qPCR with SYBR^® Premix Ex Taq™ II kit (Takara). β-actin is used as an internal reference gene to calculate the relative expression level of genes using 2^−ΔΔCt method. Primers used in this study is shown in Table S1 .

Tumors tissues and normal tissues samples were homogenated in RIPA lysis buffer with PMSF on ice. Proteins were extracted and then quantified with BCA protein assay kit (Beyotime, Shanghai, China). Then 20 µg proteins were separated on SDS–polyacrylamide gel electrophoresis and transferred onto polyvinylidene fluoride membranes. After blocking with 5% skim milk at room temperature for 2 h, membranes were incubated with primary antibodies (anti-ADH1B, anti-ALDOB, anti-ADH1A, anti-FBP1, anti-ADH6, and anti-β-actin; ABclonal Technology, Wuhan, China) at 4°C for overnight, respectively. Subsequently, membranes were incubated with HRP-linked secondary antibodies and detected with an ECL chemiluminescence kit (Beyotime). Quantification of proteins was performed by normalized to β-actin using ImageJ software.

For immunohistochemical (IHC) staining, tissues were embedded in paraffin after fixed in 4% paraformaldehyde and sectioned at 4 μM. After antigen retrieval at high temperature, sections were incubated with sheep serum albumin for blocking antigen. Sections were then incubated with primary antibodies (anti-ADH1B, anti-ALDOB, anti-ADH1A, anti-FBP1, and anti-ADH6; ABclonal Technology) for overnight. Secondary antibody was applied for 30 min. After added the diaminobenzidine solution, sections were stained with hematoxylin. The images were visualized by XSP-C204 biomicroscope.

Additionally, sections were stained with anti-mouse CD14 monoclonal antibody (Proteintech, Wuhan, China) for 15 h, then incubated with FITC-conjugated goat anti-mouse IgG (Proteintech) for 1 h. After washing with solution for 4 times, sections were incubated with sheep serum albumin for blocking antigen. Then sections were stained with anti-rat FOXP3 polyclonal antibody (Absin, Shanghai, China) for 15 h, then incubated with Cy3-conjugated goat anti-rat IgG (Proteintech) for 1 h. Nuclei were counterstained with DAPI. The images were visualized by XSP-C204 biomicroscope.

Peripheral blood samples of HCC patients and normal controls were collected and incubated with antibodies (BD Biosciences, CA, USA), including anti-CD14-ECD, anti-CD4-APC, anti-CD45RO-PE, anti-CD25-FITC, anti-CD127-APC, anti-CXCR5-PC5.5, and anti-CD68-PC7. After incubation for 18 min, cells were incubated with red blood cell lysate (BD Bioscience). After washing with PBS, cell population was gated and analyzed by BD FACSFortessa (BD Biosciences). Data were analyzed by Kaluza (v2.0).

Single‐cell‐scaled time‐of‐flight (CyTOF) mass cytometry data was collected from Mendeley Database (https://doi.org/10.17632/jxsz3hdsyg.2) with accession numbers: CRA001276 (36). Which including thirteen groups of tumor (T), and normal (N) specimens from HCC patients. Cluster analysis was performed on CD45 positive cells using t-distributed stochastic neighbor embedding (tSNE), and subsequent manual gating was performed using FlowJo v10.5.3. Treg was identified using CD4+CD3+CD25+CD127low, monocyte was identified using CD3-CD16-CD4+HLA-DR+CD45RO+CD11a+CD49d+.

R version 3.6.1 was used for bioinformatics analysis. GraphPad Prism 7.0 was used for statistical analysis and graphing. Data of experiment are expressed as mean ± SD for triplicate experiments. Statistical test was performed using Student’s t test. P < 0.05 was considered statistically significant.

To evaluate the metabolic pathways in HCC, GSEA was performed. The results showed that citrate cycle TCA cycle, fatty acid metabolism, glycerolipid metabolism, and glycolysis/gluconeogenesis were all activated in HCC in TCGA ( Figure 2A ), GSE14520 ( Figure 2B ), GSE76427 ( Figure 2C ), GSE174570 ( Figure 2D ). K-M curves showed that patients with high-score of glycolysis/gluconeogenesis had better OS compared to low-score of glycolysis/gluconeogenesis ( Figures 2E, F, G ). Unfortunately, the other three metabolic pathway scores did not significantly affect the OS of HCC patients.

To identify differentially expressed glycolysis/gluconeogenesis‐related genes, the DEGs between HCC and controls were identified. A total of 2493 DEGs in TCGA ( Figure 3A ), 2958 DEGs in GSE14520 ( Figure 3B ), 3927 DEGs in GSE76427 ( Figure 3C ), and 583 DEGs in GSE174570 ( Figure 3D ). There were 62 glycolysis/gluconeogenesis‐related genes were obtained and intersection analysis revealed 13 differentially expressed glycolysis/gluconeogenesis‐related genes in HCC to be considered as candidate genes ( Figure 3E ). The expression of candidate genes in paired tumors tissues and adjacent non-tumor tissues of GSE174570 showed that BPGM was significantly higher expressed in HCC than in controls, and HK3, ENO3, ALDH1B1, ALDH9A1, ADH6, ADH1A, ADH1B, PCK1, ALDOB, FBP1, ALDH2, and PCK2 were significantly lower expressed in HCC ( Figure 3F ). Importantly, the aberrant expression of candidate genes was validated in additional data ( Figure S1 ).

Based on the expression of 13 candidate genes, consensus cluster analysis was performed to discriminate HCC patients in TCGA. The greatest increase in the area under the cumulative distribution function (CDF) curves at k = 2 resulted in the best clustering ( Figures 4A–C ). Therefore, we obtained two clusters: C1 group and C2 group. The expression of 13 candidate genes in C1 and C2 was shown in Figure 4D . Patients in C1 group had a significantly worse OS than those in C2 group ( Figure 4E ). Interestingly, the same clustering results were also obtained in the GSE14520 and GSE141198 datasets, with patients in C1 group had a worse prognosis than those in C2 ( Figure S2 ). Additionally, TP53 was found to be mutated most frequently in C1 and CTNNB1 in C2 among somatic mutations ( Figure S3 ).

Next, the differentially expressed genes between C1 and C2 groups were identified. There were 9999 differentially expressed genes in TCGA ( Figure S4A ), 5870 differentially expressed genes in GSE14520 ( Figure S4B ), 9885 differentially expressed genes in GSE141198 ( Figure S4C ). A total of 914 differentially expressed genes were the intersection of three datasets ( Figure 5A ). Network topology analysis of soft threshold power reveals β=6 was the optimal value to construct the co-expression network ( Figure 5B ). Then seven modules were obtained ( Figure 5C ). Correlation analysis showed that brown module was greatest positive correlation was with C2 and negative correlation with C1 ( Figure 5D ). The enrichment analysis found that brown module genes were mainly involved in organic acid metabolic process, small molecule metabolic process, and carboxylic acid metabolic process of biological processes ( Figure 5E ). In the KEGG pathways, metabolic pathways, valine, leucine and isoleucine degradation, and tryptophan metabolism were mainly enriched by brown module genes ( Figure 5F ).

To investigate whether glycolysis/gluconeogenesis was associated with immune microenvironment, we explored the stromal score, immune score, ESTIMATE score, and tumor purity in C1 and C2 in TCGA. Results showed that immune score, and tumor purity were significantly higher in C1 than C2, while glycolysis/gluconeogenesis was lower in C1 than C2 ( Figure 6A ). Then the abundances of immune cells in each sample were measured using the CIBERSORT ( Figure 6B ). Macrophage M2, and CD4+ memory resting T cells were the most abundance in HCC. By comparing immune cell infiltration between HCC and controls, we found that CD4+ naive T cells, CD4+ memory activated T cells, follicular helper T cells (Tfh), regulatory T cells (Tregs), macrophage M0, myeloid dendritic cell resting, myeloid dendritic cell activated, and mast cell activated were higher infiltrated in HCC, while plasma B cells, gamma delta T cells, monocyte, macrophage M2, mast cell resting, eosinophil, and neutrophil were lower infiltrated in HCC ( Figure 6C ). Besides, memory B cell, CD4+ memory activated T cell, Tfh, Tregs, macrophage M0, myeloid dendritic cell resting, and neutrophil were higher infiltrated in C1, while monocyte, and mast cell activated were lower infiltrated in C1 ( Figure 6D ). CD4+ memory activated T cell, Tfh, Tregs, and macrophage M0 were all higher infiltration in HCC and C1 group.

Furthermore, differences in effects for patients in C1 and C2 receiving checkpoint inhibitors were predicted and found that the C1 group had a higher proportion of potential responders than the C2 group ( Figure 6E ). The responses to immunotherapies of receiving anti-PD-1 or anti-CTLA-4 were compared in C1 and C2 with the SubMap analysis ( Figure 6F ). It was found that HCC patients in C1 may be more sensitive for responses to anti-PD-1 and anti-CTLA-4.

Five intersecting DEGs were identified from candidate genes significantly affected OS in TCGA ( Figure 7A ), GSE14520 ( Figure 7B ), and GSE141198 ( Figure 7C ) by univariate Cox regression analyses. ADH1A, ADH1B, ADH6, ALDOB, and FBP1 as prognostic genes all had a protective role in HCC patient prognosis (HR<1). The nomogram was constructed for OS in HCC patients, and showed promising accuracy in predicting prognoses ( Figure 7D ). The calibration curves showed a robust calibration of nomogram ( Figure 7E ). The results of correlation analysis showed that prognostic genes and Treg or macrophage were negatively correlated, prognostic genes and monocyte was positively correlated ( Figure 7F ).

We then obtained two subtypes (high- and low-risk groups) in TCGA based on median risk score ( Figure 7G ). Prognostic genes were all lowly expressed in the high-risk group and highly expressed in the low-risk group. Area under the ROC curve (AUC) values of median risk score in 12-, 36- and 60-month were 0.66, 0.66 and 0.63, respectively ( Figure 7H ). Importantly, HCC patients in the high-risk group had a worse prognosis than those in the low-risk group ( Figure 7I ).

To determine the prognostic role of candidate genes in protein levels, the proteomic data in PXD006512 was analyzed. There were 5714 DEPs between HCC and controls ( Figure 8A ). Interestingly, the protein levels of the candidate genes were all lower expressed in HCC compared to controls ( Figure 8B ). Univariate Cox regulation analysis confirmed that ADH1A, ADH1B, ADH6, and ALDOB had protective effects on HCC ( Figure 8C ). A risk prognostic model was also established by median risk score prognostic genes to divide HCC samples into high - and low-risk groups ( Figure 8D ). Protein levels of prognostic genes were also lower expressed in high-risk group compared to low-risk group. AUC values of median risk score in 12-, 36- and 60-month were 0.79, 0.78 and 0.81, respectively ( Figure 8E ). HCC patients in the high-risk group had a worse prognosis than those in the low-risk group ( Figure 8F ).

Unsupervised cluster analysis was performed 34 clusters were established through UMAP from cells of all samples ( Figure 9A ). According to specific marker genes, we identified the following 11 major cell types: B cells, CD8+ T cells, CD45-LYZ+ cells, endothelial, epithelial, fibroblast, macrophages, monocytes, natural killer (NK), NK T cells, and Treg ( Figures 9B, C ). Interestingly, NK T cells, and CD45-LYZ+ cells were more expressed in HCC samples, NK more expressed in normal samples ( Figure 9D ). ADH1A, ADH1B, ADH6, and ALDOB were mainly expressed in monocytes, FBP1 was mainly expressed in monocytes and macrophages ( Figure 9E ).

RT-qPCR detection showed that decreased mRNA levels of ADH1B, ALDOB, ADH1A, ADH6, and FBP1 in HCC compared with controls ( Figure 10A ). The expression trend of prognostic genes was also decreased at the protein level in HCC ( Figure 10B ). We further confirmed that the expression of prognostic genes was lower in HCC than in the control group using IHC staining ( Figure 10C ).

The proportion of monocytes and Treg showed increased in HCC by flow cytometry detection ( Figure 10D ). Then, immune cells were visualized and analyzed based on CyTOF data ( Figure 11A ). According to the expression patterns of different immune cell surface markers ( Figure S5 ), monocyte and Treg cells were identified significantly enriched in normal ( Figure 11B ), and tumor ( Figure 11C ). It is found that monocyte and Treg cells were increased from normal to tumor according to the cell density map ( Figure 11D ). Importantly, the multiplex immunohistochemistry staining was performed to detect the expression of marker protein of monocytes (CD14) and Treg (FOXP3). As shown in Figure 11E , the abundance of monocytes and Treg in tumor were higher than that in control group.