Integrating bulk and single cell sequencing data to identify prognostic biomarkers and drug candidates in HBV associated hepatocellular carcinoma

Abstract

Hepatitis B virus (HBV) infection is a major driver of hepatocellular carcinoma (HCC), yet the mechanisms by which HBV triggers HCC and how it interacts with the immune system remain largely undefined. In this study, 53 immune-related key genes involved in HBV-associated HCC progression were identified. By analyzing the mean C-index of 101 machine learning models, the optimal model-combining stepwise Cox regression (forward) with RSF-was developed to characterize the immune risk index. Patients in the high-risk group exhibited worse survival outcomes and increased infiltration of immunosuppressive cells. Integrating PPI analysis with machine learning, SPP1, GHR, and ESR1 emerged as promising druggable targets, with SPP1 notably overexpressed in tumors and linked to adverse outcomes. ScRNA-seq analysis revealed SPP1 was predominantly expressed in angio-TAMs, which may impair anti-tumor immunity by limiting T and NK cell infiltration. It also involved in tumor progression via angiogenesis and EMT pathways. Drug prediction and molecular docking identified small molecules such as myricetin and mefloquine that can target the aforementioned key immune genes, thereby modulating the immune landscape of HBV-HCC. Repurposing these established drugs represents a novel therapeutic avenue, offering both efficacy and expedited clinical translation for HBV-HCC.

Keywords: HBV associated hepatocellular carcinoma; Immunosuppressive microenvironment; Machine learning; Multi-target Docking; Single-cell RNA sequencing.

Fig. 1

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^–. (F) PPI network constructed using 53 HBV-HCC related immune genes.

Fig. 2

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).

Fig. 3

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).

Fig. 4

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.

Fig. 5

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.

Fig. 6

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.

Fig. 7

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.

Fig. 8

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.

Fig. 9

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).