Identification of immune related genes involved in the onset and progression of HBV-HCC

This study was conducted following the procedures shown in flow chart (Supplementary Fig. S1). The volcano plot analysis revealed 3,637 DEGs in tumor samples from the TCGA-LIHC dataset and 683 DEGs from the GSE121248 dataset (Fig. 1A, B, Supplementary Table S2, 3). By intersecting differentially expressed genes (DEGs) from both datasets with immunity-related genes (IRGs) from the ImmPort database, we identified 53 genes that were shared across all three gene sets, which further defined as immune-related genes involved in HBV-HCC progression (Fig. 1C, Supplementary Table S4). To determine their biological roles, GO and KEGG pathway enrichment analyses were performed^29,30,31 (Fig. 1D, E). The key IRGs were primarily involved in cytokine-cytokine receptor interactions, chemokine signaling, and positive regulation of the MAPK cascade, which are closely linked to immune responses and tumorigenesis. A PPI network was then generated to further investigate the significance of these IRGs in the HBV-HCC disease network (Fig. 1F).

Identification of immune related genes involved in the onset and progression of HBV-HCC. (A) Volcano plot showing differential genes in TCGA-LIHC. (B) Volcano plot showing differential genes in GSE121248. (C) Intersecting genes of the DEGs of the above two datasets and immune related genes obtained from Immport. (D, E) GO (D) and KEGG (E) enrichment results for intersecting genes^29,30,31. (F) PPI network constructed using 53 HBV-HCC related immune genes.

Construction and verification of a prognostic immune index for HBV-HCC and HBV infection

To identify pivotal genes and simplified immune risk index with high diagnostic value, the 53 intersection genes identified earlier were used as candidate genes. In the TCGA-LIHC training cohort, 101 prediction models were built using the LOOCV framework, and the C-index for each model was calculated across both the training and test datasets. Among these, the best-performing model combined stepwise Cox regression (forward) and RSF, yielding the highest average C-index of 0.706 (Fig. 2A). Fourteen prognostic genes were identified via forward stepwise Cox regression, and an immune risk score was constructed using a Random Survival Forest model (ntree = 1000, nodesize = 5, splitrule=‘logrank’) (Fig. 2B). Using the optimal threshold determined from the TCGA training cohort, the two cohorts were classified into high-risk and low-risk groups (Fig. 2C). Notably, higher immune scores were associated with significantly worse overall survival (OS) in the TCGA-LIHC training cohort (Fig. 2D) and external validation cohorts ICGC (Fig. 2E), GSE14520 (Supplementary Fig. S2A), GSE76427 (Supplementary Fig. S2B), and GSE202069 (Supplementary Fig. S2C). Time-dependent ROC curves further demonstrated the predictive sensitivity and specificity of the immune index for OS in the training cohort (Fig. 3A) and validation cohorts (Fig. 3B, Supplementary Fig. S2D–F). Consistent with this, heatmap revealed the high-risk group had elevated SPP1 expression and advanced tumor stages (Fig. 3C). The immune index also showed strong predictive value for DSS, PFI and DFI, with high-risk patients exhibiting poorer outcomes across these parameters (Fig. 3D-F, p < 0.05). To assess the diagnostic accuracy of the prognostic immune index in predicting the disease, ROC analyses were conducted and the results showed that immune index effectively distinguished between HBV-HCC (Supplementary Fig. S2G, H), HBV (Supplementary Fig. S2I) infection patients and their corresponding controls (AUC > 0.65).

Integrated machine learning framework develops immune risk index based on HBV-HCC related immune genes for HBV-HCC. (A) 101 different combinations of machine learning algorithms and each model’s c-index were calculated for training TCGA-LIHC and test ICGC-LIRI-JP datasets. (B) The number of trees for determining with minimal error and the importance of the 14 most valuable HBV-HCC related immune genes based on the RSF algorithm. (C) The optimal cutoff of training TCGA cohort. (D, E) Kaplan-Meier survival curve of OS between high- and low-risk group in the TCGA-LIHC (D) and ICGC-LIRI-JP datasets (E).

Validation the predictive performance of the immune risk index. (A, B) ROC curves of 1-year, 3-year, and 5-year OS in the training TCGA-LIHC (A) and test ICGC-LIRI-JP datasets (B). (C) Heatmap is presented to display the relationship of immune risk groups and clinical features, as well as the expression of the 14 most valuable HBV-HCC related immune genes in patients with HCC. (D-F) Kaplan-Meier survival curve predicted the survival probability for DSS (D), PFI (E) and DFI (F).

Establishment and validation of a nomogram combined immune index and clinical characteristics

To evaluate whether the immune index serves as an independent prognostic factor for HBV-HCC, we performed both univariate and multivariate Cox regression analyses using the TCGA-LIHC dataset. Our results revealed that the immune risk score is an independent prognostic factor in both analyses (HR > 1, p < 0.05), highlighting its strong prognostic value in HBV-HCC patients (Fig. 4A, B). To enhance the clinical utility and diagnostic accuracy, we developed a nomogram that integrates the immune risk score with clinical characteristics (Fig. 4C). The calibration curves indicated well consistency between the nomogram’s predicted outcomes and the actual observations (Fig. 4D).

Independent Clinical predictive value of the immune risk index. (A, B) Forest plot of univariable (A) and multivariable (B) Cox regression results for immune risk index and clinical parameters. (C) Nomogram to predict 1-, 3-, and 5-year survival. (D) Nomogram calibration curves for 1-, 3-, and 5-year OS. (E-G) Stromal (E), Estimate (F) and immune (G) score in the two immune risk groups.

Differentiating immune characteristics based on immune index

The immune infiltration status for different immune risk subgroups was assessed using the ESTIMATE algorithm. Notably, the high-risk group displayed significantly lower stromal (Fig. 4E), ESTIMATE (Fig. 4F), and immune scores (Fig. 4G). To further investigate differences in immune cell infiltration between high- and low-risk groups, we quantified the immune cell composition of each sample using the CIBERSORT algorithm. We observed that the high-risk group had a higher infiltration of M0 and M2 macrophages, but a lower presence of CD8 T cells (Fig. 5A). In addition, the immune score was positively correlated with pro-tumor M0 and M2 macrophages, and negatively correlated with anti-tumor M1 macrophages (Fig. 5B). Moreover, samples with higher levels of M2 macrophages exhibited lower T cell infiltration (p < 0.05, R = −0.8) (Fig. 5C). The comparison of immune risk scores revealed lower scores in the responder group, which further supports the notion that the low-risk group may be more responsive to immunotherapy (Fig. 5D). Immune checkpoint genes such as CD274, PDCD1LG2, BTLA, and CD96 were expressed at significantly lower levels in the high-risk group, indicating a poorer response to immune checkpoint inhibitors (ICI) (Fig. 5E). Furthermore, we found that 14 prognostic related IRGs included in the immune index were positively associated with the infiltration of immunosuppressive M2 macrophages (Fig. 5F).

The immune landscape associated with immune index score in HBV-HCC. (A) The abundance of immune infiltrated cell between high- and low-risk groups, quantified by the CIBESORT algorithm. (B) Correlation analysis between TME infiltrated cells and immune index score. (C) The correlation between the T cell proportion and M2 macrophages proportion. (D) Box plot illustrating the differences in estimated risk scores between non-responders and responders in the GSE202069 cohort. (E) The expression of immune checkpoints in high- and low-risk groups. (F) The correlation between the expression of the 14 key IRGs and immune cell abundance.

Heterogeneity of immune index scores in single-cell dataset

After conducting quality control and removing doublets, we opting for a resolution of 0.2 to partition all cells into 12 distinct clusters (Fig. 6A). Subsequently, we manually annotated 9 cell types based on established literature (Fig. 6B, C). Differentially expressed genes between HBV-HCC and normal tissues were also explored (Fig. 6D), and up-regulated DEGs in HBV-HCC exhibited enriched activity in MYC-targets, G2M checkpoint, mTORC1 pathways indicating a hyperproliferative state (Fig. 6E). Interestingly, cells from tumor tissues generally exhibit lower immune scores than those from normal tissues, suggesting that immunosuppression is prevalent in HBV-HCC (Fig. 6F, G).

Heterogeneity of immune index scores in single-cell dataset. (A) UMAP distribution of 12 clusters at a resolution of 0.2. (B) Expression of top markers corresponding to the 12 cell clusters. (C) UMAP distribution of 9 annotated cell types. (D) Volcano plot presented differential genes between HBV-HCC and normal tissues. (E) Bar chart displaying the up and down regulated hallmarker pathways in HBV-HCC. (F) Box plot of immune index scores of samples derived from HBV-HCC and normal tissues. (G) Box plot of immune index scores of different cells derived from HBV-HCC and normal tissues. (H, I) SPP1 expression levels in 9 annotated cell types.

SPP1 macrophages drive the formation of an immunosuppressive microenvironment

SPP1, identified as the most significant gene in the immune index score, is primarily expressed in macrophages (Fig. 6H, I). To further explore macrophage heterogeneity, we isolated the macrophage and re-clustered it into 5 subclusters using a resolution of 0.1 (Fig. 7A). Referring to the study on macrophage diversity by Ma, Ruo-Yu et al.³² and top differential marker genes (Fig. 7B), we annotated the five macrophage subtypes (Fig. 7C) and characterized their functional properties (Fig. 7E). SPP1 is predominantly expressed in angio-TAMs, which show higher levels of infiltration in HBV-HCC tissues compared to normal tissues (Fig. 7D). Functional enrichment analysis revealed that the SPP1⁺ angio-TAMs exhibited enhanced metastatic potential and immune suppressive characteristics, including activation of inflammatory response pathways, NF-κB signaling, and apoptotic regulation (Fig. 7E). Moreover, SPP1 + angio-TAM infiltration was negatively correlated with T and NK cell infiltration, but positively correlated with epithelial/cancer cell abundance (Fig. 7F). Patients with higher infiltration of SPP1 + angio-TAMs showed significantly worse overall survival, further supporting the prognostic relevance of SPP1 expression (Fig. 7G, H). This finding suggests that the increased infiltration of SPP1 + TAMs likely contribute to shaping an immunosuppressive microenvironment and promoting tumorigenesis in HBV-HCC.

Functional analysis of macrophage subsets. (A) Umap showed the distribution of macrophages subclusters. (B) Top5 marker genes of macrophages subclusters. (C). Annotated macrophages types based on top-ranked marker genes. (D) Proportion of main macrophages subtypes in tumor and normal tissue. (E) Hallmarker pathway enriched in different macrophages subsets. (F) Correlation analysis between angio-TAM infiltration and other immune cell populations. (G, H) Kaplan–Meier survival curves among patients with high and low Angio-TAM infiltration in the TCGA-LIHC and GSE14520 cohorts.

Cell-cell communication network between macrophage subgroups and other cells

Cell-cell communication is intricate and plays a pivotal role in HBV-HCC pathogenesis. Compared with the normal group, the HBV-HCC group exhibited more active intercellular communication (Fig. 8A), particularly with increased signaling interactions between SPP1⁺ angio-TAMs and other cell types (Fig. 8B). Among these interactions, the SPP1 and ApoA signaling pathways contributed the most to the overall communication strength of SPP1⁺ angio-TAMs (Fig. 8B, C). Therefore, we further focused on comparing SPP1 and ApoA signaling between HBV-HCC and normal tissues. We observed that in HBV-HCC tissues, SPP1 signaling was primarily transmitted from angio-TAMs to other cell types, whereas in normal tissues, it was mainly mediated by infla-TAMs (Fig. 8D). Additionally, ApoA signaling from epithelial/tumor cells to SPP1⁺ angio-TAMs was uniquely present in tumor tissues (Fig. 8E, F).

Cell-cell Communications analysis based on HBV-HCC single-cell RNA-sequencing. (A) The number and strength of cell interaction mediated by individual. signal pathways in normal and HBV-HCC groups. (B) Circle plots illustrate the changes in the number and strength of intercellular communications. Red indicates increased signaling in HBV-HCC compared to the normal group, while blue represents decreased signaling. (C) Heatmaps displayed the overall (both outgoing and incoming) signal flows of each cell population. (D) SPP1 signal pathway in normal and HBV-HCC groups. (E) ApoA signal pathway in normal and HBV-HCC groups. (F) Bubble map of cell–cell communication mediated by individual signaling axes, with the horizontal axis showing the cell class that initiates and receives the signal, and the vertical axis showing receptor-ligand pairs of the signaling pathway.

Identification of candidate drugs for the treatment of HBV-HCC based on the hub genes

Combining machine learning-selected prognostic immune-related genes (IRGs) with hub genes from the PPI network, we identified four promising druggable targets: SPP1, ESR1, GHR, and IL33. Among these, SPP1 was a risk factor, showing high expression in tumors, while ESR1, GHR, and IL33 were protective factors, exhibiting lower expression in the tumor microenvironment (Fig. 9A, B). To identify potential drug candidates, we used the DSigDB database on the Enrichr website to investigate drugs that target these genes (Supplementary Table S1). Following this, we performed molecular docking analysis to assess the binding strength between the candidate drugs and the four core target genes. Mefloquine (Fig. 9C), myricetin (Fig. 9D), withaferin (Fig. 9E), estramustine (Fig. 9F), formononetin (Fig. 9G) and apigenin (Fig. 9H) exhibited binding energies lower than − 5 kcal/mol, indicating strong binding affinities to their corresponding targets (Supplementary Table S5). Visible hydrogen bonding and strong electrostatic interactions enable each medication candidate to bind to its protein target.

Small molecule drug prediction and molecular docking verification. (A) Forest plot of univariable Cox regression results for the 14 key IRGs. (B) Immunohistochemical analysis of promising targets in HCC tumor and corresponding adjacent normal tissues. (C-H) Visualization of docking models using PyMOL (v2.6.0).