PDF

J Cancer 2026; 17(7):1295-1317. doi:10.7150/jca.135837 This issue Cite

Research Paper

A Prognostic Risk Model for Hepatocellular Carcinoma Integrating Ferroptosis and Metabolic Reprogramming Signatures

Hao Ling , Yanzhu Hu

Department of Surgery, Klinikum rechts der Isar, TUM School of Medicine and Health, Technical University of Munich, Munich, Germany.

✉ Corresponding authors: Hao Ling (Email: hao.lingde) and Yanzhu Hu (Email: yanzhu.hude)

Received 2026-4-8; Accepted 2026-6-17; Published 2026-7-10

Abstract

Background: Hepatocellular carcinoma (HCC) continues to impose a heavy global health burden, with high incidence and mortality. The disease is highly heterogeneous and is commonly detected at late stages, which compromises treatment outcomes. Ferroptosis and metabolic reprogramming are increasingly recognized as key processes in HCC development; however, their roles in disease progression and therapeutic response remain incompletely understood. This research aimed to identify genes related to ferroptosis and metabolic reprogramming (FPMRRGs) that may serve as putative biomarkers and therapeutic targets in HCC.

Methods: The Cancer Genome Atlas (TCGA), including 369 HCC specimens and 50 normal controls, along with two Gene Expression Omnibus (GEO) datasets (GSE10143 and GSE76427), were analyzed using R (v4.3.3). From a curated list of 451 FPMRRGs, differentially expressed genes (DEGs) between tumor and normal tissues were identified. Univariate Cox regression analysis was then conducted to explore their prognostic relevance and to define molecular subtypes of HCC. Specimens were categorized into 2 subtypes using ConsensusClusterPlus, and overall survival differences were evaluated via survival analysis. Functional and pathway enrichment analyses were conducted to investigate the functional roles of these genes. Immune-related features were evaluated using the Mann-Whitney U test. A prognostic risk model was constructed using least absolute shrinkage and selection operator (LASSO) regression followed by multivariate Cox analysis. Model performance was assessed using receiver operating characteristic (ROC) curves and calibration plots. Immune cell infiltration was estimated by single-sample GSEA, and pathway activity differences were examined using gene set variation analysis (GSVA).

Results: HCC specimens were divided into 2 molecular subtypes, which demonstrated obvious differences in overall survival and immune-related features, including immune checkpoint gene expression and tumor immune dysfunction and exclusion (TIDE) scores. A prognostic model based on 12 key FPMRRGs demonstrated good predictive performance for 1- and 3-year overall survival, with moderate performance for 5-year survival. The prognostic value and expression patterns of these genes were further validated across independent datasets. In addition, these genes were mainly enriched in pathways linked to fatty acid metabolism and HIF-1 signaling, and were closely associated with patterns of immune cell infiltration.

Conclusions: This study identified numerous key genes linked to ferroptosis and metabolic reprogramming in HCC and developed a robust prognostic risk model. Our results offer new insight into the molecular basis of HCC and highlight potential biomarkers for more individualized treatment approaches. Further studies, particularly those combining clinical validation with functional experiments, are required to verify these findings and examine the therapeutic potential of targeting these pathways.

Keywords: hepatocellular carcinoma, ferroptosis, metabolic reprogramming, prognostic model, immune checkpoint genes

Introduction

Hepatocellular carcinoma (HCC) is the most common type of primary liver cancer and remains a major global health concern due to its high incidence and mortality. It accounts for nearly 90% of all liver cancer cases worldwide, with increasing rates driven largely by chronic liver diseases such as viral hepatitis, alcohol-related liver disease, and, more recently, metabolic dysfunction-associated steatohepatitis [1]. Managing HCC remains difficult, as the disease is often aggressive, diagnosed at an advanced stage, and frequently resistant to conventional therapies [2, 3]. Although progress has been made in surgical resection, liver transplantation, and locoregional treatments, outcomes for patients with advanced disease are still poor. This is partly due to the complex and multifactorial nature of HCC, which involves genetic, environmental, and metabolic factors. Together, these challenges highlight the need to identify reliable biomarkers and develop more effective targeted therapies [4, 5].

Recently, increasing attention has been paid to the molecular and metabolic changes that drive HCC development. Ferroptosis, an iron-dependent form of regulated cell death driven by lipid peroxidation, has gained increasing attention in cancer research [6, 7]. In HCC, ferroptosis has been linked to tumor progression, shaping of the tumor microenvironment, and responses to therapy, including immunotherapy [8-10]. At the same time, metabolic reprogramming allows tumor cells to adapt to changing conditions by reshaping energy production and biosynthetic pathways, thereby supporting tumor growth and survival. This phenomenon is well recognized across cancers and has also been described in liver malignancies [11, 12]. Several recent studies have focused specifically on metabolic alterations in HCC, particularly their interaction with the immune system [13-18]. In addition, disruptions in iron metabolism and lipid processing have been shown to contribute to tumor progression and treatment resistance [19, 20]. Despite these advances, the relationship between ferroptosis and metabolic reprogramming in HCC is still not fully understood, and their combined impact across different HCC subtypes remains unclear.

Herein, we sought to identify key genes related to ferroptosis and metabolic reprogramming that are associated with HCC heterogeneity and progression. Using transcriptomic data from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO), we carried out a comprehensive bioinformatics analysis to identify relevant genes. In accordance with these findings, we constructed a prognostic risk model using least absolute shrinkage and selection operator (LASSO) regression combined with Cox analysis to stratify patients according to survival risk. We also examined immune checkpoint gene expression to better characterize the interplay between tumor metabolism, ferroptosis, and immune regulation. A set of 12 key genes was identified, showing strong associations with both patient survival and the immune landscape of HCC. These genes may act as useful biomarkers for patient stratification and offer potential targets for more precise therapeutic strategies aimed at improving clinical outcomes.

Materials and Methods

Dataset Acquisition and Preprocessing

We obtained HCC transcriptomic data from TCGA via TCGAbiolinks (v2.30.0) [21] (https://portal.gdc.cancer.gov/). Following exclusion of patients with missing clinical information, 369 HCC tumor specimens and 50 normal control specimens were retained for analysis. Additional datasets, GSE10143 and GSE76427, comprising 195 HCC specimens and 359 control specimens, were downloaded from the GEO using the GEOquery package [22] (https://www.ncbi.nlm.nih.gov/geo/) (Table 1). All datasets were normalized and annotated using the limma package (v3.58.1) [23]. Ferroptosis-related genes were collected from multiple sources. First, 2,023 genes were obtained from GeneCards. In addition, a literature-based search on PubMed using the keyword “ferroptosis” [24] yielded 144 genes. After merging these lists and removing duplicates, 2,059 ferroptosis-related genes were identified. Similarly, metabolic reprogramming-related genes (MRRGs) were collected by searching “metabolic reprogramming” and retaining only protein-coding genes, resulting in 1,951 genes. An additional literature-based search [25] identified 9 more MRRGs, and after merging and removing duplicates, 1,954 MRRGs were obtained (Table S1). Comparing the ferroptosis and metabolic reprogramming lists, we identified 451 genes associated with both processes, referred to as ferroptosis- and metabolic reprogramming-related genes (FPMRRGs). Intersecting these with the TCGA-LIHC DEGs yielded 333 ferroptosis- and metabolic reprogramming-related differentially expressed genes (FPMRRDEGs).

Identification of Differentially Expressed FPMRRDEGs in HCC

The TCGA-LIHC dataset was assigned to HCC and control groups. Differentially expressed genes (DEGs) between tumor and normal specimens were identified using DESeq2 [26] with thresholds of |log2 fold-change| > 0.25 and p < 0.05. Genes with log2 fold-change > 0.25 were considered upregulated, whereas genes with log2 fold-change < -0.25 were considered downregulated. The results were visualized using volcano plots produced with ggplot2. FPMRRDEGs associated with HCC were determined by overlapping DEGs with the FPMRRGs. Univariate Cox regression analysis (p < 0.05) was then applied to identify genes relevant to HCC prognosis, which were used for subsequent subtype classification. A sensitivity analysis varying the |log2 fold-change| cutoff from 0.25 to 1.5 confirmed that the 12-gene panel is preserved at 0.25; at stricter cuts of 0.5, 1.0, and 1.5, 11, 6, and 5 genes are retained respectively, with refit Cox C-indices of 0.712, 0.698, and 0.694 versus 0.710 for the published 12-gene model (Table S6, Table S6b, Figure S6).

Identification of HCC Subtypes

Consensus clustering [27] was performed on the TCGA-LIHC dataset using expression profiles of FPMRRGs to identify HCC subtypes. Two clusters were defined, with the process repeated 50 times using 80% of the samples to ensure stability. Gene expression differences between subtypes were visualized with heatmaps, and overall survival differences were examined via Kaplan-Meier (KM) analysis.

Differential Expression Analysis and GSEA of HCC Subtypes

Using the TCGA-LIHC dataset, DESeq2 was applied to compare gene expression levels between HCC subtypes and normal liver tissue. DEGs were defined using thresholds of |log2 fold-change| > 1 and p < 0.05. Volcano plots of DEGs were produced with ggplot2 (v3.4.4). Overlapping genes with FPMRRGs were identified via a Venn diagram, yielding the set of cluster_FPMRRDEGs for subsequent analyses. GSEA was performed using clusterProfiler (v4.10.0) to evaluate ranked DEGs based on log2 fold change. Parameters included 2022 seeds, 1,000 permutations, and gene sets containing 10-500 genes, using the c2 collection from MSigDB [28]. Significance was defined as adjusted p < 0.05 and false discovery rate (FDR) < 0.05. Functional enrichment of FPMRRDEGs was further assessed using GO and KEGG analyses in clusterProfiler, with p < 0.05 deemed statistically significant.

Analysis of Immune Activity in HCC Subtypes

Immune checkpoint genes (ICGs) play key roles in modulating immune function and preserving immune balance. A total of 47 ICGs were obtained from the published literature [29] (Table S2). Expression differences of these genes among HCC subtypes were examined using the Mann-Whitney U test in the TCGA-LIHC dataset. The TIDE framework (http://tide.dfci.harvard.edu) was used to evaluate immune function and exclusion scores in each HCC subtype. Differences in TIDE scores between subtypes were analyzed with Mann-Whitney U [30] (https://www.cbioportal.org/), and differences in TMB and MSI across subtypes were similarly tested with Mann-Whitney U. Immunogenicity scores were calculated using the Immuno-Oncology Biological Research (IOBR) R package (v0.99.9), and subtype-specific differences were examined with Mann-Whitney U.

Establishment of a Prognostic Model for HCC

A prognostic model was established in accordance with FPMRRG expression. Candidate genes were first screened by univariable analysis and then refined using LASSO for variable selection. Patients were stratified into risk groups on the basis of the resulting RiskScore, and overall survival differences were examined via KM analysis. The predictive performance of the model was assessed with time-dependent ROC curves to compute the AUC values for 1/3/5-year survival. Calibration curves and decision curve analysis (DCA) were further employed to assess model precision and clinical applicability. Internal model stability was further assessed by 10-fold cross-validation × 10 repeats and by 1000-iteration bootstrap optimism correction (Harrell's enhanced bootstrap).

KM Curve Validation and Expression Analysis of Key Genes

To confirm the robustness of the identified HCC subtypes and their associated overall survival differences, KM survival analyses were conducted with survival (v3.5-7) on the TCGA, GSE10143, and GSE76427 datasets. In addition, DEG expression levels were evaluated between HCC patients and healthy subjects in the GSE10143 and GSE76427 datasets using group comparison plots to visualize differential expression.

Immune activity Analysis in HCC and Normal Samples

ssGSEA was implemented on the TCGA dataset to quantify the abundance of immune cell populations, generating an immune infiltration matrix for HCC. Expression differences between HCC and healthy subjects were visualized with ggplot2, allowing identification of immune cell populations with significant variation. Spearman correlation analysis was applied to measure relationships among immune cells, as well as between key genes and immune cell populations. Correlation results were visualized using ggplot2 (v3.4.4) and pheatmap (v1.0.12).

Functional immune signature scoring and single-cell validation

To evaluate functional immune states beyond cell-type abundance, we applied the singscore package to compute rank-based enrichment scores for 11 curated functional immune gene sets (including effector function, antigen presentation, and suppressive states) in TCGA-LIHC samples. Gene sets were derived from published literature and the MSigDB database. Differences in signature scores between low- and high-risk groups (defined by median RiskScore) were assessed using the Wilcoxon rank-sum test with Benjamini-Hochberg correction for multiple testing.

For single-cell validation, we utilized the GSE140228 dataset (15,319 CD45+ tumour-infiltrating immune cells from four HCC patients) retrieved via the TISCH2 database. Cell-type annotation was obtained directly from TISCH2. Expression of the 12 signature genes was visualized across major immune lineages and 28 finer subtypes using dot plots. No additional data processing or re-annotation was performed beyond the published TISCH2 annotations.

GSVA of Low- and High-risk HCC Subtypes

GSVA, a nonparametric approach, was applied to evaluate pathway enrichment in TCGA and GEO gene expression matrices. Using the H.all.v7.4.Symbols.gmt gene set from MSigDB (v1.50.0), HCC specimens were stratified into low- and high-risk groups. GSVA was used to identify pathways and genes enriched in the high-risk group relative to the low-risk group, with significance defined as p < 0.05.

Statistical Analysis

All statistical tests were implemented through R (v4.3.3). Continuous data are expressed as mean±standard deviation (SD). Two-group comparisons were conducted with Wilcoxon test, whereas comparisons among three or more groups used the Kruskal-Wallis test. Categorical data were evaluated using Fisher's exact or chi-square tests, as appropriate. Correlations between variables were calculated using Spearman's method unless otherwise specified. P<0.05 was deemed statistically significant. The proportional hazards assumption of the multivariate Cox model was assessed using Schoenfeld residuals. Although several individual covariates showed violation of the proportional hazards assumption, the global test was non-significant (p = 0.171), supporting the overall validity of the model (Table S7b).

Results

Technology Roadmap

The overall study design and analysis workflow are summarized in Figure 1.

 Figure 1 

Technology roadmap.

FPMRRGs Are Differentially Expressed in HCC

To identify HCC-associated FPMRRGs, we first examined differences in gene expression between tumor and healthy liver tissues. Overall, 9,881 DEGs were identified, involving 1,219 downregulated and 8,662 upregulated genes (Figure 2A). The top 20 DEGs are shown in the heatmap (Figure 2B). We then compiled a list of 451 genes related to both ferroptosis and metabolic reprogramming based on data from GeneCards and the published literature. Intersecting this list with the TCGA-LIHC DEGs yielded 333 overlapping genes (FPMRRDEGs) (Figure 2C). Among these, 15 genes with p < 0.001 (Table S3) were further analyzed using univariate Cox regression (Figure S1). The hazard ratios for these genes ranged from 1.22 to 1.82. Overall, these results indicate that a subset of genes involved in ferroptosis and metabolic reprogramming is significantly dysregulated in HCC compared with normal liver tissue.

 Figure 2 

Differential expression of FPMRRGs in HCC. (A) Volcano plot showing gene expression differences between HCC and control samples in the TCGA-LIHC dataset. (B) Heatmap of the top 20 DEGs. (C) Venn diagram revealing the overlap between TCGA-LIHC DEGs and FPMRRGs. TCGA: The Cancer Genome Atlas; LIHC: liver cancer; FPMRRGs: ferroptosis- and metabolic reprogramming-related genes.

Distinct HCC Subtypes Defined by Ferroptosis- and Metabolic Reprogramming-related Genes Show Different Clinical Features and Outcomes

We next examined whether expression patterns of FPMRRGs varied across HCC samples. Using the 15 prognostically significant FPMRRDEGs (p < 0.001), consensus clustering was implemented on the TCGA-LIHC cohort. This analysis identified two distinct subtypes: subtype 1 (n = 113) and subtype 2 (n = 256) (Figure 3A-C). The clustering results are further illustrated by a three-dimensional t-SNE plot (Figure 3D). A heatmap highlighting differences in gene expression between the 2 subtypes is shown in Figure 3E. Subsequent analysis confirmed that the 15 genes used for clustering were differentially expressed between the two subtypes (Figure 3G). To assess the clinical relevance of these subtypes, KM survival analysis was performed, revealing a significant difference in overall survival between the two groups (Figure 3F). We also compared clinical features between the subtypes (Table 2). Significant differences were observed in pathologic stage and T stage distribution, whereas no obvious differences were noted for N stage, M stage, gender or age (Figure 4). Collectively, these findings imply that HCC can be stratified into biologically distinct subtypes based on ferroptosis- and metabolic reprogramming-related gene expression, and that these subtypes differ in both clinical features and patient outcomes.

 Figure 3 

Consensus clustering identifies two molecular subtypes of HCC. (A) Consensus clustering of HCC samples. (B) Cumulative distribution function (CDF) plot of clustering results. (C) Delta area plot displaying the relative changes in CDF. (D) Three-dimensional t-SNE plot illustrating the distribution of HCC subtypes. (E) Heatmap revealing the expression of 15 FPMRRDEGs across the two subtypes. (F) KM analysis of overall survival for patients in the two subtypes. (G) Expression levels of the 15 FPMRRDEGs in each subtype. t-SNE: t-distributed stochastic neighbor embedding; CDF: cumulative distribution function; HR: hazard ratio; HCC: hepatocellular carcinoma. ***p < 0.001.

 Figure 4 

Clinical characteristics differ between the two HCC subtypes. Distribution of patients with subtype 1 and subtype 2 according to pathologic stage (A), T stage (B), N stage (C), M stage (D), gender (E), and age (F).

Differentially Expressed FPMRRGs between HCC Subtypes are Enriched in Key Biological Functions and Pathways

Given the clinical differences observed between the two HCC subtypes, we next explored the underlying functional differences that might explain these findings. Differential expression analysis between the 2 subtypes in TCGA-LIHC identified a set of DEGs (Table S4), which were subsequently analyzed using GSEA. GSEA indicated that these genes were markedly enriched in pathways and gene sets associated with liver cancer proliferation, periportal HCC, patient survival, and tumor recurrence (Figure 5, Table 3). This suggests that the molecular divergences between the subtypes are closely linked to key features of tumor progression and clinical outcome.

 Figure 5 

Enrichment of key biological pathways in HCC subtypes. (A) GSEA enrichment plot highlighting significantly enriched biological processes. (B-E) Enrichment plots for representative gene sets: CHIANG_LIVER_CANCER_SUBCLASS_PROLIFERATION_DN (B), DESERT_PERIPORTAL_HEPATOCELLULAR_CARCINOMA_SUBCLASS_UP (C), LEE_LIVER_CANCER_SURVIVAL_UP (D), and WOO_LIVER_CANCER_RECURRENCE_DN (E). Significance criteria: adjusted p-value < 0.05 and FDR < 0.05 (Benjamini-Hochberg correction). NES: Normalized Enrichment Score; BH: Benjamini-Hochberg.

To further focus on ferroptosis- and metabolic reprogramming-related genes, we intersected the subtype-specific DEGs with the previously defined list of 451 FPMRRGs, yielding 52 overlapping genes (Figure S2). Functional annotation of these genes showed enrichment across pathways related to metabolism, transport, and cellular signaling (Table 4, Figure 6A). Network analysis provided additional insight into these functional categories. Biological process terms were mainly associated with stress responses, while cellular component terms were enriched in nuclear and basolateral cell regions. Molecular function terms were largely related to fatty acid metabolism and enzyme activity, and several key metabolic and signaling pathways were identified in the KEGG analysis (Figure 6B-E). Together, these findings indicate that differences in the expression of FPMRRGs between HCC subtypes are functionally relevant and may contribute to the distinct biological behavior and clinical outcomes observed in these tumors.

 Figure 6 

GO and KEGG enrichment analyses of subtype-specific FPMRRDEGs. (A) Bubble plot showing enrichment data for the 52 differentially expressed FPMRRDEGs. (B-E) Network plots of enriched biological processes (B), cellular components (C), molecular functions (D), and KEGG (E) pathways. Node size represents the number of genes, and edges indicate functional associations. Significance threshold: p < 0.05.

The Two HCC Subtypes Show Distinct Immune Profiles and Mutation Characteristics

Given the key role of immune regulation and genomic alterations in HCC progression, we next compared immune activity and mutation-related features between the two subtypes. Analysis of ICG expression revealed significant differences between the subtypes, including ADORA2A, BTLA, and IDO1 (all p < 0.01), as well as CD40LG (p < 0.05) (Figure 7A). To assess potential responses to immunotherapy, we applied the TIDE algorithm and found an obvious difference in TIDE scores among the groups (p < 0.01; Figure 7B), suggesting distinct immune evasion profiles. Immunophenoscore (IPS) analysis further showed significant differences in effector cell activity between the subtypes (p < 0.01; Figure 7E). We also examined mutation-related features. Both tumor mutational burden (TMB) and microsatellite instability (MSI) differed remarkably across the subtypes (MSI, p < 0.05; Figure 7C; TMB, Figure 7D). Overall, these results indicate that the two HCC subtypes are characterized by distinct immune landscapes and mutation patterns, which may contribute to their different clinical behaviors and treatment responses.

 Figure 7 

Immune activity and mutation characteristics differ between HCC subtypes. (A) Comparison of ICG expression among subtypes in TCGA-LIHC. (B-D) Comparison of TIDE (B), MSI (C), and TMB (D) scores across the two subtypes. (E) Immunophenoscore (IPS) analysis showing differences in immune-related components between subtypes. ICG: immune checkpoint gene; TIDE: tumor immune dysfunction and exclusion; MSI: microsatellite instability; TMB: tumor mutation burden; IPS: immuno pheno score; EC: effector cells; MHC: major histocompatibility complex; SC: stimulatory cells; CP: checkpoint molecules. *p < 0.05; **p < 0.01; *** p< 0.001.

Development and Verification of a Prognostic Risk Model Based on FPMRRGs

Building on the observed roles of ferroptosis and metabolic reprogramming in HCC, we established a prognostic model in accordance with gene expression profiles. Univariate Cox regression (Figure 8A) followed by LASSO regression (Figure 8B, C) identified a 12-gene signature: G6PD, SLC1A5, NAP1L1, SLC16A3, KPNA2, CCNB1, SLC2A1, YWHAZ, TPX2, ANXA5, CDK4, and ETV4. Time-dependent ROC analysis indicated strong predictive performance of the model for 1/3/5-year overall survival (Figure 9A). Calibration curves supported good model calibration, with the highest accuracy at 5 years (Figure 9D-F). DCA further supported the clinical applicability of the model, again showing the strongest performance at 5 years (Figure 9G-I). Internal validation by 10-fold cross-validation × 10 repeats yielded a mean C-index of 0.676 (95% CI 0.491-0.823), and 1000-iteration bootstrap optimism correction yielded a corrected C-index of 0.679 (apparent C = 0.710, optimism = 0.031) (Figure 9J-K, Table S7). The proportional hazards assumption of the multivariate Cox model was assessed using Schoenfeld residuals. Although several individual covariates violated the proportional hazards assumption, the global test was non-significant (p = 0.171), supporting the overall validity of the model (Table S7b). A sensitivity analysis confirmed model robustness across DEG fold-change thresholds (Table S6, Table S6b, Figure S6).

 Figure 8 

A prognostic model based on FPMRRGs accurately predicts survival in HCC patients. (A) Forest plot of genes identified by univariate Cox regression analysis. (B) Construction of the prognostic risk model using LASSO regression. (C) LASSO coefficient profiles showing variable trajectories. LASSO, least absolute shrinkage and selection operator; HR: hazard ratio.

 Figure 9 

Prognostic model validation. (A) Time-dependent ROC curves for predicting overall survival in HCC. (B) Kaplan-Meier survival analysis comparing low- and high-risk groups. (C) Nomogram integrating risk score and clinical variables based on Cox regression analysis. (D-F) Calibration curves for 1-, 3-, and 5-year survival, respectively. (G-H) Decision curve analysis (DCA) for 3- and 5-year overall survival in TCGA-LIHC, comparing RiskScore alone, AJCC TNM alone, and combined RiskScore + TNM strategies against treat-all and treat-none defaults. (I) DCA for 3- and 5-year overall survival in GSE76427, comparing Local-refit RiskScore alone, BCLC stage alone, and combined Local-refit RiskScore + BCLC strategies. (J) Boxplot of 10-fold cross-validation C-index and 1/3/5-year time-dependent AUC across 100 fold-level estimates (10 folds × 10 repeats) in TCGA-LIHC; the dashed line marks the apparent C-index of 0.710. (K) Apparent C-index versus bootstrap optimism-corrected C-index (B = 1000), showing minimal optimism (0.031) and confirming low overfitting.

Patients were then stratified into low- and high-risk groups in accordance with the median RiskScore. KM analysis indicated an obvious difference in overall survival between these groups (Figure 9B). GSVA identified multiple pathways that differed between the low- and high-risk groups (Table 5, Figure S3A), and these differences were confirmed through Mann-Whitney U. Notably, pathways such as HALLMARK_BILE_ACID_METABOLISM were significantly enriched (Figure S3B). Quantitative Decision Curve Analysis at the 3-year horizon in TCGA-LIHC showed that the combined RiskScore + TNM model adds a net benefit of approximately +0.06 over TNM alone at threshold 0.50; in GSE76427 at the 5-year horizon, the Local-refit RiskScore adds a net benefit of approximately +0.13 over BCLC alone at threshold 0.20 (Figure 9G-I, Table S9).

A nomogram comprising the LASSO-derived risk scores and clinical variables was established to improve prognostic prediction. The risk score showed stronger predictive value than other clinical factors, whereas age contributed the least (Figure 9C). When compared against AJCC TNM staging in TCGA-LIHC (complete-stage sub-cohort n = 339, 116 events; 24 patients with missing AJCC pathologic stage were excluded to permit paired model comparison), the combined RiskScore + TNM model achieved a C-index of 0.717 (95% CI 0.668-0.765) versus 0.660 (95% CI 0.586-0.734) for TNM alone (paired p = 0.039), with 5-year continuous net reclassification improvement (NRI) = 0.247 (95% CI 0.035-0.442) and integrated discrimination improvement(IDI) = 0.099 (95% CI 0.034-0.170). In GSE76427 (n = 115, events = 23) against BCLC, the Local-refit RiskScore + BCLC model achieved a C-index of 0.812 (95% CI 0.740-0.883) versus 0.719 (95% CI 0.556-0.881) for BCLC alone, with 5-year continuous NRI = 0.621 (95% CI 0.004-1.003) and IDI = 0.374 (95% CI 0.050-0.632) (Table 6).

To validate the multi-gene risk model in external cohorts, we computed RiskScore in GSE10143 and GSE76427 using two complementary strategies. Direct transfer of TCGA-derived Cox coefficients (TCGA-β) failed to separate risk groups in either cohort (GSE10143 OS HR = 0.69, 95% CI 0.34-1.40, p = 0.30; GSE76427 OS HR = 0.72, 95% CI 0.32-1.65, p = 0.43), indicating that the TCGA-fit weights are not portable across microarray platforms. Local refitting of the 12-gene Cox model within each cohort, however, recovered strong prognostic separation (GSE10143 OS HR = 3.09, 95% CI 1.47-6.51, p = 0.002, C-index = 0.65; GSE76427 OS HR = 13.12, 95% CI 3.70-46.53, p < 0.001, C-index = 0.81; GSE76427 RFS HR = 3.05, 95% CI 1.65-5.67, p < 0.001, C-index = 0.73), confirming that the 12-gene panel carries genuine prognostic signal but requires cohort-specific recalibration (Figure 10, Table S8). Single-gene Kaplan-Meier curves for each of the 12 individual genes in GSE10143 and GSE76427 are provided in Supplementary Figures S4 and S5, respectively.

 Figure 10 

Multi-cohort validation of the 12-gene risk model. (A) Kaplan-Meier curve for the TCGA-LIHC training cohort (n = 363, events = 130) comparing low- and high-risk groups stratified by the median RiskScore. (B-D) External validation in GEO cohorts using cohort-specific local-refit Cox models: (B) GSE10143 overall survival (n = 80, events = 32), (C) GSE76427 overall survival (n = 115, events = 23), and (D) GSE76427 recurrence-free survival (n = 108, events = 48). Hazard ratios with 95% confidence intervals and log-rank p-values are shown in each panel. (E) Time-dependent AUC curves across the three cohorts over 5-year follow-up, demonstrating consistent discriminative performance after cohort-specific recalibration.

Taken together, these findings suggest that this 12-gene signature provides a robust prognostic model for HCC and may be useful for risk stratification and guiding clinical decision-making.

Prognostic Model Genes are Differentially Expressed in HCC and Related to Immune Cell Infiltration

We next examined whether the genes included in the prognostic model could also distinguish tumor tissue from normal liver tissue. Expression levels of the 12 key genes (G6PD, SLC1A5, NAP1L1, SLC16A3, KPNA2, CCNB1, SLC2A1, YWHAZ, TPX2, ANXA5, CDK4, and ETV4) were compared between HCC and control samples in the GSE10143 and GSE76427 datasets. In the GSE76427 dataset, significant differences were observed for NAP1L1, KPNA2, CCNB1, YWHAZ, TPX2, CDK4, and ETV4 (Figure 11A). Similarly, in the GSE10143 dataset, G6PD, NAP1L1, KPNA2, CCNB1, TPX2, ANXA5, CDK4, ETV4, and SLC16A3 showed significant differential expression between tumor and normal tissues (Figure 11B). These findings indicate that a substantial proportion of the model genes are consistently dysregulated in HCC.

 Figure 11 

Differential expression of prognostic model genes in HCC and healthy specimens. (A) Gene expression differences between HCC and healthy specimens in GSE76427. (B) Gene expression differences between HCC and healthy specimens in GSE10143. *p < 0.05, **p < 0.01, ***p < 0.001.

To explore whether these expression changes were associated with immune activity, immune cell infiltration in TCGA-LIHC was analyzed using ssGSEA. Significant differences were observed in the abundance of 25 immune cell types between HCC and healthy controls (Figure 12A). Correlation analysis revealed complex relationships among these immune cell populations (Figure 12B). Further analysis using a bubble plot demonstrated remarkable associations between key gene expression and immune cell infiltration (Figure 12C). Notably, SLC1A5 expression was positively correlated with activated CD4⁺ T cells (r = 0.564), while CCNB1 expression had a negative correlation with eosinophils (r = -0.567). This suggests that the genes incorporated in the prognostic model are not only differentially expressed in HCC but are also closely linked to the tumor immune microenvironment.

 Figure 12 

Immune cell infiltration analysis. (A) Comparison of immune cell populations between HCC and control samples using ssGSEA. (B) Correlation heatmap of immune cell infiltration in TCGA-LIHC. (C) Bubble plot showing correlations between key gene expression and immune cell infiltration. Significance levels: ns (p ≥ 0.05), *p < 0.05, **p < 0.01, ***p < 0.001.

Functional Immune Profiling and Single-Cell Validation of the 12-Gene Signature

To extend the immune analysis beyond cell-type abundance, we scored 11 functional immune signatures in TCGA-LIHC by singscore[31] and stratified samples by the median RiskScore. High-risk tumors exhibited significantly elevated MHC class I, Treg, T-cell exhaustion, and M2 TAM signatures, together with significantly reduced cytotoxicity after Benjamini-Hochberg correction (Figure 13A and Table S10). In contrast, MDSC, CD8 T cell, and IFN-γ response signatures showed the same directional changes but did not reach statistical significance after multiple-testing correction. These results indicate that high-risk tumors are characterized by a more exhausted and suppressive immune phenotype, despite relatively preserved MHC class I expression (Wilcoxon BH-adjusted p < 0.05; Figure 13A, Table S10). The bubble plot in Figure 13B summarises the magnitude and direction of these differences, with suppressive signatures positively associated with the risk score and effector signatures negatively associated.

 Figure 13 

Functional immune profiling of the 12-gene risk model in TCGA-LIHC. (A) Single-sample gene set scoring (singscore) of 11 functional immune signatures (CD8 T cell, cytotoxicity, IFN-γ response, M2 TAM, MDSC, MHC-I, NK cell, T-cell exhaustion, T-cell inflamed GEP, tertiary lymphoid structure, Treg) compared between low- and high-risk groups; boxplots show median and IQR, Wilcoxon p-values are annotated and BH-adjusted significance is indicated (*p < 0.05, **p < 0.01, ***p < 0.001). (B) Bubble plot summarising effect direction and significance of risk-group differences; the x-axis shows the delta score (high minus low), dot size encodes -log10 (adjusted p), and colour groups signatures by functional category (effector/IFN, antigen presentation/TLS, suppressive).

To resolve which immune-cell compartments express the signature genes, we examined the HCC tumour-infiltrating immune-cell atlas GSE140228[32] (15,319 cells from four donors, retrieved through TISCH2[33]; Figure 14, Table S11). The cell-cycle genes CCNB1, TPX2, KPNA2 and CDK4 were predominantly expressed in cycling T/NK cells, the metabolic/redox genes G6PD, SLC2A1 and SLC1A5 in NK and myeloid cells, and the transcription factor ETV4 was virtually absent from all immune compartments, consistent with its tumour-cell-intrinsic role. A more granular annotation across 28 finer subsets is provided in Figure S7 and Table S11.

 Figure 14 

Single-cell validation of the 12 signature genes in the HCC immune microenvironment (GSE140228, retrieved through TISCH2). Dot plot showing the expression of the 12 risk-model genes (G6PD, SLC2A1, SLC1A5, SLC16A3, CCNB1, CDK4, TPX2, KPNA2, NAP1L1, ETV4, YWHAZ, ANXA5) across seven major immune-cell compartments (T cells, NK cells, cycling T/NK, B cells, plasma cells, myeloid cells, mast cells). Dot size encodes the percentage of cells expressing each gene; dot colour encodes the mean log2(CP10K + 1) expression. Genes are grouped by functional layer (metabolism/transport, cell cycle, ferroptosis/redox). ETV4 shows near-absent expression in immune cells, consistent with its tumour-cell-intrinsic role.

Discussion

Herein, we show that genes related to ferroptosis and metabolic reprogramming have clear prognostic value in HCC. These genes were differentially expressed between tumor and healthy liver tissues and also distinguished two molecular subtypes of HCC with distinct clinical features and survival outcomes. The subtypes differed not only in prognostic outcomes but also in immune activity and mutation-associated features. According to these findings, we developed a prognostic risk model using FPMRRGs, which demonstrated excellent predictive performance for 1/3/5-year overall survival. Together, these results suggest that this group of genes may function as biomarkers to stratify risk and as potential therapeutic targets.

The 12 genes included in our model are associated with numerous biological processes known to contribute to HCC progression. For example, SLC1A5, SLC2A1, and G6PD influence cellular metabolism and redox balance, thereby indirectly regulating susceptibility to ferroptosis through effects on amino acid and glucose metabolism, as well as NADPH and glutathione levels [34-37]. ANXA5 has been linked to apoptosis, epigenetic regulation, and membrane transport, as well as tumor invasion, angiogenesis, and immune cell infiltration in HCC [38-42]. NAP1L1 promotes tumor progression through multiple mechanisms, including ubiquitination of BIRC2, activation of Wnt signaling, and facilitation of cell cycle progression, as well as recruitment of transcriptional regulators such as HDGF/c-Jun [43-45]. Other genes in the model further highlight key aspects of tumor biology. SLC16A3 has been associated with immune suppression and resistance to therapy [46, 47], while KPNA2 plays a vital role in telomere maintenance and tumor progression [48]. CCNB1 contributes to cell cycle regulation and lipid metabolism through the FOXM1-CCNB1 axis [49, 50], and TPX2 has been linked to tumor proliferation, immune evasion, and resistance to sorafenib [51, 52]. The CDK4/6-DUB3 axis stabilizes YAP1, a key effector of the Hippo pathway that promotes tumor growth [53]. ETV4 is involved in hepatic inflammation and angiogenesis [54, 55], and YWHAZ participates in multiple signaling pathways that are dysregulated in HCC [56, 57]. Compared with traditional classification systems based primarily on proliferative status or mutation profiles, this gene set captures additional dimensions of tumor biology, including metabolic reprogramming, ferroptosis sensitivity, and immune regulation. This provides a more integrated view of HCC heterogeneity and may offer new opportunities for both diagnosis and treatment.

Several gene expression-based models for HCC have been reported. For instance, a 10-gene LASSO-based model validated across TCGA and ICGC cohorts showed enrichment of cell cycle-related pathways in high-risk patients and metabolic pathways in low-risk patients [58]. Tang et al. [59] developed a four-gene model focused on ferroptosis and iron metabolism, in which the high-risk group demonstrated greater tumor mutational burden, increased immune infiltration, and elevated immune checkpoint expression. Similarly, Wu et al. [60] proposed a five-gene ferroptosis-related model that showed strong predictive performance across multiple datasets, with high-risk patients displaying features of immune escape and higher TIDE scores. Most existing models, however, are derived from broad or relatively heterogeneous gene sets. In contrast, our model specifically integrates ferroptosis and metabolic reprogramming, two closely linked processes in tumor biology. Based on this framework, HCC can be broadly divided into subtypes characterized by features such as ferroptosis tolerance with high glycolytic activity and limited immune infiltration, versus a more ferroptosis-sensitive state with lower metabolic activity and greater immune engagement. From a biological perspective, this model emphasizes pathways related to nutrient transport, glycolysis and the pentose phosphate pathway, ferroptosis regulation, and epigenetic control.

Our findings suggest that interconnected metabolic, hypoxic, and lipid signaling pathways may collectively determine whether tumor cells adopt a ferroptosis-sensitive or ferroptosis-tolerant state. This aligns with prior reports showing that these pathways play central roles in HCC development and progression. In particular, we identified four key signaling axes: the HIF-1-driven glycolysis and central carbon metabolism pathway, the PPAR-regulated fatty acid and cholesterol metabolism pathway, the bile acid-cholesterol-inflammation network, and the ROS/redox-ferroptosis regulatory system. These pathways appear to coordinate ferroptosis resistance, metabolic reprogramming, and tumor growth. For example, activation of HIF-1 promotes glycolysis through upregulation of glucose transporters and key glycolytic enzymes such as HK2 [61]. The SOX4-ChREBP-SCD1 axis drives the synthesis of monounsaturated fatty acids, which can suppress ferroptosis [62], while the SLC7A11-GSH-GPX4 axis plays a central role in controlling cellular sensitivity to ferroptosis [63, 64]. In addition, these pathways can reshape the tumor microenvironment by increasing the production of metabolites such as lactate and bile acids, as well as inflammatory mediators, thereby promoting immune modulation and angiogenesis [62, 65, 66]. Importantly, alterations in these pathways are closely associated with clinical features of HCC. Enhanced activation of immune-metabolic pathways is often associated with an immunosuppressive microenvironment, characterized by enrichment of regulatory T cell (Treg) and myeloid-derived suppressor cell (MDSC), as well as upregulated expression of immune checkpoints (e.g., PD-1/PD-L1). These changes are also linked to resistance to commonly used therapies, including sorafenib and lenvatinib, through mechanisms involving metabolic rewiring and fatty acid metabolism [67-69]. Taken together, these findings provide a rationale for exploring stratified therapeutic strategies in HCC patients defined by the 12-gene ferroptosis-metabolic reprogramming signature. High-risk tumors exhibited coordinated upregulation of glycolytic and HIF-1 signaling pathways together with an immunosuppressive immune microenvironment, characterized by elevated T-cell exhaustion, Treg, and M2 TAM signatures and reduced cytotoxicity (Figure 13 and Table S10). These features suggest that high-risk patients may be particularly vulnerable to therapeutic approaches that simultaneously target metabolic vulnerabilities and relieve immune suppression. Biologically plausible strategies for this subgroup include combining immune checkpoint inhibitors with agents that modulate glycolysis or redox balance (such as glycolysis or pentose-phosphate pathway inhibitors) or that enhance ferroptosis sensitivity through lipid metabolism regulation (such as SCD1 inhibitors). In contrast, low-risk patients, who displayed lower glycolytic activity and less pronounced immunosuppressive features, may derive sufficient benefit from standard immune checkpoint inhibitor-based regimens without the added complexity of multi-drug combinations. These hypotheses are grounded in the pathway enrichment and functional immune profiling results of the present study; however, they remain to be tested in prospective clinical trials and functional experiments.

Our immune-related analyses further highlight the clinical relevance of HCC subtype classification. Immune checkpoint activity and immune cell infiltration are closely linked and play a central role in shaping treatment response. ICGs (e.g., PD-1/PD-L1), are well-established therapeutic targets, and their inhibition can restore anti-tumor immunity [70, 71]. Moreover, the immune infiltration profile of the tumor microenvironment is an important determinant of prognosis and response to immunotherapy [72, 73]. Consistent with these concepts, we observed significant differences in ICG expression between the two HCC subtypes. Key checkpoint molecules, including LAG3, CTLA4, and PD-L1 (CD274), were more highly expressed in Subtype 1, which also showed higher TIDE scores and lower effector cell activity. This pattern suggests a more immunosuppressive phenotype, potentially associated with increased immune evasion and reduced responsiveness to immunotherapy. In contrast, Subtype 2 displayed a less immunosuppressive profile, indicating a greater likelihood of benefiting from immune checkpoint inhibitors (ICIs). These findings suggest that integrating our prognostic model with immune-related indicators, such as ICG expression and TIDE scores, could help guide individualized treatment strategies. In practice, this approach may support decisions on whether to use ICIs alone, combine ICIs with tyrosine kinase inhibitors, incorporate metabolic-targeted therapies, or avoid immunotherapy in patients unlikely to benefit.

Beyond cell-type abundance, the functional immune profile of the tumour microenvironment reinforces the immune-evasive phenotype of high-risk patients: their elevated exhaustion/Treg/MDSC signatures and reduced effector activity converge on a "cold" or immunosuppressive tumour state that may limit responses to single-agent immune checkpoint inhibitors, while the same patients could potentially benefit from combinations that target the suppressive compartment (e.g., MDSC- or Treg-depleting agents) alongside ICIs. The single-cell mapping further suggests that the cell-cycle genes in our model partly reflect tumour-infiltrating proliferating lymphocytes rather than tumour-cell proliferation alone, while the redox gene G6PD shows expression in NK cells linked to interferon response — an observation that warrants experimental follow-up in HCC immune biology.

This study has several limitations. Firstly, our findings are based entirely on in silico analyses and have not yet been validated in experimental or clinical settings. Although we included relatively large cohorts from the TCGA and GEO databases, residual batch effects and heterogeneity in data processing may limit the broader applicability of the results. In addition, we did not account for other factors, such as environmental exposures, lifestyle variables, or additional genetic alterations, that may influence HCC progression and could confound the observed associations. These limitations highlight the need for further validation in independent cohorts and in clinical tissue samples. Future studies should also incorporate functional experiments to better define the functional role of the identified genes. Additionally, the permissive |log2FC| > 0.25 threshold used for initial DEG selection may have included weakly-changing transcripts; a sensitivity analysis (Table S6, Figure S6) confirmed that the headline 12-gene signature is recoverable at this cutoff but that the candidate gene pool shrinks substantially at |log2FC| > 1. Direct comparison against BCLC staging in TCGA-LIHC was not possible because BCLC stage is not provided in the TCGA clinical files; the head-to-head comparison was therefore performed against AJCC TNM in TCGA and against BCLC in GSE76427. Finally, several important limitations of the immune and single-cell analyses should be noted. The functional immune signature scoring (Figure 13) and single-cell expression mapping (Figure 14) are entirely in silico. GSE140228 is a CD45+-enriched dataset that does not include malignant cells and comprises only four donors; therefore, the cellular localization findings require validation in larger, unsorted HCC cohorts with paired tumour and immune compartments. Additionally, the head-to-head comparison against BCLC staging in GSE76427 was limited by the small number of overall survival events (n = 23), resulting in wide confidence intervals for NRI and IDI estimates; these results should be interpreted with caution and ideally validated in larger external cohorts.

Conclusion

In summary, we identified a set of key genes associated with ferroptosis and metabolic reprogramming in HCC and established a prognostic model with good predictive performance. This model provides new insight into the molecular basis of HCC and highlights potential biomarkers for patient stratification and personalized treatment. Further work integrating clinical validation with mechanistic studies will be important to confirm these findings and to examine the therapeutic potential of targeting these pathways.

Supplementary Material

Supplementary figures and table legends.

Supplementary tables.

Abbreviations

AUC: Area under the curve; BH: Benjamini-Hochberg; CDF: Cumulative distribution function; CI: Confidence interval; CDK4: Cyclin-dependent kinase 4; CTLA4: Cytotoxic T-lymphocyte associated protein 4; DCA: Decision curve analysis; DEGs: Differentially expressed genes; DSS: Disease-specific survival; CD274: Cluster of differentiation 274; FDR: False discovery rate; FOXM1: Forkhead box M1; FPMRRDEGs: Ferroptosis- and metabolic reprogramming-related differentially expressed genes; FPMRRGs: Ferroptosis- and metabolic reprogramming-related genes; GEO: Gene Expression Omnibus; GO: Gene ontology; GPX4: Glutathione peroxidase 4; GSH: Glutathione; GSVA: Gene set variation analysis; HCC: Hepatocellular carcinoma; HIF-1: Hypoxia-inducible factor-1; Hippo: Hippo signaling pathway; HR: Hazard ratio; ICIs: Immune checkpoint inhibitors; ICGs: Immune checkpoint genes; IDI: Integrated discrimination improvement; IPS: Immunophenoscore; IQR: Interquartile range; KEGG: Kyoto Encyclopedia of Genes and Genomes; KM: Kaplan-Meier; LAG3: Lymphocyte-activation gene 3; LASSO: Least absolute shrinkage and selection operator; LIHC: Liver hepatocellular carcinoma; MDSC: Myeloid-derived suppressor cell; MSI: Microsatellite instability; MSigDB: Molecular Signatures Database; NADPH: Nicotinamide adenine dinucleotide phosphate; NES: Normalized enrichment score; NRI: Net reclassification improvement; OS: Overall survival; PD-1: Programmed cell death protein 1; PD-L1: Programmed cell death ligand 1; PFI: Progression-free interval; PPAR: Peroxisome proliferator-activated receptor; ROC: Receiver operating characteristic; ROS: Reactive oxygen species; SD: Standard deviation; SLC: Solute carrier; ssGSEA: Single-sample gene set enrichment analysis; TCGA: The Cancer Genome Atlas; TIDE: Tumor immune dysfunction and exclusion; TMB: Tumor mutational burden; Treg: Regulatory T cell; t-SNE: t-distributed stochastic neighbor embedding; YAP1: Yes-associated protein 1.

Data Availability Statement

All data used in this research were derived from publicly available databases. Details of the bioinformatics analysis are available from the corresponding author upon reasonable request.

Ethics Committee Approval and Patient Consent

This research was based entirely on publicly available datasets and did not involve the acquisition of new human/animal data. Ethical approval was therefore not required.

Competing Interests

The authors have declared that no competing interest exists.

References

1. Llovet JM, Kelley RK, Villanueva A, Singal AG, Pikarsky E, Roayaie S. et al. Hepatocellular carcinoma. Nat Rev Dis Primers. 2021;7:6

2. Yu Y, Jiang Y, Yang Y, Glandorff C, Fang W, Sun M. Synergistic strategies in systemic therapy for advanced hepatocellular carcinoma. Glob Chall. 2026;10:e00247

3. Carter JK, Cameron DC, Villanueva A. Clinical applications of immunogenomics in hepatocellular carcinoma. Clin Mol Hepatol. 2026

4. Lee JS, Choi HW, Kim JS, Lee TY, Yoon YC. Update on resection strategies for hepatocellular carcinoma: A narrative review. Cancers (Basel). 2024 16

5. Chen X, Liu HP, Li M, Qiao L. Advances in non-surgical management of primary liver cancer. World J Gastroenterol. 2014;20:16630-8

6. Liu S, Li Z, Lan S, Hao H, Baz AA, Yan X. et al. The dual roles of activating transcription factor 3 (ATF3) in inflammation, apoptosis, ferroptosis, and pathogen infection responses. Int J Mol Sci. 2024;25:824

7. Hadian K, Stockwell BR. SnapShot: Ferroptosis. Cell. 2020;181:1188-e1

8. Tu S, Zou Y, Yang M, Zhou X, Zheng X, Jiang Y. et al. Ferroptosis in hepatocellular carcinoma: Mechanisms and therapeutic implications. Biomed Pharmacother. 2025;182:117769

9. Tang D, Kroemer G, Kang R. Ferroptosis in hepatocellular carcinoma: From bench to bedside. Hepatology. 2024;80:721-39

10. Mo Y, Zou Z, Chen E. Targeting ferroptosis in hepatocellular carcinoma. Hepatol Int. 2024;18:32-49

11. Faubert B, Solmonson A, DeBerardinis RJ. Metabolic reprogramming and cancer progression. Science. 2020;368:eaaw5473

12. Raggi C, Taddei ML, Rae C, Braconi C, Marra F. Metabolic reprogramming in cholangiocarcinoma. J Hepatol. 2022;77:849-64

13. Dai P, Feng J, Dong Y, Zhang S, Cao J, Cui X. et al. Metabolic reprogramming in hepatocellular carcinoma: an integrated omics study of lipid pathways and their diagnostic potential. J Transl Med. 2025;23:644

14. Park S, Hall MN. Metabolic reprogramming in hepatocellular carcinoma: Mechanisms and therapeutic implications. Exp Mol Med. 2025;57:515-23

15. Wang H, Lun Y, Xu D, Jiang H, Yan Y, Yang X. Research progress and therapeutic strategies in hepatocellular carcinoma metabolic reprogramming. J Adv Res. 2025

16. An L, Li Z. Molecular network of metabolic reprogramming and precision diagnosis and treatment of hepatocellular carcinoma. Biomark Res. 2025;13:124

17. Gao B, Lu Y, Lai X, Xu X, Gou S, Yang Z. et al. Metabolic reprogramming in hepatocellular carcinoma: Mechanisms of immune evasion and therapeutic implications. Front Immunol. 2025;16:1592837

18. Hu N, Li H, Tao C, Xiao T, Rong W. The role of metabolic reprogramming in the tumor immune microenvironment: mechanisms and opportunities for immunotherapy in hepatocellular carcinoma. Int J Mol Sci. 2024;25:5584

19. Bernard NJ. Immunosuppressive ferroptosis. Nat Immunol. 2022;23:1653

20. Sun Y, Xue Z, Huang T, Che X, Wu G. Lipid metabolism in ferroptosis and ferroptosis-based cancer therapy. Front Oncol. 2022;12:941618

21. Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D. et al. TCGAbiolinks: An R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. 2016;44:e71

22. Davis S, Meltzer PS. GEOquery: A bridge between the gene expression omnibus (GEO) and bioconductor. Bioinformatics. 2007;23:1846-7

23. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W. et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43:e47

24. Cai Y, Deng L, Yao J. Analysis and identification of ferroptosis-related diagnostic markers in rheumatoid arthritis. Ann Med. 2024;56:2397572

25. Huang Y, Li N, Jiang J, Pei Y, Gao S, Qian Y. et al. Metabolic reprogramming-related gene classifier distinguishes malignant from the benign pulmonary nodules. Heliyon. 2024;10:e37214

26. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550

27. Lock EF, Dunson DB. Bayesian consensus clustering. Bioinformatics. 2013;29:2610-6

28. Liberzon A, Subramanian A, Pinchback R, Thorvaldsdóttir H, Tamayo P, Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinformatics. 2011;27:1739-40

29. Xu D, Liu X, Wang Y, Zhou K, Wu J, Chen JC. et al. Identification of immune subtypes and prognosis of hepatocellular carcinoma based on immune checkpoint gene expression profile. Biomed Pharmacother. 2020;126:109903

30. Gao J, Aksoy BA, Dogrusoz U, Dresdner G, Gross B, Sumer SO. et al. Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci Signal. 2013;6:pl1

31. Foroutan M, Bhuva DD, Lyu R, Horan K, Cursons J, Davis MJ. Single sample scoring of molecular phenotypes. BMC Bioinformatics. 2018;19:404

32. Zhang Q, He Y, Luo N, Patel SJ, Han Y, Gao R. et al. Landscape and Dynamics of Single Immune Cells in Hepatocellular Carcinoma. Cell. 2019 179

33. Han Y, Wang Y, Dong X, Sun D, Liu Z, Yue J. et al. TISCH2: expanded datasets and new tools for single-cell transcriptome analyses of the tumor microenvironment. Nucleic Acids Res. 2023;51:D1425-D31

34. Sun LL, He HY, Li W, Jin WL, Wei YJ. The solute carrier transporters (SLCs) family in nutrient metabolism and ferroptosis. Biomark Res. 2024;12:94

35. Kobayashi H, Yoshimoto C, Matsubara S, Shigetomi H, Imanaka S. A comprehensive overview of recent developments on the mechanisms and pathways of ferroptosis in cancer: The potential implications for therapeutic strategies in ovarian cancer. Cancer Drug Resist. 2023;6:547-66

36. Chen X, Li J, Kang R, Klionsky DJ, Tang D. Ferroptosis: Machinery and regulation. Autophagy. 2021;17:2054-81

37. Kajarabille N, Latunde-Dada GO. Programmed cell-death by ferroptosis: Antioxidants as mitigators. Int J Mol Sci. 2019;20:4968

38. Jin M, Zhang J, Sun Y, Liu G, Wei X. ANXA5: Related mechanisms of osteogenesis and additional biological functions. Front Cell Dev Biol. 2025;13:1553683

39. Oflaz FE, Bondarenko AI, Trenker M, Waldeck-Weiermair M, Gottschalk B, Bernhart E. et al. Annexin A5 controls VDAC1-dependent mitochondrial Ca²⁺ homeostasis and determines cellular susceptibility to apoptosis. EMBO J. 2025;44:3413-47

40. Zhao L, Cao H, Yuan Y, Liao C, Huang D, Li X. et al. Annexin A5 knockdown inhibits cardiomyocyte apoptosis and alleviates cardiac hypertrophy via activating the PI3K/AKT/Bcl-2 signaling pathway. Sci Rep. 2024;14:31915

41. Zheng J, Wang Y, Zhou Y, Li Z, Yang L, Gao J. et al. Augmentation of hepatocellular carcinoma malignancy by annexin A5 through modulation of invasion and angiogenesis. Scand J Gastroenterol. 2024;59:939-53

42. Wang W, Liu D, Yao J, Yuan Z, Yan L, Cao B. ANXA5: A key regulator of immune cell infiltration in hepatocellular carcinoma. Med Sci Monit. 2024;30:e943523

43. Zhang SL, Zhang SJ, Li L, Zhang YW, Wang Z, Wang L. et al. NAP1L1 regulates BIRC2 ubiquitination modification via E3 ubiquitin ligase UBR4 and hence determines hepatocellular carcinoma progression. Cell Death Discov. 2024;10:154

44. Shen B, Zhu W, Liu X, Jiang J. NAP1L1 functions as a novel prognostic biomarker associated with macrophages and promotes tumor progression by influencing the Wnt/β-catenin pathway in hepatocellular carcinoma. Front Genet. 2022;13:876253

45. Zhang YW, Chen Q, Li B, Li HY, Zhao XK, Xiao YY. et al. NAP1L1 functions as a tumor promoter via recruiting hepatoma-derived growth Factor/c-Jun signal in hepatocellular carcinoma. Front Cell Dev Biol. 2021;9:659680

46. Peng X, He Z, Yuan D, Li J, Zhu Y, Zhang X. et al. The role of solute carrier family 16 member 3 protein in hepatocellular carcinoma and sorafenib resistance. Int J Biol Macromol. 2025;330:148223

47. Kang Q, Yin X, Wu Z, Zheng A, Feng L, Ma X. et al. Integrated single-cell and spatial transcriptome reveal metabolic gene SLC16A3 as a key regulator of immune suppression in hepatocellular carcinoma. J Cell Mol Med. 2024;28:e70272

48. Ding K, Liu L, Yong W, Sun B, Zhang W. Bioinformatics analysis and experimental studies reveal KPNA2 as a novel biomarker of hepatocellular carcinoma progression and telomere maintenance. Eur J Med Res. 2025;30:628

49. Jin X, He X, Huang R, Liu Q, Wang L, Bai X. et al. SNRPB/CCNB1 axis promotes hepatocellular carcinoma progression and cisplatin resistance through enhancing lipid metabolism reprogramming. J Exp Clin Cancer Res. 2025;44:211

50. Rong MH, Li JD, Zhong LY, Huang YZ, Chen J, Xie LY. et al. CCNB1 promotes the development of hepatocellular carcinoma by mediating DNA replication in the cell cycle. Exp Biol Med (Maywood). 2022;247:395-408

51. Peng H, Wei YP, Liu XB, Wang Y, Yuan JY. Integrated multi-omics analysis identifies TPX2 as a potential prognostic and immunological biomarker in hepatocellular carcinoma. iLIVER. 2025;4:100184

52. Wang H, Chu F, Zhang XF, Zhang P, Li LX, Zhuang YL. et al. TPX2 enhances the transcription factor activation of PXR and enhances the resistance of hepatocellular carcinoma cells to antitumor drugs. Cell Death Dis. 2023;14:64

53. Huang L, Yuan W, Li X, Liu Y, Wan R, Ma X. et al. CDK4/6-mediated phosphorylation of DUB3 promotes YAP1 stability and hepatocellular carcinoma progression. Cell Death Discov. 2025;11:212

54. Su H, Shu S, Tang W, Zheng C, Zhao L, Fan H. ETV4 facilitates angiogenesis in hepatocellular carcinoma by upregulating MMP14 expression. Biochem Biophys Res Commun. 2023;684:149137

55. Qi D, Lu M, Xu P, Yao X, Chen Y, Gan L. et al. Transcription factor ETV4 promotes the development of hepatocellular carcinoma by driving hepatic TNF-α signaling. Cancer Commun (Lond). 2023;43:1354-72

56. Zheng W, Shen GL, Xu KY, Yin QQ, Hui TC, Zhou ZW. et al. Lnc524369 promotes hepatocellular carcinoma progression and predicts poor survival by activating YWHAZ-RAF1 signaling. World J Gastrointest Oncol. 2022;14:253-64

57. Xie J, Guo T, Zhong Z, Wang N, Liang Y, Zeng W. et al. ITGB1 drives hepatocellular carcinoma progression by modulating cell cycle process through PXN/YWHAZ/AKT pathways. Front Cell Dev Biol. 2021;9:711149

58. Zhou T, Cai Z, Ma N, Xie W, Gao C, Huang M. et al. A novel ten-gene signature predicting prognosis in hepatocellular carcinoma. Front Cell Dev Biol. 2020;8:629

59. Tang B, Zhu J, Li J, Fan K, Gao Y, Cheng S. et al. The ferroptosis and iron-metabolism signature robustly predicts clinical diagnosis, prognosis and immune microenvironment for hepatocellular carcinoma. Cell Commun Signal. 2020;18:174

60. Wu C, Chen X, Zhang Y, Du Y, Xu J, Peng Y. et al. A novel ferroptosis-related signature for predicting prognosis, immune characteristics, and treatment prediction in hepatocellular carcinoma. PLoS One. 2025;20:e0322158

61. Talubo NDD, Dela Cruz EWB, Fowler P, Tsai PW, Tayo LL. QSAR-based drug repurposing and RNA-Seq metabolic networks highlight treatment opportunities for hepatocellular carcinoma through pyrimidine starvation. Cancers (Basel). 2025;17:903

62. Zhang F, Wu Z, Xiang Y, He Q, Li W, Yang K. et al. SOX4 reprograms fatty acid metabolism through the CHREBP to inhibit ferroptosis in hepatocellular carcinoma. Cell Death Discov. 2025;11:246

63. Jiang Y, Yu Y, Pan Z, Glandorff C, Sun M. Ferroptosis: A new hunter of hepatocellular carcinoma. Cell Death Discov. 2024;10:136

64. Ajoolabady A, Tang D, Kroemer G, Ren J. Ferroptosis in hepatocellular carcinoma: Mechanisms and targeted therapy. Br J Cancer. 2023;128:190-205

65. Fang W, Liu J, Zhang F, Pang C, Li X. A novel cholesterol metabolism-related ferroptosis pathway in hepatocellular carcinoma. Discov Oncol. 2024;15:7

66. Zhu X, Sha X, Zang Y, Ren Q, Zhang S, Ma D. et al. Current progress of ferroptosis study in hepatocellular carcinoma. Int J Biol Sci. 2024;20:3621-37

67. Zhang X, Liu X, Zhu K, Zhang X, Li N, Sun T. et al. CD5L-associated gene analyses highlight the dysregulations, prognostic effects, immune associations, and drug-sensitivity predicative potentials of LCAT and CDC20 in hepatocellular carcinoma. Cancer Cell Int. 2022;22:393

68. Nishida N. Role of oncogenic pathways on the cancer immunosuppressive microenvironment and its clinical implications in hepatocellular carcinoma. Cancers (Basel). 2021;13:3666

69. Guo Y, Yang J, Ren K, Tian X, Gao H, Tian X. et al. The heterogeneity of immune cell infiltration landscape and its immunotherapeutic implications in hepatocellular carcinoma. Front Immunol. 2022;13:861525

70. Li B, Chan HL, Chen P. Immune checkpoint inhibitors: Basics and challenges. Curr Med Chem. 2019;26:3009-25

71. Dyck L, Mills KHG. Immune checkpoints and their inhibition in cancer and infectious diseases. Eur J Immunol. 2017;47:765-79

72. Li M, Yang Y, Xiong L, Jiang P, Wang J, Li C. Metabolism, metabolites, and macrophages in cancer. J Hematol Oncol. 2023;16:80

73. Toledo B, Zhu Chen L, Paniagua-Sancho M, Marchal JA, Perán M, Giovannetti E. Deciphering the performance of macrophages in tumour microenvironment: A call for precision immunotherapy. J Hematol Oncol. 2024;17:44

Author contact

Corresponding authors: Hao Ling (Email: hao.lingde) and Yanzhu Hu (Email: yanzhu.hude)

Citation styles

©2026 Ivyspring International Publisher. All rights reserved.
Open access content is subject to applicable licences. By accessing or using this platform, you agree to our Terms of Use.