Metabolism-related signatures is correlated with poor prognosis and immune infiltration in hepatocellular carcinoma via multi-omics analysis and basic experiments

Background: Metabolism is an ordered series of biological processes that occur in an organism. Altered cellular metabolism is often closely associated with the development of cancer. The aim of this research was to construct a model by multiple metabolism-related molecules to diagnose and assess the prognosis of patients.

Method: WGCNA analysis was used to screen out differential genes. GO, KEGG are used to explore potential pathways and mechanisms. The lasso regression model was used to filter out the best indicators to construct the model. Single-sample GSEA (ssGSEA) assess immune cells abundance, immune terms in different Metabolism Index (MBI) groups. Human tissues and cells were used to verify the expression of key genes.

Result: WGCNA clustering grouped genes into 5 modules, of which 90 genes from the MEbrown module were selected for subsequent analysis. GO analysis was found that BP mainly has mitotic nuclear division, while KEGG pathway is enriched to Cell cycle, Cellular senescence. Mutation analysis revealed that the frequency of TP53 mutations was much higher in samples from the high MBI group than in the low MBI group. Immunoassay revealed that patients with higher MBI have higher macrophage and Regulatory T cells (Treg) abundance, while NK cells were lowly expressed in the high MBI group. RT-qPCR and immunohistochemistry (IHC) revealed that the hub genes expression is higher in cancer tissues. The expression in hepatocellular carcinoma cells was also much higher than that in normal hepatocytes.

Conclusion: In conclusion, a metabolism-related model was constructed that can be used to estimate the prognosis of hepatocellular carcinoma, and the clinical treatment of different hepatocellular carcinoma patients with medications was guided.

1 Introduction

Hepatocellular carcinoma is one of the most commonly diagnosed cancers and is the 2nd major cause of cancer-related deaths (1). According to statistics, Age-standardized incidence and mortality rates (ASRs) of liver cancer worldwide are 9.5 and 8.7 per 100,000 people, respectively. The number of new liver cancer patients will rise by 55% each year in the next 20 years, with 1.4 million people likely to be diagnosed in 2040 (2). Of these, 75-85% of primary liver cancers are Hepatocellular carcinoma (3). Hepatocellular carcinoma (HCC) is mainly caused by chronic HBV or HCV infection, Excessive drinking, rare genetic disorders, and various metabolic diseases (4, 5). Due to the complex etiology of hepatocellular carcinoma, the treatment methods of hepatocellular carcinoma have great differences (6). Therefore, new biomarkers and prognostic models are needed for the precise management of individuals.

Metabolism is an orderly sequence of biological changes that take place in an organism to sustain life (7, 8). Several studies have shown that alterations in cellular metabolism are often closely associated with the development of cancer (9, 10). One study found hepatocellular carcinoma was promoted with fat by inducing glucose metabolism in unconverted hepatocytes (11). And HDAC11 can regulate the glycolytic process through LKB1/AMPK signaling pathway, which can maintain hepatocyte cancer stemness (12). In addition, PRMT6-ERK-PKM2 plays an important role in tumorigenesis. It also provides an important connection between tumor and glucose metabolism (13).. However, current studies on liver cancer and metabolism are often limited to a single molecule, while studies on multiple metabolism-related genes and liver cancer are still scarce.

The aim of this study is to integrate multiple metabolism-related key genes to construct an excellent model to diagnose and assess patient prognosis. The correlation of the model with immunotherapy, drug sensitivity allows to give different treatment regimens for different populations. In addition, we validated these results by multi-omics analysis and basic experiments.

2 Materials and methods

2.1 Data source

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

2.2 Weighted gene co-expression network analysis

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

2.3 Functional enrichment analysis

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.

2.4 Model construction

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.

2.5 Immunoassay and drug sensitivity analysis

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

2.6 Immunohistochemistry

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.

2.7 Sample collection

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.

2.8 Real-time quantitative PCR

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.

2.9 Data statistics

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.

3 Result

3.1 WGCNA screening for differential genes

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

FIGURE 1

Figure 1 A flow chart of the manuscript.

FIGURE 2

Figure 2 WGCNA screening for key genes. (A) The volcano plot clearly shows the differential genes, with a total of 362 differential genes. (B) the optimal soft threshold 8 is used to construct the co-expression network. (C, D) The genes were divided into 5 modules by co-expression network, and the P-value and correlation coefficient of each module were clearly marked. (E) VENN plot showing the intersection of differential genes and the MEbrown module of WGCNA with 52 genes.

3.2 Construction of lasso regression model

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

FIGURE 3

Figure 3 Lasso regression model and KM prognostic analysis. (A) Univariate COX prognostic analysis revealed genes with prognostic value. (B) The LASSO regression model calculated the best lamda value of 2 and selected the best modeled genes. (C, D) KM analysis shows poorer prognosis when both ATIC and G6PD are genes with high expression. (E, F) In the TCGA and ICGC databases, KM analysis identified a poorer prognosis in the high MBI group.

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

FIGURE 4

Figure 4 GO and KEGG enrichment analysis. (A-B) Gene Ontology functional enrichment analyses for differentially expressed genes. (C-D) KEGG pathway enrichment analyses for differentially expressed genes.

3.3 Prognostic value of MBI

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

FIGURE 5

Figure 5 Prognostic Value and Clinical Value of MBI. (A) Univariate COX analysis identified MBI as a risk factor. (B) multivariate COX analysis identified MBI as a risk factor. (C) The nomogram constructed to predict the probability of patient mortality. (D–F) ROC curves show that both MBI and nomogram scores are good predictors. (G) The calibration curve results found that the model works well. ***P < 0.001.

3.4 Mutation landscape of different MBI groups

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

FIGURE 6

Figure 6 Mutation landscape of different MBI groups. (A) Mutation types, mutated genes and mutation frequencies in patients in the high MBI group. (B) Mutation types, mutated genes and mutation frequencies in patients in the low MBI group.

3.5 Immunoassay and drug sensitivity analysis

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

FIGURE 7

Figure 7 Analysis of immune cells, immune terminology and immunotherapy. (A) Immunoassay revealed the expression of aDCs, iDCs, Macrophages, Treg in the high MBI group. Mast_cells and NK_cell are mainly expressed in the low MBI group. (B) APC_co_stimulation and HLA are mainly expressed in the high MBI group. While Type_Ⅰ_IFN_Response and Type_Ⅱ_IFN_Response are mainly expressed in the low MBI group. (C) Lower Stromalscore in the high MBI group in the tumor microenvironment. (D) MMR is correlated with MBI, with MSH2 having the highest correlation with MBI. (E) Immune checkpoint analysis revealed a significant correlation between MBI and PDL1, PD1, and CTLA4 etc. (F) The TIDE score is lower in high MBI patients. (G) Patients with better immunotherapy outcomes had higher MBI in IMvigor210 dataset. ***P < 0.001, **P < 0.01, *P < 0.05.

FIGURE 8

Figure 8 Sensitivity analysis of multiple chemotherapeutic agents. (A-K) Sorafenib, Cisplatin, Cytarabine and Fludarabine, Ibrutinib, Dihydrorotenone, Gemcitabine, Irinotecan, Mitoxantrone, Oxaliplatin and Ribociclib had higher IC50 values in the high MBI group. (O-S) 5-Fluorouracil, Dasatinib, Osimertinib, Lapatinib, and Gefitinib had higher IC50 values in the low MBI group. ***P < 0.001, **P < 0.01.

3.6 External and experimental validation of key gene expression

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

FIGURE 9

Figure 9 Validation of ATIC and G6PD expression in Tissues and Cell Lines. (A, B) GEO database (GSE45267) demonstrates that ATIC and G6PD are significantly more expressed in tumor tissues than in normal tissues. (C, D) qPCR demonstrates that ATIC and G6PD are significantly more expressed in tumor tissues than in normal tissues. (E, F) qPCR demonstrated that the expression of ATIC and G6PD was significantly higher in hepatocellular carcinoma cells than in normal hepatocytes. (G, H) IHC demonstrated that the expression of ATIC and G6PD was significantly higher in tumor tissues than in normal tissues.

4 Discussion

The treatment of hepatocellular carcinoma is still mainly surgical, but the recurrence rate of tumors in patients with intermediate and advanced stages is extremely high (19). With the advent of the era of tumor immunotherapy, immune checkpoint inhibitors have been used as one of the new effective approaches for tumor treatment (20). Therefore, it is difficult to know whether immunotherapy is effective for a particular patient, and how to screen this group of responding patients is a clinical problem that needs to be solved. We used metabolism-related genes to construct the model and also studied the correlation between the model and immunity to screen for this population of patients.

In this study, AITC and G6PD were screened as key genes for the model by WGCNA as well as lasso-cox based on TCGA transcriptome data and clinical information. It was shown that G6PD is activated by TSP50 along with the development of hepatocellular carcinoma (21). In addition, Wang and Zhao et al. found that G6PD could be used for building models of hepatocellular carcinoma (22, 23). All of the above studies indicate that G6PD and the development of hepatocellular carcinoma are closely related. The bifunctional enzyme ATIC can promote the proliferation of liver cancer (24). Moreover, ATIC can also be used as a prognostic biomarker for hepatocellular carcinoma (25, 26). These studies further confirmed that these two genes are good indicators for the construction of prognostic models for liver cancer.

We further performed functional enrichment analysis. Where BP is enriched to mitotic nuclear division, it has been found that mitotic transmission disruption usually leads to aneuploid offspring production. Aneuploidy is a general characteristic of tumor cells (27, 28). Interestingly, KEGG is enriched to Cell cycle, Chemical carcinogenesis, Cellular senescence and other pathways. Matthews’ study found that impairing cell exit from the cell cycle to allow continuous cell division can lead to cancer development (29–31). Cellular senescence is a condition in which the cell cycle is stalled. The senescent state leads to the development of cancer by keeping the cells alive (32, 33). These studies further confirm our results and show that MBI is a good predictor. Further exploration of the correlation of MBI and prognosis, we performed univariate and multivariate cox analyses and found MBI to be a risk factor. In different years, the ROC curve showed that the prediction efficiency of the model was better. Further nomograms can help to score clinical patients so that different treatment plans can be given to patients with different scores. Subsequently, to explore the differential profile in different MBI groups, we analyzed the mutation types and mutation frequencies. The TP53 mutation frequency was much higher in samples of the high MBI patients. It was shown that Survival of liver cancer cells increases with increasing mitochondrial fission wtih coordinated regulation of ROS-regulated NF-κB and TP53 pathways (34). And TP53 are significantly mutated in liver cancer (35, 36). Several studies above have demonstrated that mutations in are strongly correlated with hepatocellular carcinoma from multiple clinical and basic perspectives. So, the high MBI group has a poorer prognosis for this reason.

To explore the relevance of MBI and immunity. We studied 16 immune cells and 13 immune-related terms. Previous studies have demonstrated that Macrophages polarize under the influence of tumor microenvironment and form tumor-associated macrophages (TAM). And the abundance of TAM in tumors is closely associated with poor prognosis (37). Our study revealed that macrophage abundance was higher in the high MBI patients. This also explains that high MBI patients had a worse prognosis. Togashi’s study showed that Treg cells expressing FOXP3 suppress aberrant immune responses against autoantigens and also suppress anti-tumor immune responses. Large numbers of Treg cells infiltrating into the tumor tissue usually have a poor prognosis (38). This is broadly in line with our results that the higher MBI patients had higher Treg. In addition, we found that NK cells were lowly expressed in patients in the high MBI group. According to others, if cell surface markers are correlated with carcinogenic transformation, NK cells can quickly eliminate them. This property is unique among immune cells, and their ability to enhance antibody and T cell responses supports the role of NK cells as anti-cancer agents (39). This further confirms our results. After that, we observed that the interferon response was mainly expressed in the group with low MBI. Interestingly, Boukhaled et al. found that IFN contributes to antitumor immune quality and immunotherapeutic response (40, 41). Subsequently, we predicted the effect of immunotherapy by analyzing the TIDE scores of the high and low MBI groups. In addition, we found that MBI was closely associated with immunotherapy by MMR and immune checkpoint analysis. In particular, the correlation between CTLA4 and MBI was extremely high, which indicates that MBI is of great value in guiding clinical immunotherapy. Our validation using the dataset revealed that patients with high MBI had better outcomes after immunotherapy. Our study also found significant differences between the high and low MBI groups across multiple chemotherapeutic agents. These results can further guide our clinical use of medications. Finally, we verified the expression of key genes by qPCR in tumor and normal tissues and in hepatocytes and hepatoma cells. This was confirmed by further IHC results.

This study has important clinical applications. Our metabolism-related genes screened by WGCNA combined with differential analysis and prognostic analysis were more credible. Exploration of GO and KEGG helps us to understand the reasons for the difference between high and low risk groups. It helps to point the way for our next study. Moreover, MBI calculated by the model constructed by lasso-cox is a credible and independent biological marker to predict the prognosis of patients of hepatocellular carcinoma. Studying the correlation between MBI and immune and drug name susceptibility can be a useful indicator to assess the efficacy of immunotherapy and chemotherapy in patients.

There are several limitations to our study that need to be acknowledged. First, this study is an analysis using a public database and lacks validation of our own cohort. We will further investigate these hub genes in our own cohort of hepatocellular carcinoma data. Secondly, the downstream target genes of these two genes were not explored further, which also needs our further study afterwards.

5 Conclusion

In brief, we constructed a metabolism-related model. We hope that this model can be used as a reference for predicting patient survival and guiding related treatments for patients with liver cancer.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics statement

The studies involving human participants were reviewed and approved by The ethics committee of the First Affiliated Hospital of Anhui Medical University. The patients/participants provided their written informed consent to participate in this study.

Author contributions

YG have constructed and devised the research. JS and WS performed data analysis and wrote the manuscript. JYL and YL analyzed the data. JLL acquired samples and performed the experiments. 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.

Supplementary material

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

Supplementary Figure 1 | Risk curves for different data sets. (A, B) TCGA database analysis finds higher mortality in patients in high MBI group. (C, D) ICGC database analysis finds higher mortality in patients in high MBI group.

Abbreviations

WGCNA, Weighted Gene Co-expression Network Analysis; ssGSEA, Single-sample GSEA; MBI, Metabolism Index; Treg, Regulatory T cells; IHC, Immunohistochemistry; ASRs, Age-standardised Incidence and Mortality Rates; HCC, Hepatocellular Carcinoma; TCGA, The Cancer Genome Atlas; NCI, National Cancer Institute; NHGRI, National Human Genome Research Institute; GO, Gene Ontology; CC, Cellular Component; MF, Molecular Function; BP, Biological Process; HPA, Human Protein Atlas; RT-qPCR, Real-Time quantitative PCR; HR, Hazard Ratio; CI, Confidence Interval; TAM, Tumor-associated Macrophages.

References

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

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Rumgay H, Arnold M, Ferlay J, Lesi O, Cabasag CJ, Vignat J, et al. Global burden of primary liver cancer in 2020 and predictions to 2040. J Hepatol (2022) 77(6):1598–606. doi: 10.1016/j.jhep.2022.08.021

PubMed Abstract | CrossRef Full Text | Google Scholar

3. 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(3):209–49. doi: 10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Welzel TM, Graubard BI, Quraishi S, Zeuzem S, Davila JA, El-Serag HB, et al. Population-attributable fractions of risk factors for hepatocellular carcinoma in the united states. Am J Gastroenterol (2013) 108(8):1314–21. doi: 10.1038/ajg.2013.160

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Kulik L, El-Serag HB. Epidemiology and management of hepatocellular carcinoma. Gastroenterology (2019) 156(2):477–91 e1. doi: 10.1053/j.gastro.2018.08.065

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Llovet JM, De Baere T, Kulik L, Haber PK, Greten TF, Meyer T, et al. Locoregional therapies in the era of molecular and immune treatments for hepatocellular carcinoma. Nat Rev Gastroenterol Hepatol (2021) 18(5):293–313. doi: 10.1038/s41575-020-00395-0

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Judge A, Dodd MS. Metabolism. Essays Biochem (2020) 64(4):607–47. doi: 10.1042/EBC20190041

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Andrejeva G, Rathmell JC. Similarities and distinctions of cancer and immune metabolism in inflammation and tumors. Cell Metab (2017) 26(1):49–70. doi: 10.1016/j.cmet.2017.06.004

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Kroemer G, Pouyssegur J. Tumor cell metabolism: Cancer's achilles' heel. Cancer Cell (2008) 13(6):472–82. doi: 10.1016/j.ccr.2008.05.005

PubMed Abstract | CrossRef Full Text | Google Scholar

10. DeBerardinis RJ, Chandel NS. Fundamentals of cancer metabolism. Sci Adv (2016) 2(5):e1600200. doi: 10.1126/sciadv.1600200

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Broadfield LA, Duarte JAG, Schmieder R, Broekaert D, Veys K, Planque M, et al. Fat induces glucose metabolism in nontransformed liver cells and promotes liver tumorigenesis. Cancer Res (2021) 81(8):1988–2001. doi: 10.1158/0008-5472.CAN-20-1954

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Bi L, Ren Y, Feng M, Meng P, Wang Q, Chen W, et al. HDAC11 regulates glycolysis through the LKB1/AMPK signaling pathway to maintain hepatocellular carcinoma stemness. Cancer Res (2021) 81(8):2015–28. doi: 10.1158/0008-5472.CAN-20-3044

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Gu L, Zhu Y, Lin X, Tan X, Lu B, Li Y. Stabilization of FASN by ACAT1-mediated GNPAT acetylation promotes lipid metabolism and hepatocarcinogenesis. Oncogene (2020) 39(11):2437–49. doi: 10.1038/s41388-020-1156-0

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Hutter C, Zenklusen JC. The cancer genome atlas: Creating lasting value beyond its data. Cell (2018) 173(2):283–5. doi: 10.1016/j.cell.2018.03.042

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Langfelder P, Horvath S. WGCNA: An r package for weighted correlation network analysis. BMC Bioinf (2008) 9:559. doi: 10.1186/1471-2105-9-559

CrossRef Full Text | Google Scholar

16. Maeser D, Gruener RF, Huang RS. oncoPredict: An r package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform (2021) 22(6):bbac468. doi: 10.1093/bib/bbab260

CrossRef Full Text | Google Scholar

17. Uhlen M, Zhang C, Lee S, Sjostedt E, Fagerberg L, Bidkhori G, et al. A pathology atlas of the human cancer transcriptome. Science (2017) 357(6352):eaan2507. doi: 10.1126/science.aan2507

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Sun W, Shen J, Liu J, Han K, Liang L, Gao Y. Gene signature and prognostic value of ubiquitin-specific proteases members in hepatocellular carcinoma and explored the immunological role of USP36. Front Biosci (Landmark Ed) (2022) 27(6):190. doi: 10.31083/j.fbl2706190

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Couri T, Pillai A. Goals and targets for personalized therapy for HCC. Hepatol Int (2019) 13(2):125–37. doi: 10.1007/s12072-018-9919-1

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Llovet JM, Montal R, Sia D, Finn RS. Molecular therapies and precision medicine for hepatocellular carcinoma. Nat Rev Clin Oncol (2018) 15(10):599–616. doi: 10.1038/s41571-018-0073-4

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Zhang X, Gao F, Ai H, Wang S, Song Z, Zheng L, et al. TSP50 promotes hepatocyte proliferation and tumour formation by activating glucose-6-phosphate dehydrogenase (G6PD). Cell Prolif (2021) 54(4):e13015. doi: 10.1111/cpr.13015

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Wang Z, Fu Y, Xia A, Chen C, Qu J, Xu G, et al. Prognostic and predictive role of a metabolic rate-limiting enzyme signature in hepatocellular carcinoma. Cell Prolif (2021) 54(10):e13117. doi: 10.1111/cpr.13117

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Zhao C, Zhang Z, Jing T. A novel signature of combing cuproptosis- with ferroptosis-related genes for prediction of prognosis, immunologic therapy responses and drug sensitivity in hepatocellular carcinoma. Front Oncol (2022) 12:1000993. doi: 10.3389/fonc.2022.1000993

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Li M, Jin C, Xu M, Zhou L, Li D, Yin Y. Bifunctional enzyme ATIC promotes propagation of hepatocellular carcinoma by regulating AMPK-mTOR-S6 K1 signaling. Cell Commun Signal (2017) 15(1):52. doi: 10.1186/s12964-017-0208-8

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Su WJ, Lu PZ, Wu Y, Kalpana K, Yang CK, Lu GD. Identification of key genes in purine metabolism as prognostic biomarker for hepatocellular carcinoma. Front Oncol (2020) 10:583053. doi: 10.3389/fonc.2020.583053

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Ye W, Shi Z, Zhou Y, Zhang Z, Zhou Y, Chen B, et al. Autophagy-related signatures as prognostic indicators for hepatocellular carcinoma. Front Oncol (2022) 12:654449. doi: 10.3389/fonc.2022.654449

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Weaver BA, Cleveland DW. Decoding the links between mitosis, cancer, and chemotherapy: The mitotic checkpoint, adaptation, and cell death. Cancer Cell (2005) 8(1):7–12. doi: 10.1016/j.ccr.2005.06.011

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Sarkar S, Sahoo PK, Mahata S, Pal R, Ghosh D, Mistry T, et al. Mitotic checkpoint defects: En route to cancer and drug resistance. Chromosome Res (2021) 29(2):131–44. doi: 10.1007/s10577-020-09646-x

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Matthews HK, Bertoli C, de Bruin RAM. Cell cycle control in cancer. Nat Rev Mol Cell Biol (2022) 23(1):74–88. doi: 10.1038/s41580-021-00404-3

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Sherr CJ. Cancer cell cycles. Science (1996) 274(5293):1672–7. doi: 10.1126/science.274.5293.1672

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Suski JM, Braun M, Strmiska V, Sicinski P. Targeting cell-cycle machinery in cancer. Cancer Cell (2021) 39(6):759–78. doi: 10.1016/j.ccell.2021.03.010

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Calcinotto A, Kohli J, Zagato E, Pellegrini L, Demaria M, Alimonti A. Cellular senescence: Aging, cancer, and injury. Physiol Rev (2019) 99(2):1047–78. doi: 10.1152/physrev.00020.2018

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Campisi J. Aging, cellular senescence, and cancer. Annu Rev Physiol (2013) 75:685–705. doi: 10.1146/annurev-physiol-030212-183653

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Huang Q, Zhan L, Cao H, Li J, Lyu Y, Guo X, et al. Increased mitochondrial fission promotes autophagy and hepatocellular carcinoma cell survival through the ROS-modulated coordinated regulation of the NFKB and TP53 pathways. Autophagy (2016) 12(6):999–1014. doi: 10.1080/15548627.2016.1166318

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Tornesello ML, Buonaguro L, Tatangelo F, Botti G, Izzo F, Buonaguro FM. Mutations in TP53, CTNNB1 and PIK3CA genes in hepatocellular carcinoma associated with hepatitis b and hepatitis c virus infections. Genomics (2013) 102(2):74–83. doi: 10.1016/j.ygeno.2013.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Huo J, Wu L, Zang Y. Construction and validation of a reliable six-gene prognostic signature based on the TP53 alteration for hepatocellular carcinoma. Front Oncol (2021) 11:618976. doi: 10.3389/fonc.2021.618976

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Li C, Xu X, Wei S, Jiang P, Xue L, Wang J, et al. Tumor-associated macrophages: Potential therapeutic strategies and future prospects in cancer. J Immunother Cancer (2021) 9(1):e001341. doi: 10.1136/jitc-2020-001341

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Togashi Y, Shitara K, Nishikawa H. Regulatory T cells in cancer immunosuppression - implications for anticancer therapy. Nat Rev Clin Oncol (2019) 16(6):356–71. doi: 10.1038/s41571-019-0175-7

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Shimasaki N, Jain A, Campana D. NK cells for cancer immunotherapy. Nat Rev Drug Discovery (2020) 19(3):200–18. doi: 10.1038/s41573-019-0052-1

CrossRef Full Text | Google Scholar

40. Boukhaled GM, Harding S, Brooks DG. Opposing roles of type I interferons in cancer immunity. Annu Rev Pathol (2021) 16:167–98. doi: 10.1146/annurev-pathol-031920-093932

PubMed Abstract | CrossRef Full Text | Google Scholar

41. von Locquenghien M, Rozalen C, Celia-Terrassa T. Interferons in cancer immunoediting: Sculpting metastasis and immunotherapy response. J Clin Invest (2021) 131(1):e143296. doi: 10.1172/JCI143296

PubMed Abstract | CrossRef Full Text | Google Scholar