Abstract
Background. Metabolic reprogramming is associated with the carcinogenesis of hepatocellular carcinoma (HCC). The effects of metabolism-related genes on predicting survival and immune status in HCC remain unclear.
Objectives. To develop and validate metabolic models for predicting the survival and immune status of HCC patients.
Materials and methods. The metabolic core genes for overall survival (OS) and disease-free survival (DFS) were retrieved. Then, glycolysis and fatty acid metabolism prognostic models were constructed and validated using The Cancer Genome Atlas (TCGA) and the International Cancer Genome Consortium (ICGC) data. Decision trees based on machine learning were developed for classifying the prognostic risks of HCC patients. The associations between the metabolic signatures, immunotherapy and immune cell infiltration were investigated. Experimental validations were performed using reverse transcription-quantitative polymerase chain reaction (RT-qPCR) and immunohistochemistry (IHC).
Results. We identified 30 prognostic core genes for glycolysis metabolism and 12 prognostic core genes for fatty acid metabolism. Subsequently, 2 glycolysis models and 2 fatty acid metabolism models were developed to predict the OS and DFS of HCC patients, respectively. Two decision trees were constructed to classify the low-, intermediate- and high-risk groups of HCC patients for OS and DFS. Moreover, the patients in the high-risk groups of glycolysis and fatty acid metabolic models tended to have higher expression of programmed cell death ligand-1 (PD-L1 or CD274), programmed cell death 1 (PDCD1), cytotoxic T-lymphocyte-associated protein-4 (CTLA-4), and lymphocyte activating 3 (LAG3). Most of the metabolic core genes were significantly associated with immune cell infiltration. In addition, ATP-binding cassette subfamily B member 6 (ABCB6), peptidylprolyl isomerase A (PPIA), uroporphyrinogen decarboxylase (UROD), and non-SMC condensin II complex subunit H2 (NCAPH2) were positively correlated with both tumor mutational burden (TMB) and microsatellite instability (MSI) scores. The expression of ABCB6, PPIA, UROD, and NCAPH2 was validated using RT-qPCR and IHC.
Conclusions. We established novel prognostic models based on metabolism-related genes to better predict the outcome and immune status of HCC patients.
Key words: hepatocellular carcinoma, glycolysis, fatty acid metabolism, prognostic models, immune status
Background
Hepatocellular carcinoma (HCC) is characterized by a high malignancy and poor prognosis.1, 2 In China, the number of patients with hepatitis B virus-related HCC has been decreasing, while the incidence of fatty liver-related HCC is increasing due to accelerated aging and changes in diet composition.3, 4 Hepatocellular carcinoma remains a serious public health threat and increases the social and economic burdens on the government.5
The appropriate treatment depends on the physician’s assessment of the patient’s prognosis. For patients with a poor prognosis, doctors tend to use additional treatment methods, such as transcatheter arterial chemoembolization (TACE), targeted therapy and immunotherapy.6 In the past, the prognostic evaluation of HCC mainly depended on the pathological stage. Several studies have shown that in addition to the American Joint Committee on Cancer (AJCC) tumor-node-metastasis (TNM) staging system, there are other clinicopathological features and specific gene expressions that can be used to predict the prognosis of HCC patients.7, 8
Cancer development depends on the metabolic reprogramming of tumor cells, which is a hallmark of cancer and a direct or indirect result of oncogene mutations.9 The abnormal expression of metabolic genes is closely related to the prognosis of cancer patients.10 A prognostic model based on metabolic genes can be used to predict the prognosis of cancer patient.11 The liver is the primary organ involved in metabolizing glucose and lipids. Abnormal gene regulation in the liver leads to abnormal glucose and lipid metabolism, thus inducing tumor progression.12 In addition, abnormal metabolic regulation is closely associated with changes in the tumor microenvironment and immune cell infiltration.13, 14
In recent years, there has been a breakthrough in the treatment of advanced HCC patients with immune checkpoint inhibitors (ICIs).15 However, many patients do not respond to immunotherapy. Previous studies suggest that immune markers, such as programmed cell death ligand-1 (PD-L1 or CD274), cytotoxic T-lymphocyte-associated protein-4 (CTLA-4), tumor mutational burden (TMB), and microsatellite instability (MSI), can be used to predict the response to immunotherapy.16, 17 The expression of core metabolic genes is also considered a potential marker for predicting the effectiveness of immunotherapy.18, 19
Objectives
This study aims to screen for prognostic core genes of glycolysis and fatty acid metabolism and build a prediction model combining TNM features with gene risk scores to better predict overall survival (OS), disease-free survival (DFS) and the immune status of HCC patients.
Materials and methods
Data collection
RNA-sequencing expression profiles and their corresponding clinical data were obtained from The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov) and the International Cancer Genome Consortium (ICGC; https://icgc.org/). Using the bioinformatic analysis website Assistant for Clinical Bioinformatics (ACBI; https://www. aclbi.com/), RNA sequencing and survival data were obtained for 371 HCC patients from TCGA and 240 patients from the ICGC. Gene expression patterns varied significantly between histological subtypes; thus, intrahepatic cholangiocarcinoma and mixed HCC were excluded and only patients with HCC were selected. Clinical parameters, including age, gender, pathological stage (T, N and M stages), tumor grade, vital status, recurrence status, and follow-up time [years] were collected. The RNA count data were converted to transcripts per million (TPM). Hepatocellular carcinoma patients were divided into high- and low-expression groups using the median RNA expression as the cutoff value.
Acquisition of metabolic hallmark gene sets
The core metabolic genes were retrieved from the Molecular Signatures Database v. 7.5.1 (MSigDB, http://www.gsea-msigdb.org/gsea/msigdb). For further analysis, 200 glycolysis-related genes were collected from “HALLMARK_GLYCOLYSIS” and 158 fatty acid metabolism-related genes were collected from “HALLMARK_FATTY_ACID_METABOLISM” (Supplementary Table 1).
Identification of prognostic genes
for OS and DFS
The univariate Cox regression analysis was used to assess the prognostic value of the glycolysis-related and fatty acid metabolism-related genes. Assistant for Clinical Bioinformatics was used for survival analyses of batches (https://www.aclbi.com/static/index.html#/). Hazard ratios (HRs), OS and DFS were calculated for each selected gene individually with 95% confidence intervals (95% CIs). Genes significantly associated with both OS and DFS were considered core metabolic genes and entered the next stage of model building. Differences in the expression levels of core metabolic genes in tumors and normal tissues were compared using the Wilcoxon test.
Construction of metabolic risk score prediction models
With 10-fold cross-validation, the Least Absolute Shrinkage and Selection Operator (LASSO) regression algorithm was used to further filter core metabolic genes to develop the prognostic glycolysis and fatty acid metabolism signatures using the “glmnet” R package (R Foundation for Statistical Computing, Vienna, Austria).20 Risk scores for HCC patients were calculated based on the normalized gene expression levels (log2(TPM+1)) and corresponding regression coefficients in our models. The LASSO regression analysis was performed using the ACBI website. Hepatocellular carcinoma patients were then divided into high-risk and low-risk groups based on the median value of the risk scores. The Kaplan–Meier survival curves and log-rank tests were used to compare the OS and DFS between the 2 groups. To evaluate the efficacy of the glycolysis and fatty acid metabolic signatures, receiver operating characteristic (ROC) curves for 1, 3 and 5 years with the area under the curve (AUC) values were plotted using the “survivalROC” R package.
Integrated analysis of clinical parameters and metabolic risk score prediction models
Clinical information (age, gender, pathological stage: T, N and M, tumor grade) and glycolysis and fatty acid metabolism risk scores of HCC patients were integrated. To further build nomograms based on the pathological TNM data, patients with missing data were excluded. The univariate and multivariate Cox regression analyses were then performed to identify independent prognostic factors.
Development and validation of the nomograms for patients with HCC
To predict OS and DFS at 3 and 5 years in individual HCC patients, nomograms were developed by combining metabolic risks and clinical information based on the training groups, according to the results of the multivariate Cox regression (metabolic models). Corresponding models based on the pathologic TNM stage (TNM models) were constructed for head-to-head comparison using the metabolic models.
The discrimination of a predictive model was validated using concordance statistics. The concordance index (C-index) was analogous to the AUC value. In addition, the time-dependent ROC curve was used to summarize the predictive capacity of the models according to the continuous change over time.21, 22 Model calibration was validated using the calibration plot. We also developed a decision curve analysis (DCA) graph to test the clinical value of the metabolic models through visualizing the potential net benefits at each threshold probability between the metabolic and TNM models.23
Establishment and optimization of decision tree models
The Classification and Regression Trees (CART) algorithm is a non-parametric machine learning method,24, 25 and is not dependent on any type of distribution. The CART algorithm generates a top-down decision tree until the stopping criteria are met. The decision tree model has the advantage of being flexible and powerful, and can be presented in graphical form. We used the CART algorithm to establish decision tree models to identify low-, intermediate- and high-risk patient groups for OS and DFS, respectively. The glycolysis model (low risk, high risk), fatty acid metabolism (low risk, high risk) and pathological T (T1, T2, T3, and T4) were used as the covariates. The complexity parameter was used for tree pruning to optimize the decision tree. We used cross-validation to evaluate the relative error at different complexity parameter values. The “rpart” R package was used to generate the decision tree,25 and the R codes are presented in the Supplementary material.
Conjoint analysis of TMB, MSI, immune cell infiltration, and metabolic core genes
The ACBI website was used to calculate TMB and MSI scores. Spearman’s correlation analysis was then used to determine the relationship between TMB, MSI and the metabolic core genes. In addition, the relationship between immune cell infiltration and metabolic core gene expression was evaluated using the Tumor Immune Estimation Resource (TIMER) algorithm and Spearman’s correlation analysis.
Experimental validation
Ten pairs of specimens were collected from postoperative patients with HCC from the Ningbo Medical Center Lihuili Hospital, China, between in 2021 and 2022. Total RNA was isolated using FastPure Cell/Tissue Total RNA Isolation Kit V2 (Vazyme Biotech, Nanjing, China) and cDNA was synthesized with HiScript II Q RT SuperMix for qPCR (Vazyme Biotech). The method of reverse transcription-quantitative polymerase chain reaction (RT-qPCR) was described in our previous study.2 All primers were purchased from Tsingke Biotechnology (Beijing, China) and are listed in Supplementary Table 2. The RNA expression of the identified genes in different HCC tumor cells was extracted from the Cancer Cell Line Encyclopedia (CCLE, https://sites.broadinstitute.org/ccle),26 and the results were visualized using Expression Atlas (https://www.ebi.ac.uk/gxa/experiments). Immunohistochemistry (IHC) results of prognostic genes were identified from The Human Protein Atlas (https://www.proteinatlas.org/). The Institutional Ethics Committee of Ningbo Medical Center Lihuili Hospital approved the study (approval No. KY2020PJ186), and written informed consent was acquired from the patients in accordance with the Declaration of Helsinki.
Statistical analyses
The differences in gene expression between different groups were compared using the Wilcoxon test. Spearman’s correlation analysis was performed to analyze the correlations. The Kaplan–Meier survival analysis was used to compare the survival differences with the log-rank test. The univariate and multivariate Cox regression analyses were used to estimate the HR for each potential risk factor. The two-tailed Student’s t-test was used for quantitative data.
The Cox regression model met 2 assumptions: 1) the proportional hazards assumption and 2) that the relationship between each log’s hazard and continuous covariate is linear. The proportional hazards assumption was tested graphically for categorical variables and with the scaled Schoenfeld residual test for continuous covariates. The Martingale residual plot was used to test the linear relationship between continuous covariates. If the data followed a normal distribution, the parametric t-test was used, and Welch’s correction for variance heterogeneity was applied. Otherwise, the Mann–Whitney method for nonparametric t-tests was used (the results of assumption testing are shown in the Supplementary material).
A two-tailed p-value ≤0.05 was considered statistically significant. Statistical analyses were performed using R software v. 4.0.3 (R Foundation for Statistical Computing) and STATA v. 16.0 (StataCorp LLC, College Station, USA).
Results
Screening prognostic genes
for OS and DFS
The study procedure is presented in the flowchart (Figure 1). RNA sequencing and survival data for the 371 HCC patients from TCGA and 240 patients from the ICGC were acquired.
For the glycolic metabolism hallmark gene set, 85 genes were significantly associated with OS (Supplementary Fig. 1) and 35 genes were significantly associated with DFS (Supplementary Fig. 2). A Venn diagram was generated, and 30 genes were found to be statistically significant for both OS and DFS (Figure 2A). The prognostic values of these 30 genes are shown using forest plots in Figure 2B,C. Twenty-nine prognostic genes were significantly differentially expressed between tumors and normal tissues in HCC patients, except for the lactase (LCT) gene (Supplementary Fig. 3).
For the fatty acid metabolism hallmark gene set, 44 genes were significantly associated with OS (Supplementary Fig. 4), and 19 genes were significantly associated with DFS (Supplementary Fig. 5). The Venn diagram showed that there were 12 genes exhibiting statistical significance in both OS and DFS (Figure 2D). The forest plots in Figure 2E,F show the prognostic values of these 12 genes. Furthermore, these 12 prognostic genes were significantly differentially expressed between tumors and normal tissues in HCC patients (Supplementary Fig. 6).
Construction of prognostic models based on selected glycolysis-related genes
The 30 filtered genes were included in the LASSO regression algorithm to screen for glycolic prognostic signatures for OS and DFS (Figure 3A−L), and cross-validation was performed.
We selected 8 prognostic signatures for predicting an individual’s prognostic risk for OS (Figure 3B). Figure 3C shows the heatmap of the 8 genes. The risk score for OS was calculated using the following formula (the gene names in the formula indicate the normalized gene expression levels (log2(TPM+1))): risk score = (0.054 × CENPA) + (0.2044 × ABCB6) + (0.0038 × SAP30) + (0.1457 × B3GAT3) + (0.0383 × KIF20A) + (0.1171 × G6PD) + (0.0465 × HMMR) + (0.0972 × SRD5A3).
The HCC patients were divided into high-risk and low-risk groups according to the median risk score for OS (Figure 3D). The OS of the high-risk group was significantly worse than that of the low-risk group (p < 0.001, Figure 3E). The ROC curves were constructed, and the AUCs at 1, 3 and 5 years were 0.803, 0.724 and 0.654, respectively (Figure 3F). This model was validated using the independent data from the ICGC, where the AUCs at 1, 3 and 5 years were 0.706, 0.702 and 0.615, respectively.
We selected 11 prognostic signatures to predict an individual’s prognostic risk for DFS (Figure 3H). Figure 3I shows the heatmap of the 11 genes. The risk score for DFS was calculated using the following formula: risk score = (0.2161 × TXN) + (0.1311 × UGP2) + (0.1189 × ABCB6) + (0.1781 × SAP30) + (0.0303 × PPIA) + (0.0784 × PHKA2) + (0.0214 × KIF20A) + (0.022 × VEGFA) + (−0.6217 × NDST3) + (0.1059 × GLRX) + (0.157 × SRD5A3).
The HCC patients were again divided into high-risk and low-risk groups according to median risk scores for DFS (Figure 3J). The high-risk group had a significantly worse DFS (p < 0.001, Figure 3K). The ROC curves were constructed, and the AUCs at 1, 3 and 5 years were 0.720, 0.689 and 0.728, respectively (Figure 3L). As there were no follow-up data for DFS in the ICGC, we did not perform the external validation for the DFS model.
Construction of prognostic models based on selected fatty acid metabolism-related genes
The 12 filtered genes were included in the LASSO regression algorithm to screen for prognostic signatures of fatty acid metabolism for OS and DFS (Figure 4A–L), and cross-validation was performed.
We selected 7 genes to predict an individual’s prognostic risk for OS (Figure 4B). Figure 4C shows the heatmap of the 7 genes. The risk score for OS was calculated using the following formula: risk score = (−0.0793 × ACADL) + (0.1202 × APEX1) + (0.0852 × UGDH) + (0.342 × CCDC58) + (0.4259 × UROD) + (0.0578 × S100A10) + (0.1847 × SLC22A5).
The HCC patients were separated into high-risk and low-risk groups according to the median risk scores for OS (Figure 4D). The Kaplan–Meier survival analysis showed that the OS of the high-risk group was significantly worse than that of the low-risk group (p < 0.001, Figure 4E). The ROC curves were constructed, and the AUCs at 1, 3 and 5 years were 0.790, 0.684 and 0.679, respectively (Figure 4F). This model was validated using independent data from the ICGC, and the AUCs at 1, 3 and 5 years were 0.607, 0.632 and 0.671, respectively.
We selected 8 genes to predict an individual’s prognostic risk for DFS (Figure 4H). Figure 4I shows the heatmap of the 8 genes. The risk score for DFS was calculated using the following formula: risk score = (0.0554 × APEX1) + (0.0481 × UGDH) + (0.1326 × CCDC58) + (0.1365 × UROD) + (0.0106 × S100A10) + (0.0422 × NCAPH2) + (0.1631 × SLC22A5) + (−0.0451 × HAO2).
The HCC patients were again divided into high-risk and low-risk groups according to the median risk scores for DFS (Figure 4J). The high-risk group had a significantly worse DFS (p < 0.001, Figure 4K). The ROC curves were constructed, and the AUCs at 1, 3 and 5 years were 0.704, 0.616 and 0.671, respectively (Figure 4L). As there was no follow-up data for DFS in the ICGC, we did not perform the external validation for the DFS model.
Integrated analysis of clinical parameters and metabolic risk scores
After excluding patients with missing data for pathological T, N and M stages and follow-up time, 346 HCC patients with clinical data and risk scores were selected for the integrated analysis of OS, and 298 HCC patients were identified for the integrated analysis of DFS.
First, the univariate Cox regression analysis showed that T-stage, M-stage, glycolysis risk score (HR = 2.46; 95% CI: 1.68–3.59), and fatty acid metabolism risk score (HR = 2.40; 95% CI: 1.64–3.51) were significantly associated with the OS of HCC patients (Table 1). In addition, T-stage glycolysis risk score (HR = 2.50; 95% CI: 1.76–3.57) and fatty acid metabolism risk score (HR = 1.90; 95% CI: 1.34–2.69) were significantly correlated with the DFS of HCC patients (Table 1).
The multivariate Cox regression analysis was then performed, which showed that T stage, glycolysis risk score (HR = 1.70; 95% CI: 1.05–2.75) and fatty acid metabolism risk score (HR = 1.75; 95% CI: 1.09–2.78) were independent prognostic factors for OS. The T stage and glycolysis risk score (HR = 2.01; 95% CI: 1.35–2.99) were independent prognostic factors for DFS (Table 1).
Development and validation of the nomograms
The HCC patients were randomly assigned in a 2:1 ratio, considering 2/3 of the patients as the training group and the remaining patients as the validation group. Based on these results, 2 prognostic (metabolic) models were developed through combining clinical parameters and metabolic risk scores. In addition, graphical nomograms were generated to predict the probability of 3-year and 5-year survival for OS and DFS (Figure 5A,B). The calibration plots showed good agreement between the observed outcome and the predicted probability (Supplementary Fig. 7A–F). In addition, 2 prediction models (TNM models) were constructed based on the TNM status of HCC patients for head-to-head comparison (Supplementary Fig. 8).
The predictive accuracy of the metabolic nomogram calculated using AUC was 0.773 for 3-year OS and 0.701 for 5-year OS. Furthermore, the AUC of the TNM model was 0.686 for 3-year OS and 0.663 for 5-year OS (Supplementary Fig. 7C,D). The time-dependent ROC curves revealed that the predictive accuracy of the metabolic model was significantly higher than that of the TNM model in both the training and validation groups (Figure 5C,D).
The AUC of the metabolic nomogram was 0.680 for 3-year DFS and 0.676 for 5-year DFS. In addition, the AUC of the TNM model was 0.645 for 3-year DFS and 0.634 for 5-year DFS (Supplementary Fig. 7E,F). The time-dependent ROC curves showed that the predictive accuracy of the metabolic model was also significantly higher than that of the TNM model in both the training and validation groups (Figure 5E,F).
In addition, the clinical value of the metabolic models was evaluated and compared with the TNM models using DCA curves. For almost all threshold probabilities, the metabolic models outperformed the TNM models in terms of net benefit for OS and DFS prediction (Figure 5G,H).
Development and optimization of the decision tree
We generated 2 decision trees without pruning for OS or DFS. Then, the results of cross-validation showed that when the complexity parameter value was adjusted to 0.018 for OS and 0.014 for DFS, the models had the lowest relative error and minimum tree size (Supplementary Fig. 9). Finally, 2 variables (glycolysis model and pathological T) were used in the converged optimal models (Figure 6A,B).
Using the decision trees, HCC patients were classified into the low-, intermediate- and high-risk groups for OS and DFS. We found that the prognosis of the 3 groups was statistically different (Figure 6C,D). The results indicated that the decision trees were useful for the risk stratification of HCC patients.
The association between the metabolic signatures, immunotherapy and immune cell infiltration
To assess the role of the metabolic signatures in immunotherapy, we calculated the association between ICI-related genes (CD274 (PD-L1), PDCD1 (PD-1), CTLA-4, and LAG3) and metabolic risk scores. The expression of CD274, PDCD1, CTLA-4, and LAG3 was significantly higher in the high-risk groups for both glycolysis (Figure 7A–D) and fatty acid metabolic (Figure 7E–H) models for OS.
The CD274, CTLA-4 and LAG3 showed significantly higher expression levels in the high-risk group in the glycolysis model for DFS, whereas there was no difference in PDCD1 expression between high- and low-risk groups (Figure 7I–L). The high-risk group in the fatty acid metabolism model for DFS also had a higher expression of CD274, PDCD1, CTLA-4, and LAG3 (Figure 7M–P). These results suggest that high-risk patients may have a better response to immunotherapy.
We also analyzed the association between core metabolic genes and immune cell infiltration. As expected, most of the core metabolic genes were significantly associated with immune cell infiltration, including B cells, macrophages, myeloid dendritic cells, neutrophils, and CD4+ and CD8+ cells (Figure 8). These results demonstrate the potential impact of the core metabolic genes on the immune microenvironment and immune response in HCC tumors.
We analyzed the association between the core metabolic genes (listed in Figure 8) and TMB and MSI scores. The results of Spearman’s correlation analyses revealed that ABCB6 and PPIA from glycolysis models and UROD and NCAPH2 from fatty acid metabolism models were positively correlated with both TMB and MSI (Table 2, Table 3).
The prognostic values of ABCB6, PPIA, UROD, and NCAPH2 were validated using the ICGC data (Figure 9A). The RNA expression of ABCB6, PPIA, UROD, and NCAPH2 in different HCC tumor cells was visualized in Figure 9B. The RNA expression of ABCB6, PPIA, UROD, and NCAPH2 in HCC was validated using RT-qPCR (Figure 9C). The protein expression of ABCB6, PPIA, UROD, and NCAPH2 in HCC was validated with IHC (Figure 9D–G).
These results suggest that these genes are potential biomarkers for immunotherapy prediction.
Discussion
In this study, the core genes of glycolysis and fatty acid metabolism were screened from the hallmark gene sets in the MSigDB. The prognostic genes associated with OS and DFS were further identified based on the TCGA data. Then, 2 glycolysis metabolism models and 2 fatty acid metabolism models were constructed to predict the OS and DFS. Next, patients were divided into high-risk and low-risk groups according to the risk scores. The multivariate Cox regression analysis was performed using clinicopathological factors to analyze the independent risk factors. Then, metabolic models combining metabolic risk scores with pathological features were constructed based on the obtained independent risk factors, and nomograms were drawn to predict OS and DFS of HCC patients at 3 and 5 years. Compared with the TNM models, the metabolic models were found to improve the predictive accuracy of the TNM system. We also developed 2 decision trees based on machine learning to classify the prognostic risks of HCC patients. Finally, to provide ideas for basic research and clinical application of immunotherapy, the relationship between metabolic core genes and ICI-related markers, as well as the effects on immune cell infiltration, were analyzed.
Several studies presented prediction models related to glucose metabolism or lipid metabolism to predict the OS of HCC patients. For example, Chen et al. built a good prognostic signature with 4 glycolysis-related genes to predict OS in HCC patients.11 Zhang et al. constructed a glycolysis-related gene prognostic model with internal and external validation.27 Xu et al. reported a novel glycolysis-related gene signature to predict the survival and immune status of HCC patients.28 Weng et al. developed a prognostic model based on 10 metabolic genes to reflect the prognosis of HCC patients and the metabolic characteristics of tumors.10 He et al. and Hu et al. built lipid metabolism-related prognostic signatures to improve OS predictions in HCC patients.29, 30 However, few studies have focused on the prediction of DFS. Overall survival is affected by many confounding factors, such as the patient’s age, underlying diseases and social factors. Disease-free survival can better describe the malignant grade of the carcinoma in terms of the biological behavior of the tumor. The present study is a valuable addition to previous studies because it focused on the effects of glycolysis and fatty acid metabolic genes on DFS predictions.
Most of the previous studies have focused on glucose or lipid metabolism.11, 28, 29, 30 However, these are 2 interdependent processes. They communicate with and influence each other through the tricarboxylic acid cycle. Therefore, this study focused on the effects of glycolysis and fatty acid metabolism together on predicting the prognosis. The results indicate that glycolysis and fatty acid metabolism are independent risk factors for OS. Glycolysis alone was an independent risk factor for DFS in HCC patients, and fatty acid metabolism did not statistically correlate with DFS. The above results suggest that the glycolysis core genes may play a more important role in the progression of HCC.
Some studies mainly focused on gene expressions but ignored the effects of clinicopathological features on the prognosis.20, 31 In this study, we tried to eliminate these confounding factors by combining gene expressions with clinicopathological features. The TNM staging system is the most widely used prognostic evaluation system in the clinic. However, it still has some limitations and deficiencies, such as an inability to consider gene expressions of individual patients.32, 33 Gene expression differences are one of the important factors causing tumor heterogeneity. In this study, TNM staging and metabolic risk scores were innovatively combined to predict the OS and DFS of HCC patients, taking clinicopathological features and gene expression differences into account. In addition, the head-to-head comparison with the TNM models showed that metabolic models could improve the predictive accuracy of the TNM models, which increases the clinical application value of this study.
Immunotherapy is a recent breakthrough in the treatment of advanced HCC.34, 35 However, due to tumor heterogeneity and immune escape, many patients do not respond to immunotherapy. Microsatellite instability status, TMB and PD-L1 expression can be used to screen for beneficiaries of immunotherapy.36 However, the positive rate of these indicators was too low to meet the patients’ needs. In this study, we observed that the expression of ICI-related markers was significantly increased in patients from the high metabolic risk groups, and the expression of metabolic core genes was correlated with immune cell infiltration. These results suggest that metabolic-related risk scores have potential feasibility as biomarkers for predicting the effects of immunotherapy. Furthermore, the abnormal expression of metabolic core genes may accelerate the development of HCC by affecting the components of the tumor immune microenvironment.
Limitations
This study has several limitations. First, external validation of independent data played a large role in model evaluation. However, because the follow-up data for DFS was missing in the ICGC and we did not find Gene Expression Omnibus (GEO) data with sufficient RNA-seq and DFS data, we could not perform the external validation for the DFS models. Second, only PCR and IHC were used for experimental validation. It is necessary to further explore the biological function and mechanism of core metabolic genes in HCC, especially their effects on the immune microenvironment and immunotherapy. Finally, we mainly focused on genes related to glycolysis and fatty acid metabolism, but not amino acid metabolism. Previous studies have found that glutamine metabolism is also closely involved in the progression of HCC.12, 37
Conclusions
The novel prognostic models based on glycolysis- and fatty acid metabolic-related genes were developed to better predict the outcome and response of HCC patients to immunotherapy. The new metabolic models were found to improve the predictive accuracy of the TNM system. Further biological experiments are required to explore the mechanisms core metabolic genes play in HCC, especially their effects on the immune microenvironment and immunotherapy.
Supplementary material
The supplementary materials are available at https://doi.org/10.5281/zenodo.7863602. The package contains the following files:
Supplementary Fig. 1. Forest plots of the prognostic genes for OS based on glycolic hallmark gene set.
Supplementary Fig. 2. Forest plots of the prognostic genes for DFS based on glycolic hallmark gene set.
Supplementary Fig. 3. Expression of overlapping genes between tumor and normal tissues in HCC for glycolysis (p < 0.001).
Supplementary Fig. 4. Forest plots of the prognostic genes for OS based on fatty acid metabolic hallmark gene set.
Supplementary Fig. 5. Forest plots of the prognostic genes for DFS based on fatty acid metabolic hallmark gene set.
Supplementary Fig. 6. Expression of overlapping genes between tumor and normal tissues in HCC for fatty acid metabolism (p < 0.001).
Supplementary Fig. 7. A. Calibration plot of the metabolic prognostic model for OS; B. Calibration plot of the metabolic prognostic model for DFS; C. ROC curves of the metabolic and TNM models for predicting 3-year OS; D. ROC curves of the metabolic and TNM model for predicting 5-year OS; E. ROC curves of the metabolic and TNM models for predicting 3-year DFS; F. ROC curves of the metabolic and TNM model for predicting 5-year DFS.
Supplementary Fig. 8. Construction of the TNM nomograms to predict OS and DFS. A. TNM model for OS; B. TNM model for DFS.
Supplementary Fig. 9. Cost-complexity plot for the decision-tree models. A. Plot for OS; B. Plot for DFS.
Supplementary Table 1. Hallmark gene sets of glycolysis and fatty acid metabolism.
Supplementary Table 2. Primers used in this study.
