B-cell-specific signatures reveal novel immunophenotyping and therapeutic targets for hepatocellular carcinoma

Abstract

BACKGROUND

Immunotherapy presents both promises and challenges in treating hepatocellular carcinoma (HCC) due to its complex immunological microenvironment. The role of B cells, a key part of the immune system, remains uncertain in HCC.

AIM

To identify B-cell-specific signatures and reveal novel immunophenotyping and therapeutic targets for HCC.

METHODS

Using the Tumor Immune Single-cell Hub 2 database, we identified B-cell-related genes (BRGs) in HCC. Gene enrichment analysis was performed to explore the possible collaboration between B cells and T cells in HCC. We conducted univariate Cox regression analysis using The Cancer Genome Atlas liver HCC collection dataset to find BRGs linked to HCC prognosis. Subsequently, least absolute shrinkage and selection operator regression was utilized to develop a prognostic model with 11 BRGs. The model was validated using the International Cancer Genome Consortium dataset and GSE76427.

RESULTS

The risk score derived from the prognostic model emerged as an independent prognostic factor for HCC. Analysis of the immune microenvironment and cell infiltration revealed the immune status of various risk groups, supporting the cooperation of B and T cells in suppressing HCC. The BRGs model identified new molecular subtypes of HCC, each with distinct immune characteristics. Drug sensitivity analysis identified targeted drugs effective for each HCC subtype, enabling precision therapy and guiding clinical decisions.

CONCLUSION

We clarified the role of B cells in HCC and propose that the BRGs model offers promising targets for personalized immunotherapy.

Key Words: B cell; Hepatocellular carcinoma; Immune microenvironment; Immunotherapy; Molecular subtype

Core tip: We have established a reliable B-cell-related genes (BRGs) prognostic model and novel molecular subtypes in hepatocellular carcinoma (HCC). The BRGs model revealed the immune status and personalized treatment options of HCC molecular subtypes. B cells may play a role in tumor cytotoxicity by activating CD4⁺ and CD8⁺ T cells in HCC.

INTRODUCTION

Primary liver cancer is the eighth most common cause of cancer-related mortality, ranking third among all cancer-related deaths[1]. Hepatocellular carcinoma (HCC) constitutes 80% of all primary liver cancer cases. In 2019, approximately 747000 cases of HCC were reported globally, marking a 70% increase since 1990, with 480000 deaths attributed to HCC[1]. Given the malignancy of HCC, compounded by its central location within the vascular system, early systemic dissemination and metastasis are common, resulting in high recurrence rates even with curative surgery and combined radiochemotherapy[1,2]. Immunotherapy, purportedly capable of activating T cells or rejuvenating immunosurveillance against cancer[3], is regarded as the most promising approach for curing HCC[4,5]. Several multikinase inhibitors and antivascular endothelial growth factor receptor 2 antibodies have gained approval for advanced HCC treatment[6]. Nonetheless, immunotherapy exhibits only a modest response rate of 15%–30% among HCC patients[7]. This could be attributed to the current research predominantly focusing on T cells, leaving our understanding of the overall tumor immune microenvironment in HCC limited, thereby lacking specificity in identifying immunotherapeutic targets tailored to individual HCC patients. Consequently, there is an urgent need to elucidate the roles of other crucial immune cells such as B cells in HCC, facilitating improved immune stratification of patients for precision therapy, thus elevating the cure rate of HCC.

Currently, immunotherapy for tumors predominantly focuses on T-cell-mediated adaptive immunity. While B cells consistently represent a rich cellular component within tumors, little is known about their activation status and biological functions in human tumors, thus leading to a significant oversight of B cells in tumor immunotherapy[8]. Traditionally associated with humoral immune responses against viral and bacterial infections, the role of B cells in tumor immunity has long been debated and may exhibit cancer-specific characteristics[9]. The role of B cells in HCC remains obscure, yet mounting evidence from both preclinical and clinical studies suggests that B cells can induce potent anticancer immunity through humoral and cellular immune responses[8]; a hypothesis awaiting formal validation through analysis of HCC molecular features. B cells infiltrating tumors can be identified at various stages of HCC development, with their presence varying according to stage and histological subtype[10]. Given their influence on both humoral and cellular immunity, a deeper understanding of B-cell biology can offer significant opportunities for HCC immunotherapy[11]. Therefore, a thorough analysis of the specific role of B cells in HCC, along with their interplay with other crucial immune cells in HCC, is crucial for reshaping the immune landscape of the tumor microenvironment and treating HCC effectively.

In this study, we identified differentially expressed B-cell-related genes (BRGs) in HCC using the Tumor Immune Single-Cell Hub 2 (TISCH2) database. Gene enrichment analysis suggested potential synergy between B cells and T cells in HCC. Subsequently, using The Cancer Genome Atlas Liver Hepatocellular Carcinoma Collection (TCGA-LIHC) dataset, we conducted single-factor Cox regression analysis to screen BRGs associated with HCC prognosis. We used least absolute shrinkage and selection operator (LASSO) regression analysis to construct a prognostic evaluation model composed of 11 BRGs. The risk score derived from the BRG prognostic model was validated as an independent prognostic factor for HCC, confirmed by the International Cancer Genome Consortium (ICGC) dataset and GSE76427. Extensive clinical correlation analysis demonstrated the association between the prognostic model and clinicopathological factors. Immune microenvironment and immune cell infiltration analyses revealed the immune status of different risk score groups, supporting the notion of B cells and T cells synergistically inhibiting HCC. Finally, through the BRGs model, we characterized three novel molecular subtypes of HCC, delineating their immune characteristics and identifying targeted drugs effective for each subtype, thus facilitating precision therapy for HCC and aiding clinical decision-making.

MATERIALS AND METHODS

Selection of differentially expressed BRGs in HCC

We downloaded single-cell RNA data of HCC, comprising 62 530 cells, from the Gene Expression Omnibus website GSE140228 dataset. We analyzed the single-cell data using the Transcriptomic-Immunohistochemical Single-Cell Hybridization (TISH2) website[12]. We used principal component analysis (PCA) for dimensionality reduction, K-nearest neighbors, and Louvain algorithms for cluster identification, and annotated cell types using marker genes[13,14]. By applying the Wilcoxon test, we identified differentially expressed genes in B-cell clusters compared to all other cells based on the logarithmic fold change (|fold change| ≥ 1.5) and false discovery rate (FDR) < 0.05[15].

Cell–cell communication analysis

We conducted Cell Chat analysis[16] using the TISH2 website[12] to assess cell–cell communication based on the expression of known ligand–receptor (L–R) pairs across different clusters. We utilized the netVisual_circle function from the pheatmap R package and the CellChat R package to compute and display the quantity of significant L–R interaction pairs and the communication probability between two clusters. Significant L–R interaction pairs were identified for each cluster, designating them as either source or target cells, with a P value threshold set at 0.05[12].

Functional enrichment analysis of different cell clusters

To characterize the functionality of distinct cell type populations, we conducted gene set enrichment analysis (GSEA) using the TISH2 website, ranking genes based on fold changes derived from differential analysis[12]. We used the Kyoto Encyclopedia of Genes and Genomes (KEGG), Gene Ontology (GO), and GSEA analyses to identify and visualize significantly upregulated and downregulated pathways within each cell cluster (FDR ≤ 0.05), facilitating functional enrichment analysis across different clusters[12].

Collection of HCC transcriptomic and clinicopathological data

We obtained transcriptomic data from the TCGA-LIHC dataset, comprising 375 HCC samples and 50 paired adjacent normal tissues. Among these, 365 patients had complete survival data, clinicopathological information, and somatic mutation data. Additionally, we collected transcriptomic and survival data from 232 HCC patients from the ICGC database. Transcriptomic and survival data from 115 HCC patients were retrieved from the GSE76427 dataset.

Construction and validation of prognosis-related BRGs model in HCC

Using the TCGA-LIHC database as the training set, we conducted univariate Cox regression analysis (FDR < 0.05) to identify BRGs associated with HCC prognosis. LASSO regression analysis was performed using the R glm Sparse Net package[17] to mitigate model overfitting, resulting in the derivation of an HCC prognosis model composed of 11 BRGs. The BRGs incorporated into the prognosis model and their corresponding gene coefficients are presented in Supplementary Table 1. Subsequently, we computed the risk scores for each HCC sample in TCGA-LIHC. Based on a 1:1 cutoff value, TCGA-LIHC samples were split into low- and high-risk groups. We used Kaplan–Meier survival and receiver operating characteristic (ROC)[17] analysis in R survival and R time ROC to evaluate the predictive accuracy of the model for HCC prognosis. The reliability of the BRGs model was validated using the ICGC dataset and GSE76427 as validation datasets.

Functional enrichment analysis of different risk score groups

Differential gene analysis between high-risk and low-risk groups in TCGA-LIHC was conducted using the limma package in R (|log fold change| > 1, FDR < 0.05). We used cluster Profiler[17] in R for KEGG and GO analysis to explore molecular mechanisms among HCC risk score groups. GSEA assessed biological functional changes between low- and high-risk subgroups (|normalized enrichment score| > 1, FDR < 0.05).

Assessing the clinical relevance of BRGs model

To demonstrate that the risk score was an independent prognostic factor for HCC, we conducted univariate Cox regression analysis (FDR < 0.05) and multivariate Cox regression analysis (FDR < 0.05) between the risk score and other clinicopathological factors of HCC using the survival package in R. We utilized the survminer[17] and time ROC[17] packages in R to generate ROC curves, illustrating the comparative efficacy of the risk score and other clinicopathological factors for HCC prognosis. Visualization of the correlation between the risk score and other clinicopathological factors of HCC was achieved using the Complex Heatmap and reshape2 packages in R.

Evaluation of the ability of the BRGs model to characterize the immune microenvironment of HCC

Various algorithms, including tumor immune estimation resource, single-sample GSEA (ssGSEA), microenvironment cell population counter, QUANTitative immune single-cell expression, estimation, cell-type identification by estimating relative subsets of RNA transcripts (CIBERSORT), and evaluating the proportion of immune cells[18], were used to calculate the relative levels of tumor-infiltrating immune cells in TCGA-LIHC samples. ssGSEA and CIBERSORT algorithms were utilized to compute the abundance of cell subpopulations in TCGA-LIHC samples. The ssGSEA algorithm was applied to assess the immune functional status in TCGA-LIHC samples.

Assessment of response to immunotherapy

We used the R packages limma and reshape 2[18] to identify differential expression of immune checkpoint genes across various risk score groups. Elevated expression levels of these genes may correlate with improved efficacy of immune checkpoint inhibitors. For the evaluation of immunotherapy effectiveness, a prognostic model based on Biomarker Response to Immunotherapy Group was applied to 348 malignant tumor patients undergoing immunotherapy within the IMvigor210 dataset.

Identification of novel molecular subtypes of HCC

The R package consensus cluster plus[18] was used for unsupervised consensus clustering to identify novel molecular subtypes of HCC. This package generated several key outputs: the consensus matrix (CM), cumulative distribution function (CDF) plot, and consensus heatmap were essential tools for identifying the ideal number of clusters for classifying HCC subtypes. The CDF plot provided a glimpse into the stability of clustering outcomes across various cluster numbers. The CM, displayed as a matrix, indicated the frequency with which sample pairs were grouped together across iterations, providing a quantitative measure of clustering robustness. The consensus heatmap visually represented the CM, facilitating a clearer interpretation of the clustering results.

Statistical analysis

All statistical analyses were performed using R version 4.3.1. The Kruskal–Wallis test was utilized to examine variances in immune scores, immune checkpoint gene expression, and drug sensitivity across different clusters. To compare the differences in patient survival rates between two risk groups, a log-rank test from the R survival package was used for Kaplan–Meier analysis. The statistical tests were two-sided, where P < 0.05 was considered significant. Statistical significance levels are denoted with asterisks: ^aP < 0.05, ^bP < 0.01, and ^cP < 0.001.

RESULTS

B-cell-related immune microenvironmental crosstalk in HCC

The entire study followed the screening and analysis process outlined in Figure 1. We analyzed scRNA data from GSE140228 using TISCH2, and Figure 2A clearly shows the distinct segregation of B-cell subgroups. Both the pie chart and bar graph depict the number and proportion of B cells, demonstrating that they fell within the normal range (Figure 2B and C). In the tumor microenvironment, cellular interactions regulated cell function and influenced the surrounding immune environment and tumor progression. By using the cell chat algorithm in TISCH2, we predicted interactions between different cell types. The analysis showed strong interactions among B cells, most T cells, macrophages, and dendritic cells (Figure 2D and E), suggesting that B cells play a key role in the crosstalk within the HCC immune microenvironment. We examined the gene pairs involved in the interactions of B cells as donors (Figure 2F) or recipients (Figure 2G) with other cells. This revealed the probability of interactions among specific gene pairs, demonstrating that B cells significantly regulate various crucial genes of CD8 T cells such as CD8A and CD8D. To identify BRGs compared with genes of other cells, we utilized Wilcoxon tests available in TISCH2. We applied a log transformation fold change (|fold change|≥ 1.5) and FDR < 0.05 to recognize BRGs. The analysis identified 74 upregulated and 171 downregulated BRGs.

Figure 1 Detailed flowchart illustrates the identification and validation of the prognostic model associated with B-cell-related genes. HCC: Hepatocellular carcinoma; TCGA: The Cancer Genome Atlas; BRGs: B-cell-related genes; TISCH2: Tumor Immune Single-Cell Hub 2; ROC: Receiver operating characteristic; KEGG: Kyoto Encyclopedia of Genes and Genomes; GO: Gene Ontology; GSEA: Gene set enrichment analysis; LASSO: Least absolute shrinkage and selection operator; ICGC: International Cancer Genome Consortium; PCA: Principal component analysis; tSNE: T-distributed stochastic neighbor embedding; GSE: Gene Expression Omnibus Series.

Figure 2 Cell-to-cell communication in hepatocellular carcinoma. A: Cellular map illustrates the distribution and abundance of various cell subgroups within hepatocellular carcinoma (HCC); B: Pie chart displays the accurate cell counts of each cell subgroup in HCC; C: Bar graph demonstrates the proportional distribution of each cell subgroup among individual HCC patients; D: Heatmap exhibits the number of gene pairs interacting between cell groups; E: Using cell chat, the interaction probability between specific cell groups and other cell groups is depicted; F: Using cell chat, the interaction probability is delineated between B cells as donors and specific gene pairs of other cells; G: Using cell chat, the interaction probability is delineated between B cells as recipients and specific gene pairs of other cells; H: Volcano plot showcases the differentially expressed genes within B cells. Red indicates fold change > 1.5 and FDR < 0.05; green indicates fold change < 1.5 and FDR < 0.05. LIHC: Liver hepatocellular carcinoma; NK: Nature kill; DC: Dendritic cells; HLA: Human leukocyte antigen.

Functional characteristics of B cells in HCC

Various transcription factors play critical roles in the development and progression of HCC. We used the local indicators of spatial association algorithm within TISCH2 to infer the transcription factors regulating gene expression in each cell cluster. Using heatmaps, we demonstrated the expression patterns of key transcription factors across all cell clusters in the dataset. The expression patterns of transcription factors in B cells and T cells were almost identical (Figure 3A), supporting the potential synergistic immune role of B cells and T cells in HCC. Subsequently, we sought to identify cell-type-specific transcription factors. In B cells, transcription factors such as myelocytomatosis proto-oncogene, signal transducer and activator of transcription 5B, plant homeodomain finger protein 8, core-binding factor subunit, and E26 transformation-specific-related genes, NRF1, FLI1, EGR3, TFAP4, and RUNX1 showed specific expression (Figure 3B). This specificity highlighted the unique regulatory mechanisms at play within B cells in the HCC immune environment.

Figure 3 Analysis of B-cell-related immune signatures in hepatocellular carcinoma. A: Heatmap illustrates key transcription factors differentially expressed in various cells within hepatocellular carcinoma (HCC); B: Dot plot displays significantly expressed transcription factors in B cells of HCC; C: Heatmap shows immunological gene sets positively regulated by each cell subgroup; D: Heatmap displays immunological gene sets negatively regulated by each cell subgroup; E: Heatmap exhibits Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment pathways positively regulated by each cell subgroup; F: Heatmap exhibits KEGG enrichment pathways negatively regulated by each cell subgroup; G: Heatmap displays gene ontology biological process (GOBP) enrichment pathways positively regulated by each cell subgroup; H: Heatmap displays GOBP enrichment pathways negatively regulated by each cell subgroup; I: Heatmap displays gene ontology molecular function (GOMF) enrichment pathways positively regulated by each cell subgroup; J: Heatmap displays GOMF enrichment pathways negatively regulated by each cell subgroup; K: Heatmap displays gene ontology cellular component (GOCC) enrichment pathways positively regulated by each cell subgroup; L: Heatmap displays GOCC enrichment pathways negatively regulated by each cell subgroup.

To elucidate the specific functions that B cells may regulate in HCC and their synergistic interaction with T cells, we utilized TISCH2 GSEA along with a collection of 16 626 gene sets from the molecular signatures database to characterize the unique functions of B cells. These gene sets encompassed various domains including KEGG pathways, hallmark gene sets, GO terms, immunological features, oncogenic signatures, and transcription factor targets. We observed a strong correlation between B cells and multiple T-cell-related immunological gene sets (Figure 3C and D). In terms of KEGG enrichment analysis, B cells and T cells exhibited highly consistent regulation of pathways such as T-cell receptor signaling, natural-killer-cell-mediated cytotoxicity, NOD-like receptor signaling, Toll-like receptor signaling, and ribosome pathways (Figure 3E and F). Regarding GO analysis, both B cells and T cells showed significant enrichment in pathways including nuclear transcribed mRNA catabolic processes involved in nonsense-mediated decay, translational elongation, RNA catabolic processes, presynaptic membranes, and cytosolic ribosomes (Figure 3G–L). These findings collectively suggest that B cells and T cells jointly regulate mRNA metabolism, mRNA translation, and subsequent immune responses in HCC.

Construction of BRGs prognostic model of HCC

Based on the BRGs selected from Figure 1H, we extracted the expression levels of these BRGs from the TCGA-LIHC database, and merged them with clinical, pathological, and prognostic data, resulting in an integrated table containing complete information for 365 patients. Through single-factor Cox regression analysis, we identified 41 BRGs associated with prognosis of HCC patients (Figure 4A), displaying their expression levels in HCC compared to paired adjacent cancer samples (Figure 4B). We used LASSO regression analysis to mitigate model overfitting, and derived an HCC prognostic model composed of 11 BRGs (Figure 4C and D, Supplementary Table 1). Utilizing Kaplan–Meier survival estimates, we evaluated the relationship between these BRGs and HCC prognosis, revealing that high expression levels of BLNK, FYN, and KLRB1 were associated with favorable prognosis, whereas elevated expression of CALM1, LEPROTL1, S100A10, MARCKSL1, LDHA, LYAR, CKLF, and TXN correlated with poor prognosis (Figure 4E). Additionally, we presented the expression of these BRGs in B cells of HCC, noting high expression of BLNK and MARCKSL1 in B cells, while FYN, KLRB1, CALM1, LEPROTL1, S100A10, LDHA, LYAR, CKLF, and TXN exhibited low expression (Figure 4F).

Figure 4 Establishment of prognostic model of B-cell-related genes in hepatocellular carcinoma. A: Forest plot displays B-cell-related genes (BRGs) associated with hepatocellular carcinoma (HCC) prognosis, green indicates hazard ratio < 1, red indicates hazard ratio > 1; B: Heatmap illustrates the expression levels of BRGs associated with HCC prognosis in both HCC and paired adjacent non-tumor samples; C: Least absolute shrinkage and selection operator (LASSO) regression analysis cross-validation curve demonstrates the fitting effect of the model when incorporating different numbers of BRGs; D: LASSO coefficient path plot displays the lambda values when the model incorporates different numbers of BRGs; E: Survival analysis depicts the relationship between the expression levels of 11 BRGs included in the prognostic model and the prognosis of HCC patients; F: Single-cell analysis reveals the expression levels of the 11 BRGs included in the prognostic model within B cells of HCC. ^aP < 0.05; ^bP < 0.01; ^cP < 0.001.

Validation of BRGs prognostic model for HCC

To evaluate the prognostic model of BRGs, we utilized TCGA-LIHC data as the training set, with ICGC and GSE76427 serving as validation sets. The expression patterns of the 11 BRGs in the prognostic model were consistent between the training and validation sets (Figure 5A–C). The distribution of HCC patients across different risk score groups (Figure 5D–F) and their survival status (Figure 5G–I) were largely consistent across various risk score groups. We depicted the prognostic performance of patients in different risk score groups within the training and validation sets. Consistently, patients in the high-risk score group exhibited poorer prognosis compared with those in the low-risk score group (Figure 5J–L). We constructed ROC curves to validate the effectiveness of the model. The area under the ROC curve (AUC) for 1, 2, and 3 years in the TCGA-LIHC dataset was 0.758, 0.768, and 0.787 respectively (Figure 5M); 0.787, 0.707, and 0.701 in the ICGC database (Figure 5N); and 0.599, 0.655, and 0.637 in GSE76427 (Figure 5O). These results collectively affirm the efficacy of our BRG prognostic model in predicting the prognosis of HCC patients.

Figure 5 Validation of B-cell-related genes prognostic model in hepatocellular carcinoma. A–C: Heatmap displays the expression levels of the 11 B-cell-related genes (BRGs) included in the prognostic model in The Cancer Genome Atlas (TCGA) (A); International Cancer Genome Consortium (ICGC) (B); Gene set enrichment (GSE) 76427 (C); D–F: Risk curves depict the distribution of hepatocellular carcinoma (HCC) patients with different risk scores in TCGA (D); ICGC (E); GSE76427 (F); G–I: Dot plots show the distribution of the survival status of HCC patients with different risk scores in TCGA (G); ICGC (H); GSE76427 (I); J–L: Survival curves illustrate the prognosis of HCC patients with different risk scores in TCGA (J); ICGC (K); GSE76427 (L); M–O: Receiver operating characteristic curves demonstrate the predictive efficacy of the BRGs prognostic model in TCGA (M); ICGC (N); GSE76427 (O); AUC: Areas under the curves.

Clinical correlation analysis of BRGs prognostic model for HCC

To assess the guiding value of the BRGs prognostic model in clinical practice, we conducted Cox regression analysis to compare its predictive superiority regarding HCC patient prognosis with other clinicopathological factors. Single-factor Cox regression analysis revealed that the stage, T stage, M stage, and risk score of HCC were significant factors associated with poor prognosis, with risk score exhibiting the highest hazard ratio (Figure 6A). Multifactor Cox regression analysis indicated that risk score was the sole independent prognostic factor for HCC, suggesting its potential for accurate stand-alone diagnosis of HCC (Figure 6B). Compared with clinically common factors such as age, gender, grade, and stage in assessing HCC prognosis, risk score demonstrated significantly higher AUC values, underscoring the substantial advantages and efficacy of the BRGs prognostic model in evaluating HCC prognosis (Figure 6C). We also assessed the correlation between risk score and clinicopathological factors, revealing that HCC with high risk score had higher stages and T stages, while risk score showed no significant correlation with age, gender, grade, M stage, or N stage of HCC patients (Figure 6D and E). We further evaluated the value of the BRGs prognostic model in different clinical subgroups by stratifying HCC patients into different clinical subsets based on age, gender, grade, and stage. The results demonstrated that patients in the high-risk score group exhibited poorer prognosis across all clinical subsets (Figure 6F–I). These findings underscore the significant value and potential of the BRGs prognostic model in clinical diagnosis and predicting the prognosis of HCC patients.

Figure 6 Clinical correlation analysis of B-cell-related genes prognostic model in hepatocellular carcinoma. A: Single-factor Cox regression analysis shows the association between risk score and other clinicopathological factors of hepatocellular carcinoma (HCC) with prognosis; B: Multifactor Cox regression analysis demonstrates whether risk score and other clinicopathological factors of HCC can serve as independent prognostic factors for HCC; C: Receiver operating characteristic curve illustrates the effectiveness of risk score and other clinicopathological factors of HCC in predicting the prognosis of HCC patients; D: Heatmap displays the distribution of various clinical and pathological statuses of HCC patients in different risk score groups; E: Bubble plot demonstrates the proportions of various clinical and pathological statuses among HCC patients in different risk score groups and their correlation with risk score; F–I: Survival analysis shows the ability of the B-cell-related genes prognostic model to predict the prognosis of HCC patients in various clinical subgroups. AUC: Areas under the curves.

Functional enrichment in different risk score groups for HCC

To delve into the distinct features of various risk score groups and elucidate the potential role of B cells in regulating HCC, we conducted differential gene expression analysis on different risk score groups within the TCGA-LIHC dataset. We identified 420 upregulated and 28 downregulated genes (Figure 7A), with the top 50 differentially expressed genes in TCGA-LIHC displayed (Figure 7B). Subsequently, KEGG enrichment analysis was performed to identify differential pathways within different risk score groups. Pathways such as extracellular matrix (ECM)–receptor interaction, central carbon metabolism in cancer, protein digestion and absorption, hypoxia-inducible factor (HIF)-1 signaling pathway, biosynthesis of amino acids, and fructose and mannose metabolism were significantly enriched with differential genes (Figure 7C). Further GO analysis suggested that B cells were involved in crucial cellular molecular functions and biological processes in HCC, including ECM organization, apical part of cells, and G-protein-coupled receptor binding (Figure 7D and E). Lastly, GSEA was conducted to validate the significance of the enriched pathways. Results indicated that in the high-risk score group, KEGG cell cycle and KEGG ECM–receptor interaction were potentially important enriched pathways (Figure 7F), whereas in the low-risk score group, it was KEGG butanoate metabolism and KEGG fatty acid metabolism (Figure 7G).

Figure 7 Functional enrichment analysis of different risk score groups. A: Heatmap displays the expression levels of the top 50 differentially expressed genes across different risk score groups in hepatocellular carcinoma (HCC) samples; B: Volcano plot illustrates the differentially expressed genes across different risk score groups. Red indicates fold change > 2, false discovery rate (FDR) < 0.05; green indicates fold change < 2, FDR < 0.05; C: Bubble plot presents the Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis results across different risk score groups; D: Bubble plot demonstrates gene ontology (GO) enrichment analysis pathways and proportions across different risk score groups; E: Bubble plot exhibits GO enrichment analysis pathways, number of genes, and their statistical significance across different risk score groups; F: Gene set enrichment analysis (GSEA) analysis reveals the key pathways enriched in HCC patients with high risk scores; G: GSEA analysis unveils the key pathways enriched in HCC patients with low risk scores; H: Somatic mutation data analysis showcases the mutation rates of key genes in HCC patients with high risk scores; I: Somatic mutation data analysis displays the mutation rates of key genes in HCC patients with low risk scores; J: Survival analysis illustrates the prognosis of HCC patients in different tumor mutation burden groups; K: Survival analysis demonstrates the impact of tumor mutation burden combined with risk score groups on the prognosis of HCC patients. FC: Fold change; ECM: Extracellular matrix; HIF: Hypoxia-inducible factor; KEGG: Kyoto encyclopedia of genes and genomes.

The mutation of key genes in tumor cells is closely associated with the progression and prognosis of HCC[19]. We sought to determine whether the functionality of B cells in HCC correlates with this. Patients with HCC in the high-risk score group exhibited significantly higher somatic mutation rates in key genes, such as TP53 (33% vs 19%) and TTN (27% vs 20%), compared with those in the low-risk score group (Figure 7H and I). HCC patients with a high tumor mutation burden (TMB) had significantly poorer prognosis than those with a low TMB (Figure 7J). When evaluating the prognosis of HCC patients using both TMB and risk score, those with high TMB and high risk scores had the worst prognosis. Conversely, those with low TMB and low risk scores had the best prognosis. Patients with either high TMB or high risk scores alone had prognoses that fell between these extremes (Figure 7K). These results suggest that the functionality of B cells in HCC is associated with somatic mutations and that the BRGs prognosis model, when combined with TMB, can enhance the prediction of prognosis in HCC patients.

Assessment of immune microenvironment in different risk score groups of HCC

B cells, as crucial immune cells, are likely to play a significant role in the regulation of the immune microenvironment in HCC. Regarding immune microenvironment analysis, patients with HCC in the high-risk score group showed a trend towards decreased stromal score (P = 0.063), immune score (P = 0.056), and estimate score (P = 0.051), although the lack of significant differences may be attributed to the small sample size (Figure 8A–C). Concerning immune cell infiltration, there was a significant negative correlation between risk score and B-cell infiltration, with various types of T cells showing a similar trend to B cells (Figure 8D). ssGSEA algorithm or CIBERSORT algorithm both revealed significantly suppressed B-cell infiltration in patients with high risk scores, with T cells exhibiting a similar infiltration pattern (Figure 8E–G). These findings align with those shown in Figures 2 and 3, indicating a probable collaboration between B cells and T cells in HCC. Based on the ssGSEA algorithm, we also observed significant suppression of functions related to cytolytic activity, inflammation promotion, T-cell co-inhibition, and T-cell co-stimulation in high-risk patients (Figure 8H), further corroborating our standpoint.

Figure 8 Assessment of the immune microenvironment in different risk score groups. A: Box plot illustrates the stromal score across different risk score groups; B: Box plot displays the immune score across different risk score groups; C: Box plot showcases the estimate score across different risk score groups; D: Bubble plot demonstrates the correlation between risk score and immune cell infiltration within hepatocellular carcinoma samples; E: Box plot shows the abundance of immune cells across different risk score groups calculated based on the single-sample gene set enrichment analysis (ssGSEA) algorithm; F: Box plot exhibits the abundance of immune cells across different risk score groups calculated based on the cell-type identification by estimating relative subsets of RNA transcripts (CIBERSORT) algorithm; G: Heatmap displays the distribution and proportions of immune cells across different risk score groups calculated based on the CIBERSORT algorithm; H: Box plot demonstrates the immune functional status across different risk score groups calculated based on the ssGSEA algorithm; I: Association between The Cancer Genome Atlas immune subtypes and BRGs prognostic subtypes; J: Box plot showcases the expression of immune checkpoint genes across different risk score groups; K: Survival analysis depicts the immune therapy response and efficacy in different risk score groups of the IMvigor210 immune therapy cohort. ^aP < 0.05; ^bP < 0.01; ^cP < 0.001. LIHC: Liver hepatocellular carcinoma; CD: Cluster of differentiation; NK: Nature kill; DC: Dendritic cells; HLA: Human leukocyte antigen; TCGA: The cancer genome atlas.

BRGs prognostic model and TCGA official immune typing

The TCGA official team has classified HCC into four subtypes: (1) Wound Healing (Immune C1); (2) Interferon- dominant (Immune C2); (3) Inflammatory (Immune C3); and (4) Lymphocyte depleted (Immune C4), which has become a significant reference for classification[20]. Our findings indicate that the BRGs prognosis model can differentiate immune subtypes, particularly with a notable increase in the number of Immune C1 patients with high-risk scores and a marked decrease in the number of Immune C3 patients (Figure 8I).

BRGs prognostic model evaluates immunotherapy efficacy

Our findings indicate that the BRGs prognosis model can differentiate immune subtypes, particularly with a notable increase in the number of Immune C1 patients with high-risk scores and a marked decrease in the number of Immune C3 patients (Figure 8I). Immune checkpoints, crucial factors expressed on immune cells capable of modulating immune activation levels, play a pivotal role. These genes, under normal circumstances, can inhibit T-cell function to prevent excessive immune system activation. Simultaneously, in tumor tissues, they might be exploited by tumors to facilitate immune evasion. The effectiveness of immune checkpoint inhibitors in treating HCC has been demonstrated[21,22]. Our analysis suggests that in the high-risk score group, expression of most immune checkpoint genes (17 of 25) significantly increased (Figure 8J), indicating that patients in the high-risk score group may benefit from these immune checkpoint inhibitors. ENTPD1, NT5E, and HAVCR2 are also immune checkpoint genes coexpressed at high levels in both B cells and T cells[23], further indicating potential synergistic action between B cells and T cells in HCC. We evaluated the efficacy of immunotherapy in the IMvigor210 immune treatment cohort using the BRGs prognosis model and found that high-risk patients showed significantly poor response to immunotherapy (Figure 8K). This suggests that high-risk group patients may require alternative therapeutic agents to improve efficacy. Therefore, based on drug sensitivity analysis, we have identified potential targeted drugs such as osimertinib, gefitinib, GDC0810, and paclitaxel (Supplementary Figure 1) that may be beneficial for high-risk group patients. These drugs, in combination with immunotherapy, could serve as future references for further validation in clinical experiments.

Novel HCC molecular classification based on 11 BRGs

Current immunotherapy for HCC patients yields only a 15%–30% remission rate, possibly due to inadequate differentiation among patients with distinct immune microenvironments, thus failing to achieve personalized precision treatment. Therefore, identifying novel molecular subtyping for HCC is imperative[7]. Leveraging all 11 BRGs in the prognosis model, we classified HCC patients into molecular subtypes using consensus cluster plus. Based on delta area and CDF curves, we observed good stability in sample clustering when k = 3 (Figure 9A and B). We depicted sample distribution under different k values (Figure 9C). A heatmap of CM demonstrated consistent blue shading when k = 3 (Figure 9D). Subsequently, the C1 group had the poorest prognosis, while those in the C2 and C3 groups exhibited better prognosis based on survival analysis (Figure 9E). A Sankey diagram illustrated the correspondence between novel molecular subtyping and prognosis model grouping (Figure 9F). PCA and t-distributed stochastic neighbor embedding analyses both distinguished HCC patients based on the novel molecular subtyping (Figure 9G and H). We attempted to ascertain whether the novel molecular subtyping could differentiate HCC patients with distinct immune features. Analysis of the immune microenvironment indicated that patients in the C1 group had high estimate scores, stromal scores, and immune scores, whereas those in the C2 group had low scores, with C3 group patients falling between the two (Figure 9I–K). Regarding immune cell infiltration, C1 group patients exhibited markedly elevated levels across various immune cell types, possibly contributing to their poorer prognosis due to immune resource depletion, whereas C2 and C3 group patients showed low levels of immune cell infiltration (Figure 9L). Immune checkpoint analysis revealed significant upregulation of most immune checkpoint genes (42 of 46) in C1 group patients, indicating their potential benefit from immune checkpoint inhibitors. C2 patients might benefit from CD40 inhibitors, while C3 patients might benefit from IDO2, KIR3DL1, and CD160 inhibitors (Figure 9M). We analyzed potential targeted drugs suitable for HCC patients in each immune subtype. C1 group patients were most likely to benefit from staurosporine, pevonedistat, and midkine (MK)-8776, while C2 group patients might benefit from podophyllotoxin bromide and elephantin, and C3 group patients might benefit from epirubicin, extracellular regulated protein kinases (ERK)-2440, and Janus kinase (JAK)-8517 (Supplementary Figure 2).

Figure 9 Novel molecular subtyping for identification of hepatocellular carcinoma based on B-cell-related genes. A: Cumulative distribution function curves for different numbers of classifications, where different curves represent the stability of models at different k values; B: Distribution of samples across different numbers of classifications; C: Heatmap displays distribution of hepatocellular carcinoma (HCC) patients under different k values; D: Consistency clustering plot displays the clustering results when k = 3; E: Survival curves illustrate the prognosis of HCC patients with different molecular subtypes; F: Sankey diagram shows the correspondence between different molecular subtypes of HCC patients and different risk score groups; G: Principal component analysis shows the distribution of samples across different HCC molecular subtypes; H: T-distributed stochastic neighbor embedding analysis shows the distribution of samples across different HCC molecular subtypes; I: Box plot demonstrates the estimate score across different HCC molecular subtypes; J: Box plot showcases the stromal score across different HCC molecular subtypes; K: Box plot displays the immune score across different HCC molecular subtypes; L: Heatmap displays the immune cell infiltration status across different HCC molecular subtypes based on different algorithms; M: Box plot demonstrates the expression of immune checkpoint genes across different HCC molecular subtypes. ^aP < 0.05; ^bP < 0.01; ^cP < 0.001.

DISCUSSION

Although the role of T cells in antitumor immunity has been extensively studied, the role of B cells, including in HCC, remains poorly understood, hindering efforts to utilize B-cell responses for cancer immunotherapy[24]. B cells, as a crucial subset of antigen-presenting cells, play a significant role, which is often overlooked in tumor immunotherapy, including HCC. Unlike dendritic cells and macrophages, which primarily present antigens through phagocytosis, B cells predominantly present antigens through the high-affinity binding of the B-cell receptor (BCR) to the antigen, thereby effectively recognizing low levels of antigens and enhancing antigen presentation efficiency. Upon binding of antigens to specific BCRs, the process of antigen internalization initiates, and the antigens are displayed on the surface of B cells via MHC class II molecules, activating cognate CD4⁺ helper T cells, thus inducing antibody production[25]. Based on gene enrichment analysis of BRGs, B cells and T cells (including CD4⁺ and CD8⁺ T cells) significantly enrich pathways associated with HCC occurrence and development, such as T-cell receptor signaling pathway and natural-killer-cell mediated cytotoxicity. B cells and T cells jointly regulate important biological processes such as GO nuclear transcribed mRNA catabolic process nonsense mediated decay and GO translational elongation. Our analysis reveals a potential synergistic interaction between B cells and T cells in HCC. Our findings suggest that B cells play a crucial role in activating CD4⁺ T cells and may act in conjunction with CD8⁺ T cells. This discovery represents a novel breakthrough in tumor therapy.

The role of B cells in most tumors, including HCC, remains a subject of considerable controversy[25]. Among the 11 key genes in our BRGs prognostic model, BLNK is associated with a favorable prognosis in HCC, exhibiting high expression in B cells within HCC. Conversely, CALM1, LEPROTL1, S100A10, LDHA, LYAR, CKLF, and TXN are linked to an adverse prognosis in HCC, displaying low expression in B cells within HCC. The alignment of these BRGs with the prognostic outcomes of HCC, as reflected in their expression within B cells, suggests that B-cell activation is a favorable prognostic factor for HCC. However, certain BRGs exhibit discrepancies between their overall expression in HCC and their expression specifically within B cells. For instance, FYN and KLRB1 are associated with a favorable prognosis in HCC but demonstrate low expression in B cells within HCC, while MARCKSL1 is linked to an adverse prognosis in HCC yet exhibits high expression in B cells within HCC. These findings suggest potential complexities in the roles of these BRGs, indicating that their functions may extend beyond B cells, warranting further investigation. The high risk score group identified by the BRGs prognostic model shows significant enrichment in pathways such as ECM–receptor interaction, central carbon metabolism in cancer, protein digestion and absorption, HIF-1 signaling pathway, biosynthesis of amino acids, and fructose and mannose metabolism. This implies that B cells may regulate the occurrence and progression of HCC through these pathways, necessitating further validation in future studies.

Similar to regulatory T cells, B cells can also exhibit regulatory phenotypes to suppress immune responses. They express inhibitory receptors on their surface, akin to immune checkpoint genes. These include a variety of immune checkpoint molecules like programmed cell death protein 1 (PD-1) and its ligands PD ligand 1 (PD-L1)/PD ligand 2, which play diverse regulatory roles in both humoral and cellular immunity[26]. Recent findings in HCC reveal that B cells with elevated PD-1 expression, when encountering PD-L1-expressing cells, induce T-cell suppression by releasing interleukin-10. Blocking PD-1 on B cells with checkpoint blockade antibodies has shown promise in enhancing cancer immunotherapy[27]. Anti-PD-1 therapy has become the preferred immunotherapy for treating HCC. However, anti-PD1 therapy targeting T cells has only achieved a 15% response rate in HCC patients[7]. In this context, anti-PD-1 therapy targeting B cells may potentially become an important approach to improving the efficacy of immunotherapy for HCC patients. Our results indicate significant upregulation in expression of most T-cell immune checkpoint genes in the high risk scoring group of HCC patients, with ENTPD1, NT5E, and HAVCR2 being among the immune checkpoint genes highly expressed in B cells[23]. Targeting inhibitors against these shared immune checkpoint genes in both B cells and T cells may effectively harness the synergistic action of these cells in HCC therapy; a hypothesis deserving validation in future multicenter large-scale clinical trials.

Our prognostic model for HCC, termed BRGs, offers guidance in diagnosis and prognosis prediction. Through univariate Cox regression and LASSO regression analysis, we constructed the BRGs prognostic model. Our analysis indicates a significant reduction in overall survival for patients in the high-risk group compared to those in the low-risk group. Furthermore, through visualizations such as risk curves, risk heatmaps, ROC curves, and Kaplan–Meier survival curves, we validated the robust predictive power of our risk model, which outperformed other significant risk factors in HCC such as age, grade, and stage. Importantly, similar robust results were replicated in two external validation sets for HCC (ICGC and GSE76427), bolstering the reliability and reproducibility of our research findings.

Beyond elucidating the role of B cells in HCC, identifying precise molecular subtypes is crucial for enhancing personalized treatment, given the widespread resistance of HCC to targeted therapies[28]. Leveraging 11 BRGs from our prognostic model, we classified HCC patients into three molecular subtypes. Our results indicate that C1 patients exhibited the most extensive immune cell infiltration, potentially indicating immune exhaustion and consequently the poorest prognosis. However, these patients may benefit from targeted drugs such as staurosporine, pevonedistat, and MK-8776. Conversely, C2 and C3 patients showed better prognosis. C2 patients may improve their prognosis with drugs such as podophyllotoxin bromide and elephantin, while C3 patients may benefit from drugs such as epirubicin, ERK-2440, and JAK-8517. Our proposed novel HCC molecular subtypes hold significant implications for precise HCC treatment.

CONCLUSION

Our research suggests that in HCC, B cells play a role in tumor cytotoxicity by activating CD4⁺ and CD8⁺ T cells. We have established a reliable BRGs prognostic model and novel molecular subtypes, offering pivotal and valuable references for the diagnosis, prognosis, and personalized treatment strategies of HCC.