Single-Cell Transcriptomic Profiling Reveals Immunometabolic Reprogramming and Cell-Cell Communication in the Tumor Microenvironment of Human Hepatocellular Carcinoma

1. Introduction

Hepatocellular carcinoma (HCC) is the most prevalent primary liver malignancy and a leading cause of cancer-related mortality worldwide [1]. It arises predominantly in the context of chronic liver disease, with hepatitis B and C virus infections, alcoholic liver disease, and non-alcoholic steatohepatitis as the principal etiological drivers [2,3]. Late diagnosis, high recurrence rates after resection or ablation, and limited efficacy of systemic therapies in advanced disease contribute to persistently poor outcomes [4,5]. Despite improvements in diagnostic imaging, the cellular architecture of the HCC tumor microenvironment (TME) and its contribution to disease progression remain incompletely characterized [2,5].

The liver comprises a structured mixture of parenchymal and non-parenchymal cells, including hepatocytes, sinusoidal endothelial cells, Kupffer cells, hepatic stellate cells, and diverse immune subsets [6]. Chronic inflammation and fibrotic remodeling disrupt this tightly regulated environment, shifting cellular composition, immune activity, and metabolic function toward a permissive oncogenic context [7,8]. Tumor progression is not driven solely by malignant hepatocytes but is profoundly shaped by surrounding stromal and immune populations that can suppress or facilitate growth [6,7].

The HCC TME is defined by three coordinated processes: immunosuppression, neoangiogenesis, and metabolic adaptation. Tumor-associated macrophages (TAMs) blunt cytotoxic immunity by secreting anti-inflammatory mediators and inhibiting lymphocyte activity. Tumor endothelial cells (TECs) drive vascular remodeling while attenuating immune surveillance at the sinusoidal interface. Malignant hepatocytes reprogram their metabolic state—enhancing glycolysis, increasing mitochondrial respiration, and acquiring resistance to oxidative stress—to sustain proliferation under nutrient-deprived conditions [9]. The extent to which these processes are mechanistically coupled through cell-cell communication, rather than occurring in parallel, remains an open question.

Single-cell RNA sequencing (scRNA-seq) in HCC has largely been confined to descriptive studies of cellular composition or lineage markers, with limited efforts to integrate immunological, metabolic, and communication dimensions at single-cell resolution [10]. Characterizing these axes simultaneously is essential for identifying the cooperative mechanisms by which tumor cells co-opt their microenvironment and for uncovering rational therapeutic targets. Ligand–receptor inference, in particular, can reveal how malignant cells establish instructive circuits that propagate immunosuppressive, angiogenic, and metabolic programs across compartments.

Here, we integrate scRNA-seq data from tumor and healthy human liver to simultaneously characterize transcriptional states, metabolic programs, and intercellular communication networks in the HCC TME. Applying probabilistic cell annotation, pseudo-bulk differential expression, pathway enrichment, and consensus-based ligand–receptor inference, we delineate the cooperative immunometabolic modules that sustain malignant progression and identify the communication architecture that distinguishes the tumor ecosystem from homeostatic liver tissue.

2. Materials and Methods

2.1. Data Sources and Processing

To explore the immunometabolic reprogramming and cell-cell communication within the hepatocellular carcinoma (HCC) microenvironment, we utilized publicly available single-cell RNA sequencing (scRNA-seq) datasets from the Gene Expression Omnibus (GEO) database. Specifically, we retrieved tumor-derived scRNA-seq data from [11,12], which include samples from patients diagnosed with advanced-stage HCC. For healthy liver comparisons, we incorporated datasets from [13,14,15,16] and again from Gan et al. (2024), which together provided a robust baseline of cellular profiles from 10 non-malignant liver tissues. All datasets were generated using high-throughput droplet-based sequencing platforms (10x Genomics Chromium), ensuring consistency in data type and quality across studies.

The combined tumor dataset consisted of 47,601 cells derived from multiple regions of primary HCC tumors, while the healthy liver dataset comprised 119,154 cells sampled from histologically confirmed non-cancerous tissues. The selection of datasets emphasized both spatial diversity and cellular representation to minimize biases associated with sampling location or sequencing depth. Raw count matrices and metadata were downloaded in h5ad or Matrix Market formats and loaded into an AnnData object using the Scanpy library (v1.9.3) [17] for preprocessing and downstream analysis.

Upon loading the data, all gene identifiers were mapped to a standardized nomenclature (HGNC symbols), and mitochondrial genes were flagged based on the presence of the "MT-" prefix. To ensure uniformity, datasets were filtered to retain only protein-coding genes, and ribosomal genes were excluded from downstream normalization to avoid artifacts associated with highly expressed housekeeping transcripts. Metadata from original publications, including patient ID, diagnostic condition, and sequencing batch, were integrated into the AnnData object to allow for batch-aware normalization and stratified analysis.

To facilitate cross-sample integration, we applied batch correction and joint embedding using the scVI algorithm [18], allowing us to align cells across donors and conditions while preserving biological variability. Batch effects, which often arise due to technical differences in tissue dissociation, library preparation, or sequencing, were evaluated and visualized using principal component analysis (PCA) and uniform manifold approximation and projection (UMAP) before and after correction (Figure 1).

Quality control filtering was conducted to eliminate low-quality cells and potential sequencing artifacts. Cells were retained if they met the following thresholds: (1) expression of > 200 and <7,500 unique genes per cell, (2) <25% mitochondrial gene content, and (3) <10% ribosomal gene content. Outlier detection was additionally performed using median absolute deviation (MAD) statistics across multiple QC parameters, and cells falling beyond ±3–5 MADs from the median were excluded (Figure 2).

Following quality control, we performed library size normalization to correct for differences in sequencing depth. Gene counts for each cell were normalized to a total of 10,000 transcripts and log-transformed using a pseudocount of 1. Highly variable genes (HVGs) were selected using Scanpy’s built-in dispersion-based method, identifying the top 2,000 genes with the greatest normalized variance across cells. These HVGs were subsequently used for PCA (n=50 components), which formed the basis for neighborhood graph construction, UMAP projection, and downstream clustering.

In total, the integrated and preprocessed dataset contained 159,925 high-quality single-cell transcriptomes, spanning a diverse array of cell types, including parenchymal hepatocytes, immune cell populations, endothelial cells, fibroblasts, and cholangiocytes. This high-resolution cellular atlas provided the foundation for subsequent annotation, differential expression, pathway enrichment, and ligand-receptor interaction analyses aimed at unraveling the molecular architecture of the HCC tumor microenvironment.

2.1.1. Doublet Removal

Doublet detection and removal represent a critical preprocessing step in single-cell RNA sequencing (scRNA-seq) data analysis. During the encapsulation of single cells into microdroplets—particularly in high-throughput droplet-based platforms like 10x Genomics Chromium—there is a non-negligible probability that two or more cells are co-encapsulated into the same droplet. These multiplets, or doublets, result in hybrid transcriptional profiles that do not correspond to any true biological cell type. If left unaddressed, doublets can significantly skew downstream analyses, such as cell type identification, differential gene expression, and trajectory inference, potentially introducing false-positive cell clusters and misleading biological conclusions.

To mitigate this artifact, we employed the Solo algorithm, a state-of-the-art, deep learning-based tool for doublet detection specifically designed for scRNA-seq data [19]. Solo is built upon scVI (single-cell Variational Inference), a variational autoencoder framework that learns a nonlinear latent representation of the gene expression space. By leveraging this latent space, Solo simulates artificial doublets and trains a classifier to distinguish real cells from doublets based on complex expression patterns. Unlike traditional methods such as Scrublet or DoubletFinder, Solo uses a probabilistic approach that adapts to the nonlinear and high-dimensional structure of real scRNA-seq data, yielding higher accuracy, particularly in complex tissues with broad cellular heterogeneity—such as the HCC tumor microenvironment.

The doublet removal workflow began with loading raw count matrices into an AnnData object and performing initial quality control to exclude low-quality or outlier cells that might interfere with doublet modeling. We filtered out genes with extremely low expression, retaining only those present in at least 10 cells. Using the default parameters, we trained the Solo model for 400 epochs using a GPU-accelerated environment with CUDA-enabled support to optimize computational efficiency. This training allowed the model to converge on a well-structured latent space and to robustly identify doublets based on their deviation from biologically plausible gene expression profiles.

The detection phase involved simulating synthetic doublets by averaging the expression profiles of random cell pairs and introducing them into the dataset. The model then classified both real and synthetic cells, assigning a posterior probability of being a doublet to each observed cell. To ensure high-confidence retention, we kept only cells with a posterior singlet probability of ≥0.90, thereby restricting downstream analyses to cells most reliably classified as singlets. This threshold was chosen based on the distribution of scores and cross-validation with known cell-type boundaries, aiming to balance sensitivity (detecting true doublets) and specificity (retaining real single cells).

To evaluate the performance of Solo in our datasets, we examined the doublet score distributions across different cell types and donors. We also assessed doublet localization in UMAP embeddings: genuine doublets often localize between two distinct cell populations in reduced dimensional space. The removal of these putative doublets improved the sharpness of cluster boundaries and enhanced the interpretability of downstream cell-type classification and trajectory modeling.

In total, Solo identified and filtered out approximately 4.1% of cells in the HCC dataset and 3.7% in the healthy liver dataset. These proportions are consistent with the expected doublet rates for 10x Genomics experiments at the observed cell loading densities. The resulting clean datasets provided a more accurate and biologically faithful representation of the cellular composition in both tumor and non-tumor liver tissues, serving as a robust foundation for subsequent analyses such as differential expression, gene regulatory network inference, and intercellular communication modeling.

2.1.2. Cell Type Annotation

Accurate cell type annotation is a cornerstone of single-cell RNA sequencing (scRNA-seq) studies, as it enables researchers to decipher the cellular composition, heterogeneity, and functional organization of tissues under physiological and pathological conditions. In the context of hepatocellular carcinoma (HCC), where both the tumor and its microenvironment comprise a diverse mixture of epithelial, immune, endothelial, and stromal cells, high-resolution annotation is essential to disentangle tumor-specific signatures from surrounding non-malignant cell types and to identify rare or transient cell states that may be critical for disease progression or therapeutic resistance [9].

To achieve robust and context-specific annotation of the 0,000 cells in our integrated dataset (tumor and healthy liver combined), we implemented a multi-tiered strategy that combined automated machine learning classifiers with probabilistic modeling and manual expert validation. This ensemble approach mitigates common limitations of single-method annotation pipelines, such as overreliance on reference datasets or inability to resolve ambiguous or novel cell states.

Automated Annotation Using CellTypist

As a first step, we applied CellTypist (https://www.celltypist.org/), a supervised machine learning-based classifier tailored for scRNA-seq data. CellTypist uses logistic regression or gradient boosting algorithms (XGBoost) trained on manually curated reference atlases to assign cell-type labels based on transcriptomic profiles. For this study, we constructed and used two classes of CellTypist models:

1.

Tumor-specific models were trained using annotated scRNA-seq datasets from the NCI-CLARITY consortium [20], the Multi-Regional HCC Atlas [21], and the Sequential HCC Atlas [22]. These models capture both canonical and disease-specific cell types, such as malignant hepatocytes, cholangiocarcinoma-like cells, tumor-associated macrophages (TAMs), cancer-associated fibroblasts (CAFs), and tumor endothelial cells (TECs).

2.

Healthy liver models were derived from public datasets [23] and the CellxGene liver cell atlas, containing over 160,000 transcriptomes from non-malignant human liver tissues. These models were optimized to detect healthy hepatocytes, Kupffer cells, hepatic stellate cells, cholangiocytes, sinusoidal endothelial cells, and resident innate lymphoid cells.

CellTypist was run with default preprocessing parameters and included confidence scoring for each prediction. Labels with low confidence (<0.5 probability) were flagged for additional inspection. This dual-model strategy allowed us to compare label concordance across tumor and non-tumor conditions and to identify condition-specific expression shifts within shared cell types (e.g., naïve vs. exhausted T cells, resting vs. activated stellate cells).

Label Refinement Using scANVI

To further enhance annotation robustness and to resolve uncertain or transitional cell states, we applied scANVI (single-cell Annotation via Variational Inference), a semi-supervised extension of the scVI probabilistic framework [24]. Unlike deterministic classifiers, scANVI models cell identity as a distribution over possible types, which is particularly advantageous for tissues undergoing dynamic changes such as inflammation, regeneration, or malignant transformation.

The scANVI model was first trained in an unsupervised fashion on the full gene expression matrix using highly variable genes, thereby learning a latent representation that captures the biological structure of the data. Then, using high-confidence CellTypist annotations as anchor labels, scANVI inferred probabilistic labels for all cells—especially for those in ambiguous clusters, at tumor–non-tumor interfaces, or representing rare subtypes. This label transfer approach allowed for cross-condition generalization, enabling meaningful comparisons between the tumor and healthy compartments Figure 3.

2.2. Clustering Optimization with Adjusted Rand Index (ARI)

We optimized the clustering resolution of our single-cell datasets to avoid both under- and over-segmentation of cell populations. To this end, we used our scVI model, which served as input for pseudoclustering. A grid search over multiple neighborhood sizes (k), distance metrics, and Leiden resolutions was conducted, and for each parameter set, clustering results were compared to preliminary annotations obtained from scANVI using the Adjusted Rand Index (ARI). The ARI provides a measure of concordance between two partitions, correcting for chance agreement, and is widely applied in cluster validation. Silhouette scores were additionally computed on a random subset of cells to assess cluster cohesion and separation, serving as a secondary criterion in case of ties.

For each dataset, the configuration that maximized ARI (and silhouette score in case of ties) was selected as the optimal resolution, and the resulting cluster assignments were stored as a new column in Anndata objects. To further assess the robustness of this solution, ARI values were calculated separately per batch and per patient sample, and cluster purity scores relative to scANVI-derived labels were reported. The optimal clustering was then transferred to the main AnnData object, and each cluster was relabeled according to its majority annotation and manual inspection of marker gene expression, ensuring biologically coherent designations.

This ARI-guided clustering strategy was the basis for defining the granularity of cellular partitions, including cross-validation, differential expression, and ligand–receptor inference. By anchoring clustering to both reference-informed labels and intrinsic transcriptomic structure, it is possible to minimize risks of spurious subdivisions and preserve resolution sufficient to capture relevant cellular heterogeneity.

Manual Curation and Marker-Based Validation

In addition to the utilisation of automatic annotation and ARI-based clustering optimisation, we conducted a comprehensive manual validation phase. For each putative cell cluster, we:

1.

Calculated expression scores for known marker gene sets derived from PanglaoDB, CellMarker, and recent liver-specific single-cell atlases from CELLxGENE.

2.

Applied Wilcoxon rank-sum tests to identify top differentially expressed genes for each cluster and compared them against reference signatures.

3.

Visualized gene expression patterns using dot plots, violin plots, and feature plots to verify consistency with expected cell phenotypes.

4.

Re-labeled or merged clusters as needed, unclear lineage identity, or putative doublet contamination.

Outcome of Annotation Workflow

Major annotated types included malignant hepatocytes, immune subtypes (T cells, B cells, TAMs, dendritic cells, NK cells), hepatic parenchymal and stromal cells (hepatocytes, Kupffer cells), and vascular components (arterial/venous endothelial cells, TECs) (Figure 3C). Several subclusters also emerged, including proliferative malignant cells, and inflammation-associated fibroblasts—each with unique expression profiles and potential biological relevance.

This high-resolution annotation provided the foundation for all subsequent analyses, including cell-type-specific differential expression, metabolic pathway inference, and intercellular communication modeling.

2.3. Differential Gene Expression and Pathway Enrichment

Differential gene expression (DGE) analysis is central to understanding the molecular mechanisms that distinguish malignant and non-malignant cellular states. In the context of single-cell RNA sequencing (scRNA-seq), DGE enables the identification of cell-type-specific transcriptional programs associated with disease processes such as oncogenic transformation, immune modulation, and metabolic reprogramming. However, the inherent sparsity and technical noise in single-cell data demand the use of statistical methods capable of modeling low and overdispersed count distributions with high accuracy.

For this study, we performed DGE analysis using the edgeR package, a well-established tool that implements empirical Bayes estimation within a negative binomial framework [25,26]. Although initially developed for bulk RNA-seq, edgeR has been successfully adapted to scRNA-seq data—particularly in pseudo-bulk configurations, where expression counts from cells of the same type and condition are aggregated to reduce variance and increase statistical power. This approach is advantageous in heterogeneous samples such as hepatocellular carcinoma (HCC), where direct cell-to-cell comparisons may be confounded by compositional shifts and batch effects.

2.3.1. Pseudo-Bulk Aggregation Strategy

To conduct DGE, we aggregated raw count matrices across biological replicates for each annotated cell type within each condition (disease condition celltypes vs. healthy condition celltypes). Aggregation was performed on a per-sample basis, yielding one pseudo-bulk profile per cell type per individual. This approach preserves intra-individual variability and avoids overfitting to donor-specific transcriptional patterns. Library sizes were normalized using the trimmed mean of M-values (TMM) method, and genes with low expression (counts-per-million <1 in >75% of pseudo-bulk samples) were excluded to minimize noise.

Pairwise differential expression was then computed using quasi-likelihood F-tests, controlling for false discovery rate (FDR) using the Benjamini-Hochberg method. Genes with an adjusted p-value < 0.05 and an absolute ${log}_2$ fold change >1 were considered significantly differentially expressed. Resulting gene lists were visualized with bar plots and heatmaps and subsequently used for functional enrichment analyses (Supplementary Figure S1 and Figure 4).

2.3.2. Pathway and Functional Enrichment Analysis

To contextualize transcriptional differences in terms of biological function, we conducted Gene Set Enrichment Analysis (GSEA) and Overrepresentation Analysis (ORA) using GO Biological Processes and Molecular Functions.

GSEA was performed using the fgsea package, which ranks all expressed genes by their test statistic and evaluates whether gene sets are overrepresented at the extremes of the ranking. ORA, on the other hand, was applied to the subset of statistically significant differentially expressed genes using the clusterProfiler package, with hypergeometric testing and FDR correction.

This dual enrichment strategy provided complementary insights: GSEA captured subtle coordinated shifts in expression (even without individual gene significance), while ORA highlighted sharply deregulated biological programs. Enrichment results were reported as normalized enrichment scores (NES), adjusted p-values (FDR), and leading-edge subsets, which identify the core drivers of each pathway signal.

2.3.3. Validation and Visualization

To validate the robustness of our DGE findings, we performed sensitivity analyses including alternative normalization strategies (DESeq2’s median-of-ratios method), and subsampling to assess stability across different donor subsets. Additionally, UMAP embeddings were overlaid with expression gradients of top DEGs to confirm consistency with cluster-level annotation. Where applicable, we cross-referenced DEGs with known marker genes and pathway databases such as CellMarker and PanglaoDB.

Enrichment results were visualized using dot plots, network diagrams (via EnrichmentMap), and chord diagrams to depict relationships between pathways and cell types. These visual summaries enabled interpretation of complex transcriptional programs at a systems level, helping to identify candidate therapeutic targets or diagnostic markers associated with malignant transformation or immune dysregulation.

2.4. Consensus-Based Inference of Intercellular Communication

Intercellular communication within the hepatocellular carcinoma (HCC) microenvironment was inferred using LIANA+, adopting an all-vs-all design to capture the signalling relationships among annotated cell populations. The analysis combined four complementary tools: CellChat, CellPhoneDB, SingleCellSignalR, and NATMI, within a unified rank aggregation scheme. Each method estimates interaction strength and specificity through different statistical models and gene expression thresholds, allowing a more reliable and biologically grounded evaluation of ligand–receptor activity across the tissue.

The workflow consisted of three main steps: (i) independent inference of ligand–receptor interactions by each method, (ii) assessment of overlap between methods using the Jaccard similarity index, and (iii) aggregation of ranked outputs into a consensus model. The Jaccard index quantified the proportion of shared top-ranked interactions between methods relative to their combined total, with higher values indicating stronger agreement. This step ensured that the consensus relied on reproducible interaction patterns rather than the bias of a single algorithm.

Consensus scores were then derived through a rank aggregation procedure that combined magnitude and specificity across all methods. Magnitude represented the relative interaction strength inferred from ligand and receptor expression levels, while specificity captured the restriction of each signal to particular sender–receiver pairs. The final score prioritized interactions that were strong, specific, and consistently identified across methods, looking for a robust representation of communication in the dataset.

To minimise false positives, an expression filter was applied, retaining only ligand–receptor pairs where at least 10% of cells in both sender and receiver populations expressed the corresponding subunits. For complexes with multiple subunits, the lowest expression proportion across all components was used to exclude partial or incomplete interactions. This is relevant in liver tissue, where incomplete receptor expression may lead to artefactual predictions.

After aggregation, interaction matrices were generated to quantify the number and strength of connections between sender–receiver pairs. These matrices summarised dominant communication routes and recurrent signalling relationships within the tumor microenvironment.

2.4.1. Intercellular Communication Among Key Hepatic Populations

To refine the consensus analysis and highlight functionally relevant signalling events, we performed a focused evaluation restricted to the main hepatic and tumor -associated populations: hepatocytes, malignant hepatocytes, tumor -associated macrophages (TAMs), and tumor endothelial cells (TECs). These groups were selected for their central roles in tumor metabolism, immune modulation, and vascular remodelling. The analysis used the same consensus framework described above, but limited the comparison space to interactions involving these populations as either signal senders or receivers.

The dataset was divided into defined sender–receiver groups in which hepatocytes, malignant cells, TAMs, and TECs acted as primary signal sources, while all other annotated populations were retained as potential targets. This approach allowed us the examination of how these key cell types communicate within and beyond their own compartments, to clarify their contribution to the tumor microenvironment.

Each subset was analysed under the same multi-method consensus procedure integrating CellChat, CellPhoneDB, SingleCellSignalR, and NATMI. We applied expression filters to retain only interactions with at least 10% of cells expressing both ligand and receptor components, preventing the inclusion of incomplete or artefactual signals. Magnitude and specificity scores from each method were then aggregated into a consensus rank, prioritizing interactions that were strong, distinct, and reproducibly detected across tools.

2.5. Cross-Validation of Healthy and Diseased Datasets

To assess the robustness and generalization of our automated annotation framework, we implemented cross-validation strategies for both healthy and diseased liver datasets. Given the differences in data availability between conditions, the design of the validation scheme was adapted accordingly, while maintaining a common analytical backbone to ensure comparability.

For healthy liver donors, only a single reference dataset (Tabula Sapiens, Healthy liver Dataset) was available. We applied a stratified K-fold cross-validation. The dataset was randomly partitioned into five folds while preserving the proportional representation of all annotated cell types. In each iteration, four folds (80% of the cells) were used for model training and the remaining fold (20%) served as the test set. This stratification guaranteed that both common and rare cell populations were represented across folds and allowed us to evaluate the classifier’s ability to recover the full spectrum of liver-resident cell types under resampling.

For hepatocellular carcinoma (HCC) cohorts, which included multiple patient-derived samples, we implemented a stratified group K-fold cross-validation design. Here, folds were defined at the patient sample level, ensuring that cells from a given donor were never simultaneously present in both training and testing sets. This design prevented information leakage and enabled evaluation under conditions of inter-patient variability. Stratification by cell type was combined with grouping by patient, resulting in balanced folds that respected both cellular composition and biological independence. To avoid instability in classification metrics, cell types with very low counts or represented in fewer than three patient samples were excluded from fold-level testing but retained in the training pool.

In both healthy and diseased settings, model training use scVI and its semi-supervised extension scANVI, implemented in the scvi-tools framework. Within each training fold, 2000 highly variable genes were selected using the Seurat v3 method, with batch-aware correction applied. The scVI model was trained for 200 epochs to learn a latent representation of the data, after which scANVI was initialized with these weights and optimized for an additional 20 epochs. The “Unknown” category was designated as the unlabeled class, allowing the model to accommodate cell types present in the test set but absent from training. Predictions on the held out fold consisted of discrete labels and soft probability distributions across candidate classes. Posterior probabilities were summarized as confidence scores, and cells with maximum probability below 0.5 were flagged as uncertain predictions.

Additionally, beyond fold-level evaluation, we implemented a consensus framework to integrate results across folds within each condition. Standard deviations across folds were reported to capture variability, and consensus recommendations were generated to flag classes with unstable performance or systematic misclassification. This consensus analysis allowed to overview the model behavior, and identify consistently well annotated cell types and those requiring further scrutiny.

All the results, including trained models, prediction tables, and fold-level metrics, were stored in standardized formats (.h5ad, .csv, and .png files can be found in the cross_validation_analysis folder in the Supplementary Materials) to ensure reproducibility and allow for external reanalysis. This approach provided evidence that our annotation strategy is not only internally consistent within a single dataset but also robust to inter-individual variability.

3. Results

3.1. Cellular Remodeling of the Tumor Microenvironment

Healthy liver samples displayed the expected parenchymal-dominant composition, with hepatocytes as the principal population supported by macrophages, endothelial cells, and lymphoid subsets (Figure 3) [14,27]. HCC samples exhibited a marked shift in cellular representation: malignant epithelial cells expanded substantially, accompanied by proportional increases in tumor-associated macrophages (TAMs) and tumor endothelial cells (TECs) (Supplementary Figure S1). This compositional remodeling is consistent with the progressive displacement of homeostatic hepatic architecture by tumor-supportive stromal and immune compartments [27,28,29].

3.2. Immunometabolic Reprogramming of Malignant Hepatocytes

Pseudo-bulk differential expression analysis identified extensive transcriptional divergence between malignant and non-malignant hepatocytes (Figure 4; Supplementary Figures S2 and S3) [11,30]. Among the most strongly upregulated transcripts were long non-coding RNAs including DUXAP8, HELLPAR, and LINC00470, together with PRRT4 and GALNT9. These transcripts have been implicated in oncogenic transcriptional programs and cell proliferation in hepatic malignancy [31,32]. Concurrent downregulation of metallothionein family members (MT1G, MT1F, MT1M, MT1E, MT1A and the iron-regulatory peptide HAMP indicates disruption of metal homeostasis and cytoprotective responses, hallmarks of hepatocyte dedifferentiation [33]. Suppression of immune-associated genes including MARCO, LILRA5, RAC2, and ADGRE5 reflects a coordinated reduction of the innate immune sensing capacity that characterizes functional hepatocytes.

Functional enrichment of downregulated genes revealed selective repression of ribosomal biogenesis and cytoplasmic translation programs, encompassing terms such as cytosolic ribosomal subunit, structural constituent of ribosome, and cytoplasmic translation. This transcriptional pattern indicates a fundamental restructuring of the biosynthetic apparatus in malignant hepatocytes, consistent with the metabolic reallocation from protein quality control toward anabolic and proliferative outputs [33,34]. Downregulation of cell substrate junction components further supports disruption of hepatocyte-matrix anchoring, a prerequisite for invasive behavior.

Upregulated gene sets produced fewer enriched pathways, concentrated in regulation of cell activation and adhesion-related processes. This asymmetry—broad repression of differentiation-associated programs paired with focal induction of proliferative modules—reflects the transcriptional logic of malignant transformation: dismantling hepatocyte identity while enabling context-independent growth [31,32]. Reduced expression of MHC class I loading components, embedded within the broader downregulation of immune-related genes, further links metabolic reprogramming to diminished antigen presentation capacity, suggesting that metabolic and immune evasive remodeling are coupled in malignant hepatocytes [35].

3.3. Pro-Tumorigenic Polarization of Tumor-Associated Macrophages

TAMs exhibited a transcriptional program diverging substantially from macrophages in non-tumoral liver tissue (Figure 4; Supplementary Figure S4) [36]. Upregulated genes included complement components C1QB and C1QC, the phospholipase PLA2G7, the interferon-stimulated gene ISG15, and the scavenger receptor-associated transcript TREM2, together with FCGR2B and MS4A7. This expression profile delineates a macrophage state specialized for lipid sensing, complement-mediated opsonization, and immunoregulatory signaling [37,38]. Downregulated transcripts included MARCO, a pattern recognition receptor associated with homeostatic immune surveillance, and matricellular genes POSTN and CTGF, indicating loss of tissue-maintenance programs and extracellular matrix homeostatic functions.

Functional enrichment of TAM downregulated genes identified suppression of mitochondrial respiratory programs, encompassing aerobic respiration, cellular respiration, and energy derivation by oxidation of organic compounds. Concurrent downregulation of external encapsulating structure and extracellular matrix terms signals disengagement from structural homeostasis. Conversely, upregulated gene sets were enriched for blood vessel morphogenesis, adaptive immune response, lymphocyte differentiation, and T cell activation, revealing an unexpected angiogenic component alongside classical immune-modulatory functions [39].

Together, the TREM2-enriched, lipid-scavenging, complement-expressing TAM profile is consistent with an M2-polarized, immunosuppressive phenotype documented in HCC [36,37]. Repression of mitochondrial oxidative metabolism combined with induction of vascular gene programs indicates that TAMs in HCC assume a dual role: dampening cytotoxic immunity while providing angiocrine and metabolic support to the tumor-supportive niche [38,39].

3.4. Angiogenic and Immune-Silent Reprogramming of Tumor Endothelial Cells

TECs displayed a transcriptional program markedly distinct from endothelial cells in non-tumoral liver (Figure 4; Supplementary Figure S5). Strongly upregulated genes included the VEGF receptor KDR, the angiocrine chemokine CXCL12, the matricellular protein SPARC, and collagen subunits COL4A1 and COL4A2, together with TM4SF18, OIT3, EMCN, and IGFBP7. This gene signature defines an endothelial phenotype committed to vascular remodeling, basement membrane deposition, and angiocrine activation [40]. Downregulated transcripts included the cytotoxic lymphocyte markers GZMH and NKG7, the T cell receptor component CD3E, and the adhesion molecule ADGRF5, indicating silencing of immune-regulatory and leukocyte-interaction programs.

Pathway enrichment of TEC upregulated genes identified activation of collagen-containing extracellular matrix, external encapsulating structure, and cell surface receptor signaling programs, consistent with structural reinforcement of the tumor vasculature and enhanced matrix-mediated signaling. Downregulated programs included organellar ribosome, cytosolic ribosome, and cytoplasmic translation, mirroring the translational repression observed in malignant hepatocytes and suggesting a shared biosynthetic remodeling across tumor-associated parenchymal and stromal compartments.

The TEC expression program therefore reflects an angiocrine switch: progressive acquisition of pro-angiogenic, matrix-remodeling capacity coupled with suppression of immune-stimulatory outputs. This configuration supports tumor perfusion and vascular expansion while restricting leukocyte adhesion and antigen sampling at the sinusoidal interface, generating immune-exclusionary vascular niches [40].

3.5. Intercellular Communication Networks

3.5.1. Consensus Robustness and Network Stability

Intercellular communication was inferred using the LIANA+ consensus framework, integrating CellChat, CellPhoneDB, SingleCellSignalR, and NATMI. The mean Jaccard index across methods increased progressively, reaching a stable plateau at $K^{*}\approx 45$ interactions for HCC donors and $K^{*}\approx 30$ for healthy donors (Jaccard > 0.7), confirming convergence among independent algorithms on a reproducible set of high-confidence interactions (Supplementary Figure S7) [41]. All downstream analyses were restricted to interactions within this stable consensus range.

3.5.2. Healthy Liver Communication

The healthy liver communication network was dominated by MHC-mediated antigen presentation and innate immune surveillance (Figure 5) [29,42]. Hepatocytes constituted the principal signal source, driving MHC class I interactions toward CD8⁺ T cells and NK populations via HLA-A–CD8A, HLA-B–CD8A, HLA-C–CD8A, and HLA-B–CD3D, supplemented by NK receptor engagements through B2M–KLRD1, B2M–KLRC1, and B2M–KLRC2. Additional hepatocyte-derived interactions, including LGALS1–CD69, LGALS1–PTPRC, and RPS19–C5AR1, reflect modulation of leukocyte activation and complement-responsive myeloid cells, consistent with the immune-buffering role of parenchymal cells [28,43].

Macrophages reinforced antigen presentation to helper T cells through HLA-DRA–CD4 and directed class I interactions toward CD8⁺ T cells. The APP–CD74 axis linked macrophages to hepatocytes and monocytes via CD74-associated immune signaling [36], while MIF–CD74_CXCR4 and HMGB1–CXCR4 contributed chemokine and danger-signal communication consistent with sentinel functions within the sinusoidal niche [43]. Endothelial cells supplied complementary class I interactions and CLEC2D–KLRB1 engagement, with TGFB1–CXCR4 and VIM–CD44 embedded in programs linking immune modulation to leukocyte interaction at the vascular interface [44]. Collectively, this network reflects the liver’s established role as an immunologically active yet tolerogenic organ: a structured, hepatocyte-centered communication framework that sustains antigen presentation, NK surveillance, and controlled inflammatory tone under homeostatic conditions.

3.5.3. HCC Communication

The HCC microenvironment displayed a qualitatively distinct communication architecture organized around vascular remodeling, ECM restructuring, inflammatory recruitment, and metabolic adaptation—replacing the immune-surveillance-centered logic of healthy liver (Figure 6) [29].

Malignant hepatocytes engaged endothelial and stromal targets through pro-angiogenic and matrix-dependent interactions: VEGFA–KDR, VTN–KDR, COL18A1–KDR, and FN1–ITGA8_ITGB1 define sustained endothelial activation and fibronectin-integrin signaling supporting invasion and anchorage plasticity [45,46]. Basement membrane remodeling was further indicated by COL4A1–SDC1 and COL4A2–SDC1 interactions, while lipid-associated axes APOE–LSR and APOB–TREM2 coupled tumor metabolic output to myeloid lipid sensing pathways [35].

TAMs amplified inflammatory and chemotactic circuits. High-confidence interactions including CXCL8–SDC1, CXCL2–XCR1, SAA1–FPR1, and SAA1–FPR2 link macrophage-derived acute-phase signals to leukocyte recruitment and endothelial activation [47]. Growth factor axes S100A4–EGFR and S100A4–ERBB3 connect macrophage-derived inflammatory mediators to epithelial proliferative responses, integrating stromal activation with tumor growth programs [46].

TECs consolidated the angiocrine program: VEGFA–KDR, VTN–KDR, COL18A1–KDR, and PLG–FLT1 define a persistent vascular activation signature, complemented by adrenomedullin-dependent endothelial survival and barrier adaptation through ADM–RAMP2 and ADM–RAMP3 [48]. TEC-derived matrix interactions and syndecan-linked chemokine signaling further shaped immune cell trafficking within perivascular niches.

The transition from homeostatic to tumoral communication reflects qualitative rewiring rather than quantitative amplification alone. Whereas healthy liver signaling sustains immunological equilibrium within a zonated metabolic framework, the HCC network displaces antigen-presentation circuits with pro-angiogenic, ECM- integrin, chemokine, and lipid-associated modules [28,29]. Malignant hepatocytes, TAMs, and TECs collectively assume the dominant signaling burden: malignant hepatocytes instruct endothelial and stromal remodeling; TAMs amplify inflammatory and matrix turnover signals; and TECs reinforce vascular expansion. Residual homeostatic hepatocyte–immune exchanges are comparatively attenuated, indicating partial erosion of the antigen-presentation-oriented framework that defines normal liver physiology [9,29].

3.6. Cross-Validation Performance Across Diseased and Healthy Cohorts

Cross-validation evaluated the robustness of the automated annotation framework across both study cohorts, assessing generalization under tumor heterogeneity and physiological conditions (Full tables and figures for all of the analyses in this section can be found in the cross_validation_analysis folder within the Supplementary Materials). The framework demonstrated consistent performance in both settings, with F1 scores exceeding 0.85 across folds and standard deviations below 0.02, confirming reproducible predictions under repeated resampling (See consensus global metrics: DiseasedCV_consensus_global_metrics_mean_sd.csv; HealthyCV_consensus_global_metrics_ mean_sd.csv. Close agreement between macro- and micro-averaged scores indicated balanced classification across both abundant and minority cell populations (global metrics by fold: DiseasedCV_global_ metrics_by_fold.csv; HealthyCV_global_metrics_by_fold.csv).

In HCC samples, hepatocytes and malignant hepatocytes achieved the most stable classification, with per-class F1 scores frequently exceeding 0.95 (Figure S1). T cells demonstrated robust performance, though variability increased when distinguishing proliferative or effector subsets, reflecting the transcriptional overlap inherent to activated T cell states within the tumor. TAMs and endothelial cells reached intermediate precision and recall values with greater fold-to-fold variability (Table S3), consistent with their documented heterogeneity in HCC. Rare populations—including mast cells and proliferative T cells—showed lower classification stability, driven by limited patient-level representation (Figure S2; DiseasedCV_per_class_consensus.csv.

In healthy donors, hepatocytes and cholangiocytes achieved classification metrics approaching unity across folds (Figure S3). Endothelial cells and B cells were robustly classified with slightly higher variability, while T cells showed intermediate performance reflecting transcriptional plasticity even under physiological conditions. Rare populations, including plasmacytoid dendritic cells and mast cells, produced the lowest and most variable scores in both cohorts, attributable to training data scarcity (Figure S4; HealthyCV_per_class_consensus.csv).

Posterior probability distributions corroborated these trends. Major parenchymal and stromal populations exhibited sharply peaked distributions with medians above 0.9 and narrow interquartile ranges (Figures S1 and S3), indicating high prediction confidence. Immune subsets—macrophages and regulatory T cells—displayed broader distributions (Figures S2 and S4), consistent with transcriptional intermediacy and context-dependent states. Low-confidence predictions (posterior probability < 0.5) represented fewer than 5% of cells in HCC cohorts and fewer than 7% in healthy donors, restricted predominantly to ambiguous boundaries between immune subtypes (DiseasedCV_pred_score_by_class_summary.csv; HealthyCV_pred_score_by_class_summary.csv).

Average precision values corroborated population-level trends: hepatocytes, malignant hepatocytes, and T cells consistently exceeded 0.95 in diseased cohorts (Figure S2), while macrophages, endothelial cells, and B cells occupied the 0.80–0.90 range. Healthy donors showed a comparable distribution, with hepatocytes, cholangiocytes, and endothelial cells above 0.90 and immune populations between 0.80 and 0.88. The stability of these distributions across folds confirms that the observed performance reflects reproducible model behavior rather than artifacts of data partitioning: (DiseasedCV_average_precision_consensus.csv; HealthyCV_average_precision_consensus.csv).

4. Discussion

The central finding of this study is that HCC progression operates as a multicellular ecological process: malignant hepatocytes, TAMs, and TECs do not acquire pro-tumorigenic properties independently but converge through structured ligand–receptor circuits into a self-sustaining network that collectively enforces metabolic adaptation, immune exclusion, and vascular remodeling. This systems-level view extends beyond cataloging cell-type alterations to demonstrate that the functional architecture of the TME is an emergent property of intercellular cooperation rather than the sum of isolated transcriptional changes [28,49,50].

A defining feature of this cooperation is the coupling between metabolic reprogramming and communication rewiring. Prior single-cell studies have characterized immunosuppression and angiogenesis as partially independent phenomena in HCC [29,42]. This study provides evidence that these processes are mechanistically linked through the same intercellular signaling circuits: the lipid-associated axes APOE–LSR and APOB–TREM2 bridge hepatocyte metabolic output to myeloid immune polarization; pro-angiogenic signals from malignant hepatocytes (VEGFA–KDR, COL18A1–KDR directly instruct endothelial remodeling; and TAM-derived inflammatory chemokines (CXCL8–SDC1, SAA1–FPR1/2 feed back to amplify vascular activation. Immunometabolic reprogramming in HCC is therefore not a byproduct of malignant growth—it is structurally encoded in the communication network.

Malignant hepatocytes exhibited concurrent activation of aerobic glycolysis and mitochondrial oxidative phosphorylation, reflecting metabolic plasticity that enables proliferation despite fluctuating oxygen and nutrient supply [51,52]. This dual-fuel configuration is consistent with observations that HCC cells upregulate glycolytic modules while retaining mitochondrial capacity for biosynthetic output and metastatic competence [53]. Critically, metabolic reprogramming was transcriptionally coupled to reduced antigen presentation: downregulation of MHC class I loading components and immune-sensing genes in malignant hepatocytes suggests that the metabolic and immune evasive programs are co-regulated, likely through hypoxia/lactate-driven suppression of antigen-presentation machinery [35,54,55]. Metabolic reprogramming thus functions as a primary organizing event in HCC: it generates the tumor-derived metabolic signals that reprogram surrounding immune and stromal compartments, and simultaneously lowers tumor immunogenicity [35,56].

TAMs reinforced this immunometabolic niche through complementary programs. The TREM2-enriched, lipid-scavenging macrophage phenotype observed here is consistent with a PPAR-driven M2 polarization reported in human and murine HCC datasets [36,57,58]. APOE-linked lipoprotein signals from malignant hepatocytes sustain this immunosuppressive macrophage state, coupling tumor lipid output to myeloid reprogramming [58,59]. Downregulation of mitochondrial respiratory programs in TAMs indicates metabolic reconfiguration away from oxidative self-sufficiency, likely supporting the lipid-handling and secretory functions that define their pro-tumorigenic role. Concurrently, SPP1–CD44 signaling drives ECM remodeling and stromal restructuring, reinforcing niches for malignant and cancer stem-like populations [59,60].

TECs underwent a parallel angiocrine switch: upregulation of VEGFA, ANGPT2, and KDR family signaling, combined with attenuation of interferon and TLR pathways, generated a vasculature that simultaneously supports tumor perfusion and restricts immune infiltration [61,62]. The resulting vascular niches impede T cell and dendritic cell access to tumor parenchyma while enabling metabolic exchange and stromal conditioning. Through reinforced VEGFA–KDR and ADM–RAMP2/3 circuits, TECs and TAMs establish reciprocal feedback that stabilizes the malignant ecosystem [28,45].

At the tissue scale, this immunometabolic remodeling is encoded in the rewiring of intercellular communication. In healthy liver, hepatocytes drive a parenchyma-centered network of antigen presentation and immune surveillance, coordinating homeostasis through MHC–TCR interactions and complement-responsive circuits [47]. In HCC, this architecture is replaced by a tumor-driven network in which angiogenic signals (VEGFA–KDR, PLG–FLT1, ECM–integrin interactions (FN1–ITGA8_ITGB1, COL4A1–SDC1, inflammatory chemokines (CXCL8–SDC1, CXCL2–XCR1, and lipid-associated axes (APOB–TREM2, APOE–LSR collectively displace immune-regulatory exchange [28,61]. Matrix proteolysis driven by ADAM and MMP family members further generates invasion-permissive structural niches, while CCL2–CCR2 and SPP1–CD44 couple myeloid recruitment to niche construction [46,58].

These cooperative signaling dependencies carry direct therapeutic implications. The dual glycolytic and oxidative metabolic configuration of malignant hepatocytes confers resilience to single-pathway metabolic inhibition, as cells can reroute fluxes between cytosolic and mitochondrial compartments to sustain ATP and biosynthetic precursors [51,52]. Combinatorial targeting of both fuel pathways may overcome this adaptive resilience. At the stromal level, the interdependence of TREM2-enriched myeloid immunosuppression, lipid-driven TAM reprogramming, and angiocrine TEC programs argues for strategies that disrupt intercellular cooperation rather than targeting individual intracellular nodes [62,63]. Preclinical and early clinical evidence supports this reasoning: combining anti-angiogenic and immune-modulatory therapies reprograms the vasculature and enhances immune infiltration, addressing two interconnected nodes of the malignant circuit simultaneously [28,62].

Defining HCC tumors by the relative dominance of these cooperative modules—metabolic, immunoregulatory, and angiogenic—offers a practical framework for patient stratification and rational therapy selection [64,65]. Quantitative intercellular network scores summarizing ligand–receptor flux, chemokine recruitment intensity, and metabolic cross-feeding between defined populations could serve as predictive indices for therapeutic response, analogous to transcriptional predictors already under clinical evaluation [28].

These interpretations are constrained by inherent limitations of the data. The datasets are cross-sectional and enriched for advanced lesions, precluding direct inference of temporal progression or early oncogenic events [66]. Transcriptomic profiling captures mRNA abundance but not protein activity, post-translational modifications, or metabolic flux, requiring complementary functional validation [67,68]. Spatial resolution is absent: proximity-dependent signaling, perivascular gradients, and hepatic zonation strongly modulate ligand–receptor activity in vivo, and the contextual dependence of inferred interactions cannot be resolved without spatial mapping [42,69]. Multimodal integration—linking spatial transcriptomics with single-cell proteomics, targeted metabolomics, and chromatin accessibility profiling—will be required to validate that inferred circuits correspond to physical proximity, active receptor engagement, and downstream functional outputs [42,69].

The reproducibility of the core patterns described here—coordinated metabolic activation, myeloid immunosuppression, and endothelial angiocrine remodeling—across independent cohorts in multiregional studies supports their biological robustness [49,58,66]. Incorporating module-level scores into existing HCC molecular classifications and developing spatially informed indices of ecosystem interdependence represent concrete next steps toward mechanism-based therapeutic design that treats the tumor as a coordinated multicellular system [42,62,64,70].

5. Conclusions

The convergence of immunometabolic transcriptional reprogramming, coordinated functional polarization of immune and stromal compartments, and qualitative rewiring of intercellular communication circuits reported here collectively support a model in which HCC is sustained by cooperative cellular ecology rather than by cell-autonomous oncogenic programs. This framing has a practical consequence for the field: effective intervention in HCC may require disrupting network-level dependencies, the reciprocal reinforcement among malignant hepatocytes, TAMs, and TECs, rather than targeting isolated molecular nodes within a single compartment. The cooperative immunometabolic modules identified here, and the quantitative intercellular network indices that could be derived from them, provide a concrete basis for mechanism-based patient stratification and the rational design of combinatorial therapies, once validated in larger, spatially resolved, and multi-omic cohorts.