Construction of a Prognostic Model and Subgroup Characteristics Related to Hypoxia-Immune Evasion in T Cells of Hepatocellular Carcinoma Based on Single-Cell and Bulk RNA Analysis

Introduction

Hepatocellular carcinoma (HCC) remains a leading hepatic malignancy in both incidence and mortality on a global scale.^1,2 The pathogenesis and progression of this neoplasm is predominantly driven by infection with hepatitis B or C virus,^3,4 non-alcoholic fatty liver disease, and alcohol-related liver injury.^5,6 Notwithstanding incremental therapeutic progress, the enduring obstacles of frequent recurrence, pronounced chemoresistance, and the biologically complex tumor microenvironment (TME) continue to undermine curative efficacy and prognosis.^7–9 Thus, comprehensive dissection of HCC’s molecular landscape and the discovery of innovative biomarkers are indispensable for timely diagnosis, precision-targeted interventions, and prognostic prediction.

Immune escape is a cardinal feature of HCC orchestrated by multiple mechanisms in the TME.¹⁰ These include the infiltration of immunosuppressive cells (eg, regulatory T cells, M2 macrophages),^11–13 aberrant activation of immune checkpoints such as PD-1/PD-L1,¹⁴ and secretion of inhibitory factors like TGF-β and VEGFA.¹⁵ Hypoxia, a common driver of HCC in the TME, arises from the imbalance between rapid tumor growth and dysfunctional angiogenesis.¹⁶ Beyond reshaping tumor metabolism, hypoxia profoundly impairs immune cell functions.^17,18 Cytotoxic CD8⁺ T cells, for instance, frequently manifest functional impairment when confined to hypoxic tumor territories.¹⁹ Mechanistically, hypoxia-inducible transcription factors, such as HIF-1α, alter T cell metabolism and effector programs, driving exhaustion, senescence, and immune escape.^20,21 In HCC, hypoxia can induce tumor cells to secrete autophagy-related proteins or circRNAs, thereby inhibiting T-cell infiltration and enhancing immunosuppression.^22,23 Collectively, T cells, hypoxia, and immune escape form a complex, interactive network within the HCC microenvironment. Yet, how these elements precisely intertwine remains incompletely understood, and clinically actionable predictive models are still lacking. Therefore, constructing a risk scoring model based on T cells, hypoxia, and immune escape is of great significance for revealing the complexity of the TME in HCC, achieving timely diagnosis, customized treatment, and prognosis prediction. It is worth noting that there is diversity in the functions of T cell subgroups. Different subgroups may play different roles in HCC, such as promoting or inhibiting immune escape. Based on this point, a more refined risk assessment model can be constructed.

To address this gap, we integrated single-cell and bulk RNA-sequencing data to construct a T cell hypoxia and immune escape risk score model for HCC patients. We further delineated the prognostic attributes of distinct T cell subpopulations and their potential associations with hypoxia and immune escape. Our findings not only provide novel biomarkers for HCC prognostication but also offer fresh perspectives on modulating the HCC microenvironment.

Materials and Methods

Data Acquisition and Preprocessing

Gene expression profiles and clinical annotations (age, sex, tumor grade, TNM stage) of HCC samples were downloaded from TCGA (https://portal.gdc.cancer.gov/), including 50 normal and 374 tumor samples. Samples with complete overall survival (OS) information and an OS time of more than 30 days were retained as the training cohort. Additionally, the HCC-related single-cell RNA sequencing datasets GSE149614 and GSE282701 from the GEO (Gene Expression Omnibus) database (https://www.ncbi.nlm.nih.gov/) were downloaded, along with the corresponding bulk RNA sequencing data for subsequent analysis. The GSE149614 dataset contained 10 HCC samples and 8 healthy control samples. The GSE282701 dataset consisted of a total of 12 samples, among which 6 were HBV positive (HBV+) (3 tumor tissues and 3 paired adjacent normal tissues), and 6 were HBV negative (HBV-) (3 tumor tissues and 3 paired adjacent normal tissues). The bulkRNA data included GSE76427-GPL10558, GSE14520-GPL571, and GSE14520-GPL3921, totaling 357 HCC sample data from 3 validation cohorts, and were used as the validation set for survival analysis. Among them, the data from the two platforms of GSE14520 both adopted the RMA standardized expression matrix provided by the original study, and were subjected to background correction, quantile normalization, log2 transformation, and probe aggregation processing. The data of GSE76427-GPL10558 were standardized using the RSN method based on the R package “lumi” in the original study. This study directly uses the above standardized expression data for risk score calculation and validation analysis. The gene expression values used for estimating the risk score in each data set are shown in Supplementary Tables S1–S4. Immune escape-related genes (relevance score > 10) were extracted from GeneCards (https://www.genecards.org/) (Supplementary Table S5). Hypoxia-related genes were taken from a previous study²⁴ (Supplementary Table S6). The two gene lists were merged to generate the hypoxia-immune escape-related genes (HIERGs) (Supplementary Table S7). Table 1 lists the information on samples from GEO. This study was based on national legislation guidelines, item 1 and 2 of Article 32 of the Measures for Ethical Review of Life Science and Medical Research Involving Human Subjects dated February 18, 2023, China.

Table 1 Informaiton of HCC Samples From GEO

Processing of HCC Single-Cell RNA-Seq Data

The preprocessing of the single-cell RNA sequencing dataset GSE149614 was accomplished based on the Seurat R package (version 5.1.0).²⁵ Firstly, the original expression matrix was constructed into a Seurat object. Low-quality cells, ie, cells with genes that can only be detected as being expressed in 3 or fewer cells, or those cells with fewer than 200 genes, were then removed (min.cells = 3, min.features = 200). Filtering thresholds for retained cells were: nFeature_RNA > 200, nCount_RNA < 40,000, percent.mt < 5%. Expression matrices were log-normalized with the “NormalizeData” function, the top 2000 variable genes were identified with the “FindVariableFeatures” function, and data were scaled with the “ScaleData” function. Principal component analysis (PCA) was run with the number of principal component (PC) set to 20. PCA was conducted with the “RunPCA” function, followed by Harmony for batch correction. A k-nearest-neighbor graph was built with the “FindNeighbors” function, and clusters were called at resolution 0.5 with the “FindClusters” function. Then, the “RunTSNE” function was used to visualize cell clusters. Cell-type annotation relied on the “singleR” package²⁶ and canonical markers from CellMarker 2.0 (http://117.50.127.228/CellMarker/). Proportions of cell types in each sample were calculated with the “prop.table” function and compared between tumor and normal tissues using t-tests, whose results were visualized by the “ggplot2” package.²⁷ Differential expression analysis (“FindMarkers” function, logfc.threshold = 0.25, min.pct = 0.1) defined differentially expressed genes in T cells (T-DEGs). Lastly, GO and KEGG enrichment analyses were performed with the “clusterProfiler” package,²⁸ with results visualized by “ggplot2”.

The single-cell RNA sequencing dataset GSE282701 was processed using the same analysis procedure as described above. A Seurat object was constructed using the “Seurat” package (version 5.1.0), and genes expressed in ≤ 3 cells and cells with fewer than 200 detected genes (min.cells = 3, min.features = 200) were filtered during the quality control stage.²⁹ The further screening criteria were set as 200 < nFeature_RNA < 6000 and percent.mt < 20%. Data normalization, high-variable gene screening, and standardization processing were completed through the NormalizeData, FindVariableFeatures (the top 2000 variable genes), and ScaleData functions. Subsequently, principal component analysis (PCA) was performed using RunPCA, and the first 35 PCs were selected for subsequent analysis. Batch effect correction was carried out using the Harmony algorithm. Based on the dimensionality reduction integration, an adjacency graph was constructed using FindNeighbors, cell clustering was performed using FindClusters (resolution = 0.7), and visualization was conducted using RunTSNE. Cell type annotation was also determined based on the SingleR package in combination with the marker genes from the CellMarker 2.0 database.

Identification of Model Genes

T-DEGs (Supplementary Table S8) were intersected with HIERGs to obtain T cell hypoxia-immune-escape-related differentially expressed genes (T-HIERDEGs). Univariate Cox regression was conducted with the “survival” package³⁰ to determine the relationship between T-HIERDEG expression and OS. LASSO regression was conducted on the significant genes with P < 0.05 from the univariate Cox regression using the “glmnet” package.³¹ The optimal cut-off was calculated by the “survminer” package³⁰ to run multivariate Cox regression on the screened results by LASSO, yielding the final prognostic model genes and risk-score formula:

This risk score was calculated using a multivariate Cox proportional hazards regression model. The coefficient of each gene represents the regression coefficient (β value) obtained from the multivariate Cox analysis. The gene expression value represents the standardized transcriptome expression level of each gene. The final risk score for each patient is the sum of the product of each gene’s expression value and its corresponding regression coefficient.

Validation and Independent Prognostic Value of the Risk Model

Patients were dichotomized into high- and low-risk groups based on the median risk score. Kaplan-Meier (K-M) survival curves were drawn by the “survival” package and time-dependent ROC curves were generated by the “timeROC”³² and “survival” packages to evaluate the sensitivity and specificity of the risk model on all the bulk RNA cohort (including TCGA training cohort, GSE76427, and GSE14520). Additionally, univariate and multivariate Cox analyses evaluated the independent prognostic value of the model.

Nomogram Construction

Using the “rms” R package,³³ we constructed a nomogram that combined the risk score and clinical variables of HCC samples to predict the 1-, 3-, and 5-year OS of patients. Calibration curves evaluated the accuracy of the nomogram’s predictions.

Immune Infiltration and Immunotherapy Analyses

The “estimate” package was used to compute immune scores (reflecting the degree of immune cell infiltration), the stromal scores (reflecting the content of stroma), and the ESTIMATE scores (reflecting the combined content of stroma and immune markers). The scores of two groups were compared using Wilcoxon tests. Immune-cell fractions and pathway activities were quantified by ssGSEA and CIBERSORT, shown in boxplots. Tumor immune dysfunction and exclusion (TIDE) scores were calculated to infer immune-escape likelihood in high- and low-risk groups. Immunophenoscores (IPS) were obtained from TCIA to further compare the likelihood of immunotherapy response between the two risk groups.

Mutation Profiling and Drug Sensitivity Analysis

Top 10 mutated genes in the two risk groups were listed using the “maftools”³⁴ package. Mutation profiling of model genes in the two risk groups was visualized with waterfall plots. Tumor mutational burden (TMB) was compared through violin plots and a Wilcoxon test was applied to compare TMB scores of the two risk groups. Drug sensitivity was projected with the “oncoPredict” package³⁵ by leveraging chemotherapy response data from the GDSC2 database (https://www.cancerrxgene.org/). Wilcoxon tests were applied to pinpoint differential drug sensitivities between high- and low-risk groups.

Clustering and Characterization of T Cell Subpopulations

T cells were extracted from the single-cell RNA sequencing datasets GSE149614 and GSE282701, respectively, for secondary dimensionality reduction and re-clustering analysis. DEGs among clusters were identified with the “FindAllMarkers” function (min.pct = 0.25, logfc.threshold = 0.25, |log2FC| > 0.25, adjusted P < 0.05). The marker of each DEG was used to annotate each subpopulation. To evaluate the overall expression activity of model characteristic genes at the single-cell level, we calculated the signature score for each T cell based on the predefined set of model characteristic genes. This was done using the AddModuleScore function from the “Seurat” package. First, the average normalized expression level of the characteristic gene set in each cell was calculated. Then, genes with similar expression levels from the entire gene set were selected as the control gene set and their average expression value was calculated. The difference between the two was used as the characteristic score for this cell to reflect the overall activity level of the characteristic genes relative to the background expression. Subsequently, the Wilcoxon test was used to compare the score differences between the normal group and the tumor patient group in different T cell subpopulations. Pathway alterations were examined by GSVA using the Hallmark gene sets from MSigDB with the “GSVA” package in R.

In the GSE149614 dataset, further cell communication analysis and network visualization were conducted using the “CellChat” R package³⁶ to explore the interactions among T cell subsets. The plotGeneExpression function was utilized to plot the gene expression distribution of signal transduction genes related to signaling pathways.

Cell Culture

Human HCC cell line HepG2 (IM-H038) and human normal liver cell line THLE-2 (IM-H473) were purchased from IMMOCELL (China). HepG2 cells were maintained in 90% DMEM+10% FBS+ penicillin/streptomycin. THLE-2 cells were cultured in their dedicated medium (IM-H473-1). All cells were incubated in a thermostatic incubator at 37 °C with 5% CO₂ and 95% air.

Cell Transfection

sh-LGALS3 and its negative control (sh-NC) vectors were purchased from GenePharma (Shanghai, China). Transfection was performed using Lipofectamine™ 3000 transfection reagent (Thermo Fisher, USA) for transient transfection of sh-LGALS3 and sh-NC, respectively. A portion of cells was harvested to detect LGALS3 expression levels via qRT-PCR and Western blot analysis to validate knockdown efficiency.

Dasatinib Treatment and Experimental Grouping

Dasatinib (HY-10181) was purchased from MedChemExpress (MCE, USA). Dasatinib was dissolved in DMSO (D2650, Sigma-Aldrich, USA) at a concentration of 10 mM to prepare a stock solution (stored at −20°C). To investigate the effect of LGALS3 knockdown on Dasatinib sensitivity, four experimental groups were established: sh-NC + DMSO; sh-LGALS3 + DMSO; sh-NC + Dasatinib; sh-LGALS3 + Dasatinib. All functional assays followed a standardized protocol. After 24 hours of transfection, the original medium was removed, and fresh complete medium containing the respective treatment (DMSO or Dasatinib) was added to all cell groups. This time point was designated as “0 hours” for drug treatment.

qRT-PCR

Total RNA from each cell line was isolated using TRIzol reagent (Invitrogen, USA), and RNA concentration and purity were determined with a NanoDrop 2000 spectrophotometer. cDNA was synthesized using HiScript II Q RT SuperMix (Vazyme, China). qRT-PCR was performed on an ABI 7500 system with ChamQ Universal SYBR qPCR Master Mix (Vazyme). GAPDH served as the internal reference, and all reactions were run in triplicate. Relative gene expression was calculated using the 2^−ΔΔCt method. (Table 2).

Table 2 Primer Sequences

Western Blot (WB)

Total protein was extracted using RIPA lysis buffer (P0013B, Beyotime, China) containing protease and phosphatase inhibitors. Protein concentration was determined using the BCA Protein Concentration Assay Kit (P0009, Beyotime, China) via the BCA method. Equal amounts of protein were subjected to SDS-PAGE electrophoresis and transferred to PVDF membranes. Membranes were blocked with 5% skim milk and incubated overnight at 4°C with primary antibodies against LGALS3 (ab53082; Rabbit; Abcam; UK) and GAPDH (ab8245; Rabbit; Abcam; UK) overnight at 4°C. Membranes were then incubated with the corresponding HRP-labeled secondary antibodies: Goat anti-Rabbit IgG (HRP) (ab6721; Goat; Abcam; UK), and visualized using ECL chemiluminescence.

Cell Counting Kit-8 (CCK-8)

At 0 h post-drug treatment, HepG2 cells from different treatment groups (sh-NC, sh-LGALS3, sh-NC + Dasatinib, and sh-LGALS3 + Dasatinib) were digested and counted. They were then replated at a density of approximately 2 × 10³ cells per well in a 96-well plate and immediately treated with the corresponding agent (DMSO or Dasatinib). Each group included three replicate wells, with one well containing only medium (no cells) serving as a blank control. Timing commenced from the start of seeding and drug addition. At 0, 24, 48, and 72 h post-treatment, CCK-8 reagent (CK04, Dojindo, Japan) (10 μL/well) was added to each well and incubated at 37°C for 2 h. Absorbance values (OD450) were subsequently measured at 450 nm using a microplate reader to assess cell proliferation capacity in each group.

Colony Formation Assay

After 72 h of drug treatment, HepG2 cells from different treatment groups (sh-NC, sh-LGALS3, sh-NC + Dasatinib, and sh-LGALS3 + Dasatinib) were trypsinized, prepared as a single-cell suspension, and counted. They were then seeded at a density of approximately 500 cells per well in a 6-well plate. After seeding, fresh medium containing the corresponding concentration of DMSO or Dasatinib was replaced every 3 days. Cells were cultured at 37°C with 5% CO₂ for approximately 10–14 days. Once visible colonies formed, the medium was discarded. We gently washed them twice with PBS, fixed them with 4% PFA for 15 min, and stained them with 0.1% crystal violet (C0121, Beyotime, China) for 20–30 min. Subsequently, after thorough washing with PBS, air-drying of the colonies was performed at room temperature, which was photographed and counted. The number of colonies per well was recorded to assess cellular clonogenic potential.

Transwell Invasion Assay

Cell invasion assays were conducted using Transwell chambers (Corning, USA, 8 μm pore size). Prior to the experiment, the upper chamber membrane surface was pre-coated with Matrigel (356234, BD Biosciences, USA). 24 h after drug treatment, HepG2 cells from different treatment groups (sh-NC, sh-LGALS3, sh-LGALS3 + Dasatinib) were digested, prepared into a single-cell suspension using FBS-free medium containing the respective treatment (DMSO or Dasatinib), and 100 μL (at a density of 1 × 10⁵ cells/mL) were seeded into the Matrigel-coated upper chamber. The lower chamber was filled with 600 μL of medium containing 10% FBS and the same concentration of DMSO or Dasatinib. After 24 h of incubation at 37°C and 5% CO₂, non-invading cells were removed. Cells that invaded the sub-membrane surface were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet solution. Three random fields per well were counted to determine the number of cells that had traversed the membrane.

Apoptosis Analysis

48 h after drug treatment, cells from different treatment groups (sh-NC, sh-LGALS3, sh-NC + Dasatinib, and sh-LGALS3 + Dasatinib) were collected. Following PBS washing, cells were stained and analyzed by flow cytometry according to the Annexin V-FITC/PI dual staining kit (560931, BD Biosciences, USA) protocol. Specifically, 500 μL cell suspension was mixed with 5 μL Annexin V-FITC and 5 μL PI, then incubated in the dark for 15 min. Samples were subsequently analyzed on a BD FACSCanto™ II flow cytometer (Agilent, USA), with data analyzed using FlowJo software. The proportions of early and late apoptotic cells were used to assess apoptosis levels in each group.

Statistical Analysis

Data analysis was performed using R software (version 4.1.2), GraphPad Prism 10.1.2, and Microsoft Excel (version 12.1.0.24031). The data were presented as mean ± standard deviation (SD). The differences between groups were analyzed using one-way analysis of variance (ANOVA) or t-test. P < 0.05 was considered statistically significant.

Results

Identification of HCC Cell Subpopulations

After stringent quality control of scRNA-seq data from 10 HCC samples and 8 normal control samples (GSE149614), 46,537 high-quality cells were retained for downstream analysis (Figure 1A and B). Variance analysis of the core cells’ genes identified 2000 highly variable genes (Figure 1C). Top 20 PCs with P < 0.05 were generated from PCA for subsequent analysis (Figure 1D). t-SNE visualization resolved 24 distinct clusters (Figure 1E). Manual annotation classified these clusters into 12 major cell types: T cells, NK cells, Kupffer cells, hepatocytes, monocytes, endothelial cells, plasma cells, dendritic cells, fibroblasts, B cells, epithelial cells, and minor populations (Figure 1F). Bubble plots revealed canonical marker-gene expression for each cell type (Figure 1G). Proportional analyses showed that all 12 cell types were present across tumors, yet with differential proportions indicating substantial inter-patient heterogeneity (Figure 1H). Notably, T cells were the most abundant population but displayed markedly altered fractions between tumor and normal tissues, suggesting a central regulatory role within the HCC immune microenvironment. We therefore focused subsequent analyses on T cells.

Figure 1 Quality control of HCC scRNA-seq data. (A) Distribution of detected genes (nFeature_RNA), total RNA molecules (nCount_RNA), and mitochondrial percentage (percent.mt). (B) Correlation heatmaps among nCount_RNA, nFeature_RNA, and percent.mt. Data reliability was confirmed by the strong positive correlation between nFeature_RNA and nCount_RNA, and the Independence of both from percent.mt. (C) Variance analysis with Top 10 highly variable genes. (D) PCA JackStraw plot. (E) Sample-wise t-SNE embedding. (F) Cell-type annotation t-SNE plot. (G) Bubble plots of marker-gene expression across clusters. (H) Cell-type proportions in each sample.

In the independent GSE282701 single-cell dataset, we conducted a confirmatory analysis using the same analytical procedure. After dimensionality reduction and clustering, a total of 30 cell clusters were identified (Figure S1A), and they were annotated as 9 major cell types, including T cells, endothelial cells, NK cells, liver cells, monocytes/macrophages, dendritic cells, B cells, fibroblasts, and plasma cells (Figure S1B). The expression patterns of marker genes for each cell type are shown in Figure S1C. In this cohort, T cells were also observed to be the most abundant cell population, which was consistent with the results of the main cohort.

GO and KEGG enrichment of T-DEGs between tumor and normal tissues revealed involvement in biological processes including mitochondrial oxidative phosphorylation, protein stabilization, etc. (Figure 2A); cellular components such as ribosomes, mitochondrial protein-containing complex, mitochondrial inner membranes, vesicle lumen, and secretory granule lumens (Figure 2B); and molecular functions including NADH dehydrogenase activity, cadherin binding, primary active transmembrane transporter activity (Figure 2C). KEGG pathways enriched included apoptosis, antigen processing and presentation, RNA degradation, and Th17 cell differentiation (Figure 2D).

Figure 2 GO and KEGG enrichment analyses of T-DEGs between normal and tumor tissues. (A) Biological processes (BP) enriched by T-DEGs. (B) Cellular components (CC) enriched by T-DEGs. (C) Molecular functions (MF) enriched by T-DEGs. (D) KEGG pathways enriched by T-DEGs.

Development and Validation of a Prognostic Model with 13 Signature Genes

Intersecting T-DEGs with HIERGs yielded 266 T-HIERDEGs (Figure 3A). Univariate Cox regression in the TCGA cohort identified 41 T-HIERDEGs significantly associated with OS (P < 0.01; Supplementary Table S9). Minimizing the cross-validation error, LASSO regression selected 21 candidate genes (Figure 3B and Supplementary Table S10), and multivariate Cox regression finalized a prognostic model with 13 signature genes (Figure 3C). HSP90AA1, CFH, BCL11B, HLA-E, LGALS3, S100A9, and ADA were shared between T cell and immune-escape signatures. TUBA1C, UGP2, ANXA2, UBB, and PSMD1 intersected hypoxia and immune-escape signatures. SPP1 was a common gene among T cell, immune-escape, and hypoxia signatures (Figure 3D). The formula for the risk model was:

Figure 3 Risk model construction. (A) Venn diagram showing the intersection between T-DEGs and HIERGs. (B) Coefficient distribution plot generated for the log(λ) sequence and LASSO coefficient profile from the LASSO Cox regression. (C) Forest plot of the multivariate Cox regression. (D) Venn diagram illustrating overlaps among T cell, immune-escape, and hypoxia signatures. *P < 0.05; **P < 0.01; ***P < 0.001.

Patients were stratified into high- and low-risk groups by the median risk score (Figure 4A). High-risk patients exhibited significantly shorter OS and worse prognosis (Figure 4A and B). ROC curves demonstrated robust predictive accuracy with AUCs > 0.78 at 1, 3, and 5 years (Figure 4C). A heatmap showed that SPP1, LGALS3, PSMD1, TUBA1C, S100A9, HSP90AA1, ADA, and ANXA2 were up-regulated in the high-risk group, whereas CFH, HLA-E, UBB, UGP2, and BCL11B were down-regulated (Figure 4D). The results of the box plot analysis showed that, compared with normal tissues, the expression levels of SPP1, LGALS3, HLA-E, UBB, PSMD1, TUBA1C, HSP90AA1, ADA and ANXA2 in tumor tissues were significantly increased (all P < 0.05), while the expression levels of UGP2 and S100A9 were significantly decreased (P < 0.05). Additionally, no significant differences were observed between tumor tissues and normal tissues for CFH and BCL11B (P > 0.05) (Figure 4E). Consistent performance was validated in 3 independent external cohorts (GSE76427-GPL10558, GSE14520-GPL571, and GSE14520-GPL3921) (Figure 5A–I). Overall, the prognostic model built in the study demonstrated relatively excellent efficiency in predicting prognosis in HCC.

Figure 4 Continued.

Figure 4 Internal validation of the risk model. (A) Distribution of risk scores and survival status in the TCGA cohort. (B) K-M survival curves comparing high- and low-risk groups within the TCGA cohort. (C) ROC curve for the TCGA cohort. (D) Expression heatmap of model genes in high- and low-risk groups from the TCGA cohort. (E) Differential expression analysis of model genes between the tumor group and normal group. NS: >0.05; *: P<0.05; **: P<0.01; ***: P<0.001.

Figure 5 External validation of the risk model. (A–C) Distribution of risk scores and survival status in the GSE76427-GPL10558 cohort (A), the GSE14520-GPL571 cohort (B), and the GSE14520-GPL3921 cohort (C). (D–F) K-M survival curves comparing the high- and low-risk groups in the GSE76427-GPL10558 cohort (D), the GSE14520-GPL571 cohort (E), and the GSE14520-GPL3921 cohort (F). (G-I) ROC curves for the GSE76427-GPL10558 cohort (G), the GSE14520-GPL571 cohort (H), and the GSE14520-GPL3921 cohort (I).

Independent Prognostic Value and Nomogram Construction

Univariate and multivariate Cox analyses of risk scores and clinical variables confirmed the risk score as an independent predictor of OS (P < 0.001) (Figure 6A and B). A nomogram integrating risk score and clinical variables was thus established (Figure 6C). According to calibration curves, the nomogram effectively predicted 1-, 3-, and 5-year OS with high accuracy (Figure 6D–F).

Figure 6 Assessment of the model’s independent prognostic value and construction and validation of the nomogram. (A) Univariate Cox regression of clinical variables and risk scores. (B) Multivariate Cox regression of clinical variables and risk scores. (C) Nomogram integrating clinical variables and risk score. (D–F) Calibration curves for 1-year (D), 3-year (E), and 5-year (F) OS prediction.

Immune Landscape and Immunotherapy Potential Between Risk Groups

To systematically evaluate the immune microenvironment characteristics of different risk groups and their potential benefits from immunotherapy, we conducted immune infiltration and immunotherapy sensitivity prediction analyses in the TCGA, GSE14520, and GSE76427 cohorts. In the TCGA cohort, the ESTIMATE algorithm results showed that the immune score, stroma score, and ESTIMATE comprehensive score of the low-risk group were significantly higher than those of the high-risk group (p < 0.05), suggesting a more active immune state (Figure 7A). CIBERSORT demonstrated enrichment of naive B and CD8⁺ T cells in the high-risk group, while M0 and M2 macrophages accumulated in the low-risk group (P < 0.05) (Figure 7B). ssGSEA further showed greater infiltration of 11 immune cell types (eg, B cells, dendritic cells, NK cells) and elevated activity of six immune-related functions (eg, pro-inflammatory response, T cell co-inhibition, interferon response) in the low-risk group (P < 0.05) (Figure 7C and D). These results revealed a more complex immune microenvironment and more active immune activities in the low-risk group. IPS and TIDE analyses indicated that the low-risk group that harbored higher IPS and lower TIDE scores was more likely to benefit from immunotherapy (P < 0.05) (Figure 7E and F).

Figure 7 Immune landscape and immunotherapeutic potential between high- and low-risk groups. (A) Immune-related scores (immune score, stromal score, and ESTIMATE score) of two risk groups. (B) CIBERSORT assessment of the two risk groups’ immune microenvironment. (C–D) ssGSEA assessment of the two risk groups’ immune microenvironment. (E) IPS scores of the two risk groups. (F) TIDE scores of the two risk groups. ns: P > 0.05; *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001.

A consistent trend was observed in the external validation cohort. In the GSE14520 cohort, the ssGSEA analysis revealed that aDCs, macrophages, neutrophils, Th2 cells, and Treg cells had higher infiltration levels in the high-risk group, while iDCs and T helper cells were more abundant in the low-risk group (p < 0.05) (Figure S1A and B). In the GSE76427 cohort, most immune cell types and immune function pathways were significantly enriched in the low-risk group (p < 0.05) (Figure S1C and D).

Further analysis of immune treatment response prediction showed that in the TCGA cohort, the low-risk group had higher IPS scores and antigen presentation scores (AZ), and a lower TIDE score, suggesting that it was more likely to benefit from immunotherapy (p < 0.05) (Figure 7E–F and Figure S1E). In the GSE14520 cohort, the effector cells (EC) and immune checkpoint-related scores (CP) of the low-risk group were higher (p < 0.05) (Figure S1F). In the GSE76427 cohort, the low-risk group also showed higher IPS, AZ, and EC scores (p < 0.05) (Figure S1G). These results consistently supported that the low-risk group had more favorable immune treatment response characteristics in multiple independent cohorts.

Mutation Profiling and Drug Sensitivity Analysis Between Risk Groups

Mutational profiling revealed higher mutation frequencies and elevated TMB in the high-risk group (Figure 8A–C). Waterfall plots showed low-level mutations of PSMD1, HSP90AA1, BCL11B, CFH, UGP2, and ANXA2 in the high-risk group, and PSMD1, HSP90AA1, BCL11B, CFH, UGP2 in the low-risk group (Figure 8D and E). Drug sensitivity predicted that the low-risk group was more sensitive to AZD8055 1059, Dactolisib 1057, NU7441 1038, and R0-3306 1052 (with lower IC₅₀ values), while the high-risk group was more sensitive to GNE-317 1926, WIKI4 1940, Afatinib 1032, Dasatinib 1079, Erlotinib 1168, Osimertinib 1919, and Pictilisib 1058 (P < 0.05) (Figure 9A and B).

Figure 8 Mutation profiles of high- and low-risk groups. (A) Mutation landscape in the high-risk group. (B) Mutation landscape in the low-risk group. (C) TMB comparison. (D) Waterfall plot of model-gene mutations in the high-risk group. (E) Waterfall plot of model-gene mutations in the low-risk group.

Figure 9 Predicted drug sensitivities in high- and low-risk groups. (A) Drugs with higher predicted sensitivity in the high-risk group. (B) Drugs with higher predicted sensitivity in the low-risk group. (C) Prediction of drug sensitivity for LGALS3. The vertical axis represents the IC50 value of the drug, and the horizontal axis represents the expression level of LGALS3.

Given the potential role of LGALS3 in HCC, we further conducted a drug sensitivity prediction analysis for LGALS3. The results of the correlation analysis showed that high expression of LGALS3 was significantly negatively correlated with the IC₅₀ values of the drugs Dasatinib and Pictilisib, indicating that high expression of LGALS3 was more sensitive to these drugs (P < 0.05) (Figure 9C).

Association Between T Cell Subpopulations and the Model

Re-dimensional reduction and clustering annotation analysis of T cell clusters in HCC identified 4 major subpopulations: ANXA1⁺, CXCL13⁺, FCRL5⁺, and TIGIT⁺ T cells (Figure 10A–C). Bubble plots confirmed expression of model genes across all subpopulations, with ANXA2, HSP90AA1, UGP2, UBB, HLA-E, and SPP1 most highly expressed (Figure 10D). Signature scores of the 4 subpopulations differed significantly. The highest signature score was in TIGIT⁺ T cells, followed by ANXA1⁺ and CXCL13⁺ subpopulations, and the lowest was in FCRL5⁺ T cells (P < 0.05) (Figure 10E). These findings suggested that these T cell subpopulations might be closely related to hypoxia and immune escape mechanisms.

Figure 10 Continued.

Figure 10 Relationship between T cell subpopulations and the model. (A) UMAP of T cell subpopulations from tumor samples. (B) Heatmap of differentially expressed genes among tumor T cell subpopulations. (C) Annotated UMAP of tumor T cell subpopulations. (D) Bubble plot of model-gene expression across tumor T cell subpopulations. (E) Violin plots of “AddModuleScore” score based on model genes in tumor T cell subpopulations (HIE_signature indicates model genes; H = hypoxia, IE = immune escape). (F) GSVA of tumor T cell subpopulations based on Hallmark gene sets. (G) Number and strength of cell communication among tumor T cell subpopulations. (H) Pseudotime trajectory of distinct T cell subpopulations. (I) Model-gene expression along the pseudotime trajectories of distinct T cell subpopulations. (J) Heatmap of model-gene expression patterns across T cell subpopulations.

To further explore its biological significance, related genes of each subpopulation were collected and analyzed by GSVA. Results revealed that TIGIT⁺ T cell-related genes were enriched in hypoxia, PI3K-Akt-mTOR, Notch, Hedgehog, Wnt-β-catenin, glycolysis, oxidative phosphorylation, interferon response, and DNA-repair pathways, whereas ANXA1⁺ T cells were enriched in immune-related signaling pathways such as PI3K-Akt-mTOR, Notch, Hedgehog, and Wnt-β-catenin, and additionally showed significant involvement in apoptosis, inflammation, and epithelial-mesenchymal transition. CXCL13⁺ T cell-related genes were mainly associated with signal pathways such as TGF-β signaling, angiogenesis, and interferon-γ response. The FCRL5⁺ T cell subpopulation is related to KRAS signal inhibition (Figure 10F).

Furthermore, CellChat demonstrated strong communication among T cell subpopulations, most notably between TIGIT⁺ and ANXA1⁺ T cells (Figure 10G). Trajectory analysis positioned FCRL5⁺ T cells at the earliest stage of T cell differentiation, with subsequent bifurcation toward TIGIT⁺ and ANXA1⁺ T cells along pseudotime (Figure 10H). Further analysis revealed that the model genes ANXA2, HSP90AA1, LGALS3, TUBA1C, S100A9, SPP1, and TUBA1C exhibited the highest expression levels in the TIGIT⁺ T cells, predominantly during the mid-stage of T cell differentiation, whereas HLA-E, ADA, BCL11B, and CFH were most highly expressed in the ANXA1⁺ T cells, mainly in the mid-to-late stages of T cell differentiation (Figure 10I and J). Notably, the TIGIT⁺ T cell subpopulation was enriched with more risk genes from the model (eg, HSP90AA1, LGALS3, TUBA1C, S100A9, SPP1, TUBA1C), while the ANXA1⁺ T cell subpopulation was enriched with more protective factors from the model (eg, HLA-E, BCL11B). Additionally, combining the signature scores of T cell subpopulations with findings from GSEA showed that the gene set of TIGIT⁺ T cells with the highest signature score was significantly enriched in pathways related to hypoxia, oxidative phosphorylation, and glycolysis, whereas the ANXA1⁺ T cells were downregulated in these pathways and more enriched in pathways related to apoptosis and inflammatory responses. Thus, we speculated that TIGIT⁺ T cells in the mid-differentiation stage might be more adapted to hypoxic environments and promote immune escape by modulating metabolic pathways and immune responses. In contrast, the ANXA1⁺ T cells in the late-differentiation stage might have a weaker role in immune escape but could play other functions in the TME, such as regulating apoptosis and inflammatory responses to influence tumor cell behaviors.

Overall, the model genes we identified might participate in T cell immune escape mechanisms and hypoxic responses in the TME by regulating the T cell differentiation, and functions and interactions of T cell subpopulations, thereby affecting immune escape and tumor progression in HCC.

After conducting the same analysis process on the GSE282701 dataset, two main T cell subpopulations were identified: the ANXA1⁺ T cell subpopulation and the TIGIT⁺ T cell subpopulation (Figure S3A–C). Since the expression levels of CXCL13 and FCRL5 were relatively low in each cluster, they were classified as other T cell populations. The bubble plot results showed that most model genes were expressed in both of the two main T cell subpopulations, among which ANXA2, HSP90AA1, UBB, and HLA-E had relatively higher expression levels (Figure S3D). The feature score analysis indicated significant differences between the two T cell subpopulations (p < 0.05), with the ANXA1+ T cell subpopulation having the highest score (Figure S3E). The GSVA analysis revealed that the TIGIT+ T cell subpopulation was mainly enriched in pathways such as Glycolysis, mTORC1 signaling, PI3K-Akt-mTOR signaling, Oxidative phosphorylation, Reactive oxygen species pathway, Interferon response, and DNA repair; while the ANXA1+ T cell subpopulation was enriched in pathways such as IL2-STAT5 signaling, Complement, P53 pathway, Inflammatory response, Hypoxia, Angiogenesis, and Epithelial-mesenchymal transition (Figure S3F). These results supported the functional differences of different T cell subpopulations in metabolic regulation and immune-related pathway activities from an indirect perspective in the independent dataset.

Experimental Validation of Key Genes in the Model

To validate the in vitro expression characteristics of key genes in the 13-gene prognostic model constructed based on bioinformatics analysis, we performed qRT-PCR detection of mRNA expression levels for model-related genes in the HCC cell line HepG2 and the normal hepatocyte line THLE-2. Results showed that compared with normal hepatocytes THLE-2, genes SPP1, LGALS3, PSMD1, TUBA1C, S100A9, HSP90AA1, ADA, HLA-E, UBB, and ANXA2 were significantly upregulated (P < 0.05), while UGP2 was significantly downregulated (P < 0.05). No significant differences were observed for CFH and BCL11B between tumor and normal groups (P > 0.05) (Figure 11). Overall, these in vitro expression results largely align with the expression trends of model genes in high- and low-risk groups observed in bioinformatics analysis, further supporting the stability and reliability of this 13-gene prognostic model in HCC.

Figure 11 qRT-PCR validation of key genes expression in HCC cells. * P < 0.05.

Functional Investigation of LGALS3

Based on prior bioinformatics analysis, LGALS3 was identified as a key gene in the prognostic model and exhibited significantly elevated expression in HCC. Single-cell transcriptomics further revealed that LGALS3 was expressed across multiple cell populations within the tumor microenvironment, particularly in TIGIT+ T cell subpopulations, where it showed significant association with high model scores and pathways related to hypoxia and metabolic reprogramming. Given that LGALS3 participates in regulating intrinsic malignant behavior of tumor cells while potentially influencing immune escape through tumor-immune interactions, this study selected LGALS3 as a representative risk gene to focus on investigating its functional role in tumor cells.

Based on drug sensitivity prediction results, we further validated the relationship between LGALS3 expression levels and Dasatinib treatment response in an in vitro model.

qRT-PCR and WB analyses confirmed the successful establishment of a HepG2 cell model with LGALS3 knockdown (P < 0.05) (Figure 12A and B). Subsequent in vitro functional assays demonstrated that LGALS3 knockdown significantly inhibited HepG2 cell proliferation and clonogenic capacity, markedly reduced invasive ability, and markedly promoted apoptosis (P < 0.05) (Figure 12C–F). Notably, Dasatinib failed to exhibit additional antitumor effects in LGALS3-knockdown cells (Figure 12C–F), suggesting that the antitumor activity of Dasatinib is partially dependent on LGALS3 expression levels. These findings indicate that LGALS3 plays a key regulatory role in the malignant biological behavior of HCC cells and may influence tumor cell response to Dasatinib treatment.

Figure 12 In vitro validation of LGALS3 function and its association with Dasatinib treatment response. (A and B) qRT-PCR (A) and WB (B) detection of LGALS3 mRNA and protein expression levels in HepG2 cells to assess transfection efficiency. (C-D) CCK-8 assay and colony formation assay evaluated the effects of LGALS3 overexpression and Dasatinib treatment on HepG2 cell proliferation. (E) Transwell assay assessed the impact of LGALS3 overexpression and Dasatinib treatment on HepG2 cell invasion capacity. (F) Flow cytometry analysis of the effects of LGALS3 overexpression and Dasatinib treatment on apoptosis in HepG2 cells. * P < 0.05.

Discussion

HCC ranks among the most prevalent and deadly malignancies globally. However, its significant heterogeneity and the complexity of the TME present substantial hurdles for early detection and effective treatment.³⁷ As a result, the identification of reliable biomarkers that can enhance personalized treatment strategies and improve prognostic precision remains a critical unmet clinical need. In contrast to conventional bulk RNA sequencing, scRNA-seq provides a high-resolution view of the TME at the single-cell level, and allows for the detailed characterization of individual cell phenotypes, functional states, and dynamic changes throughout tumor progression, thus serving as a robust tool for uncovering tumor-specific biomarkers.³⁸ By integrating bulk and single-cell transcriptomic data from HCC patients, we developed a risk model of T cell hypoxia and immune escape, and explored the relationship of the model and T cell differentiation within the TME. Our findings provided novel insights that could potentially improve HCC prognosis and therapy.

First, we identified 12 core cell types from the scRNA-sq data of HCC, among which T cells constituted the largest fraction. We therefore centered our subsequent analyses on T cells and integrated T cell, hypoxia, and immune-escape gene sets to develop a prognostic model with 13 T-HIERDEGs: SPP1, LGALS3, PSMD1, TUBA1C, S100A9, HSP90AA1, ADA, CFH, HLA-E, UBB, UGP2, ANXA2, and BCL11B. The first 7 genes emerged as dominant risk drivers, and in-depth mechanistic exploration was conducted on these genes’ functions and impacts on HCC progression. SPP1 is a multifunctional glycoprotein implicated in tumorigenesis and progression across multiple cancers. In lung cancer, SPP1 correlates with poor prognosis and chemoresistance,³⁹ while in melanoma it fosters proliferation, migration, and invasion.⁴⁰ In this study, the SPP1 gene was not only an intersecting gene in the model but also a common gene among T cell-related genes, hypoxia genes, and immune escape genes. Therefore, it might serve as a key gene linking T cells, hypoxia genes, and immune escape. The specific mechanisms of this gene in HCC have remained incompletely understood. Nevertheless, a recent study on HCC cancer stem cells reported that the colocalization of cancer stem cells and SPP1-type macrophages in hypoxic areas may lead to poor prognosis of HCC.⁴¹ PSMD1 is a key component of the 26S proteasome and plays a significant role in the genesis and progression of HCC. PSMD1 is closely related to copy number variations, TMB, and methylation of HCC cells and significantly affects patient survival.⁴² Further research indicates that HCC patients with high PSMD1 expression have a poor prognosis.⁴³ Additionally, depletion of PSMD1 induces G2/M cell cycle arrest in HCC cells, leading to DNA damage and promoting cancer cell apoptosis, a process independent of p53 status.⁴⁴ TUBA1C is a microtubule component involved in various cancers,⁴⁵ and its high expression in HCC can promote migration and proliferation of HCC cells and predict poor prognosis.⁴⁶ S100A9 is a member of the S100 protein family, mainly expressed in immune cells, and regulates inflammatory responses, proliferation, and migration.⁴⁷ In HCC, expression of S100A9 is induced by HIF-1α. High expression of S100A9 can cause mitochondrial fission and reactive oxygen production, thereby promoting HCC growth and metastasis and resulting in a poor prognosis.⁴⁸ HSP90AA1 is a molecular chaperone protein involved in tumor cell metastasis and drug resistance.^49–51 In HCC, abnormal overexpression of HSP90AA1 is associated with enhanced drug resistance, proliferation, and metastatic ability of tumor cells.^52,53 ADA is involved in adenosine metabolism and regulates immune responses. As an immune suppressive signal, ADA deficiency can cause T cell immunodeficiency.^54–56 Patients with ADA deficiency may develop liver dysfunction or HCC even after successful hematopoietic stem cell transplantation due to immunodeficiency.⁵⁷ LGALS3 (Galectin 3, also known as Gal-3) is a glycoprotein that participates in various processes such as cell-cell interactions, immune regulation, apoptosis, proliferation, and migration.^58–61 In HCC, the high expression of LGALS3 is associated with the malignancy and poor prognosis of HCC, and it promotes bone metastasis of HCC.^62,63 This study found that LGALS3, as a key risk gene in the model, not only shows high expression in the TIGIT⁺ subset in single-cell data and is related to hypoxia and metabolic reprogramming pathways, but also has a direct functional role in the malignant behavior of HCC cells in vitro experiments. Knockdown of LGALS3 significantly inhibits the proliferation, colony formation, and invasion of tumor cells, promotes apoptosis, and affects the drug sensitivity of Dasatinib. This result suggests that LGALS3 may play a core regulatory role in tumor-immune interaction, and its expression level is not only related to prognosis but may also affect treatment response, providing a potential target for the precise treatment of HCC. Compared to other risk genes, the functional validation of LGALS3 provides strong support for the biological rationality of the model and strengthens its key position in the T-cell hypoxia-immune evasion process. Overall, these findings indicate that the prognostic genes identified in this study play important roles in the occurrence, development, and immune escape of HCC. They affect the prognosis, drug resistance, and metastasis of HCC by influencing cell proliferation, migration, immune responses, and TME. Therefore, further research on these genes in the future may provide new ideas and potential targets for the early diagnosis, targeted therapy, and immunotherapy of HCC.

Furthermore, to understand the relationship between T cells and hypoxia and immune escape, in addition to establishing a T cell hypoxia and immune escape-related model, we also compared the signature scores of different T cell subpopulations and explored the expression of model genes. Through our analysis, we uncovered significant disparities in signature scores across various T cell subpopulations, suggesting that these subpopulations exhibit distinct functional characteristics and immune regulatory roles within the hypoxia and immune escape-related model. We then delved deeper by conducting enrichment, trajectory, and cell communication analyses for each T cell subpopulation. Our results revealed robust interactions among the 4 T cell subpopulations. Notably, the interaction between the TIGIT⁺ and ANXA1⁺ T cell subpopulations was particularly close, indicating that these two subpopulations may collaborate within the TME. Further investigation showed that the genes of the TIGIT⁺ and ANXA1⁺ T cell subpopulations, which were in the mid-to-late stages of differentiation, were commonly enriched in several key signaling pathways, including PI3K-Akt-mTOR, Notch, Hedgehog, and Wnt-β-catenin. Enrichment in these pathways has been previously demonstrated to facilitate immune escape in HCC.^64–68 This suggested that both TIGIT⁺ and ANXA1⁺ T cells were associated with strong immune suppression. The TIGIT⁺ T cell subpopulation, which boasted the highest signature score, appeared to play a more significant role in immune escape. In addition to this, pathways related to hypoxia, oxidative phosphorylation, and glycolysis were significantly enriched in this subpopulation, while these pathways were suppressed in the ANXA1⁺ T cell subpopulation. Hypoxia-related genes or signatures can significantly promote immune escape in HCC.^23,69,70 Moreover, hypoxia can induce a metabolic shift in cancer cells from oxidative phosphorylation to aerobic glycolysis, also known as the Warburg effect, which can promote the malignancy of the tumor.⁷¹ Our study also found that the ANXA1⁺ T cell subpopulation in the late stage of differentiation was enriched in pathways related not only to immune escape but also to apoptosis and inflammatory responses. This suggested that, although this subpopulation’s role in immune escape is relatively weaker than that of the TIGIT⁺ T cell subpopulation in the mid-stage of T cell differentiation, it may play a significant role in apoptosis and tissue remodeling. In summary, T cell differentiation and their enrichment in hypoxia and immune escape-related signaling pathways may play a pivotal role in the immune escape, tumor progression, and prognosis of HCC. Abnormal activation of PI3K-Akt-mTOR, Notch, Hedgehog, and Wnt-β-catenin signaling pathways may be crucial mechanisms driving immune escape. In addition, we found that the model genes ANXA2, HSP90AA1, LGALS3, TUBA1C, S100A9, SPP1, and TUBA1C exhibited the highest expression levels in TIGIT⁺ T cells, predominantly during the mid-stage of T cell differentiation. In contrast, HLA-E, ADA, BCL11B, and CFH had the highest expression levels in ANXA1⁺ T cells, mainly in the mid-to-late stages. According to the mentioned findings, TIGIT⁺ and ANXA1⁺ T cell subpopulations had higher signature scores and were significantly enriched in signaling pathways that might regulate immune escape. Therefore, we speculated that T cells in the mid-to-late stages of differentiation are more susceptible to regulation by hypoxia and immune escape-related signals, and that the TIGIT⁺ and ANXA1⁺ T cell subpopulations are likely the primary T cell subpopulations suppressed in the TME.

In summary, the study built a T cell hypoxia-immune-escape prognostic model by integrating bulk RNA and scRNA-sq data, and elucidated the central role of T cells in mediating hypoxia-driven immune escape in HCC. Our research results indicate that the T-cell hypoxia and immune escape characteristics associated with HCC are closely related to prognosis and the immune microenvironment, providing a useful framework for risk stratification and biological interpretation.

However, this study still has several limitations that need to be clarified. Firstly, the prognostic model constructed in this study was mainly trained and validated using retrospective transcriptome data from public databases. Although it showed good robustness in multiple independent datasets, it lacks further validation through prospective, multicenter clinical cohorts. Therefore, the current findings should be interpreted as providing a biology-based prognostic assessment method and insights into hypothesis generation, but not a tool for immediate clinical decision making. Its generalization ability in real clinical scenarios still needs to be further evaluated. Secondly, although we conducted multi-level functional analysis of T-cell subpopulation characteristics, hypoxia status, and immune escape-related pathways using single-cell transcriptome data, the related conclusions were mainly based on bioinformatics inference. The causal relationships at the mechanistic level still need to be verified through more systematic in vitro and in vivo functional experiments. Moreover, the analysis of immune infiltration in this study originated from bioinformatics inference, mainly providing evidence of the bioinformatics association between high-risk and low-risk HCC in the composition of the immune microenvironment. Although we observed significant differences such as naive B cells, CD8+ T cells, and M2-type macrophages between groups, the specific molecular mechanisms and functional interactions of these differences affecting the prognosis of liver cancer still need to be further verified and clarified through in vitro and in vivo experiments. Finally, this study only conducted representative cell function experiments and drug response validations for key genes in the model. The specific biological roles of other characteristic genes in different cell subpopulations still need further in-depth research. In the future, combining prospective clinical samples, multi-omics data integration, and more detailed immune function experiments will help further enhance the clinical translational value and mechanism explanation depth of the model.