Your new experience awaits. Try the new design now and help us make it even better

ORIGINAL RESEARCH article

Front. Immunol., 02 February 2026

Sec. Cancer Immunity and Immunotherapy

Volume 17 - 2026 | https://doi.org/10.3389/fimmu.2026.1739293

This article is part of the Research TopicImmuno-metabolic Approaches for the Treatment of Hepatobiliary and Pancreatic TumorsView all 16 articles

Single-cell RNA sequencing and spatial transcriptomic analysis reveal a distinct population of G6PD+ cells with aberrant bile acid metabolism in hepatocellular carcinoma

Xing JiangXing Jiang1Haiyan QuanHaiyan Quan2Ting YinTing Yin1Hailun YaoHailun Yao2Yajun LiYajun Li2Bin PengBin Peng1Xinye YuanXinye Yuan1Weiguang ZengWeiguang Zeng1Honghui Chen*Honghui Chen1*Rong Li*Rong Li1*
  • 1Hunan Provincial Key Laboratory of Basic and Clinical Pharmacological Research of Gastrointestinal Cancer, Department of Pharmacy, Institute of Pharmacy and Pharmacology, the Second Affiliated Hospital, University of South China, Hengyang, Hunan, China
  • 2School of Pharmaceutical Technology, Hunan Polytechnic of Environment and Biology, Hengyang, Hunan, China

Background: Metabolic reprogramming is a hallmark of hepatocellular carcinoma (HCC). Among various metabolic pathways, bile acids act not only as crucial metabolites but also as key signaling molecules that regulate diverse physiological and pathological processes in the liver. However, the biological functions and clinical implications of bile acid metabolism in HCC progression remain largely unclear.

Methods: Single-cell transcriptomic data from 67 patients with HCC were integrated to construct a bile acid metabolism scoring system. Pseudotime trajectory analysis was employed to characterize the differentiation patterns of cells exhibiting abnormal bile acid metabolism. Spatial transcriptomics was used to explore their spatial distribution features. Furthermore, machine learning algorithms were applied to analyze transcriptomic data from HCC cohorts to develop a prognostic prediction model. The findings were complemented by immune infiltration analysis, molecular characterization, and drug sensitivity prediction using CellMiner, followed by molecular docking validation.

Results: G6PD+ malignant tumor cells with high bile acid metabolism scores exhibited enhanced bile acid metabolic activity, accompanied by activation of macrophages and endothelial cells. These cells were predominantly localized at the tumor boundary region. A prognostic prediction model based on G6PD+ expression successfully identified a high-risk subgroup with significantly poorer outcomes. In vitro experiments demonstrated that knockdown or overexpression of G6PD markedly affected the proliferative, migratory, and invasive capacities of HCC cells.

Conclusion: This study reveals that bile acid metabolism promotes HCC progression by facilitating vascular network formation and establishing an immunosuppressive tumor microenvironment. The bile acid metabolism scoring system may serve as a novel prognostic biomarker and provide a theoretical foundation for developing precision therapeutic strategies in HCC.

1 Introduction

Hepatocellular carcinoma (HCC) is a highly malignant tumor with substantial global incidence and mortality, typically arising in the setting of liver cirrhosis, chronic hepatitis virus infection, or metabolic dysfunction–associated fatty liver disease (MAFLD) (1). To sustain their rapid proliferation, HCC cells undergo profound metabolic reprogramming, which has become one of the defining hallmarks of their malignant phenotype (2). Among these alterations, a striking feature is their strong dependence on aerobic glycolysis — even under oxygen-rich conditions, tumor cells preferentially obtain energy via glycolysis, a phenomenon known as the Warburg effect. This process leads to excessive lactate production and secretion, thereby creating an acidic tumor microenvironment that promotes immune evasion and tumor growth (3, 4). Compared with normal hepatocytes, HCC cells exhibit markedly increased expression of glycolysis-related genes such as GLUT1, HK2, and PKM2, which enhances glucose uptake and utilization, accelerates tumor progression, and correlates with poor patient outcomes (5, 6). Beyond glycolysis, metabolic abnormalities in HCC also encompass the pentose phosphate pathway, fatty acid metabolism, and amino acid metabolism (2, 7). Consequently, targeting metabolic reprogramming has emerged as a critical strategy in the development of novel HCC therapies.

The liver serves as the central organ for bile acid synthesis and metabolism. In addition to their well-known role in digestion, bile acids act as pivotal signaling molecules that regulate metabolism, inflammation, and immune responses, thereby exerting widespread influence on hepatic physiological and pathological processes (8, 9). Disruption of bile acid homeostasis—whether in synthesis, transport, or secretion—can lead to intrahepatic cholestasis, persistent hepatocellular injury, and chronic inflammation, ultimately predisposing the liver to malignant transformation (10, 11). In addition, cholesterol homeostasis mediated by bile acid metabolism can also promote tumorigenesis (12).

Clinically, this connection is supported by the observation that up to 40% of HCC patients present with cholestatic jaundice at diagnosis, suggesting a strong association between bile acid dysregulation, tumor progression, and adverse prognosis (13). Mechanistically, loss of farnesoid X receptor (FXR) disrupts the feedback regulation of bile acid synthesis, leading to excessive accumulation, chronic hepatic injury, and spontaneous HCC development (14, 15). Conversely, FXR activation restores bile acid homeostasis by repressing the rate-limiting enzyme CYP7A1 and inducing the bile salt export pump BSEP, thereby alleviating hepatic inflammation and injury (16). Collectively, bile acid dysregulation represents a fundamental pathogenic driver of HCC. However, whether there exist distinct cellular subtypes with abnormal bile acid metabolism in HCC tissues, and whether these can serve as independent prognostic indicators, remains to be further investigated.

In this study, we integrated bulk, single-cell, and spatial transcriptomic datasets to develop an innovative bile acid metabolism scoring system and comprehensively characterize bile acid metabolism in HCC at the molecular and spatial levels. We identified a unique subset of G6PD+ malignant cells exhibiting aberrant bile acid metabolism, which are closely associated with angiogenesis and the establishment of an immunosuppressive microenvironment. By fulfilling these findings, we seek to establish a mechanistic link between bile acid metabolic reprogramming and HCC progression, providing a translatable framework for targeting this pathway.

2 Materials and methods

2.1 Data sources

We systematically explored multiple biological dimensions of HCC by integrating data from four complementary public databases. Single-cell transcriptomic datasets were obtained from GSE149614 (n=10 patients) and GSE156625 (n=57 patients), which reveal the cellular heterogeneity of HCC (17, 18). Spatial transcriptomic data (10x Genomics, 226334_D5_2) provided a tissue-level spatial context for molecular expression patterns. In addition, bulk transcriptomic data from TCGA-LIHC, ICGC-LIRI-JP, GSE14520, and GSE43619 cohorts were used to validate the robustness of our findings across diverse clinical populations. The bile acid metabolism-related gene set was curated from the GeneCards database (Supplementary Table 1), encompassing comprehensive pathway components to ensure biological relevance. This study involved the analysis of existing, de-identified human genomic and transcriptomic data obtained from public repositories. All data were accessed and used in accordance with the respective database policies. Therefore, no separate ethics approval was required for this secondary analysis, as confirmed by our institutional review board.

2.2 Single-cell RNA sequencing analysis

Single-cell RNA sequencing (scRNA-seq) data were processed using the Seurat package (v4.3.0) in the R environment. To ensure high data quality, genes expressed in fewer than three cells were excluded, and cells with fewer than 50 detected genes or with mitochondrial gene content exceeding 5% were filtered out. Data normalization was performed using the SCTransform method, followed by IntegrateData to correct batch effects and construct a unified analytical object. Dimensionality reduction was achieved via principal component analysis (PCA), and cells were clustered using the FindNeighbors and FindClusters functions to identify subpopulations with similar transcriptional profiles. Cluster visualization was accomplished through Uniform Manifold Approximation and Projection (UMAP). Cell types were annotated based on cluster-specific marker genes identified by FindAllMarkers, with each subpopulation defined by its unique transcriptional signature (19). The irGSEA function was applied to score bile acid metabolism activity, enabling evaluation of functional enrichment across distinct cellular subtypes. This comprehensive analytical workflow ensured robust identification and functional characterization of the cellular components within the tumor microenvironment.

2.3 Cell–cell communication analysis

To systematically dissect the complex intercellular communication network within the tumor microenvironment, we employed CellChat, an integrative computational framework containing a manually curated ligand–receptor interaction database and downstream signaling pathway annotations (20). This approach enabled the identification and quantification of significantly enriched ligand–receptor interactions based on co-expression patterns across cell types. Through this rigorous framework, we uncovered previously unrecognized signaling axes that coordinate intercellular communication in the HCC microenvironment, providing mechanistic insights into how specific cellular subpopulations cooperate to promote tumor progression and immune evasion.

2.4 Single-cell trajectory construction

To capture the dynamic evolution of cellular states during tumor progression, we performed pseudotime trajectory inference on epithelial cells using Monocle (v2.30.1) (21). The analysis focused on epithelial subsets previously integrated into the R environment. A new cell dataset was created using newCellDataSet, and genes dynamically expressed along pseudotime were identified. Dimensionality reduction was achieved using the DDRTree algorithm implemented in reduceDimension, mapping cells onto a minimum spanning tree while preserving developmental relationships inferred from transcriptomic similarity. The resulting trajectory was visualized with plot_cell_trajectory, illustrating gene expression dynamics along pseudotime. Pseudotime-dependent genes were further characterized, and their expression trends were visualized using plot_pseudotime_heatmap, thereby depicting the molecular reprogramming events accompanying epithelial phenotypic transitions in HCC.

2.5 Spatial transcriptomics analysis

Spatial transcriptomic data were processed and visualized using the Seurat R package. Data normalization was performed with SCTransform, followed by unsupervised clustering to identify spatially distinct transcriptional domains. Cell-type annotation was validated through histopathological inspection of H&E-stained sections. To achieve accurate cell identity mapping within spatial coordinates, we applied a multimodal computational framework combining Robust Cell Type Decomposition (RCTD) and Multimodal Intersection Analysis (MIA), integrating reference scRNA-seq datasets with spatial transcriptomic profiles (22). This multimodal integration preserved the native spatial context of cell–cell interactions and revealed previously unrecognized region-specific cellular distributions within the HCC microenvironment.

2.6 Machine learning algorithms

An ensemble machine learning framework was employed to construct a comprehensive prognostic model for HCC. Ten algorithms were systematically evaluated, including random survival forests, elastic net regression, Lasso regression, stepwise Cox regression, ridge regression, CoxBoost, partial least squares regression for Cox models, supervised principal component analysis, gradient boosting machines, and survival support vector machines (23). The model with the highest average concordance index (C-index) across all comparisons was selected as the final predictive model. Supervised Principal Component Analysis (SuperPC) was selected as the final modeling algorithm based on its superior average Harrell’s Concordance Index (C-index). The model’s key hyperparameters specifically, the correlation cutoff (threshold) for feature gene selection and the number of principal components were optimized via a grid search combined with 10-fold cross-validation on the training set, with the combination maximizing the cross-validated C-index chosen for the final model refit. To ensure generalizability and prevent overfitting, we relied on SuperPC’s inherent dimensionality reduction regularization, conducted all tuning via cross-validation to avoid bias, and most critically, validated the fixed model on an independent testing cohort where it achieved a C-index of 0.71, demonstrating consistent performance. Model performance was rigorously validated using time-dependent receiver operating characteristic (ROC) curves, and the area under the curve (AUC) was calculated via the timeROC package. Multivariate Cox regression analysis confirmed that the risk score derived from this model was an independent and statistically significant prognostic factor for patient outcomes.

2.7 Immune infiltration and survival analysis

To comprehensively characterize the immune landscape of HCC, we conducted gene set variation analysis (GSVA) in combination with deconvolution approaches implemented in the IOBR platform (24, 25). Enrichment patterns of various immune cell subsets were systematically assessed, and expression differences in immune checkpoint molecules and functional immune signatures were compared between predefined subgroups. The clinical relevance of the prognostic model was further evaluated by survival analyses across multiple HCC cohorts (GEO, LIHC and ICGC datasets) using the survival (v3.5-8) and Survminer (v0.4.9) packages. Patients were stratified into high- and low-risk groups based on the median risk score, and Kaplan–Meier survival curves were plotted to compare outcomes between groups.

2.8 Drug sensitivity prediction

To identify potential therapeutic compounds, we performed integrative drug sensitivity modeling using the CellMiner database. Gene expression matrices were log2-transformed and batch-corrected for normalization before integration with drug response data. The resulting computational model predicted candidate compounds with potential therapeutic relevance. Structural modeling and molecular docking were conducted using Discovery Studio 2019. Protein crystal structures were obtained from the Protein Data Bank (PDB, ID: 7UAG), and small-molecule structures from PubChem (ID: 321674-73-1). After structural optimization, molecular docking simulations were carried out with rigorously defined parameters to characterize the binding modes between candidate compounds and their target proteins.

2.9 Cell culture and functional assays

The HCC lines HCCLM3 were sourced from the Cell Bank of Shanghai Institutes for Biological Sciences, Chinese Academy of Sciences, with identity validation performed by Short Tandem Repeat (STR) analysis. HCCLM3 cells were grown in DMEM containing 10% FBS and 1% streptomycin/penicillin at 37°C in a 5% CO2 incubator. To evaluate clonogenic potential, HCCLM3 cells were plated at 6- or 24-well plates. Following 7 days of culture in 6-well plates, colonies were fixed in 4% paraformaldehyde, stained with crystal violet, and imaged under a microscope. For EdU assay, follow the operating steps of the EdU kit (C0071S; Beyotime; China). To specifically assess cell migration and minimize the contribution of proliferation in scratch assay, cells were cultured in medium containing 0.5% FBS for the duration of the assay. The wound closure was monitored by time-lapse microscopy at the same pre-marked locations at 0-, 12-, and 24-hours post-scratch. Small interfering RNA (siRNA)-mediated knockdown of G6PD was performed to validate its functional role. The target sequence for si-G6PD in our study as follows:

si-G6PD-1: 5’-GCTGAAGAAGTATGACAACATTTCAAGAGAATGTTGTCATA-3’, si-G6PD-2: 5’-GCAGTTCGTGTGGGTGATAAATTCAAGAGATTTATCACCCA-3’, si-G6PD-3: 5’-GGACAAGCCCATCATTCACAGTTCAAGAGACTGTGAATGA-3’. Cells were seeded in appropriate culture plates and grown to 60-70% confluence. Transfection was carried out using Lipofectamine 3000 according to the manufacturer’s instructions, with a final siRNA concentration of 50 nM. After 48 hours of transfection, cells were harvested for subsequent functional assays and Western blot analysis. All in vitro experiments were performed with at least three independent biological replicates. Quantitative data are presented as the mean ± standard deviation (SD).

2.10 Cell migration and invasion

Cell migration was evaluated using a wound healing assay. Briefly, cells were seeded in 6-well plates and grown to full confluence. A straight scratch was created using a sterile pipette tip, and the detached cells were gently washed away with PBS. Images of the wound were captured at the same location at 0, 12, and 24 hours under a microscope. The wound area at each time point was measured using Image J software, and the changes in area were analyzed to quantify cell migration ability.

Cell invasion was assessed using Transwell chambers (8 μm). The upper chambers were pre-coated with Matrigel matrix, and cells resuspended in serum-free medium were seeded into the upper compartment. The lower chambers were filled with medium containing 10% FBS as a chemoattractant. After 24 hours of incubation, non-invaded cells on the upper surface of the membrane were carefully removed with a cotton swab. The cells that had invaded to the lower side were fixed with 4% paraformaldehyde, stained with 0.1% crystal violet, and then counted under a microscope in randomly selected fields for statistical analysis.

2.11 Statistical analysis

All statistical analyses were performed using R software (v4.3.1) and GraphPad Prism (v8.0). Kaplan–Meier survival analysis with log-rank testing was used to compare survival differences between groups. Continuous variables were analyzed using two-tailed Student’s t-tests, and correlations were assessed using Spearman’s rank correlation. Statistical significance was defined as p < 0.05.

3 Results

3.1 Single-cell transcriptomic features of HCC patients

To comprehensively characterize the heterogeneity of the HCC tumor microenvironment, we integrated the single-cell RNA sequencing datasets GSE149614 and GSE156625. After standard QC filtering (as described in Methods), we retained a total of 94,900 high-quality cells and 34,326 genes for downstream analysis. Through unsupervised clustering, we identified nine distinct cellular populations with unique transcriptional signatures (Figures 1A, B), including T cells, B cells, NK cells, endothelial cells, fibroblasts, macrophages, plasma cells, hepatocytes, and epithelial cells. These clusters exhibited characteristic expression patterns of canonical marker genes (Figure 1C), reflecting population-specific transcriptional identities and functional states (Figure 1D). Functional enrichment analysis further revealed that these cell populations were significantly enriched in biological processes such as cell adhesion, immune activation, and metabolic pathways (Figure 1E). CellChat analysis delineated a complex network of ligand–receptor interactions within the HCC microenvironment, highlighting the extensive intercellular communication among diverse cell types (Figures 1F, G). Notably, myeloid-derived cells exhibited particularly strong signaling activity, underscoring their central role as communication hubs that orchestrate cross-talk between immune and stromal components in the tumor microenvironment.

Figure 1
(A) and (B) UMAP plots displaying cell clusters with labels such as T_cells and NK_cells. (C) Bubble plot showing feature intensity across cell types. (D) Heatmap illustrating gene expression across various cell types. (E) Network diagram of signaling pathways involved in cell interactions, color-coded by cluster. (F) and (G) Dot plots indicating the number and weight of interactions between sender and receiver cell types, with gradients representing interaction count and weight.

Figure 1. Single-cell RNA sequencing analysis of HCC patients. (A) UMAP visualization of cell clusters derived from scRNA-seq data of HCC samples. (B) UMAP projection showing annotated cell types within the HCC tumor microenvironment. (C) Bubble plot illustrating the expression of marker genes across identified cell populations. (D) Heatmap displaying the top 10 characteristic genes for each cell type. (E) Network representation of KEGG pathway enrichment analysis across cell types. (F, G) CellChat analysis depicting intercellular communication networks within the HCC microenvironment.

3.2 Enhanced bile acid metabolism in malignant HCC cells

Our single-cell bile acid metabolism profiling revealed pronounced metabolic heterogeneity across cellular populations, with malignant cells exhibiting the highest bile acid metabolism scores (Figures 2A–C). Dimensionality reduction and clustering analyses further identified five distinct malignant cell subtypes, among which G6PD+ malignant cells displayed a uniquely elevated bile acid metabolic signature (Figures 2D, E; Supplementary Figure S1). Pathway enrichment analysis demonstrated significant activation of bile acid–related signaling pathways within this metabolically distinct subpopulation (Figure 2F). Pseudotime trajectory analysis positioned the G6PD+ malignant cells at a critical transitional state within the malignant cells differentiation continuum (Figures 2G, H), showing a dynamic expression pattern closely paralleling that of the characteristic bile acid–synthesizing enzyme CYP7A1 (Figure 2I).Collectively, these findings suggest that G6PD+ malignant cells represent a pivotal subpopulation in the regulation of bile acid metabolism, providing novel mechanistic insight into the metabolic dysregulation underlying HCC progression.

Figure 2
Graphical representation of data on bile acid pathway enrichment scores across different cell types, including B cells, T cells, and hepatocytes. UMAP plots display cell clustering and scoring, with color gradients showing pseudotime and expression levels. Dot plots illustrate percent and average expression for bile acid pathways, while a network diagram categorizes gene activities like transferase and kinase activities. The gene expression over pseudotime is shown for G6PD and CYP7A1.

Figure 2. Bile acid metabolic features in scRNA-seq analysis of HCC patients. (A) Activity scores of bile acid metabolism across different cell types. (B, C) UMAP visualization of bile acid metabolism scores stratified by cell type. (D) UMAP projection depicting malignant and epithelial subclusters. (E) Comparative analysis of bile acid metabolism scores between malignant and epithelial subpopulations. (F) Network representation of enriched signaling cascades in G6PD+ malignant cells. (G) Pseudotime trajectory reconstruction of malignant cell subpopulations. (H) Mapping of differentiation states within the G6PD+ malignant cell cluster. (I) Dynamic expression patterns of G6PD+ malignant cell markers along the differentiation trajectory.

3.3 Spatial niches of bile acid metabolism

To characterize the spatial distribution patterns of bile acid metabolism, we performed spatial transcriptomic analysis using dataset 226334_D5_2. Spatial imaging and regional annotation identified seven distinct microenvironmental regions within HCC tissues (Figure 3A). Cell enrichment analysis based on the Multimodal Intersection Analysis (MIA) algorithm revealed that Region 4 represented a key area enriched with G6PD+ malignant cells (Figure 3B). Spatial pattern analysis demonstrated that G6PD+ malignant cells were non-uniformly distributed across the tumor tissue (Figure 3C), showing marked accumulation in Region 4, which displayed a complex composition of multiple cell types (Figure 3D). Notably, key genes involved in bile acid metabolism, including CYP7A1, CYP27A1, NR1H4, CYP7B1, and BSEP, exhibited significantly variable expression levels among these spatial domains (Figures 3E, F), indicating pronounced spatial heterogeneity of bile acid metabolism within the tumor microenvironment. Collectively, these findings suggest that bile acid metabolism in G6PD+ malignant cells may play a region-specific role within distinct spatial niches of the HCC microenvironment.

Figure 3
A composite image with multiple panels displaying scientific visualizations and analyses:   Panel A shows tissue sections with COL1A1 and GPC3 staining, highlighting spatial niches.   Panel B features a heatmap illustrating gene depletion and enrichment across different cell types and niches.   Panel C presents a distribution of G6PD+ Hepatocytes, color-coded by value from low to high.   Panel D includes spatial maps with pie charts identifying cell types such as B cells and Macrophages.   Panel E shows a map of bile acid metabolic activity.   Panel F displays maps for bile acid metabolic genes, including CYP3A4 and BAAT, color-coded by expression levels.

Figure 3. Spatial transcriptomic analysis of bile acid metabolic features in G6PD+ malignant cells. (A) Spatial imaging and regional annotation of HCC tissue sections. (B) Multimodal Intersection Analysis (MIA) of cell enrichment across different spatial regions. (C) Spatial distribution pattern of G6PD+ malignant cell populations within tumor tissues. (D) Cellular composition map of tumor microenvironmental domains enriched in G6PD+ malignant cells. (E, F) Spatial expression profiles of bile acid metabolism–related genes within tumor regions.

3.4 Construction of a risk prediction model using machine learning

To develop a robust prognostic model, we integrated bile acid metabolism–related features using multiple machine learning algorithms. The results demonstrated that different algorithmic combinations yielded varying concordance index (C-index) values across cohorts, with our integrated model achieving superior predictive performance in all datasets (Figure 4A). Survival analyses across four independent cohorts—TCGA, GSE14520, GSE43619, and ICGC—consistently revealed that patients in the high-risk group had significantly worse overall survival compared to those in the low-risk group (p < 0.001), with corresponding hazard ratios (HRs) of 2.42, 2.37, infinity, and 3.56, respectively (Figure 4B). When compared with existing prognostic models, our model achieved consistently higher C-index values across multiple datasets, underscoring its superior prognostic accuracy (Figure 4C). These findings indicate that the bile acid metabolism–based risk scoring model can effectively predict clinical outcomes in HCC patients, thereby providing valuable guidance for precision prognosis assessment and clinical decision-making.

Figure 4
Three panels display data related to survival analysis and C-index scores across different cohorts.  Panel A: Heatmap shows varying C-index scores for methods across cohorts GSE14520, GSE43619, ICGC, and TCGA, with a color gradient from 0.5 to 0.9.  Panel B: Four Kaplan-Meier plots compare survival probabilities in cohorts TCGA, GSE14520, GSE43619, and ICGC, highlighting hazard ratios and significance values.  Panel C: Bar charts compare C-index scores from different studies across the four cohorts, with specific values for each case.

Figure 4. Construction of a bile acid metabolism–based risk scoring model derived from G6PD+ malignant cells using integrated machine learning. (A) Risk assessment of different models across four independent cohorts evaluated by concordance index (C-index). (B) Comparative survival analysis between high-risk and low-risk patient subgroups. (C) Performance comparison between the G6PD+ malignant cell–based predictive model and previously established published prognostic models.

3.5 Performance evaluation of the G6PD+ malignant cell–based predictive model

To validate the independent prognostic value of the bile acid metabolism–derived risk score, we performed univariate and multivariate Cox regression analyses. Univariate analysis revealed that, apart from sex, race, age, and N stage, variables such as M stage, T stage, overall clinical stage, and the risk score were all significant prognostic factors (Figure 5A). Multivariate analysis further confirmed that the risk score remained an independent prognostic indicator after adjustment for other clinical covariates (p < 0.001, HR = 1.060, 95% CI: 1.052–1.068) (Figure 5B). Based on these findings, we constructed an integrated nomogram combining the risk score with key clinical parameters (Figure 5C), which demonstrated excellent calibration performance (Figure 5D). Time-dependent ROC curve analysis showed that our model outperformed other clinical variables in predicting 1-, 3-, and 5-year survival, achieving AUC values of 0.730, 0.752, and 0.689, respectively (Figure 5E). Collectively, these results establish the bile acid metabolism–based risk score as an independent and powerful tool for prognostic assessment in HCC.

Figure 5
Charts illustrating statistical analysis results:   A) Univariate Cox regression forest plot displaying hazard ratios and p-values for variables like gender, race, and age.   B) Multivariate Cox regression forest plot with similar data and additional variables.   C) A nomogram illustrating point assignments for different clinical variables to predict outcomes.   D) Calibration plot comparing observed and nomogram-predicted overall survival percentages over one, three, and five years.   E) ROC curves showing sensitivity and specificity for one, three, and five-year predictions with respective AUC values.

Figure 5. Validation of the independent prognostic value of the bile acid metabolism–based risk score. (A) Univariate Cox regression analysis of clinicopathological parameters and the G6PD+ malignant cell–derived risk score. (B) Multivariate Cox regression analysis assessing the independence of the G6PD+ malignant cell–based risk score from clinical covariates. (C) Nomogram integrating common clinical parameters with the G6PD+ malignant cell–based risk score for prognostic prediction. (D) Calibration curve evaluating the predictive accuracy of the nomogram. (E) Time-dependent ROC curves illustrating the predictive performance of the model at 1-, 3-, and 5-year survival intervals.

3.6 Distinct molecular characteristics between risk groups

Comparative analyses between risk groups revealed substantial differences in molecular and clinical features. Pathway enrichment analysis indicated that pro-tumorigenic signaling pathways were significantly activated in the high-risk group (Figure 6A). Mutation profiling demonstrated a higher TP53 mutation frequency among high-risk patients (Figure 6B), accompanied by a markedly elevated tumor mutation burden (TMB) (p < 0.001; Figure 6C). Both high-risk classification and increased TMB were associated with significantly poorer overall survival (Figure 6D), thereby validating the prognostic reliability of our risk stratification framework. Collectively, these findings provide a molecular and pathological basis for bile acid metabolism–related risk stratification in HCC, highlighting its potential utility in guiding prognosis and individualized therapeutic strategies.

Figure 6
Heatmap and cluster analysis with risk levels showing gene expressions and categories. Mutation plots compare low and high risk groups, categorized by mutation types. Scatter plots depict correlations between risk scores and variants or tumor mutational burden (TMB). Kaplan-Meier survival curves compare different risk and variant groups, with statistical significance indicated.

Figure 6. Molecular characteristics of risk groups based on G6PD+ malignant cells. (A) Comparative pathway enrichment analysis between high-risk and low-risk patient groups. (B) Mutation landscape analysis illustrating genomic alterations across high- and low-risk groups. (C) Correlation analysis between the G6PD+ malignant cell–derived risk score and tumor mutation burden (TMB). (D) Survival analysis stratified by combined risk score and TMB status.

3.7 Immune infiltration characteristics of risk subgroups

To characterize the immune microenvironment associated with G6PD+ malignant cells, we conducted a comprehensive immune analysis using single-sample gene set enrichment analysis (ssGSEA). Distinct immune landscapes were observed between the risk subgroups (Figure 7A), with the high-risk group showing significant enrichment of specific immune cell populations, including various T-cell subsets, macrophage subtypes, and other immune cells (Figure 7B). Pathway enrichment analysis further revealed pronounced differences in immune-related functional pathways among high-risk patients (Figure 7C), accompanied by differential expression of immune checkpoint molecules (Figure 7D). Collectively, these findings suggest that G6PD+ malignant cells with aberrant bile acid metabolism may contribute to HCC progression by modulating the immune microenvironment, offering mechanistic insight into the interplay between metabolic reprogramming and immune regulation in HCC.

Figure 7
This image consists of four panels related to TCGA LIHC samples. Panel A is a heatmap showing immune infiltration and various cell types across 363 samples with annotations, including sex, status, age, stage, risk score, and immune infiltration levels. Panel B displays box plots comparing cell fractions of immune cell types between high and low-risk groups. Panel C presents box plots of various scores related to risk groups. Panel D shows box plots of gene expression levels, comparing high and low-risk groups across several genes. Color codes indicate risk and statistical significance levels.

Figure 7. Immune infiltration characteristics of risk groups based on G6PD+ malignant cells. (A) Immune infiltration landscape of HCC patients. (B) Comparative analysis of immune cell infiltration between high-risk and low-risk groups. (C) Functional characterization of immune responses in high- and low-risk cohorts. (D) Differential expression patterns of immune checkpoint molecules across risk groups. *p < 0.05, **p < 0.01, ***p < 0.001, ns > 0.05.

3.8 Strong interactions between G6PD+ malignant cells and endothelial cells in aberrant bile acid metabolism

CellChat analysis revealed that G6PD+ malignant cells exhibited the strongest intercellular interactions with endothelial cells (Figure 8A), accompanied by significant activation of the VEGFA–VEGFR1 and NAMPT–INSR signaling pathways (Figure 8B). Consistently, GSEA analysis demonstrated that the high-risk signature was closely associated with specific signaling pathways involving bile acid metabolism and endothelial cell regulation (Figure 8C). Both pathways were significantly enriched in the high-risk group, suggesting a potential synergistic role in promoting HCC progression. These findings confirm that G6PD+ malignant cells with aberrant bile acid metabolism may drive HCC progression by reshaping tumor microenvironmental signaling networks, particularly through the enhancement of angiogenic activity.

Figure 8
Panel A shows network graphs indicating relationships between cell types based on parameters labeled “weight” and “count.” G6PD+Hepatocytes and Endothelial_cells are highlighted. Panel B features a strip with VEGFA-VEGFR1 and NAMPT-INSR, indicating a pathway from G6PD+Hepatocytes to Endothelial_cells with a p-value marker. Panel C contains two line graphs for “Bile Acid Metabolic Pathway” and “Endothelial Cells” showing enrichment scores (NES of 1.608 and 2.571, respectively) across risk levels, both with FDR under 0.05.

Figure 8. Roles of bile acid metabolism and specific immune cell pathways in HCC prognosis. (A) Cell–cell communication network between G6PD+ malignant cells and other cell populations within the HCC microenvironment. (B) Signaling crosstalk between G6PD+ malignant cells and endothelial cells, highlighting activation of the VEGFA–VEGFR1 and NAMPT–INSR pathways. (C) GSEA analysis demonstrating the correlation between the G6PD+ malignant cell–derived risk score and pathways related to bile acid metabolism and endothelial cell function.

3.9 BIBR-1532 as a potential therapeutic compound targeting G6PD+ malignant cells

Based on machine learning–driven gene importance analysis, G6PD was identified as the most critical regulatory factor (Figure 9A), with its elevated expression strongly correlated with poor prognosis in HCC patients (Figure 9B). Drug sensitivity analysis further revealed a significant association between G6PD expression and the anticancer response to BIBR-1532 (Figure 9C). Molecular docking simulations confirmed that BIBR-1532 forms stable interactions with the G6PD protein through key residues PHE173, TYR147, and PHE253 (Figures 9D–F), suggesting a high binding affinity and potential inhibitory effect. Collectively, these findings provide a theoretical foundation for precision therapeutics targeting G6PD, highlighting BIBR-1532 as a promising candidate for the treatment of HCC driven by aberrant bile acid metabolism.

Figure 9
Chart A shows the importance of genes, with G6PD highlighted. Graph B illustrates overall survival rates, comparing low and high G6PD groups. Graph C correlates drug sensitivity with gene expression, showing a positive trend. Diagram D depicts a molecular structure with highlighted interactions. Close-up E highlights specific amino acid interactions (PHE 253, TYR 147, PHE 173). Diagram F presents a molecular interaction model with a LibDockScore of 57.13, detailing different types of interactions such as Van der Waals and hydrogen bonds.

Figure 9. G6PD as a key regulator and potential drug target in bile acid metabolism. (A) Feature importance analysis of genes in the risk prediction model. (B) Survival outcomes stratified by G6PD expression levels. (C) Correlation analysis between G6PD expression and therapeutic drug sensitivity. (D) Molecular docking simulation of BIBR-1532 with the G6PD-encoded protein. (E) Structural characterization of the BIBR-1532–G6PD binding interactions. (F) Two-dimensional visualization of the binding interface between BIBR-1532 and G6PD.

3.10 In Vitro validation and functional characterization of G6PD

To verify the functional role of G6PD in HCC progression, we performed a series of in vitro experiments. First, qRT-PCR confirmed the efficiency of both G6PD overexpression and knockdown (Figure 10A). Western blot analysis further validated corresponding changes in protein expression—G6PD protein levels were markedly increased in the OE-G6PD group, whereas significantly reduced expression was observed in the Si-G6PD-1 and Si-G6PD-3 groups (Figure 10B). Cell proliferation assays demonstrated that G6PD overexpression promoted proliferation, reflected by stronger fluorescence intensity signals (Figure 10C), while G6PD knockdown significantly inhibited HCC cell growth (Figure 10D). Colony formation assays further confirmed that G6PD knockdown (Si-G6PD-1 and Si-G6PD-3) markedly reduced colony numbers, whereas overexpression of G6PD enhanced colony-forming ability (Figures 10E, F). Transwell invasion assays revealed that G6PD knockdown (Si-G6PD-1 and Si-G6PD-3) led to a substantial reduction in the number of invading cells, while the OE-G6PD group exhibited a marked increase in invasion potential (Figures 10G, H). Similarly, silencing G6PD notably impaired cell migratory capacity compared with controls, whereas overexpression significantly enhanced migration (Figures 10I–K). Collectively, these findings demonstrate that G6PD+ malignant cells with aberrant bile acid metabolism promote HCC progression by enhancing tumor cell proliferation, migration, and invasion capabilities, thereby underscoring G6PD as a functional driver and potential therapeutic target in HCC.

Figure 10
Graphs and images depict the effects of manipulating G6PD expression on various cellular processes. Bar graphs show relative mRNA and protein levels, colony formation, and cell migration. Fluorescent images illustrate EDU-positive cells. Colony images and migration assays further evaluate cell behavior. Statistical markers show significance.

Figure 10. Functional validation of G6PD in HCC cells. (A) qRT-PCR analysis confirming the efficiency of G6PD overexpression and knockdown in HCCLM3 cells (n = 3). (B) Western blot validation of G6PD protein expression levels following overexpression or silencing in HCCLM3 cells (n = 3). (C, D) EdU assays evaluating the proliferative capacity of HCCLM3 cells after G6PD interference (n = 3). (E, F) Colony formation assays assessing the tumorigenic potential of HCCLM3 cells after G6PD knockdown (n = 3). (G, H) Transwell invasion assays evaluating the invasive capacity of HCCLM3 cells after G6PD knockdown or overexpression (n = 3). (I–K) Wound-healing assays examining the migratory ability of HCCLM3 cells following G6PD modulation (n = 3). *p < 0.05, **p < 0.01.

4 Discussion

Tumor tissues exhibit pronounced spatial heterogeneity, most notably at the invasive front where cells interacting with the extracellular matrix and vasculature display enhanced invasiveness and metastatic capacity (26, 27). In our analysis, we identified a population of G6PD+ malignant cells with aberrant bile acid metabolism that are enriched at the tumor boundary regions of HCC. These cells exhibit high G6PD expression, indicating a hyperactivated pentose phosphate pathway (PPP) that provides sufficient NADPH to counter oxidative stress and support biosynthetic demands. This metabolic feature likely reflects adaptation to a bile acid–rich, inflammatory microenvironment, where excessive bile acids induce ROS accumulation, compelling tumor cells to upregulate G6PD and antioxidant defenses (28, 29). Notably, the enrichment of G6PD+ cells at the invasive front suggests their potential role as drivers of invasion and dissemination (30). Previous studies have shown that elevated G6PD expression drives malignant phenotypes in HCC, including enhanced migration, invasion, and EMT (31). Our findings are consistent with these reports, implying that G6PD+ tumor cells at the periphery may act as the “vanguard” of tumor expansion.

Dysregulated bile acid metabolism not only induces hepatotoxic accumulation and inflammation but also affects tumor cell behavior via multiple signaling pathways. In HCC, elevated bile acids can activate inflammation-related cascades such as NF-κB and JAK–STAT3, as well as metabolic sensing pathways including Hippo–YAP and AMPK (3234),indicating extensive crosstalk between bile acid signaling and proliferative signaling. Meanwhile, bile acid receptors such as TGR5 and FXR regulate both lipid and glucose metabolism, linking metabolic and oncogenic signaling (35, 36). Excessive bile acid activates GPBAR1 specifically in cancer-associated fibroblasts, prompting CXCL10 production that enhances EMT and metastasis in cholangiocarcinoma (37). This process also recruits neutrophils to establish an immunosuppressive microenvironment. Importantly, inhibiting bile acid synthesis in liver cancer cells reprograms this milieu, reinvigorating T-cell responses, suppressing tumor growth, and consequently sensitizing tumors to anti-PD-1 immunotherapy (38). While previous studies have confirmed the role of disordered bile acid metabolism in HCC progression at the molecular level, our research focuses on the G6PD+ malignant cellular subtype with aberrant bile acid metabolism and its functional impact, thereby further substantiating the critical role of bile acid metabolism in HCC progression at the cellular level.

Importantly, glucose-6-phosphate dehydrogenase (G6PD), the rate-limiting enzyme of the PPP, may itself be regulated by bile acid–related oxidative signaling. Oxidized bile acids can induce ROS production and activate transcription factors like NRF2, which transcriptionally upregulates G6PD to enhance redox balance (28, 29). Mice deficient in NRF2 exhibit reduced HCC incidence in carcinogenic models, partially due to decreased G6PD expression and PPP activity, leading to uncompensated oxidative stress (28). These findings imply that bile acid–induced oxidative stress promotes tumor survival through the NRF2–G6PD axis, forming a feedback loop that links bile acid and glucose metabolism. We propose that metabolic crosstalk drives HCC progression by enhancing tumor cell invasiveness and anti-apoptotic resistance, allowing them to thrive in harsh microenvironments.

HCC is a hypervascular tumor, whose growth and metastasis are tightly dependent on vascular remodeling, including both classical angiogenesis and non-classical vasculogenic mimicry (VM) formation (3941). We observed that G6PD+ malignant cells with aberrant bile acid metabolism maintain intimate interactions with endothelial cells and macrophages in the tumor microenvironment. At the invasive margin, G6PD+ malignant cells are often spatially adjacent to endothelial cells, suggesting direct cellular crosstalk. On one hand, these tumor cells may secrete pro-angiogenic factors such as VEGF and basic fibroblast growth factor (bFGF) to stimulate endothelial proliferation and neovascularization (42, 43). Conversely, endothelial cells release pro-invasive cues, enhancing tumor cell motility and invasiveness (42, 44, 45). This bidirectional feedback facilitates mutual activation between invasive-front tumor cells and neighboring endothelium, thereby driving angiogenesis and local invasion.

Macrophages likewise play pivotal roles in angiogenesis and immune remodeling of the tumor microenvironment (46, 47). We speculate that metabolic byproducts derived from G6PD+ malignant cells may drive macrophage M2 polarization, promoting secretion of pro-angiogenic and immunosuppressive mediators. Previous studies have shown that under conditions of bile acid dysregulation, tumor-derived inflammatory mediators (e.g., SAA proteins) recruit macrophages and bias them toward an M2 phenotype (30). This suggests that bile acid imbalance modulates the malignant cell–macrophage axis, reshaping both immune and vascular niches. In our results, the interaction between G6PD+ malignant cells and macrophages was particularly strong, implying that G6PD+ cells may secrete cytokines that attract macrophages to the tumor margin. The recruited macrophages, in turn, release growth factors and cytokines that further enhance angiogenesis and invasion (48). Collectively, G6PD+ malignant cells with aberrant bile acid metabolism interact dynamically with endothelial and immune cells, establishing a pro-angiogenic and immunosuppressive microenvironment that supports HCC progression. Therefore, targeting these intercellular interactions or blocking key pro-angiogenic signals may represent promising strategies to impede vascular remodeling and metastatic potential in HCC.

Our study identifies a promising bile acid metabolism-related prognostic signature, yet its clinical application faces two main challenges. First, standardization is difficult due to technical variations in RNA-seq and sample processing, as well as biological heterogeneity in ethnicity, lifestyle, and gut microbiota, which may affect signature performance. Second, the model requires broader validation, and multi-center evaluation. Moreover, since bile acid metabolism influences multiple cancers, testing this signature across tumor types will clarify whether it serves as a pan-cancer biomarker or a context-specific tool, guiding its use in precision oncology.

Nevertheless, several limitations should be acknowledged. First, the causal link between G6PD+ malignant cells and bile acid metabolic dysregulation remains correlative; functional studies are needed to confirm whether bile acid signaling directly modulates G6PD activity and its oncogenic consequences. Second, while spatial analyses revealed G6PD+ cell enrichment at the tumor margin, it remains unclear whether this is a universal phenomenon across different HCC stages and differentiation states—requiring validation in larger cohorts. Third, therapeutic targeting of metabolism must consider potential safety constraints; for example, G6PD inhibition may trigger hemolysis, particularly in individuals with underlying G6PD deficiency. Thus, achieving tumor-specific targeting while sparing normal tissues remains a key challenge. Moreover, all functional assays for proliferation, migration, and invasion were conducted using a single HCC cell line (HCCLM3). While this line was selected for its aggressive metastatic phenotype, which aligns with our research focus on progression, conclusions drawn from a single model system require caution in generalization to the highly heterogeneous landscape of HCC. Finally, the bile acid–immune–metabolic interplay is exceedingly complex; its dual (“double-edged sword”) effects warrant balanced consideration. Future studies should leverage multi-omics and spatial single-cell technologies to map the bile acid–metabolism–immunity network, identify critical nodes, and guide the discovery of precise intervention targets and specific regulatory mechanism. Future validation may involve a multicenter retrospective cohort with broader clinical features, combined with PDOX models derived from co-cultured patient derived organoids and autologous immune cells.

In summary, this study definitively identifies a G6PD+ malignant hepatocyte subset with reprogrammed bile acid metabolism, situated within immunosuppressive and angiogenic niches, and establishes its derived gene signature as a robust, independent prognostic indicator for HCC. Looking forward, integrating metabolic-targeted approaches with current treatment paradigms, guided by cross-disciplinary strategies, may open new avenues for precision therapy in HCC and improve patient survival outcomes.

Data availability statement

The data analyzed in this study are publicly available and were obtained from public repositories. Processed data and analysis results generated in this study are provided in the article and Supplementary Material. Further inquiries can be directed to the corresponding authors.

Author contributions

XJ: Writing – original draft, Visualization, Formal Analysis, Conceptualization, Data curation. HQ: Software, Writing – review & editing, Funding acquisition, Formal Analysis, Validation. TY: Validation, Resources, Investigation, Writing – review & editing, Software. HY: Writing – review & editing, Investigation, Funding acquisition. YL: Validation, Funding acquisition, Writing – review & editing. BP: Software, Writing – review & editing, Formal Analysis. XY: Resources, Writing – review & editing, Investigation. WZ: Investigation, Validation, Writing – review & editing. HC: Project administration, Resources, Supervision, Conceptualization, Methodology, Writing – review & editing, Data curation. RL: Methodology, Data curation, Supervision, Resources, Conceptualization, Writing – review & editing, Project administration.

Funding

The author(s) declared that financial support was received for this work and/or its publication. This work was funded by grants from the Hunan Provincial Natural Science Foundation Project (No. 2023JJ60206, 2025JJ80322, 2025JJ80345, 2026JJ81759).

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

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

References

1. Wang H, Lun Y, Xu D, Jiang H, Yan Y, and Yang X. Research progress and therapeutic strategies in hepatocellular carcinoma metabolic reprogramming. J Adv Res. (2025), S2090–1232. doi: 10.1016/j.jare.2025.08.023

PubMed Abstract | Crossref Full Text | Google Scholar

2. Park S and Hall MN. Metabolic reprogramming in hepatocellular carcinoma: mechanisms and therapeutic implications. Exp Mol Med. (2025) 57:515–23. doi: 10.1038/s12276-025-01415-2

PubMed Abstract | Crossref Full Text | Google Scholar

3. Wang K, Li X, Guo S, Chen J, Lv Y, Guo Z, et al. Metabolic reprogramming of glucose: the metabolic basis for the occurrence and development of hepatocellular carcinoma. Front Oncol. (2025) 15:1545086. doi: 10.3389/fonc.2025.1545086

PubMed Abstract | Crossref Full Text | Google Scholar

4. Du D, Liu C, Qin M, Zhang X, Xi T, Yuan S, et al. Metabolic dysregulation and emerging therapeutical targets for hepatocellular carcinoma. Acta Pharm Sin B. (2022) 12:558–80. doi: 10.1016/j.apsb.2021.09.019

PubMed Abstract | Crossref Full Text | Google Scholar

5. Vander Heiden MG, Cantley LC, and Thompson CB. Understanding the Warburg effect: the metabolic requirements of cell proliferation. Science. (2009) 324:1029–33. doi: 10.1126/science.1160809

PubMed Abstract | Crossref Full Text | Google Scholar

6. Huang Q, Tan Y, Yin P, Ye G, Gao P, Lu X, et al. Metabolic characterization of hepatocellular carcinoma using nontargeted tissue metabolomics. Cancer Res. (2013) 73:4992–5002. doi: 10.1158/0008-5472.Can-13-0308

PubMed Abstract | Crossref Full Text | Google Scholar

7. Liu Y, Sun Y, Guo Y, Shi X, Chen X, Feng W, et al. An overview: the diversified role of mitochondria in cancer metabolism. Int J Biol Sci. (2023) 19:897–915. doi: 10.7150/ijbs.81609

PubMed Abstract | Crossref Full Text | Google Scholar

8. Liu Y, Zhu J, Jin Y, Sun Z, Wu X, Zhou H, et al. Disrupting bile acid metabolism by suppressing Fxr causes hepatocellular carcinoma induced by YAP activation. Nat Commun. (2025) 16:3583. doi: 10.1038/s41467-025-58809-z

PubMed Abstract | Crossref Full Text | Google Scholar

9. Zhang W, Zhang Y, Wan Y, Liu Q, and Zhu X. A bile acid-related prognostic signature in hepatocellular carcinoma. Sci Rep. (2022) 12:22355. doi: 10.1038/s41598-022-26795-7

PubMed Abstract | Crossref Full Text | Google Scholar

10. Jüngst C, Berg T, Cheng J, Green RM, Jia J, Mason AL, et al. Intrahepatic cholestasis in common chronic liver diseases. Eur J Clin Invest. (2013) 43:1069–83. doi: 10.1111/eci.12128

PubMed Abstract | Crossref Full Text | Google Scholar

11. Hirschfield GM, Heathcote EJ, and Gershwin ME. Pathogenesis of cholestatic liver disease and therapeutic approaches. Gastroenterology. (2010) 139:1481–96. doi: 10.1053/j.gastro.2010.09.004

PubMed Abstract | Crossref Full Text | Google Scholar

12. Riscal R, Gardner SM, Coffey NJ, Carens M, Mesaros C, Xu JP, et al. Bile acid metabolism mediates cholesterol homeostasis and promotes tumorigenesis in clear cell renal cell carcinoma. Cancer Res. (2024) 84:1570–82. doi: 10.1158/0008-5472.Can-23-0821

PubMed Abstract | Crossref Full Text | Google Scholar

13. Qin LX and Tang ZY. Hepatocellular carcinoma with obstructive jaundice: diagnosis, treatment and prognosis. World J Gastroenterol. (2003) 9:385–91. doi: 10.3748/wjg.v9.i3.385

PubMed Abstract | Crossref Full Text | Google Scholar

14. Goodwin B, Jones SA, Price RR, Watson MA, McKee DD, Moore LB, et al. A regulatory cascade of the nuclear receptors FXR, SHP-1, and LRH-1 represses bile acid biosynthesis. Mol Cell. (2000) 6:517–26. doi: 10.1016/s1097-2765(00)00051-4

PubMed Abstract | Crossref Full Text | Google Scholar

15. Lu TT, Makishima M, Repa JJ, Schoonjans K, Kerr TA, Auwerx J, et al. Molecular basis for feedback regulation of bile acid synthesis by nuclear receptors. Mol Cell. (2000) 6:507–15. doi: 10.1016/s1097-2765(00)00050-2

PubMed Abstract | Crossref Full Text | Google Scholar

16. Pan Z, Jin X, Li Q, Zhou Y, Zeng Y, Wang X, et al. Citrus pectin supplementation alleviated hepatic lipid accumulation through gut microbiota indole lactic acid promoting hepatic bile acid synthesis and excretion. Int J Biol Sci. (2025) 21:5015–33. doi: 10.7150/ijbs.116929

PubMed Abstract | Crossref Full Text | Google Scholar

17. Tian Z, Song J, She J, He W, Guo S, and Dong B. Constructing a disulfidptosis-related prognostic signature of hepatocellular carcinoma based on single-cell sequencing and weighted co-expression network analysis. Apoptosis. (2024) 29:1632–47. doi: 10.1007/s10495-024-01968-z

PubMed Abstract | Crossref Full Text | Google Scholar

18. Li X, Li R, Miao X, Zhou X, Wu B, Cao J, et al. Integrated single cell analysis reveals an atlas of tumor associated macrophages in hepatocellular carcinoma. Inflammation. (2024) 47:2077–93. doi: 10.1007/s10753-024-02026-1

PubMed Abstract | Crossref Full Text | Google Scholar

19. Liang R, Hong W, Zhang Y, Ma D, Li J, Shi Y, et al. Deep dissection of stemness-related hierarchies in hepatocellular carcinoma. J Transl Med. (2023) 21:631. doi: 10.1186/s12967-023-04425-8

PubMed Abstract | Crossref Full Text | Google Scholar

20. Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan CH, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun. (2021) 12:1088. doi: 10.1038/s41467-021-21246-9

PubMed Abstract | Crossref Full Text | Google Scholar

21. Ye J, Gao X, Huang X, Huang S, Zeng D, Luo W, et al. Integrating single-cell and spatial transcriptomics to uncover and elucidate GP73-mediated pro-angiogenic regulatory networks in hepatocellular carcinoma. Res (Wash D C). (2024) 7:387. doi: 10.34133/research.0387

PubMed Abstract | Crossref Full Text | Google Scholar

22. Cable DM, Murray E, Zou LS, Goeva A, Macosko EZ, Chen F, et al. Robust decomposition of cell type mixtures in spatial transcriptomics. Nat Biotechnol. (2022) 40:517–26. doi: 10.1038/s41587-021-00830-w

PubMed Abstract | Crossref Full Text | Google Scholar

23. Liu H, Zhang W, Zhang Y, Adegboro AA, Fasoranti DO, Dai L, et al. Mime: A flexible machine-learning framework to construct and visualize models for clinical characteristics prediction and feature selection. Comput Struct Biotechnol J. (2024) 23:2798–810. doi: 10.1016/j.csbj.2024.06.035

PubMed Abstract | Crossref Full Text | Google Scholar

24. Zhengdong A, Xiaoying X, Shuhui F, Rui L, Zehui T, Guanbin S, et al. Identification of fatty acids synthesis and metabolism-related gene signature and prediction of prognostic model in hepatocellular carcinoma. Cancer Cell Int. (2024) 24:130. doi: 10.1186/s12935-024-03306-4

PubMed Abstract | Crossref Full Text | Google Scholar

25. Guo Y, Lin Z, Zhang W, Chen H, Chen Y, Liu Y, et al. Comprehensive multi-omics analysis of nucleotide metabolism: elucidating the role and prognostic significance of UCK2 in bladder cancer. Funct Integr Genomics. (2025) 25:133. doi: 10.1007/s10142-025-01642-w

PubMed Abstract | Crossref Full Text | Google Scholar

26. Schürch CM, Bhate SS, Barlow GL, Phillips DJ, Noti L, Zlobec I, et al. Coordinated cellular neighborhoods orchestrate antitumoral immunity at the colorectal cancer invasive front. Cell. (2020) 182:1341–1359.e1319. doi: 10.1016/j.cell.2020.07.005

PubMed Abstract | Crossref Full Text | Google Scholar

27. Grünwald BT, Devisme A, Andrieux G, Vyas F, Aliar K, McCloskey CW, et al. Spatially confined sub-tumor microenvironments in pancreatic cancer. Cell. (2021) 184:5577–5592.e5518. doi: 10.1016/j.cell.2021.09.022

PubMed Abstract | Crossref Full Text | Google Scholar

28. Evert M, Calvisi DF, Evert K, De Murtas V, Gasparetti G, Mattu S, et al. V-AKT murine thymoma viral oncogene homolog/mammalian target of rapamycin activation induces a module of metabolic changes contributing to growth in insulin-induced hepatocarcinogenesis. Hepatology. (2012) 55:1473–84. doi: 10.1002/hep.25600

PubMed Abstract | Crossref Full Text | Google Scholar

29. Ngo HKC, Kim DH, Cha YN, Na HK, and Surh YJ. Nrf2 mutagenic activation drives hepatocarcinogenesis. Cancer Res. (2017) 77:4797–808. doi: 10.1158/0008-5472.Can-16-3538

PubMed Abstract | Crossref Full Text | Google Scholar

30. Wu L, Yan J, Bai Y, Chen F, Zou X, Xu J, et al. An invasive zone in human liver cancer identified by Stereo-seq promotes hepatocyte-tumor cell crosstalk, local immunosuppression and tumor progression. Cell Res. (2023) 33:585–603. doi: 10.1038/s41422-023-00831-1

PubMed Abstract | Crossref Full Text | Google Scholar

31. Lu M, Lu L, Dong Q, Yu G, Chen J, Qin L, et al. Elevated G6PD expression contributes to migration and invasion of hepatocellular carcinoma cells by inducing epithelial-mesenchymal transition. Acta Biochim Biophys Sin (Shanghai). (2018) 50:370–80. doi: 10.1093/abbs/gmy009

PubMed Abstract | Crossref Full Text | Google Scholar

32. Attia YM, Tawfiq RA, Ali AA, and Elmazar MM. The FXR agonist, obeticholic acid, suppresses HCC proliferation & Metastasis: role of IL-6/STAT3 signalling pathway. Sci Rep. (2017) 7:12502. doi: 10.1038/s41598-017-12629-4

PubMed Abstract | Crossref Full Text | Google Scholar

33. Beyoğlu D and Idle JR. The gut microbiota - A vehicle for the prevention and treatment of hepatocellular carcinoma. Biochem Pharmacol. (2022) 204:115225. doi: 10.1016/j.bcp.2022.115225

PubMed Abstract | Crossref Full Text | Google Scholar

34. Aoki T, Nishida N, Kurebayashi Y, Sakai K, Fujiwara N, Tsurusaki M, et al. Molecular classification of hepatocellular carcinoma based on zoned metabolic feature and oncogenic signaling pathway. Clin Mol Hepatol. (2025) 31:981–1002. doi: 10.3350/cmh.2024.1088

PubMed Abstract | Crossref Full Text | Google Scholar

35. Claudel T, Staels B, and Kuipers F. The Farnesoid X receptor: a molecular link between bile acid and lipid and glucose metabolism. Arterioscler Thromb Vasc Biol. (2005) 25:2020–30. doi: 10.1161/01.ATV.0000178994.21828.a7

PubMed Abstract | Crossref Full Text | Google Scholar

36. Inagaki T, Choi M, Moschetta A, Peng L, Cummins CL, McDonald JG, et al. Fibroblast growth factor 15 functions as an enterohepatic signal to regulate bile acid homeostasis. Cell Metab. (2005) 2:217–25. doi: 10.1016/j.cmet.2005.09.001

PubMed Abstract | Crossref Full Text | Google Scholar

37. Huang F, Liu Z, Song Y, Wang G, Shi A, Chen T, et al. Bile acids activate cancer-associated fibroblasts and induce an immunosuppressive microenvironment in cholangiocarcinoma. Cancer Cell. (2025) 43:1460–1475.e1410. doi: 10.1016/j.ccell.2025.05.017

PubMed Abstract | Crossref Full Text | Google Scholar

38. Varanasi SK, Chen D, Liu Y, Johnson MA, Miller CM, Ganguly S, et al. Bile acid synthesis impedes tumor-specific T cell responses during liver cancer. Science. (2025) 387:192–201. doi: 10.1126/science.adl4100

PubMed Abstract | Crossref Full Text | Google Scholar

39. Liu M, Ruan X, Liu X, Dong W, Wang D, Yang C, et al. The mechanism of BUD13 m6A methylation mediated MBNL1-phosphorylation by CDK12 regulating the vasculogenic mimicry in glioblastoma cells. Cell Death Dis. (2022) 13:1017. doi: 10.1038/s41419-022-05426-z

PubMed Abstract | Crossref Full Text | Google Scholar

40. Qin LN, Zhang H, Li QQ, Wu T, Cheng SB, Wang KW, et al. Vitamin D binding protein (VDBP) hijacks twist1 to inhibit vasculogenic mimicry in hepatocellular carcinoma. Theranostics. (2024) 14:436–50. doi: 10.7150/thno.90322

PubMed Abstract | Crossref Full Text | Google Scholar

41. Zhu AX, Duda DG, Sahani DV, and Jain RK. HCC and angiogenesis: possible targets and future directions. Nat Rev Clin Oncol. (2011) 8:292–301. doi: 10.1038/nrclinonc.2011.30

PubMed Abstract | Crossref Full Text | Google Scholar

42. Feng T, Yu H, Xia Q, Ma Y, Yin H, Shen Y, et al. Cross-talk mechanism between endothelial cells and hepatocellular carcinoma cells via growth factors and integrin pathway promotes tumor angiogenesis and cell migration. Oncotarget. (2017) 8:69577–93. doi: 10.18632/oncotarget.18632

PubMed Abstract | Crossref Full Text | Google Scholar

43. Shi S, Zhu C, Hu Y, Jiang P, Zhao J, and Xu Q. ENG is a biomarker of prognosis and angiogenesis in liver cancer, and promotes the differentiation of tumor cells into vascular ECs. Front Biosci (Landmark Ed). (2024) 29:315. doi: 10.31083/j.fbl2909315

PubMed Abstract | Crossref Full Text | Google Scholar

44. Lamalice L, Le Boeuf F, and Huot J. Endothelial cell migration during angiogenesis. Circ Res. (2007) 100:782–94. doi: 10.1161/01.RES.0000259593.07661.1e

PubMed Abstract | Crossref Full Text | Google Scholar

45. Qiu C, Tang C, Tang Y, Su K, Chai X, Zhan Z, et al. RGS5(+) lymphatic endothelial cells facilitate metastasis and acquired drug resistance of breast cancer through oxidative stress-sensing mechanism. Drug Resist Update. (2024) 77:101149. doi: 10.1016/j.drup.2024.101149

PubMed Abstract | Crossref Full Text | Google Scholar

46. Jin X, Zhang N, Yan T, Wei J, Hao L, Sun C, et al. Lactate-mediated metabolic reprogramming of tumor-associated macrophages: implications for tumor progression and therapeutic potential. Front Immunol. (2025) 16:1573039. doi: 10.3389/fimmu.2025.1573039

PubMed Abstract | Crossref Full Text | Google Scholar

47. Zeng W, Li F, Jin S, Ho PC, Liu PS, and Xie X. Functional polarization of tumor-associated macrophages dictated by metabolic reprogramming. J Exp Clin Cancer Res. (2023) 42:245. doi: 10.1186/s13046-023-02832-9

PubMed Abstract | Crossref Full Text | Google Scholar

48. Cassetta L and Pollard JW. A timeline of tumour-associated macrophage biology. Nat Rev Cancer. (2023) 23:238–57. doi: 10.1038/s41568-022-00547-1

PubMed Abstract | Crossref Full Text | Google Scholar

Keywords: bile acids, G6PD, hepatocellular carcinoma, machine learning, tumor microenvironment

Citation: Jiang X, Quan H, Yin T, Yao H, Li Y, Peng B, Yuan X, Zeng W, Chen H and Li R (2026) Single-cell RNA sequencing and spatial transcriptomic analysis reveal a distinct population of G6PD+ cells with aberrant bile acid metabolism in hepatocellular carcinoma. Front. Immunol. 17:1739293. doi: 10.3389/fimmu.2026.1739293

Received: 04 November 2025; Accepted: 09 January 2026; Revised: 07 January 2026;
Published: 02 February 2026.

Edited by:

Jianwei Xu, Shandong University, China

Reviewed by:

Hai-Feng Zhou, Nanjing Medical University, China
Doblin Sandai, University of Science Malaysia (USM), Malaysia

Copyright © 2026 Jiang, Quan, Yin, Yao, Li, Peng, Yuan, Zeng, Chen and Li. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Honghui Chen, Y2hoaGhjY2hoQHNpbmEuY29t; Rong Li, MTk5MTAyMDAxNUB1c2MuZWR1LmNu

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.