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.

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.

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.

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.

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.

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.

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.

To identify potential therapeutic compounds, we performed integrative drug sensitivity modeling using the CellMiner database. Gene expression matrices were log₂-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.

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% CO₂ 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).

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.