The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/) program was launched by the National Cancer Institute (NCI) and the National Human Genome Research Institute (NHGRI), and currently studies 36 cancer types in total (14). 374 LIHC samples and 50 normal samples were included for study, with clinical data including gender, age, stage, grade, survival time, survival status, etc. In addition, the corresponding mutation data was downloaded of hepatocellular carcinoma for subsequent mutation analysis. The ICGC database (https://dcc.icgc.org/projects/LIRI-JP) is an international tumor genome collaborative group, a global collaborative database containing samples from different countries and regions. We downloaded data from it for liver cancer as a validation cohort. The GEO dataset (GSE45267) was used for external validation (https://www.ncbi.nlm.nih.gov/).

WGCNA is to explore whether there is co-expression between genes and to divide a certain cluster of co-expressed genes into a module based on certain values, so that different clusters of genes clustered together are divided into different modules. Metabolism-related genes were used to construct the weighted correlation network, where modules that differed between cancer and paracancer were selected for further analysis (15).

KEGG is a database for studying biological functions from genomic and molecular level information. Its PATHWAY sub-database integrates current knowledge in molecular interaction networks and allows prediction of pathways enriched by differential genes. The Gene Ontology (GO) database divides the functions of genes into three components: cellular component (CC), molecular function (MF), and biological process (BP). To explore the underlying biological processes and signaling pathways associated with the acquisition of differential genes.

A risk score model was constructed based on the integrated role of key genes in liver cancer to comprehensively evaluate the function of these molecules in patient prognosis. The model was construct by lasso-cox regression analysis. The metabolic index (MBI) for each patient with HCC was computed with the equation:

MBI = ∑ (Expression level of Gene i × coefficient i)

The TCGA dataset is the training set and the ICGC is the validation set. Kaplan-Meier curves were used to investigate survival situation between different MBI groups of patients with liver cancer in the training and validation set. Risk-survival curves were used to explore the survival and mortality of the samples in different MBI groups, and the differences in the key genes that were modeled between the two groups. The roc curve is used to determine the effectiveness of this prediction model.

The abundance scores of 16 immune cells and 13 immune-associated terms in each sample were assessed with ssGSEA method. The estimate algorithm was used to calculate the immunescore, stromalscore, and Estimatescore in the tumor microenvironment. The IMvigor210 dataset was used to assess the efficacy of immunotherapy. We used the “oncopredict” package to differentiate the sensitivity of the different groups to the drug, of which a total of 198 chemotherapeutic agents could be evaluated. We screened the most sensitive drugs from them (16).

Human Protein Atlas (HPA) (https://www.proteinatlas.org/) is a freely available public database (17). This database investigates the expression of proteins at the protein level in various human tissues and organs. It contains the IHC profile in normal and tumor tissues. The expression of key genes in normal and tumor tissues was investigated.

Ten pairs of HCC tissues and normal tissues were provided by the First Affiliated Hospital of Anhui Medical University, and the samples were stored at -80°C for a long time. The ethics committee of the First Affiliated Hospital of Anhui Medical University approved the study.

The specific experimental procedure is described previously (18). TRIZOL (#15596026, Invitrogen, USA) was used to extracted RNA according to manufacturer’s instructions. According to protocol, cDNA was derived from RNA by reverse transcription by PrimeScript™ RT Master Mix (Takara Bio, Japan). TB Green^® Premix Ex Taq™ II (Tli RNaseH Plus), Bulk (#RR820L, Takara, Japan) was used to perform RT-qPCR according to the manufacturer’s instructions. RT-qPCR reaction program: Preheat: 95°, 30; Denaturation: 95°, 30s; Protein refolding: 55°, 30s; Cycle: 72°, 60s, 40 times; Termination: 72°, 5min. The primer sequences are as follows: ATIC-F: ACCTGACCGCTCTTGGTTTG, ATIC-R: TACGAGCTAGGATTCCAGCAT; G6PD-F: CGAGGCCGTCACCAAGAAC, G6PD-R: GTAGTGGTCGATGCGGTAGA; GAPDH-F: GGAGCGAGATCCCTCCAAAAT, GAPDH-R: GGCTGTTGTCATACTTCTCATGG.

Analysis of differences was performed using the Wilcoxon test. The correlation analysis was based on the Spearman correlation test. Kaplan-Meier analysis was used for survival analysis. R package “survival” was used to perform Cox regression analysis, along with hazard ratio (HR) and 95% confidence interval (CI). P < 0.05 was regarded as statistically significant. R software (version 4.1.2) performs statistical analysis and plotting.

In order to understand the framework and logic of the article, a flowchart was drawn ( Figure 1 ). 947 metabolism-related genes were collected from the GSEA database, and a total of 362 genes were differential genes after performing differential analysis, as can be clearly seen in the volcano plot (P < 0.05, |log 2 FC| > 1) ( Figure 2A ). After calculation, the optimal soft threshold 8 is used to construct the co-expression network ( Figure 2B ). After dividing the genes into different modules, it is used to draw the gene clustering tree ( Figure 2C ). The clustering was performed by WGCNA, and we can see that the genes were divided into 5 modules, among which the MEbrown module was selected for subsequent analysis, with a total of 90 genes ( Figure 2D ). After intersecting the differential analysis and the genes selected by WGCNA, 52 genes were used for our subsequent analysis ( Figure 2E ).

As previously described, we performed a univariate COX prognostic analysis of these differential genes, 19 of which were statistically significant (P < 0.05) ( Figure 3A ). Next, a further LASSO regression model calculated the optimal lamda value of 2, where ATIC and G6PD were used to construct the model ( Figure 3B ). KM analysis was used to analyze the prognosis of the two genes, and it can be seen from the graph that patients in the high gene expression group had a poorer prognosis ( Figures 3C, D ). We used the TCGA database and applied the model to evaluate every patient’s MBI. Patients were grouped equally to two teams according to MBI, and it was found that high MBI patients had a worse prognosis ( Figure 3E ). The model was further validated using the ICGC dataset and the outcome is in agreement with the above ( Figure 3F ). Further risk curves were then plotted, and both the TCGA database and the ICGC database found higher mortality rates for patients at higher MBI ( Supplementary Figure 1A-D ).

To explore the differences between differential MBI groups, GO and KEGG analysis were performed using differential genes. Among them, BP was mainly enriched to mitotic nuclear division, CC enrichment to chromosomal regions and collagen-containing extracellular matrix, and MF is mainly enriched to lipid transporter activity ( Figures 4A, B ). KEGG analysis revealed that the enriched pathways were mainly cell cycle, chemical carcinogenesis and cellular senescence ( Figures 4C, D ).

Subsequently, we performed a COX prognostic analysis of the clinical characteristics and MBI of these patients. Among them, staging and MBI were statistically significant in both univariate and multivariate analyses ( Figures 5A, B ). Considering the clinical use of the model, we further plotted the nomogram. MBI were included in addition to the clinical information of the patients ( Figure 5C ). The ROC curves analyzed the efficacy of the model at 1, 3 and 5 years and found that all were well evaluated. Considering that the inclusion of clinical information may be more effective for the model, we calculated the Nomogram score. and the ROC curve was found to be better when it was evaluated ( Figures 5D, F ). The calibration curve results found the model to be good ( Figure 5G ).

To study the profile of mutations, we mapped the major mutated genes and frequencies in the different MBI groups. The predominant mutation type in both MBI groups was Missense Mutation. And TP53, CTNNB1, and MUC16 were mutated more frequently in high MBI patients. Among them, the TP53 mutation frequency was 36% in the high MBI group, while it was 16% in the low MBI group ( Figures 6A, B ).

To explore the relationship between MBI and immune cells, we investigated the expression of 16 immune cells. Multiple immune cells differed significantly between high and low MBI groups. Among them, aDCs, iDCs, Macrophages, and Treg have higher expression in high MBI patients. And Mast cells, NK cell were mainly expressed in the low MBI group ( Figure 7A ). Subsequently, we further evaluated the differences in immune-related functions in high and low MBI groups, and we found that APC_co_stimulation and HLA were predominantly expressed in the high MBI group. While Type_I_IFN_Response and Type_II_IFN_Response was expressed in low MBI patients ( Figure 7B ). Next, we found a lower Stromalscore in the high MBI group in the tumor microenvironment ( Figure 7C ). Considering the relationship between MBI and immunotherapy, we evaluated the relationship between MBI and MMR. MMR all correlated with MBI, with MSH2 having the highest correlation with MBI ( Figure 7D ). Interestingly, immune checkpoint analysis revealed a significant correlation between MBI and PDL1, PD1, and CTLA4 ( Figure 7E ). The model was found to be a better predictor of the effect of immunotherapy through the dataset, where the TIDE score was lower in high MBI patients ( Figure 7F ). Afterwards, we found that patients with better immunotherapy outcomes had higher MBI ( Figure 7G ). Subsequently, we investigated the sensitivity of different chemotherapeutic agents in different MBI groups, and we found IC50 values of Sorafenib, Cisplatin, Cytarabine, Fludarabine, Ibrutinib, Dihydrorotenone, Gemcitabine, Irinotecan, Mitoxantrone, Oxaliplatin and Ribociclib were higher in the high MBI group. While 5-Fluorouracil, Dasatinib, Osimertinib, Lapatinib, and Gefitinib had higher IC50 values in the low MBI group ( Figures 8A–S ).

We first found that both ATIC and G6PD were significantly higher expression in cancerous tissue through the GEO database ( Figures 9A, B ). The subsequent qPCR performed on 10 pairs of hepatocellular carcinoma and paraneoplastic tissues used both revealed higher expression of these two key genes in the tumor tissues ( Figures 9C, D ). Immediately after, we further performed qPCR validation using LO2 normal hepatocytes and three types of hepatocellular carcinoma cells, BEL-7402, HEPG2 and HCCLM3, and we found hub genes were higher in cancer cells than in normal hepatocytes ( Figures 9E, F ). Finally, we observed the immunohistochemistry of these two critical genes. The results were basically consistent with the above, that is, higher expression in HCC tissues ( Figure 9G, H ).