USP21 is involved in the development of chronic hepatitis B by modulating the immune microenvironment

Abstract

Hepatitis B virus (HBV) infection is a global public health challenge that alters the immune microenvironment of the liver and drives disease progression by triggering chronic inflammation that leads to hepatic cell death through multiple programmed cell death (PCD) modalities. Due to the persistence of covalently closed circular DNA in hepatocytes, there is a lack of curative drugs that can completely eradicate HBV. Therefore, revealing how HBV infection leads to changes in the hepatic immune microenvironment, as well as searching for specific molecular targets, is crucial for controlling the onset and progression of chronic hepatitis B (CHB). In this study, we used the single sample gene set enrichment analysis and CIBERSORT algorithms to assess immune cell infiltration in the livers of CHB patients. With three advanced machine learning algorithms, random forest, least absolute shrinkage and selection operator, and selected support vector machine recursive feature elimination, we identified the PCD signature genes associated with CHB from the candidate genes. We further validated that ubiquitin-specific peptidase 21 could differentiate CHB patients with different natural courses by receiver operating characteristic analysis. These findings enhance our understanding of the mechanisms of HBV infection.

Introduction

HBV is a double-stranded DNA virus that persistently infects its host by forming covalently closed circular DNA (cccDNA) in the nucleus [1]. According to the World Health Organization (WHO), approximately 257 million people worldwide have tested positive for hepatitis B surface antigen (HBsAg), indicating that HBV infection remains a serious global public health problem [2].

The liver, as a vital organ in the human body, is constantly exposed to various pathogens from the bloodstream and gastrointestinal tract due to its unique blood supply mechanism, thus playing a role in immune surveillance within the body [3]. Chronic HBV infection not only impairs the immune system’s surveillance function but also affects multiple metabolic pathways in hepatocytes, thereby activating immune responses and leading to changes in the liver’s immune microenvironment [4]. The current research indicates that the host’s immune system plays a crucial role in the progression of chronic HBV infection [5]. For instance, T cell exhaustion is a prominent feature during HBV infection. Exhausted T cells exhibit mitochondrial dysfunction that impairs oxidative phosphorylation (OXPHOS), thereby affecting the antiviral function of T cells [6]. However, the mechanism by which HBV regulates the liver’s immune microenvironment to facilitate the occurrence and progression of chronic infection is not yet fully understood.

Hepatocyte injury and death resulting from HBV infection are associated with various forms of PCD. In the current literature on HBV infection-related diseases, the forms of PCD that have been more extensively studied include autophagy, apoptosis, ferroptosis, and necroptosis. Autophagy is a lysosome-mediated intracellular catabolism process that plays a crucial role in maintaining cellular and metabolic homeostasis within the liver [7]. Activation of hypoxia-inducible factor-1α or high-mobility group box-1 has been reported to activate autophagy in hepatic stellate cells [8], leading to liver fibrosis. Apoptosis is a highly regulated form of cell death characterized by cell shrinkage, membrane blebbing, and the formation of apoptotic bodies [9]. Most studies suggest that HBV infection can suppress apoptosis, which facilitates viral proliferation and promotes disease progression [10, 11]. Ferroptosis is an iron-dependent cell death caused by lipid-based reactive oxygen species accumulation [12]. It has been shown that HBV-induced up-regulation of heat shock protein family A member 8 (HSPA8) promotes hepatocarcinogenesis by inhibiting ferroptosis and stimulating HBV replication [13]. Necroptosis is a novel form of programmed cell death mediated by receptor-interacting protein kinase 1 (RIPK1), phosphoinositol-3-kinase (PI3K), and the mixed lineage kinase domain-like pseudo kinase (MLKL) [14]. Cellular responses can induce apoptosis to eliminate HBV-infected hepatocytes, but sustained antiviral immune responses can exacerbate disease progression [15]. These forms of programmed cell death play a significant role in the pathogenesis of CHB. They are not only involved in the injury and death of hepatocytes but may also influence the liver’s immune microenvironment and inflammatory responses.

In this study, we leveraged the publicly available dataset GSE83148 to perform a comprehensive bioinformatic analysis of RNA sequencing data from CHB patients, employing differential expression analysis, gene enrichment profiling, and machine learning algorithms. Previous studies using GSE83148 have focused on different aspects of CHB pathogenesis, such as Cao et al.’s study of gene expression changes during TiaoGanYiPi (TGYP) therapeutic interventions [16] and Tang et al.’s description of the transcriptional regulatory network that drives hepatic inflammatory outbreaks [17]. However, our study provides new PCD-mediated immune dysregulation insights. We observed significant differences in the infiltration of liver immune cells and the expression of immune-related pathways between patients with CHB and healthy controls (HCs); such differences were also evident among different subtypes of CHB. These differences may be associated with the occurrence and progression of CHB, revealing the complex role of immune responses in the advancement of liver disease. In addition, through the initial analysis of PCD-related genes, using different machine learning methods for dimensionality reduction, we finally targeted ubiquitin-specific peptidase 21 (USP21) as a feature gene for further analysis. USP21 plays an important role in the cell cycle and tumorigenesis, and its role in CHB deserves further investigation.

Materials and methods

Data acquisition

We obtained three gene expression profile datasets related to HBV infection from the Gene Expression Omnibus (GEO) (Table 1). Standardize and normalize each dataset as it is collated. By annotating the GSE83148 dataset, we obtained gene expression matrices and clinical information for each sample and grouped the samples for subsequent analysis, with a workflow diagram as in Fig. S1.

Table 1 Details from the datasets related to patients

Identification of differentially expressed genes

We used the "limma" package in R to calculate differentially expressed genes (DEGs) between CHB and HC samples. DEGs were determined by threshold criteria of LogFC > 0.5 and p < 0.05. Volcano plots and heatmaps were visualized using the "ggplot2" and "pheatmap" R packages.

Gene ontology (GO) and Kyoto Encyclopedia of genes and genomes analysis (KEGG)

To uncover the functional roles of DEGs and their involvement in critical biological pathways, we conducted GO and KEGG pathway enrichment analyses. This analysis was performed using several R packages, including “clusterProfiler,” “ggplot2,” “dplyr,” “ggsci,”and the broader “tidyverse” suite. For the results of GO and KEGG enrichment analyses, a p < 0.05 was considered statistically significant. The GO enrichment analysis results were visualized using circular plots generated in ChiPlot (https://www.chiplot.online/). The KEGG enrichment analysis results were visualized using chord diagrams generated in SangerBox 3.0 (SangerBox—a bioinformatics data analysis toolbox).

Weighted gene co-expression network analysis (WGCNA)

WGCNA is an algorithm that constructs gene clustering modules based on similar gene expression patterns. To identify co-expressed genes associated with CHB, we used the "WGCNA" R package to construct a gene co-expression network from liver samples of HC and CHB. First, we calculated and selected the optimal soft threshold, choosing a value of 9, as it satisfied the scale-free topology condition and exhibited strong average connectivity. Second, we fixed the scale-free topology index R² at 0.85 to further enhance network connectivity. Finally, based on the similarity of gene expression, we constructed a hierarchical clustering tree and grouped highly co-expressed genes into the same module. Through the heatmap of modules and traits, a total of 28 modules were identified. We selected core genes from co-expression modules with a correlation coefficient greater than 0.5 for further analysis. To reveal a high degree of correlation with CHB, we plotted a scatter plot of module membership (MM) versus gene significance (GS) for the core co-expression modules.

Feature selection of core genes via three machine learning methods

We obtained genes associated with 13 types of PCD from the literature [18]. The intersection of PCD-related genes, genes obtained from WGCNA, and DEGs was taken to obtain the candidate genes. We further reduced the dimensionality of the candidate genes using RF, least absolute shrinkage and selection operator (LASSO), and support vector machine recursive feature elimination (SVM-RFE). RF, a randomized algorithm, was implemented using the “randomForest” R package to estimate the number of decision trees and corresponding error rates. The LASSO algorithm was executed via the “glmnet” R package, determining the optimal penalty parameter λ through tenfold cross-validation, yielding seven core genes. SVM-RFE was used for feature selection and identified four core genes using five-fold cross-validation. We intersected the core genes obtained from the three machine learning methods using a Venn diagram and identified three PCD feature genes that are closely associated with CHB.

Single-gene gene set enrichment analysis (GSEA)

Using the R package "clusterProfiler," we conducted single-gene GSEA enrichment analysis to explore the functions and biological significance of the signature genes. We then selected and visualized the five most highly enriched functional pathways, along with the five least enriched, for comparative analysis.

Single sample gene set enrichment analysis (ssGSEA)

We conducted ssGSEA was conducted on CHB and HCs samples using the “GSVA” R package, aiming to investigate immune cell infiltration and function expression. Differences between the groups were visualized using box-line plots created with the “ggpubr” R package. Additionally, we used the “ggplot2” package to explore the correlation of immune cells and immune pathways with the characterized genes and presented our findings through heatmaps.

Consensus clustering of CHB

We utilized the "ConsensusClusterPlus" R package to perform CHB subtyping based on the expression of signature genes. We generated a consistency matrix, a cumulative distribution function (CDF) curve, and a trace plot to identify the optimal number of subtypes. The selection of the optimal clusters was guided by assessing the clarity of blank regions within the blue modules and evaluating the relative changes in the CDF curve and trace plot. To confirm the reliability of the clustering results, the “ggplot2” R package was used to create scatter plots of the samples between subtypes and observe whether the scatter plots were differentiated between subtypes. Additionally, box plots were used to display the expression differences of feature genes across different subtypes.

Immune infiltration analysis of different subtypes of CHB

To investigate the differences in the immune microenvironment between different subtypes of CHB, we estimated the degree of immune cell infiltration using the CIBERSORT algorithm and visualized it with stacked bar histograms and box plots generated by the "ggplot2" R package. Additionally, the associations between signature genes and immune cells were illustrated using lollipop and scatter plots.

Validation cohorts and blood samples

Peripheral blood samples for this study were collected from Qilu Hospital of Shandong University from July 2022 to June 2024. This study comprised a total of 185 participants, including 24 HCs and 161 CHB patients. HCs had normal liver biochemistry and no history of liver disease or alcohol abuse. The inclusion criteria for the study were as follows: (1) HBsAg positivity for more than 6 months; (2) no prior treatment with interferon; and (3) age over 18 years old. The exclusion criteria were as follows: (1) co-infection with other types of hepatitis viruses or human immunodeficiency virus (HIV); (2) liver disease related to alcohol or obesity; and (3) presence of hepatocellular carcinoma.

Fasting blood samples were collected from each subject in the early morning after an 8-h fast. Peripheral blood mononuclear cells (PBMCs) were isolated from the blood using Ficoll-Paque Plus density gradient centrifugation. Total RNA was then extracted using Trizol reagent, following the manufacturer’s instructions, and further synthesized into cDNA. This study has been reviewed and approved by the Medical Ethics Committee of Qilu Hospital of Shandong University, and all participants have signed written informed consent forms. The research process strictly adheres to the principles of the Declaration of Helsinki.

Statistical analysis

All statistical analyses were conducted using R software (version 4.3.0) and GraphPad Prism (version 9.5.1). The Wilcoxon test was used to analyze differences between the two groups, while the Kruskal–Wallis test was employed to assess differences among multiple groups. Correlation analyses were performed using Pearson’s or Spearman’s tests. The diagnostic value of USP21, alanine aminotransferase (ALT), and aspartate aminotransferase (AST) was evaluated using receiver operating characteristic (ROC) curve analysis. p < 0.05 was considered statistically significant.

Results

Differential expression analysis

To explore the DEGs that play a key role in the complex landscape of CHB, we compared the gene expression profiles between CHB and HCs liver tissues using the "limma" package in the R programming language. We successfully identified 1724 DEGs with significant changes in gene expression levels, including 1031 up-regulated genes and 693 down-regulated genes. These insightful findings were visualized through a volcano plot (Fig. 1A) and a heatmap (Fig. 1B).

Fig. 1

Identification of the differential genes. A Heatmap displaying the DEGs. B Volcano plots of DEGs

Functional enrichment analysis of DEGs

To gain a deeper understanding of the biological functions of DEGs, we performed GO and KEGG enrichment analyses on the DEGs. The GO analysis allowed us to comprehensively understand the changes in biological functions caused by gene expression dysregulation in CHB, including biological processes (BP), cellular components (CC), and molecular functions (MF). Histograms (Fig. 2A) and circular plots (Fig. 2B) displayed the top-ranked terms for BP, CC, and MF. In terms of biological processes, DEGs were involved in processes such as immune response, leukocyte activation, regulation of immune system process, immune system process, cell activation, and positive regulation of immune system process, indicating their potentially important roles in regulating immune responses and maintaining immune balance. In terms of cellular components, DEGs were mainly associated with the extracellular region, endomembrane system, extracellular exosome, cell surface, extracellular space, and vesicle functions, which are necessary for cell communication and immune regulation. In terms of molecular functions, DEGs were primarily involved in glycosaminoglycan binding, peptide binding, chemokine receptor binding, amide binding, chemokine activity, and identical protein binding. These findings suggest that the dysregulated genes may play important roles in cell migration, adhesion, differentiation, and immune responses. To further reveal the pathways associated with CHB, we conducted a KEGG enrichment analysis, with key results visualized by Lollipop Plots (Fig. 2C) and chord diagrams (Fig. 2D). These key pathways included cytokine–cytokine receptor interaction, human T cell leukemia virus one infection, intestinal immune network for IgA production, and Th1 and Th2 cell differentiation. The enriched pathways suggest their potentially important roles in immune responses, cell signaling, inflammatory responses, and viral defense.

Fig. 2

Functional enrichment study based on DEGs. A Histogram of the top six terms for GO enrichment analysis. B Circular plot of the top eight terms for GO enrichment analysis. C Lollipop plot of KEGG enrichment analysis. D Chord diagram for KEGG enrichment analysis

WGCNA

In the dataset GSE83148, we employed WGCNA to investigate the connections between key genes in CHB and clinical characteristics [19]. We first screened the CHB and HCs samples, removing outliers and presented the results through a heatmap (Fig. 3A). The ideal soft threshold was determined to be 9 at which the scale-free topology fitting index R^2 exceeded 0.85 (Fig. 3B). After merging related modules, we obtained 28 modules, which were displayed by a clustering tree (Fig. 3C). In the module–clinical trait relationship heatmap, the dark turquoise module exhibited the strongest positive correlation with CHB (r = 0.78) (Fig. 3D). We then selected the dark turquoise and royal blue modules to construct scatter plots of gene-CHB correlations (Fig. 3E, F). We chose module genes with an absolute r value greater than 0.5 for further analysis and presented the correlation plots of each module with the phenotype (Fig. 3G).

Fig. 3

Construction of gene co-expression network. A Cluster tree and heatmap for all samples in GSE83148, the red heatmap corresponds to the clinical characteristics of the samples in that cluster tree. B The optimal soft threshold β = 9 is chosen for the scale-free topological fitting exponent R² = 0.85 in the left panel to obtain the best average connectivity of the co-expression network in the right panel. C Example of a clustering tree for a merged module. D Heatmap of correlations between modules and clinical features. Red indicates a positive correlation; blue indicates a negative correlation, and the darker the color, the stronger the correlation. The number in parentheses is the p-value of the correlation between the module and the feature. The number above the parentheses indicates the correlation size between the module and the feature. E Scatter plot between dark turquoise module and CHB gene significance, cor = 0.51, p < 0.001. F Scatter plot between royal blue module and CHB gene significance, cor = 0.42, p < 0.01. G Correlation between modules and phenotypes

Identification of feature genes

Using WGCNA, we filtered out genes closely associated with CHB phenotypes from a large number of genes then performed an intersection analysis with PCD genes and DEGs, successfully identifying 17 potential candidate genes (Fig. 4A). To further explore PCD signature genes related to CHB, we employed three advanced machine learning techniques, RF, LASSO regression, and SVM-RFE, to conduct an in-depth analysis and dimensionality reduction of these candidate genes. First, we used the RF algorithm (ntree = 1000) to assess the importance of these 17 candidate genes and ranking them accordingly, ultimately selecting the top 10 genes (Fig. 4B). Then, through LASSO regression analysis, we identified seven genes with the minimum binomial deviation among these candidates (Fig. 4C, D). Finally, using the SVM-RFE algorithm for risk assessment, we found that the model achieved the highest accuracy (Fig. 4E) and the lowest error rate (Fig. 4F) when considering four genes. Integrating the results from these three machine learning techniques, we ultimately identified three signature genes: AKT1, POR, and USP21 (Fig. 4G). AKT1 is involved in two important cellular processes: autophagy and apoptosis. The POR gene is associated with the process of ferroptosis. USP21 is involved in necroptosis. Among the programmed cell death genes, these are the most closely related to the occurrence and development of CHB.

Fig. 4

Screening feature genes by different machine learning methods. A Candidate genes are obtained by taking the intersection of Venn diagrams. B Random Forest calculated the top 10 genes in terms of gene importance. C, D Key genes were screened using the LASSO algorithm. E, F The SVM-RFE algorithm screens biometric genes, and the point with the lowest accuracy and error rate is used as the number of key genes screened by SVM-RFE. G The key genes obtained from the three machine learning algorithms are taken as an intersection to get the feature genes. H Demonstration of the location of feature genes in the genome by circus plot

Feature gene expression differences

In the GSE83148 dataset, we observed that the mRNA expression levels of the three signature genes were significantly higher in patients with CHB compared to the HCs group (Fig. 5A). To further validate the expression patterns of these signature genes, we obtained two additional datasets, GSE58208 and GSE65359, from the GEO database. In the GSE58208 dataset, we analyzed gene expression data from PBMCs of 12 CHB patients and 5 HCs individuals. The analysis revealed that the expression levels of these three signature genes remained elevated in CHB patients within PBMCs (Fig. 5B), consistent with our preliminary findings in the GSE83148 dataset. Furthermore, we analyzed the GSE65359 dataset, which included gene expression data from the liver tissues of 83 CHB patients at different immune stages. Our analysis revealed no significant differences in the expression levels of AKT1 and POR genes across different immune stages. However, the expression levels of the USP21 gene showed significant variation across these stages (Fig. 5C), indicating that USP21 plays a crucial role in the immune response of CHB.

Fig. 5

Feature genes expression differences. A Expression of feature genes in the GSE83148 dataset. B Expression of feature genes in the GSE58208 dataset. C Expression of feature genes in the GSE65359 dataset. IC inactive carrier state, ACH active chronic hepatitis, including immune-active hepatitis and hepatitis B e antigen-negative hepatitis, IT immune tolerant

Functional enrichment analyses for feature genes

To uncover the molecular functions and potential biological significance of the three signature genes, we conducted a GSEA. Our analysis revealed that the pathways affected by the POR gene include arginine biosynthesis, glycosaminoglycan biosynthesis—keratan sulfate, other glycan degradation, phenylalanine metabolism, proximal tubule bicarbonate reabsorption, DNA replication, glycosylphosphatidylinositol anchor biosynthesis, oxidative phosphorylation, ribosome, and terpenoid backbone biosynthesis (Fig. 6A). The regulation of these pathways may be closely related to the role of POR in cellular metabolism and signal transduction. Pathways in which the USP21 gene is involved in down-regulation include allograft rejection, asthma, graft-versus-host disease, rheumatoid arthritis, and type 1 diabetes. These results suggest that USP21 may play an inhibitory role in immune responses and inflammatory processes. In addition, up-regulated pathways in which the USP21 gene is involved include arginine biosynthesis, glycine serine and threonine metabolism, glyoxylate and dicarboxylate metabolism, histidine metabolism, and phenylalanine metabolism (Fig. 6B). These results reveal the important role of USP21 in amino acid metabolism and immune regulation. Down-regulated pathways affected by AKT1 include graft-versus-host disease, p53 signaling pathway, rheumatoid arthritis, viral life cycle of HIV-1, and viral protein interaction with cytokine and cytokine receptor. The down-regulation of these pathways may be related to the inhibitory role of AKT1 in cell growth and immune regulation. The up-regulated pathways involved in AKT1 include arginine and proline metabolism, arginine biosynthesis, other glycan degradation, thiamine metabolism, ubiquinone, and other terpenoid-quinone biosynthesis (Fig. 6C). AKT1 may play an active role in cellular metabolism and energy production. Through GSEA analysis, we not only revealed the roles of these three characterized genes in a variety of biological processes but also identified their potential importance in immune regulation and metabolic regulation.

Fig. 6

Functional enrichment analyses for feature genes. A The POR gene acquires the first five positive and negative pathways through GSEA. B The USP21 gene acquired the first five positive and negative pathways through GSEA. C The AKT1 gene acquired the first five positive and negative pathways through GSEA. D Ridge map of the first 20 pathways acquired by the AKT1 gene through GSEA

Analysis of immune cell and immune function infiltration between HCs and CHB

In this study, we utilized the ssGSEA algorithm to depict the differences in immune infiltration between HCs and CHB patients. The significant differences in the infiltration levels of 28 immune cell types were presented through box plots (Fig. 7A), revealing that all significantly differing immune cell infiltration levels were higher in the CHB group compared to the HCs group. Additionally, we analyzed the infiltration of 17 immune pathways, and the results were consistent with those related to immune cell infiltration (Fig. 7C), further confirming the significant differences in immune status between HCs and CHB. In further gene expression analysis, we examined the expression of the three signature genes in HCs and CHB patients. The results showed that the expression of these three signature genes was decreased in CD56bright natural killer cells and eosinophils, while it was significantly increased in monocytes (Fig. 7B). Lastly, among the 17 immune functions associated with the three signature genes, we found that POR was highly expressed in cytokine receptors and the BCR signaling pathway. At the same time, USP21 was significantly reduced in natural killer cell cytotoxicity, interleukins, interferons, and chemokines pathways (Fig. 7D). These findings reveal specific expression patterns in certain immune cell subsets and immune pathways.

Fig. 7

Immune infiltration analysis between HCs and CHB. A Box plots of the difference in the infiltration of 28 immune cells between HCs and CHB. B Heatmap of the correlation between feature genes and 28 immune cell infiltrates in CHB samples, with red indicating positive correlation and blue indicating negative correlation. C Box plots of the difference in expression of 17 immune functions between HC and CHB. D Heatmap of the correlation between feature genes and 17 immune functions in CHB samples, with red indicating positive correlation and blue indicating negative correlation. *p < 0.05, **p < 0.01, ***p < 0.001, ns means p > 0.05

Consensus clustering of the identification of two subtypes

To explore the heterogeneity of CHB, we performed a subtyping analysis using the expression profiles of three signature genes, AKT1, POR, and USP21, from the GSE83148 dataset. The optimal number of subtypes, k = 2, was determined by evaluating the shared matrix plot, the CDF plot, and the relative change in the area under the CDF curve. Consequently, we classified CHB samples into two subtypes: Cluster1 and Cluster2 (Fig. 8A–C). We performed principal component analysis (PCA) on these two subtypes and clearly distinguished the Cluster1 and Cluster2 subtypes through the clustering distribution on the scatter plot (Fig. 8E). The findings revealed that the expression levels of AKT1, POR, and USP21 were markedly elevated in the Cluster2 subtype compared to the Cluster1 subtype (Fig. 8F). The high expression of AKT1, POR, and USP21 in the Cluster2 subtype may suggest that patients of this subtype have distinct biological characteristics and potential therapeutic targets.

Fig. 8

Consensus clustering of CHB samples. A Consistent matrix plot, where cleaner blank areas between blue modules indicate more successful analysis. B CDF plots of the consistency clusters showing the relative change in consistency indices with CDF values from k = 2 to k = 6, with the k value of the curve with the most stable change being the optimal fractal number. C Area under the curve of the CDF curve from k = 2 to k = 6. D Trace plot of k = 2 to k = 6. E PCA plot of CHB samples. The scatterplot allows visualization of the feature genes that divide CHB into two subtypes, Cluster1 and Cluster2. F Box plots of the difference in expression of feature genes between Cluster1 and Cluster2 of two subtypes of CHB. *p < 0.05, **p < 0.01, ***p < 0.001

Immune infiltration analysis of two subtypes

To further investigate the role of the immune microenvironment in CHB subtypes, we performed an immune infiltration analysis of different subtypes of CHB samples. The CIBERSORT algorithm was used to determine the proportion of infiltration of 22 immune cell types in Cluster1 and Cluster2, and the results were visualized through stacked histograms (Fig. 9A). Our analysis showed that, compared to Cluster1, Cluster2 had significantly higher infiltration expression of T cells CD8, Plasma cells, and T cells regulatory (Tregs) (Fig. 9B). This finding suggests significant differences in the immune microenvironment between different CHB subtypes.

Fig. 9

Immune infiltration analysis between the two subtypes of CHB. A Stacked histogram of infiltration rates of 22 immune cells in CHB samples. B Box-line plot of differences in the infiltration of 22 immune cells between Cluster1 and Cluster2. *p < 0.05, **p < 0.01, ***p < 0.001, ns means p > 0.05

Next, we investigated the correlation between the signature genes (AKT1, POR, and USP21) and immune cell infiltration. Based on the magnitude and statistical significance, we depicted the correlations of these signature genes with 22 types of immune cells and visualized them through lollipop plots (Fig. 10A–C). The results indicated that both AKT1 and POR were positively correlated with Tregs and Plasma cells while negatively correlated with T cells gamma delta. USP21 was positively correlated with Tregs, NK cells resting, B cells memory, mast cells resting, macrophages M2, and neutrophils, while negatively correlated with T cells follicular helper, macrophages M1, and T cells gamma delta. To further understand the relationship between the feature genes and immune cells, we demonstrated the relationship between the expression level of the feature genes and the degree of immune cell infiltration by scatter plots (Fig. 10D–I). These analyses not only revealed the specific roles of signature genes in the immune microenvironment but also provided new perspectives for the classification of CHB subtypes and the development of immunotherapeutic strategies.

Fig. 10

Feature genes immune cell infiltration analysis. A–C Lollipop plots of the correlation of feature genes with 22 immune cell types. The size of the circle represents the strength of the correlation. D–I The linear relationship between feature gene expression levels and immune cell expression levels. p < 0.05 statistically significant

Differentiated value of USP21 for the natural courses of CHB

In our cohort, a significant increase in the expression levels of USP21 was observed in patients with CHB (Fig. 11A), a finding that is consistent with results from public databases. Through ROC curve analysis, USP21 demonstrated high accuracy in distinguishing between HCs and CHB patients, with an area under the ROC curve (AUC) of 0.834 (Fig. 11C). To further explore the association between USP21 and the natural course of HBV infection, we divided our study cohort into four distinct groups: HBeAg-positive chronic HBV infection, HBeAg-positive CHB, HBeAg-negative chronic HBV infection, and HBeAg-negative CHB. The results showed a statistically significant increase in the expression levels of USP21 during the HBeAg-negative chronic HBV infection stage, which was statistically significant (Fig. 11D). Further, we constructed a random forest (RF) model for demonstrating the ROC curves and AUC values of USP21 in distinguishing these four different natural courses. The results showed that the expression level of USP21 was effective in differentiating the different disease stages of CHB patients (Fig. 11E). To enhance the differentiated efficacy of USP21, we used its expression level in conjunction with two routine liver function indicators: ALT and AST (Fig. 11F). The application of this joint indicator significantly improved the accuracy of differentiating between the four disease stages. These results suggest that the combined use of USP21 expression levels and conventional liver function indices can be an effective strategy to differentiate and monitor the different stages of CHB disease.

Fig. 11

Differentiated value of USP21 for the natural courses of CHB. A Differential gene expression of USP21 in HCs and CHB. B Correlation of USP21 gene expression levels with various clinical indicators. Red indicates a positive correlation, and blue indicates a negative correlation. The numbers in the circles indicate the magnitude of the correlation between gene expression levels and clinical indicators. C ROC curve showing the ability of USP21 to distinguish between HCs and CHB. D Differences in USP21 gene expression levels in different natural courses of CHB. E ROC curves showing the efficacy of USP21 in differentiating patients with CHB of different natural courses. F ROC curve showing the efficacy of USP21 in combination with ALT and AST in differentiating CHB of different natural courses. G1 HBeAg-positive chronic HBV infection, G2 HBeAg-positive CHB, G3 HBeAg-negative chronic HBV infection, G4 HBeAg-negative CHB.*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, ns means p > 0.05

Discussion

Over the past two decades, significant progress has been made in the treatment and clinical management of CHB. Nevertheless, due to the persistence of cccDNA in hepatocytes, there are still no curative drugs that can completely eradicate HBV [20, 21]. The changes in the liver’s immune microenvironment caused by chronic inflammation triggered by HBV infection are key factors in disease progression [22]. Therefore, immunomodulatory therapies and gene-silencing technologies have emerged as the most promising treatment strategies for eliminating the virus [23]. Previous studies have shown that various modes of cell death, including ferroptosis and pyroptosis, are involved in the occurrence and development of HBV infection-related diseases and are closely related to disease prognosis [24]. In this study, we introduced 12 different modes of PCD and comprehensively analyzed multiple RNA sequencing datasets from CHB patients using various bioinformatics tools. We further screened feature genes related to immunity and explored in depth the impact of HBV infection on the liver’s immune microenvironment and its role in disease development. Consequently, we identified core genes that may be involved in the entire natural course of HBV infection and demonstrated their ability to distinguish between different natural stages of CHB by constructing ROC curves.

We used a variety of machine learning methods and went through a scientifically rigorous screening process to finally identify three signature genes that are closely related to CHB. Validation across multiple public datasets revealed that the expression levels of these three signature genes were significantly elevated in patients with CHB. AKT is a serine/threonine kinase that inactivates proapoptotic factors by phosphorylating them, resulting in a pro-cell survival effect [25]. This mechanism facilitates the survival of infected cells and promotes viral replication [26, 27]. AKT1 is the predominant isoform in the cytoplasm, and it has been shown that the interaction of hepatitis B virus x protein (HBx) with AKT1 may affect cell proliferation, apoptosis, and oncogenic processes [28]. POR is a redox protein involved in controlling the detoxification of exogenous substances and maintaining redox homeostasis; it plays an important role in the regulation of iron death [29]. USP21 is a member of the ubiquitin-specific protease family [30]. Numerous studies have shown that USP21 is highly expressed in hepatocellular carcinoma cells and is strongly associated with poor survival in HCC [31, 32]. The single-gene GSEA analysis of the three signature genes showed associations with immune regulation, inflammatory responses, and cellular metabolism, suggesting that these genes may be involved in the occurrence and development of HBV infection through these pathways.

In the comparative analysis of the immune microenvironment of patients with HBV infection, we found significant differences between the immune microenvironment of CHB patients and HCs. Specifically, the infiltration levels of activated CD8 T cells, natural killer cells, macrophages, dendritic cells, and B cells were significantly higher in CHB patients. These immune cells play crucial roles in the immune response to HBV infection. Activated CD8 T cells can directly recognize and kill infected hepatocytes, thereby slowing disease progression [33]. Natural killer cells are key players in antiviral defense. Studies have observed enhanced antigen processing and presentation capabilities of myeloid cells in patients with HBV-related cirrhosis [34]. Macrophages, dendritic cells, and B cells, as professional antigen-presenting cells, were also found to be significantly enriched in CHB patients. Additionally, dendritic cells can activate B cells and promote antibody production in HBV-infected patients, thereby participating in humoral immune responses [35]. Our analysis also revealed that the expression of immune pathways such as antigen processing presentation, chemokines, chemokine receptors, cytokines, cytokine receptors, and natural killer cell cytotoxicity was significantly higher in CHB patients compared to HCs, suggesting that CHB patients have a stronger antibody-dependent cellular cytotoxicity response and a tendency toward chemokine/cytokine signaling. These findings are consistent with the results of previous studies and further confirm the unique immune functions of CHB patients.

To further explore the role of the immune microenvironment in the development of CHB, we employed a signature gene-based subtyping approach to meticulously categorize CHB patients and conducted an immune infiltration analysis of the different subtypes. The results showed that the expression levels of the three characterized genes were higher in the Cluster 2 subtype compared to the Cluster 1 subtype, and the infiltration of T cells CD8, plasma cells, and regulatory T cells (Tregs) was more significant. Specifically, the high infiltration levels of T cells CD8 and Plasma cells may suggest that Cluster 2 has a stronger cytotoxic response and antibody-producing capacity, which is essential for the clearance of HBV-infected hepatocytes [36]. Tregs play a key role in maintaining immune homeostasis and preventing excessive immune responses [37]. The high infiltration of Tregs in the Cluster 2 subtype may indicate enhanced regulation and tolerance of immune responses in CHB patients in vivo. These findings not only deepen our understanding of the immune microenvironment in CHB patients but also provide potential targets for the development of therapeutic strategies against specific subtypes.

Chronic HBV infection is a complex and dynamically changing biological process that reflects the delicate interaction between HBV replication and the host immune system [38]. To gain deeper mechanistic insights into the host immune status during HBV infection, we stratified the natural course of the infection into four stages based on virological and hepatic functional parameters and subsequently collected corresponding cohorts for analysis. Through an analysis of the public database GSE65359, we identified USP21 as a key gene with the potential to differentiate among the four stages of the natural course among three feature genes related to PCD. In previous studies, USP21 is involved in DNA repair and promotion of cell proliferation and invasion, thereby promoting the development of HCC, and is considered an independent prognostic factor [31, 32, 39]. USP21 is crucial for the physiological function of Treg and plays a role in the formation of liver fibrosis and immune tolerance [40]. Our study further revealed the potential of USP21 gene expression levels in differentiating the natural course of CHB patients. More importantly, the discriminatory ability of USP21 was significantly enhanced when it was used in combination with traditional liver function indices, such as ALT and AST. In summary, USP21 is biologically important in the natural course of HBV infection.

This study advanced bioinformatics techniques to identify characterized genes associated with CHB and to explore the crucial role they may play in the disease by regulating the immune microenvironment. However, we must recognize the limitations of this study. Firstly, this study was based on cross-sectional validation of CHB patients at different natural course stages, lacking long-term follow-up data to observe the dynamic development of the disease. Secondly, although we identified USP21 as a candidate gene, its specific molecular role in the pathogenesis of CHB remains unexplored. Future studies need to overcome these limitations to understand the pathologic process of CHB better and develop more effective therapeutic approaches.

Institutional review board statement

This study was conducted by the Declaration of Helsinki and approved by the Ethics Committee of Qilu Hospital of Shandong University (KYLL-202306-021-1). Informed consent was obtained from all subjects involved in this study.