Tumor Prognostic Risk Model Related to Monocytes/Macrophages in Hepatocellular Carcinoma Based on Machine Learning and Multi-Omics

Abstract

Tumor-associated macrophages (TAMs) are crucial in hepatocellular carcinoma (HCC) development and invasion. This study explores monocyte/ macrophage-associated gene expression profiles in HCC, constructs a prognostic model based on these genes, and examines its relationship with drug resistance and immune therapy responses. Single-cell RNA sequencing(scRNA-seq) data from 10 HCC tissue biopsy samples, totaling 24,597 cells, were obtained from the GEO database to identify monocyte/macrophage-associated genes. A prognostic model was constructed and validated using external datasets and Western blot. Relationships between the model, clinical correlates, drug sensitivity, and immune therapy responses were investigated. From scRNA-seq data, 2,799 monocyte/macrophage marker genes were identified. Using the TCGA dataset, a prognostic model based on the single-gene UQCRH was constructed, stratifying patients into high-risk and low-risk groups based on overall survival rates. High-risk group patients showed reduced survival rates and higher UQCRH expression in tumor tissues. Western blot analysis further confirmed the elevated expression of UQCRH in HCC cell lines. Spatial transcriptomics analysis revealed that high UQCRH expression co-localized with malignant cells in the tumor tissue. Drug sensitivity analysis revealed that the high-risk group had lower sensitivity to sorafenib and axitinib. Immune therapy response analysis indicated poorer outcomes in the high-risk group, with more pronounced APC inhibition and a weaker IFN-II response. Clinical indicator analysis showed a positive correlation between high UQCRH expression and tumor invasion. Enrichment analysis of UQCRH and associated molecules indicated involvement in oxidative phosphorylation and mitochondrial electron transport. This study introduces a prognostic model for HCC patients based on monocyte/macrophage marker genes. The single-gene model predicts HCC patient survival and treatment outcomes, identifying high-risk individuals with varying drug sensitivities and immune suppression states.

Introduction

Hepatocellular carcinoma (HCC) is the primary malignancy of the liver and constitutes a significant health challenge globally, particularly in regions such as Africa and Asia, where there is a high incidence of the disease related to risk factors like chronic hepatitis B and C virus infections, cirrhosis, and non-alcoholic fatty liver disease (NAFLD) [1]. As the third leading cause of cancer-related deaths worldwide, there is an urgent need for improved treatment strategies and accurate diagnostic methods for HCC [2]. The incidence of HCC is significantly influenced by geographical differences, with varying prevalence and related risk factors across different regions, highlighting the complexity of HCC with environmental factors and genetic susceptibility. Advances in targeted therapies and the emergence of immunotherapy bring new hope to the field, yet treatment remains challenging due to high recurrence rates and the late stage of disease at diagnosis for many patients [2, 3]. Current treatments focus on innovative methods such as local region treatments, systemic treatments, and immunotherapy, aiming to utilize the human immune system to fight cancer [4]. Despite these advancements, the prognosis for HCC remains poor [5], emphasizing the urgent need for early detection methods and more effective treatment options. Research on molecular biomarkers and the tumor immune microenvironment remains crucial, with ongoing studies aiming to develop personalized treatment plans that could significantly improve patient outcomes [1]. These efforts are essential for managing the disease, which is often diagnosed at an advanced stage and difficult to treat curatively [5].

Against this backdrop, the role of monocytes/macrophages in the HCC immune microenvironment has become a key area of investigation. These immune cells play a critical role in shaping the immune landscape of the tumor and directly influence its progression [6]. The intrinsic characteristics of tumor cells, along with various external factors in the tumor microenvironment, contribute to significant changes in the metabolism, phenotype, and biological functions of macrophages [7, 8]. These alterations drive the reprogramming of macrophages within the tumor microenvironment to support tumor growth [9, 10].

In recent years, bioinformatics approaches have been used to identify gene clusters and regulatory networks related to HCC prognosis, immune microenvironment, and ferroptosis, providing valuable insights into HCC outcomes [11,12,13]. However, traditional RNA-seq analysis, which reflects the average gene expression across all cells in a sample, fails to capture the molecular expression changes in specific cell types within the tumor microenvironment [14].

Single-cell RNA sequencing (scRNA-seq) is an advanced technology that enables gene expression analysis at the single-cell level, revealing differences between cell types and their interactions, as well as the heterogeneity of disease gene expression [15,16,17,18,19,20,21]. Machine learning (ML), particularly in oncology genomics, has become a crucial tool in improving medical accuracy and patient prognosis. Its applications in risk assessment, early diagnosis, prognosis evaluation, and therapeutic drug selection have already demonstrated significant effectiveness [22].

This study integrates machine learning with multi-omics sequencing to perform single-cell RNA sequencing and transcriptomic analysis on HCC tumor samples, focusing on cellular and molecular changes. A single-gene prognostic risk model for HCC patients is developed and validated using transcriptomic and proteomic datasets through survival analysis and expression differentials. Experimental validation of gene expression and spatial transcriptomics is also performed to examine gene localization in malignant cells. Finally, the model’s drug sensitivity, predicted immune therapy responses, and immune function are analyzed. The study workflow is illustrated in Fig. 1.

Fig. 1

The workflow of the present study. This study integrates machine learning with multi-omics sequencing to analyze HCC. A prognostic risk model is developed and validated, with analyses of drug sensitivity, immune therapy responses, and immune function

Materials and methods

Acquisition and Preprocessing of Single-Cell Data

The single-cell sequencing data were sourced from the article by Lu et al. titled “A single-cell atlas of the multicellular ecosystem of primary and metastatic hepatocellular carcinoma”, available in the GEO database under the dataset identifier GSE149614 [23]. We extracted single-cell sequencing data from ten primary tumor samples within this dataset for analysis. The raw data for each sample were processed using the Seurat package in RStudio [24]. Quality control was performed on the analysis samples, selecting cells containing more than 300 genes and fewer than 2500 genes, totaling 24,597 cells for further bioinformatics analysis.

Data Integration and Dimension Reduction

The Seurat objects containing gene expression data for each sample were processed using the Read10x() function [25]. After quality control, data were normalized using the ‘NormalizeData’ function. Gene expression was represented by scaling each gene’s expression per sample by a factor of 10,000, followed by transformation to natural logarithm and normalization after adding one to avoid log of zero. The top 2000 highly variable genes (HVG) were identified using ‘FindVariableFeatures’ and the top 10 HVG were visualized.

Cell Clustering and Annotation

Cells were clustered using Seurat based on PCA scores, where each principal component (PC) represents a ‘meta-feature’ combining information from a related feature set. Dimensionality reduction was performed on the selected top 2000 highly variable genes using the ‘RunPCA’ function, and the effects of the reduction were visualized. Seurat employed spectral clustering to identify cell clusters and used ‘FindNeighbors’ and ‘FindClusters’ functions to distinguish all cell types. T-distributed Stochastic Neighbor Embedding (tSNE) was used for nonlinear dimension reduction visualization, as shown in Fig. 2(C) [26, 27]. Cell markers for hepatocellular carcinoma from CellMarker 2.0 (http://117.50.127.228/CellMarker/CellMarker_ help.html) [28] were referenced to annotate the discerned clusters of cell types. Subsequently, the proportion of each cell type in each sample was analyzed.

Fig. 2

Preliminary analysis of single-cell RNA sequencing samples. (A) Visualization of highly variable genes. (B) Principal component analysis (PCA) of the 10 included samples. (C) Clustering of cells from the included samples. (D) Annotation of cell types through marker genes. (E) Bubble chart showing the expression of marker genes across eight cell types. (F) Violin plots displaying the expression characteristics of marker genes in eight cell types. (G) Bar charts analyzing the average expression ratio of eight cell types across 10 single-cell samples and their individual proportions in each of the 10 samples

Marker Gene Identification

Using the ‘FindAllMarkers’ function in R software, marker genes for monocytes/macrophages in single-cell data were identified with a minimum percentage (min.pct) of 0.01 and a log2 fold change threshold (log2FC.threshold) of 0.01.

Cell-Cell Communication Analysis

Systematic analysis of signaling pathways between cells helps us understand the overall communication between them [29]. Using the R package CellChat, cell communication networks were constructed to identify key signaling pathways between macrophages and other cell types in the tumor microenvironment, such as T cells, inferring intercellular communication at the signaling pathway level.

Construction of Transcription Factor Gene Regulatory Networks

SCENIC is a computational method that reconstructs gene regulatory networks(GRN) and identifies cell states simultaneously from single-cell RNA-seq data [30]. Using SCENIC, stable cell states can be identified. Initially, co-expression networks of transcription factors and their target genes were analyzed using the SCENIC framework’s Genie3 and AUCell on single-cell RNA data. These tools were used to derive transcriptional regulatory networks for monocyte/macrophage cell types and identify binding motifs of regulators using RcisTarget, thereby determining direct target relationships in the network. Subsequently, Regulons of single cells were scored based on each regulator and its target gene set, and cell activity states were predicted using a generated binary matrix. Finally, the transcription factor gene regulatory networks were constructed using Cytoscape. This process not only aids in understanding the functional states of cells but also reveals complex interaction networks between them.

Extraction and Analysis of Transcriptomic and Proteomic Data

Data on HCC mRNA and related clinicopathological features were downloaded from the TCGA database (https://portal.gdc.cancer.gov/), which included gene expression data and clinical features from 424 samples [31]. Additionally, transcriptomic sequencing samples from 389 HCC patients and proteomic sequencing samples from 330 patients were obtained from the ICGC database (https://dcc.icgc.org/) and the CPTAC database (https://proteomic.datacommons.cancer.gov/pdc/) respectively [32, 33]. After organizing and analyzing gene expression in LIHC samples from the TCGA database, genes in the top 70% of average gene expression in the TCGA-LIHC database were intersected with monocyte/macrophage marker genes from the single-cell data, and the extracted expression matrix was normalized.

Construction and Validation of a Single-Gene Prognostic Model

To identify significant gene expression markers affecting patient survival times, univariate Cox regression analysis was first performed on the intersected genes to identify potential differentially expressed genes [34]. Then, variables with p < 0.05 were selected for regression analysis using the LASSO method, and the Cox proportional hazards model-based Lasso regression was constructed using the glmnet() function in the R package glmnet [35]. To evaluate the predictive performance of the model and prevent overfitting, the model’s C-index was calculated and implemented using the cv.glmnet() function, selecting the C-index as the performance metric. Significant non-zero coefficient genes identified from the Lasso model were extracted, and their expression data along with patient survival information were used for constructing the prognostic model. To further evaluate the role of key genes identified from Lasso regression analysis in predicting patient survival, a multivariate Cox regression model was employed, and the prognostic model’s Risk Score was calculated [36, 37]. The Risk Score is calculated using the following formula:

$${\rm{Risk}}\,{\rm{Score}}\, = \,\sum\limits_{i = 0}^k {\beta i*exp{i^1}} $$

(1)

Patients were then divided into high and low-risk groups based on the median Risk Score using the R packages survival and survminer, and Kaplan-Meier survival curves and risk plots were drawn and visualized along with survival status and differences for each sample. Additionally, the risk score’s predictive performance for 1-year, 3-year, and 5-year patient survival was evaluated using the R packages pROC and timeROC. Integrating risk and clinical data, the function of the timeROC package was applied using the Aalen additive model to compute ROC curves and the corresponding area under the curve (AUC) at various time points without competing risks, quantifying the efficacy of risk scores and other clinical variables in predicting patient survival, with the AUC value visually demonstrating the risk model’s ability to predict long-term patient survival [38]. Furthermore, HCC datasets from the ICGC database, GEO database(GSE14520) and the CPTAC dataset were selected as independent external cohorts to validate the prognostic model’s effectiveness.

Correlation of Risk Model and its Gene Expression With Clinicopathological Features

To assess the correlation of clinical variables with patient prognosis, we used a nomogram for visual analysis. Integrating risk and clinical data, a Cox regression model was built using the cph function in the R package rms, and based on this model’s results, a nomogram was generated using the regplot package. The nomogram intuitively shows the relative contribution of each clinical variable to the prognosis, with each variable assigned points according to its coefficient in the model. By summing the points of all variables, the predictive probability, i.e., the survival probability of patients at 1 year, 3 years, and 5 years, is obtained. Additionally, to verify the accuracy of the model, we also used the calibrate function to plot calibration curves, comparing the model-predicted survival probabilities with the actual observed survival rates. This method not only enhances the interpretability of the predictive model but also increases its practical value in clinical decision-making [39]. Furthermore, we used the R package ComplexHeatmap to draw heatmaps to visually assess the correlation between genes in the risk model and clinical features.

Drug Sensitivity Analysis, Predicting Immune Therapy Response, and Immune Function Scoring

Next, we conducted drug sensitivity analysis using the R package oncoPredict to assess the differences in drug responses between the high and low-risk groups of the prognostic model. By combining TCGA-LIHC data and GDSC drug sensitivity data, oncoPredict was used to calculate drug sensitivities for the test samples [40]. Subsequently, the Wilcoxon test was used to compare drug sensitivities between the high and low-risk groups, and the results were visually presented using box plots drawn by the R package ggplot2. We then analyzed immune therapy responses in tumor samples using the Tumor Immune Dysfunction and Exclusion (TIDE) method [41]. By organizing and normalizing expression files from the TCGA-LIHC dataset and merging them with TIDE prediction results and clinical risk data, we compared TIDE scores between high and low-risk groups within the prognostic model to assess the potential immune therapy responses of patients with different risk levels. The analysis results were graphically displayed using the ggviolin in the R package ggpubr, further revealing the relationship between immune therapy responses and patient prognostic risk. Lastly, we conducted immune function scoring on the TCGA-LIHC dataset using the ssGSEA method and analyzed the immune-related gene sets of each sample using the R package GSVA. By evaluating the normalized scores combined with clinical risk data, we assessed the relationship between immune function scoring and patient prognosis.

Immunohistochemical Expression of Single Genes in the Prognostic Model, Paired Sample Differential Expression, and Enrichment Analysis of Directly Associated Proteins

We obtained immunohistochemical data for single genes in the prognostic model from normal liver tissue and HCC tissue from the HPA database (https://www.proteinatlas.org/) [42]. Additionally, we used the TCGA-LIHC dataset to analyze differential expression between normal and tumor samples for single genes in the prognostic model. Finally, based on the STRING database (https://string-db.org/), proteins directly associated with the single protein in the prognostic model were identified, and their protein interaction networks were visualized using Cytoscape. Subsequent KEGG and GO enrichment analyses were performed to understand the biological processes involving the single gene in the prognostic model [43].

Cell Lines and Antibodies

Human HCC cell lines, Hep-G2 and Huh-7, as well as LX2 cells, were obtained from Guangzhou Saiku Biotechnology Co., Ltd. (China). All cells were cultured at 37 °C in a humidified atmosphere with 5% CO2 and 95% air, in RPMI-1640 medium (Life Technologies, Carlsbad, CA, USA) supplemented with 10% fetal bovine serum (FBS, Gibco, USA) and 0.5% penicillin-streptomycin sulfate (Invitrogen, Carlsbad, CA, USA). Cell counts were performed using an automated cell counter (Star, Invitrogen, Carlsbad, CA, USA). GAPDH antibody was purchased from Cell Signaling Technology Inc. (Beverly, MA, USA), and anti-UQCRH (ab154803) antibody was purchased from Abcam (Cambridge, MA, USA).

Western Blot Analysis

Cells were harvested, washed, and lysed with 1× RIPA buffer. Protein concentrations were determined using the Thermo BCA Protein Assay Kit. Equal amounts of protein from the cell lysates were mixed with 5× SDS sample buffer, separated by 12% SDS-polyacrylamide gel electrophoresis, and transferred to polyvinylidene fluoride (PVDF) membranes. Membranes were blocked with 5% non-fat milk in TBST and incubated with primary antibodies against UQCRH and GAPDH at 4 °C overnight. After washing, the membranes were incubated with a secondary antibody against rabbit IgG for 1 h, followed by further washing. Detection was performed using ECL solution (Millipore, Darmstadt, Germany) and images were captured using the Bio-Rad ChemiDoc XRS + Chemiluminescence Imaging System (Bio-Rad, Hercules, CA, USA). Protein expression levels were quantified using ImageJ software.

Spatial Transcriptomics Analysis of UQCRH Expression

We downloaded HCC spatial transcriptomics samples (Sample ID: VISDS000508, VISDS000511) from the CROST database (https://ngdc.cncb.ac.cn/crost/). Using the R package spacexr, we mapped the single-cell transcriptomic data onto the tissue space and analyzed the gene expression patterns in different tissue regions based on spatial coordinates.

Statistical Analysis

The Wilcoxon signed-rank test was used to compare the expression differences of individual genes between normal and tumor samples. LASSO regression and Cox regression were employed to establish predictive models for disease prognosis, and Kaplan-Meier survival curves were used to evaluate the survival differences between different risk groups, with the log-rank test determining statistical significance. All statistical analyses were performed in Rstudio, and a two-sided p-value of less than 0.05 was considered statistically significant.

Results

Cellular Composition of Hepatocellular Carcinoma and Gene Expression in Monocytes/Macrophages

This study included scRNA-seq samples from 10 primary hepatocellular carcinoma tumors. Each sample was analyzed by integrating the top 2000 highly variable genes, followed by PCA (Principal Component Analysis) and subsequent clustering dimension reduction (Fig. 2A, B). t-SNE clustering visualization showed that 24,597 cells from hepatocellular carcinoma tissues were grouped into 19 subgroups (Fig. 2C), and clusters were identified as different cell types based on marker genes for hepatocellular carcinoma cell types retrieved from CellMarker 2.0 (Fig. 2D). Information about the cell subgroups is provided in Table 1. To clarify the characteristics of each cell subgroup, we defined the marker genes for each cell subgroup and visualized them through bubble plots and violin plots (Fig. 2E, F). Through cell annotation of the clusters, a total of 8 cell types were identified, including T cells, monocytes/macrophages, plasma cells, hepatocytes, fibroblasts, endothelial cells, epithelial cells, and B cells, with T cells, monocytes/macrophages, and hepatocytes present in multiple cell clusters. Subsequently, we analyzed and visualized the average number of each cell type across all samples and their numbers in each individual sample (Fig. 2G), finding a higher number of monocytes/ macrophages, which suggests that monocytes/macrophages may play a significant role in the development and invasion progression of HCC. Consequently, we further identified marker genes related to monocytes/macrophages in the single-cell data, with related data presented in Supplementary Table 1.

Table 1 Information about the cell subpopulations

Cell Communication Analysis Related To Monocytes/Macrophages in scRNA-seq

Using CellChat, we constructed intercellular communication networks by calculating cell communication probabilities, while independently analyzing the communication networks between each cell type and other cell types (Fig. 3A-C). Additionally, we inferred the direct communication networks involving monocytes/macrophages with other cell types and discovered pathways such as FN1 and SPP1, which exhibited intensive communication between monocytes/ macrophages and other cell types (especially hepatocytes) (Fig. 3D, E).

Fig. 3

Cell communication analysis. (A) Construction of a cell communication network for the 10 included tumor samples. (B) Independent display of the number of cell communication events between each cell type. (C) Independent display of the communication weights between each cell type. (D) In-depth analysis of the role of the FN1 signaling pathway in cell communication using chord diagrams and heatmaps. (E) Further analysis of the SPP1 signaling pathway in cell communication through chord diagrams and heatmaps

Transcription Factor Regulatory Networks Related To Monocytes/ Macrophages in scRNA-seq

Based on SCENIC analysis, key regulatory networks of transcription factors within monocytes/macrophages were predicted. ZNF385A, NFKB2, and NONO (three core transcription factors) will serve as central hubs controlling HCC-Mono/Macro (Fig. 4A). Subsequently, we visualized the predicted transcription factors and their regulatory genes using Cytoscape 3.10.1 to construct the transcription factor-gene regulatory network (Fig. 4B), where red represents transcription factors, and green represents regulated genes.

Fig. 4

Predicting transcription factor regulatory networks associated with monocytes/macrophages in scRNA-seq samples using SCENIC. (A) ZNF385A, NFKB2, and NONO (three core transcription factors) will become the central hubs controlling the LIHC state in monocytes/macrophages. (B) In-depth analysis of the transcription factor-gene regulatory network centered on ZNF385A, NFKB2, and NONO

Construction and Validation of Prognostic Model and Gene Expression

In the TCGA-LIHC database, genes in the top 70% of average gene expression were intersected with 2799 monocyte/macrophage marker genes from single-cell data (Supplementary Table 2), and visualized using a Venn diagram (Fig. 5A). The extracted expression matrix was then normalized. Subsequently, univariate Cox regression analysis was performed on 700 intersecting genes, identifying 197 genes with potential prognostic significance (p < 0.05) (Supplementary Table 3). To prevent overfitting of the prognostic model, Lasso regression analysis was conducted, resulting in the selection of 8 genes with prognostic significance (Fig. 5B, Supplementary Table 4). Ultimately, multivariate Cox regression analysis led to the identification of one gene with significant prognostic relevance, namely, Ubiquinol-cytochrome c reductase hinge protein (UQCRH), for which a single-gene prognostic model was constructed (Supplementary Tables 5–6). First, paired sample differential expression analysis of UQCRH revealed that this gene is significantly overexpressed in tumor samples from the TCGA-LIHC dataset (Fig. 5C). Immunohistochemical results from the HPA database showed that UQCRH has higher expression in hepatocellular carcinoma tumor cells compared to normal liver tissue (Fig. 5D, E). Subsequently, based on the median Risk Score calculated from the single-gene prognostic model, patients were divided into high and low-risk groups for survival analysis and ROC analysis. Survival curve analysis indicated that the high-risk group and the single-gene high-expression group had worse overall survival (Fig. 5F-H). Additionally, in the TCGA-LIHC data, ROC analysis demonstrated that the single-gene prognostic model’s risk score effectively predicts patient overall survival, with AUC values of 0.739, 0.628, and 0.605 for 1-year, 3-year, and 5-year OS, respectively (Fig. 5I). Moreover, by analyzing the GEO database, ICGC database and CPTAC database, Kaplan-Meier survival curve analysis and differential expression analysis of the UQCRH gene in HCC at the transcriptomic and proteomic levels revealed that high expression of UQCRH is associated with poor OS prognosis and is overexpressed in tumor samples, as shown in Fig. 5J, K, L.

Fig. 5

Construction of a single-gene prognostic risk model and validation of UQCRH overexpression in HCC using external independent datasets and Western blot. (A) Intersection of monocyte/macrophage expressed genes from scRNA-seq with TCGA-LIHC expressed genes using a Venn diagram. (B) To prevent model overfitting, the prognostic risk model is further constructed based on LASSO regression analysis. (C) Differential expression analysis in paired samples from the TCGA-LIHC database revealed overexpression of UQCRH in tumor samples. (D, E) Analysis from the HPA database shows high expression of UQCRH in HCC samples. (F-H) Survival risk curves indicate that overexpression and high-risk groups in the single-gene prognostic model are associated with poor prognosis in HCC patients. (I) ROC analysis of the single-gene prognostic model risk shows good performance in predicting patient OS with an AUC score. (J) Transcriptomic analysis from the ICGC-HCC database confirms the poor prognosis associated with high UQCRH expression in patients (p = 0.0427). (K) Transcriptomic analysis from the GEO database(GSE14520) confirms overexpression of UQCRH in HCC tumor samples (p < 2.2e-16). (L) Proteomic analysis from the CPTAC database confirms overexpression of UQCRH in HCC tumor samples (p < 0.001). (M) Western blot analysis confirms the significantly elevated expression of UQCRH in HCC cell lines

UQCRH Expression Differs Significantly between Non-HCC and HCC Cells

Our Western blot results show that UQCRH is barely expressed in the non-HCC cell line LX2, but significantly upregulated in the HCC cell lines Hep-G2 and Huh-7 (Fig. 5M). The high expression in Hep-G2 and Huh-7, compared to LX2, suggests that UQCRH is specific to HCC cell lines, indicating its potential role in the development and progression of HCC.

Spatial Transcriptomic Analysis of UQCRH Expression in HCC

We performed cell annotation using scRNA-seq data of HCC from GSE166635 (Fig. 6A, B), followed by joint analysis and spatial mapping of two HCC spatial transcriptomics datasets obtained from the CROST database. Our analysis revealed that high expression of UQCRH co-localizes spatially with malignant cells within the HCC microenvironment (Fig. 6C-F). Additionally, high expression of UQCRH also co-localizes with monocyte/macrophages, particularly in HCC sample 2 (Supplementary Fig. 1).

Fig. 6

Spatial analysis of UQCRH expression and malignant cell distribution in HCC samples. (A) UMAP clustering showing different cell types in HCC, including malignant cells, T cells, and others. (B) Dot plot of marker gene expression across cell clusters, highlighting markers for malignant cells. (C) Spatial distribution of malignant cells in HCC tissue, with redder colors indicating higher abundance. (D) Spatial expression of UQCRH in HCC sample 1, with redder colors representing higher expression levels. (E) Spatial distribution of malignant cells in HCC sample 2, similar to sample 1. (F) Spatial expression of UQCRH in HCC sample 2, showing co-localization with malignant cells

Analysis of Prognostic Model and its Single Gene’s Clinical Relevance

We first extracted clinical data from the TCGA-LIHC dataset, including patient age, gender, T, N, M, and Stage, to analyze the correlation between UQCRH expression and these clinical parameters. We found that high expression of UQCRH is positively correlated with tumor invasion and infiltration (especially T and Stage) (Fig. 7A). Next, based on the algorithm, age, gender, stage, T, N, M, and risk score were included in a nomogram (Fig. 7B). Then, the points of each variable were combined into a total score to predict the probability of OS at 1 year, 3 years, and 5 years (Fig. 7C), and the calibration curve shows that the nomogram’s predictions are well aligned with the actual outcomes.

Fig. 7

Prognostic value of UQCRH expression and risk score in HCC. (A) In the TCGA-LIHC database, high expression of UQCRH is positively correlated with clinical indicators T and Stage (p < 0.001). (B) Nomogram analysis reveals that patient prognosis is primarily associated with Risk Score (p < 0.01). (C) Calibration curves show good agreement between the predicted values of the nomogram and actual values for 1-year, 3-year, and 5-year OS probabilities

UQCRH-Related Protein Interaction and Enrichment Analysis

We used the STRING database to retrieve proteins directly associated with UQCRH and used Cytoscape to construct a protein interaction network, followed by enrichment analysis. We found that UQCRH is associated with proteins such as UQCRC1, UQCRC2, and UQCR11, which are involved in biological processes and signaling pathways including mitochondrial electron transport, cellular energy metabolism, and oxidative phosphorylation (Fig. 8A-C).

Fig. 8

Functional analysis of UQCRH and its associated proteins in HCC. (A) Protein-protein interaction analysis identifies direct associations between UQCRH and proteins such as UQCRB and UQCRQ. (B, C) KEGG and GO enrichment analyses show that UQCRH and its directly associated proteins are involved in biological processes and signaling pathways such as mitochondrial electron transport, maintenance of cellular energy metabolism, and oxidative phosphorylation

Drug Sensitivity, Immune Therapy, and Immune Function Analysis of the Prognostic Model

Using machine learning methods, we conducted drug sensitivity analysis on the single-gene prognostic model. Our analysis found that the high-risk group is less sensitive to anti-cancer drugs such as sorafenib and axitinib. (Fig. 9A). Additionally, in the immune therapy prediction analysis, we found that the high-risk group has higher TIDE scores, indicating poorer outcomes with immune therapy, whereas the low-risk group shows better effects (Fig. 9B). We then used the ssGSEA algorithm to score immune functions in the high and low-risk groups of the prognostic model, finding that the high-risk group may exhibit more active suppression of antigen presentation, while the low-risk group shows better IFN-II response (Fig. 9C).

Fig. 9

Drug sensitivity, immune therapy response, and immune function analysis in high-risk and low-risk groups of HCC patients. (A) Drug sensitivity analysis indicates that the high-risk group in the prognostic model is less sensitive to commonly used anti-cancer drugs such as sorafenib and axitinib. (B) Immune therapy response prediction within the prognostic model shows that the high-risk group exhibits less pronounced responses to immune therapy. (C) ssGSEA immune function scoring reveals that the high-risk group in the prognostic model exhibits suppressed antigen presentation and poor response to IFN-II in terms of tumor immunity

Discussion

This study utilized a scRNA-seq dataset originally from Lu et al.‘s research, with no overlap in the analysis focus compared to the original study [23]. This research analyzed scRNA-seq data from 10 primary HCC samples from a public dataset, aiming to identify clinically relevant tumor immune microenvironments and tumor characteristics. The analysis revealed a high frequency of interactions, particularly between monocytes/macrophages and other cell types such as hepatocytes and fibroblasts. These interactions are not merely incidental but indicative of the pivotal role monocytes/macrophages play in sculpting the TME. Monocytes/macrophages engage in signaling processes that contribute to tumor progression, immune evasion, and the establishment of an immunosuppressive microenvironment. This observation aligns with established knowledge in HCC pathogenesis, highlighting the prominent role of immune cells in facilitating cancer invasion and progression. Further insights were gleaned from the interaction strength between these cell types. Strong communication signals were noted between immune cells and tumor cells, underscoring the deep involvement of monocytes/macrophages in altering the immune landscape. For instance, existing studies have found that FN1 is involved in cell adhesion and migration processes, and it is associated with promoting tumor invasion and poor prognosis [44]. Xu demonstrated through experiments that targeting the downregulation of FN1 significantly inhibits the invasion of HCC cells [45]. Liu considers SPP1 as an immune-related prognostic factor in HCC, mediated by SPP1-CD44 and SPP1-PTGER4 interactions between HCC cells and macrophages [46]. In vitro experiments have shown that SPP1 can trigger macrophage polarization towards M2-like tumor-associated macrophages, revealing its carcinogenic mechanism in HCC. Additionally, analysis of the transcription factor-gene regulatory network in monocytes/ macrophages of HCC found that ZNF385A, NFKB2, and NONO play significant gene regulatory roles. Peng’s research found that Zinc finger protein 385 A (ZNF385A) is overexpressed in HCC and related to poor prognosis [47]. Hepatitis B virus infection can lead to overexpression of ZNF385A, accompanied by increased apoptosis and chronic inflammation. Moreover, they are positively correlated with immune-suppressive cells, inflammatory cytokines, immune checkpoint genes, and poor immunotherapy outcomes. Knockdown of ZNF385A expression inhibited tumor cell invasion and migration. NFKB2, a classic cancer-related transcription factor, encodes the p100 protein, which transforms into p52 in the non-canonical NF-κB signaling pathway, regulating cell survival, proliferation, and apoptosis. Although its role in HCC is not clearly defined, NFKB2 mutations and fusions are observed in other cancers like bone, skin, and cervical cancers [48]. Tang’s study indicated that inhibiting the NF-κB signaling pathway, involving Akt and NFKB2, can suppress liver cancer cell proliferation, migration, and invasiveness [49]. Shen found that NONO plays a significant role in HCC, especially under hypoxic conditions, where NONO interacts with and stabilizes the HIF-1 and HIF-2 complexes, thereby activating the transcription of hypoxia-induced genes, allowing malignant proliferation of tumor cells and promoting the development of liver cancer [50]. However, it is important to note that the transcription factors identified in the SCENIC analysis primarily regulate gene expression related to monocytes/macrophages, and these factors may not directly regulate all genes involved in HCC progression. In contrast, the key prognostic markers in our study were identified through a broader multi-omics analysis, which integrated various data beyond specific networks. This broader analysis allowed us to specifically focus on the expression of monocyte/macrophage-related genes, as they play a key role in shaping the immune microenvironment of HCC.

Building on these observations, our study aims to explore the expression of monocyte/macrophage marker genes in HCC and construct a prognostic model for HCC. To identify differentially expressed genes in monocytes/macrophages in HCC, we utilized scRNA-seq data from the literature, selecting tissue biopsy samples from 10 primary HCC patients, encompassing data from 24,597 single cells. Based on these observations, our study aims to explore the expression of monocyte/macrophage marker genes in HCC and construct a prognostic model for HCC. To identify differentially expressed genes in monocytes/macrophages in HCC, we utilized scRNA-seq data from literature, selecting tissue biopsy samples from 10 primary HCC patients, encompassing data from 24,597 single cells. Based on this data, we developed a prognostic risk model for HCC targeting monocyte/macrophage marker genes on the TCGA database, further validated through data from the ICGC and CPTAC databases. Our analyses, including single-gene survival analysis, differential expression analysis, prognostic model risk scoring, and correlation analysis with clinical indicators, demonstrated that UQCRH is overexpressed in HCC tumor samples. Notably, although immunohistochemical analysis shows UQCRH expression in normal tissues, its expression is significantly higher in HCC, and its presence in normal tissues may be related to its basic function in mitochondrial metabolism. UQCRH and its prognostic model may play a role in promoting cancer progression in HCC by facilitating tumor cell invasion. Mitochondrial membrane complexes, such as Complexes I, II, III, and IV of the respiratory chain, and ATP synthase, mainly located on the inner mitochondrial membrane, are directly involved in the energy conversion process, generating ATP through oxidative phosphorylation. The function of mitochondrial membrane complexes is closely related to mitochondrial membrane dynamics and plays a central role in cell death, particularly programmed cell death, with increased permeability of the mitochondrial outer membrane allowing the release of apoptotic inducers like cytochrome c into the cytoplasm, a critical step in apoptosis [51, 52]. UQCRH, a vital component of mitochondrial respiratory chain Complex III, helps maintain the stability of the complex’s structure and function, indirectly influencing the transfer of electrons from ubiquinone through Complex III to cytochrome c. This function is crucial in the energy metabolism and cellular respiration processes of normal cells but is often altered in expression and functionality in cancer cells [53,54,55,56]. Modena et al. found that UQCRH expression is reduced or even absent in some cancer cell lines, such as ovarian and breast cancer, possibly due to structural rearrangements or epigenetic inactivation mechanisms like promoter methylation in these cancer cell lines [57]. Additionally, in clear cell renal cell carcinoma (ccRCC), UQCRH often exhibits DNA hypermethylation, leading to mRNA downregulation. This hypermethylation status is closely associated with poor clinical outcomes, indicating UQCRH’s potential tumor-suppressive role in tumor progression. In ccRCC, UQCRH downregulation is associated with shorter survival times, and its overexpression can restore mitochondrial function, increase oxygen consumption, reduce the Warburg effect, and induce apoptosis in cancer cells [53, 54]. Conversely, in lung adenocarcinoma, UQCRH can induce the production of reactive oxygen species (ROS), and its expression in patients with lung adenocarcinoma is significantly higher than in patients with pneumonia and healthy controls, showing high sensitivity and specificity for diagnosing lung adenocarcinoma [55]. Eun-Ran Park found through clinical sample analysis that UQCRH is overexpressed in HCC patient samples and is associated with tumor size, poor differentiation, and vascular invasion. Overexpression of UQCRH in HCC patients significantly shortens overall and recurrence-free survival times, especially in patients with elevated alpha-fetoprotein levels, thus indicating that UQCRH could be an indicator of poor prognosis in HCC [56]. This finding aligns with the results from our bioinformatics analysis, further confirming UQCRH’s oncogenic role in HCC, suggesting its potential as an important clinical prognostic marker.

Based on this understanding, in order to more comprehensively evaluate the guidance significance of UQCRH expression on therapeutic strategies, we conducted further analyses on drug sensitivity, immune therapy responses, and immune function scoring in our prognostic model. These analyses helped us validate the biomarker potential of UQCRH from multiple perspectives and explore its value in precision medicine. In our prognostic model, the high-risk group exhibited insensitivity to various targeted therapy drugs, including sorafenib and axitinib. This may be due to mitochondrial dysfunction caused by abnormal expression of UQCRH, which affects the energy and drug metabolism of HCC tumor cells, leading to reduced intracellular drug accumulation, weakened drug efficacy, and the promotion of drug resistance in HCC tumor cells [58]. Additionally, in our prognostic model, patients in the high-risk group showed poorer efficacy in immune therapy, which is associated with suppressed antigen presentation responses and insensitivity to Type II interferon (IFN-γ) found in immune function analysis. Interestingly, existing literature indicates that IFN-γ is a critical factor determining the success of immune therapy, primarily through activating immune cells and promoting antigen presentation in tumor cells. However, tumor cells may hijack this pathway by disrupting the IFNG/IFNGR/JAK/STAT signaling pathway or preferentially expressing certain interferon-stimulated genes (ISGs), increasing resistance to IFN-γ and thus circumventing immune system attack [59]. Understanding this mechanism not only strengthens the biological foundation of our prognostic model but also highlights the key challenges that need to be overcome in immune therapy strategies. Further research into these resistance mechanisms will aid in the development of new treatment methods, potentially including inhibitors targeting specific signaling pathways or strategies to modulate ISG expression, to enhance the responsiveness of HCC patients to immune therapy.

Furthermore, we analyzed the spatial expression pattern of UQCRH in HCC using spatial transcriptomics. The results demonstrated that high expression of UQCRH showed significant spatial co-localization with malignant cells, suggesting that this gene may play an important role in the tumor microenvironment of HCC. Through joint analysis of two different HCC samples, we found that UQCRH exhibited elevated expression in regions enriched with malignant cells, further supporting its association with tumor progression. These findings reveal the critical role of UQCRH in HCC, particularly in the spatial distribution within the malignant cell microenvironment, indicating that UQCRH may be involved in tumor initiation and progression. Future research should further explore the specific mechanisms of UQCRH in the tumor microenvironment, including its functional differences in TAMs and tumor cells, and how these differences influence HCC progression and therapeutic responses.

There are several limitations in this study. First, our study is based on data from public databases to construct the prognostic model, and further prospective clinical studies are needed to verify the authenticity and stability of this model. Second, although a prognostic model has been developed, its validation remains limited. So far, we have only performed a preliminary analysis of UQCRH expression in non-HCC and HCC cells using Western blot. Future experiments are needed to explore the mechanisms of action of this model on targeted drug resistance, immune therapy responses, and immune function. Additionally, the spatial transcriptomics analysis used publicly available datasets with a small sample size, and future studies should incorporate more spatial omics data for further validation. In constructing the prognostic model, we used Lasso regularization analysis to avoid overfitting, which is particularly suitable for high-dimensional data with small sample sizes. While Elastic Net provides a more balanced feature selection and coefficient shrinkage, Lasso was sufficient for our needs due to the limited multicollinearity in the data. Future studies may explore Elastic Net as a complementary method to further validate the robustness of the results. Moreover, although the UQCRH single-gene prognostic model demonstrated good performance for 1-year survival prediction (AUC = 0.735), its predictive performance for 3-year and 5-year survival (AUC = 0.632 and 0.607, respectively) still requires improvement. To enhance predictive performance, we integrated UQCRH risk scores with key clinical variables and constructed a comprehensive nomogram, which helps improve the predictive accuracy of the model to some extent.

Conclusion

Through scRNA-seq, we constructed a comprehensive single-cell transcriptomic atlas of HCC, revealing the heterogeneity of monocytes and macrophages in intercellular communication and transcription factor regulatory networks. By integrating scRNA-seq with transcriptomic analysis, we identified prognostic features for predicting overall survival in HCC patients. Spatial transcriptomics further revealed UQCRH’s spatial expression patterns in HCC, showing significant co-localization with malignant cells, suggesting its role in the tumor microenvironment and progression. Our findings demonstrate that the UQCRH-centered prognostic risk model provides independent prognostic value for survival prediction in HCC. This research also advances our understanding of targeted drug resistance, immune therapy responses, and immune function in HCC, warranting further validation through clinical and basic research.