Hepatocellular carcinoma (HCC) ranks sixth (4.3%) among newly diagnosed cancer cases globally and is the third (7.8%) leading cause of cancer-related deaths worldwide (Bray et al., 2024). In China, liver cancer is the third most common cancer and the second leading cause of cancer-related mortality (Zheng et al., 2024). In African and Asian populations, the primary etiological factor among hepatocellular carcinoma patients is chronic infection with hepatitis B virus (HBV) (Llovet et al., 2021). Over time, patients often progress from chronic hepatitis B to cirrhosis and eventually HBV-HCC. Due to the marked heterogeneity of hepatocellular carcinoma, the early stages of the disease often present with subtle symptoms, progress rapidly, and make prognosis challenging. Many patients are diagnosed at an advanced stage when symptoms appear, Once the opportunity for surgical intervention is missed, non-surgical treatment options generally offer limited efficacy for most patients (Yang et al., 2019). Given the high incidence, mortality, and significant heterogeneity of HBV-HCC, it is essential to identify more efficient diagnostic biomarkers and new therapeutic targets.

Liu et al. (2023) discovered a novel metabolism-related form of programmed cell death called disulfidptosis. They found that under glucose deprivation, cells with high expression of solute carrier family 7 member 11 (SLC7A11) experience accelerated depletion of cytoplasmic nicotinamide adenine dinucleotide phosphate (NADPH). The lack of NADPH prevents cells from reducing cystine to cysteine, leading to the accumulation of disulfides. This accumulation causes disulfide stress within the cell, resulting in aberrant disulfide bonding between actin cytoskeleton proteins, leading to actin network collapse and cell death. Inducing disulfidptosis by exploiting cancer metabolism weaknesses could provide new therapeutic approaches. Currently, disulfidptosis-related genes have been used to screen for prognostic markers and potential therapeutic targets in various malignancies such as colorectal adenocarcinoma (Li et al., 2024) and lung cancers. HBV-HCC, as one of the etiological types of hepatocellular carcinoma, is also the most prevalent type of this cancer in Asia and Africa (Zheng et al., 2024; Llovet et al., 2021). However, research related to HBV-HCC is still lacking. Furthermore, compared to tumors of other origins, hepatocellular carcinoma exhibits a high degree of heterogeneity (Yang et al., 2019). This indicates that we need to explore in greater depth and with more specificity the therapeutic value that disulfide death-related genes can offer to patients.

This study aims to identify subtypes of HBV-HCC patients with different prognoses, tumor immunity, somatic mutations, and clinical characteristics based on the expression of disulfidptosis-related genes (DRGs). We constructed a prognostic model for HBV-HCC that includes five genes: DLAT, STC2, POF1B, S100A9, and CPS1. We analyzed the tumor microenvironment of different risk groups and assessed the efficacy of immunotherapy and drug treatment based on risk scores. In summary, this study provides a new method for predicting the prognosis of HBV-HCC patients and offers new biomarkers for personalized treatment.

The screening criteria for our study population included: 1. Positive results for serum hepatitis B virus surface antigen (HBsAg), positive HBV-DNA results, or a history of hepatitis B; 2. Pathological results indicating primary hepatocellular carcinoma. Patients who met at least one of the criteria in point 1 and also met the criteria in point 2 were identified as HBV-HCC patients in our study. We downloaded transcriptome data for 424 liver cancer cases from The Cancer Genome Atlas database (TCGA, https://gdc.cancer.gov/). After screening, we obtained 23 samples of adjacent normal tissue and 226 samples of HBV-HCC tumor tissue. We also obtained clinical data for 377 patients and gene mutation data (Simple Nucleotide Variation, SNV) for 368 liver cancer patients from the TCGA database. Additionally, we downloaded the GSE45114 dataset from the Gene Expression Omnibus database (GEO, https://www.ncbi.nlm.nih.gov/geo/), which includes RNA transcriptome data for 19 HBV-HCC samples. And we downloaded the GSE14520 (Roessler et al., 2010) dataset along with the corresponding follow-up data, which were then processed and filtered for analysis. The GSE14520 dataset contains the results of gene expression profiling conducted by Roessler et al. on tumor and paired non-tumor samples, as well as normal liver samples from 64 patients. From the UCSC Xena server (https://xena.ucsc.edu/), we acquired copy number variation (CNV) data. We have collected 118 disulfidptosis-related genes (Supplementary Data 2) from previous studies (Liao et al., 2023; Liu et al., 2023).

First, we used the CIBERSORT algorithm (https://cibersort.stanford.edu/) (Newman et al., 2015) combined with the LM22 signature matrix (Newman et al., 2015) to analyze the infiltration of 22 immune cell types in the TME. In the “Configure Basic CIBERSORT Options” section, since the data had already been normalized, we chose to reject the normalization of the uploaded data, and all other configurations were set to default. Based on the results of the CIBERSORT algorithm, we generated visualizations illustrating the infiltration of immune cells in high and low-risk groups. Furthermore, to assess the levels of stromal and immune cells in the TME in HBV-HCC, we employed the “estimate” R package to evaluate the differences in immune scores, stromal scores, and estimate scores between risk groups.

Subsequently, we performed differential analysis of immune checkpoint gene expression in high and low risk groups. We also assessed the effectiveness of immunotherapy between high and low risk groups using the tumor immune dysfunction and exclusion (TIDE) scores (http://tide.dfci.harvard.edu/). IMvigor 210 is a clinical trial cohort that investigates the efficacy of anti-PD-L1 immunotherapy in patients with metastatic urothelial carcinoma. We obtained transcriptomic data from the IMvigor-210 and selected patient samples with complete immunotherapy data. Using data from the IMvigor-210, we evaluated the ability of our model to predict patient responses to immunotherapy. Additionally, we evaluated the TMB of patients using TCGA somatic mutation data, comparing the TMB status between different risk groups. All results were visualized using the R package.

Unless otherwise specified in the text, we used the Wilcoxon rank-sum test to compare differences in continuous variables between the two groups, and the Kaplan-Meier curve along with the log-rank test to evaluate survival outcomes among different groups. Cox proportional hazards regression analysis was employed to identify independent prognostic factors. The Kruskal–Wallis rank-sum test was utilized to assess differences in gene expression across different patient classifications. Chi-square tests were conducted to compare the associations of categorical variables. Lastly, Pearson correlation analysis was used to evaluate the linear correlation between two continuous variables. R version 4.2.0 was used for statistical analyses. Perl was used to organize the data. To minimize batch effects within the dataset, we applied the ComBat function from the “sva” R package for data adjustment (Leek et al., 2012). For statistical significance in the entire text, numbers, and figure legends, the following terms were used: ***p < 0.001, **p < 0.01, *p < 0.05.

The overall process of this study is illustrated in Figure 1. Comparing the gene expression between HBV-HCC samples from the TCGA and normal tissue samples with a history of HBV infection, we identified 8,110 DEGs, including 2,620 downregulated genes and 5,490 upregulated genes (Figure 2A). We intersected these genes to obtain 17 BHDRGs and plotted a box plot of their differential expression in HBV-HCC tumor and non-tumor tissues (Figure 2B). We also analyzed the copy number variation frequency of the BHDRGs in the samples (Figure 2C), finding that the increase in BOP1 copies was significantly higher than the deletion frequency, while the deletion frequency of WASF2 copies was significantly higher than the increase. Figure 2D illustrates the locations of BHDRGs on the chromosomes along with their copy number variations.

Figure 1. Flow chart of the study. N: cancer-adjacent normal tissue, T: tumor tissue.

Figure 2. Expression and genetic alteration of BHDRGs in HBV-HCC. (A) DEGs between HBV-HCC samples and normal tissues with a history of HBV in the TCGA. (B) The expression of 17 BHDRGs in HBV-HCC and normal tissues. (C) The copy number variation frequency of BHDRGs in HBV-HCC. (D) Variation in the copy number of the BHDRGs at various places on the chromosome. (E) The mutation frequency of BHDRGs in HBV-HCC. *p < 0.05, p < 0.01. ***p < 0.001.

Next, we examined the somatic mutation data in the samples, finding mutations in 29 out of 224 HBV-HCC, with MYH7B having the highest mutation rate (Figure 2E). We then combined the GSE45114 data with the TCGA data, obtaining 240 HBV-HCC samples. The univariate Cox analysis results of BHDRGs in HBV-HCC samples are shown in Table 1. We evaluated the association between significant genes from univariate analysis and patient prognosis using K-M survival analysis (Figures 3A–N). The results indicating that MYH3, GLUD1, and EPAS1 are protective factors, while WASF2, SQSTM1, SLC7A11, PPM1F, PPIH, NDUFS1, MYH7B, ME1, FLNC, FLNA, DBN1, BOP1, ACTN2, and ACTG2 are risk factors for the prognosis of HBV-HCC patients. As shown in Figure 3O, the prognostic network outlines the interactions, interconnections, and prognostic value of BHDRGs in HBV-HCC patients.

Table 1. UniCOX and K-M analysis of BHDRGs.

Figure 3. The K-M analysis revealed that 14 genes within the BHDRGs were significantly associated with HBV-HCC patient prognosis. (A) ACTG2. (B) ACTN2. (C) BOP1. (D) DBN1. (E) FLNA. (F) FLNC. (G) ME1. (H) MYH7B. (I) NDUFS1. (J) PPIH. (K) PPM1F. (L) SLC7A11. (M) SQSTM1. (N) WASF2. (O) Prognostic network of BHDRGs.

To clarify the relationship between BHDRGs and HBV-HCC, we performed consensus clustering analysis on HBV-HCC samples. Using the PAC algorithm in combination with clustering results (Figure 4A), we divided the samples into BHDRG cluster A and BHDRG cluster B. Principal component analysis (PCA) was then used for dimensionality reduction and visualization (Figure 4B). The PCA results showed good discrimination between the two subtypes. Subsequently, we conducted KM survival curve analysis (Figure 4C), showing that the 1 year, 3 year, and 5 year survival rates of group B were significantly better than those of group A. Combined with age, sex and T stage, we plotted a heatmap of BHDRG subtypes in HBV-HCC patients (Supplementary Figure 3A).

Figure 4. (A) The clustering analysis of HBV-HCC samples based on gene expression profiles. (B) The principal component analysis of subtypes. (C) Kaplan-Meier survival analysis shows the two clusters. (D) The differential analyse between immune cells mong the two clusters. *p < 0.05, **p < 0.01, *** k p < 0.001.

GSVA shows the differences between the top20 most significant pathways in cluster A and cluster B (Supplementary Figure 3B). The results indicate that, compared to the BHDRG cluster A, pathways such as linoleic acid metabolism, alpha-linolenic acid metabolism, maturity-onset diabetes of the young, alanine, aspartate, and glutamate metabolism, nitrogen metabolism, and glycine, serine, and threonine metabolism, as well as signaling pathways including the peroxisome proliferator-activated receptor signaling pathway, olfactory transduction, and neuroactive ligand-receptor interaction, are more active in the BHDRG cluster B. In contrast, gene sets related to lysosomes, pathogenic Escherichia coli infection, endocytosis, epithelial cell signaling in Helicobacter pylori infection, neurotrophin signaling pathway, pathways in cancer, small cell lung cancer, pancreatic cancer, and colorectal cancer show higher expression levels in the BHDRG cluster A.

And the boxplot showed the difference of immune cell infiltration in different subtypes (Figure 4F).In the BHDRG cluster A group, the infiltration levels of activated CD4 T cells, activated dendritic cells, CD56dim natural killer cells, immature B cells, immature dendritic cells, myeloid-derived suppressor cells (MDSCs), natural killer T cells, natural killer cells, plasmacytoid dendritic cells, regulatory T cells, and T follicular helper cells are significantly higher than those in the B group. Conversely, neutrophils exhibit higher infiltration levels in the BHDRG cluster B group.

To explore the potential reasons for the prognostic differences between subtypes, we performed differential analysis to obtain interGenes (Supplementary Data 1). And we performed GO and KEGG enrichment analyses to identify the functions and pathways enriched for the upregulated and downregulated genes of interGenes. The GO analysis demonstrated (Figure 5A) that the downregulated genes in interGenes were significantly enriched in collagen-containing extracellular matrix, extracellular matrix organization, extracellular structure organization, and external encapsulating structure organization. And the upregulated genes were primarily enriched in collagen-containing extracellular matrix, blood coagulation, lipid transport, and peptidase inhibitor activity. The KEGG enrichment analysis results indicated (Figure 5B) that the downregulated genes were involved in biological pathways such as Proteoglycans in cancer, Cytoskeleton in muscle cells, ECM-receptor interaction, protein digestion and absorption, and focal adhesion. Meanwhile, the upregulated genes showed significant enrichment in the complement and coagulation cascades, PPAR signaling pathway, and various metabolic pathways.

Figure 5. GO and KEGG analysis (A) Analysis of GO enrichment demonstrated the possible function of the interGene. (B) KEGG pathway enrichment analysis revealed the possible pathways.

For more accurate sample typing, we conducted univariate Cox regression analysis on interGenes to identify uniSigGenes. Based on uniSigGenes, we performed another clustering analysis (Figure 6A). Both the clustering heatmap and PAC algorithm suggested that the optimal clustering model was achieved when K = 2, resulting in geneclusterA and geneclusterB. KM survival curve analysis (Figure 6B) showed that the prognosis of group B was significantly better than group A. We plotted a heatmap of gene clusters, considering patients’ age, sex, and T stage (Figure 6C). We also evaluated the differential expression of BHDRGs between gene clusters (Figure 6D).

Figure 6. (A) Concordance matrix of subtypes. (B) Kaplan-Meier survival curve of the geneClusters. (C) A heatmap illustrated the expression patterns. (D) Expression of BHDRGs between geneCluster A and geneCluster (B, E) LASSO regression analysis. (F) Partial likelihood deviance on the prognostic genes. (G) Kaplan-Meier survival curves for patients in the high- and low-risk groups defined by the prognostic model in the all set, test set, and train set. (H) ROC curve analyses for 1 year, 3 year, and 5 year survival rates in the all set, test set, and train set. (I) The differences in risk score in the BHDRGcluster and the geneCluster. (J) The expression of BHDRGs in the high- and low-risk groups. (K) Sankey diagram shows the construction of the prognostic model.

In the test group, we conducted LASSO-Cox analysis on the expression levels of uniSigGenes (Figures 6E, F), identifying 9 genes (DLAT, SLC2A1, STC2, SLC38A1, POF1B, S100A9, AP1M2,MMP9, CPS1) with optimal λ values. We then performed multivariate regression analysis on these 9 genes to construct a new prognostic risk model, including key genes such as DLAT, STC2, POF1B, S100A9, and CPS1. The risk scoring formula for the model is as follows:

We assigned risk scores to all samples based on the risk formula, dividing patients into high-risk and low-risk groups using the median risk score and conducting Kaplan-Meier survival analysis (Figure 6G). The survival curves of the total group, training group, and test group all showed that the survival rate of the low-risk group was significantly higher than that of the high-risk group. We plotted 1 year, 3 year, and 5 year ROC curves to evaluate the sensitivity and specificity of the gene signature and calculated the area under the ROC curve (AUC), indicating that our model has good predictive power (Figure 6H). And we validated patient groups and their prognostic differences in the GSE14520 dataset (Supplementary Figure 1). Additionally, we compared our model with previous risk models (Supplementary Figure 2). Compared to the Tang signature (Tang et al., 2020), Tian signature (Tian et al., 2024), Li signature (Li et al., 2022), and Weng signatures (Weng et al., 2023),BHDRGs signature demonstrated superior performance in predicting the prognosis of HBV-HCC patients.

We also evaluated the risk score differences between BHDRG subtypes and gene clusters (Figure 6I) and plotted heatmaps related to risk grouping (Supplementary Figure 4). We found that the risk score of genecluster-A was significantly higher than that of genecluster-B, and the risk score of BHDRG cluster A was higher than that of BHDRGcluster-B, consistent with our previous survival analysis results. Finally, we evaluated the expression differences of BHDRGs between high-risk and low-risk groups (Figure 6J). The genes MYH7B, MYH3, and GLUD1 show no differential expression between the high-risk and low-risk groups. In contrast, the expression levels of genes such as EPAS1, WASF2, SQSTM1, SLC7A11, PPM1F, PPIH, NDUFS1, ME1, FLNC, FLNA, DBN1, BOP1, ACTN2, and ACTG2 are significantly higher in the high-risk group compared to the low-risk group.

Next, we used Sankey diagrams (Figure 6K) to show the distribution of HBV-HCC patients in the clustering process. Most genecluster-A samples came from the BHDRGcluster-A group, which had worse survival data. In the risk scoring, they were mainly distributed in the high-risk score group, with deceased patients primarily concentrated in the high-risk group. The Sankey results were consistent with previous analysis results.

We conducted univariate and multivariate Cox regression analyses on the clinical features of all cases combined with riskScores (Figures 7A, B). Univariate Cox regression analysis showed that riskScore, age and stage were significantly associated with patient survival probability. Multivariate Cox regression analysis also found that riskScores, age and stage were independent predictors of prognosis. We then constructed a nomogram comprising riskScores, age, and stage (Figure 7C). We verified the accuracy of the nomogram by plotting calibration curves (Figure 7D). The ROC curve of the nomogram (Figure 7E) suggested that it had good discriminative power compared to other influencing factors. The DCA curve (Figure 7F) also showed that the nomogram had good applicability in clinical practice. Next, we compared our nomogram with those of Tian (Tian et al., 2024), Tang (Tang et al., 2020), and Li (Li et al., 2022). Our findings indicated that (Supplementary Figure 5), in the ROC analysis, calibration curves, and DCA curve, our nomogram consistently demonstrated superior performance. Parameters such as specificity, sensitivity, precision, positive and negative likelihood ratios, and Youden index can be found in Supplementary Data 3.

Figure 7. (A, B) The forest plots for univariate and multivariate. (C) The nomogram of the risk and clinical parameters (age, stage) of all sets. (D) The calibration curves displayed the accuracy of the nomogram in the l year, 3 year, and 5 year. (E) The ROC curve of the nomogram. (F) The DCA curve of the nomogram.

To assess the therapeutic value of the risk model for HBV-HCC patients, we explored immune cell infiltration between high-risk and low-risk groups. Immune infiltration analysis (Figure 8A) indicated that Macrophages M0 and M2 were positively correlated with risk scores, while T cells gamma delta and T cells CD8 were negatively correlated. In the tumor microenvironment scores (Figure 8B), StromalScore, ImmuneScore, and ESTIMATEScore were all significantly higher in the high-risk group than in the low-risk group. To assess the value of risk scores in immunotherapy, we examined the differential expression of immune checkpoint genes between high-risk and low-risk groups (Figure 8C). Results showed differential expression of 40 immune checkpoint genes between the two groups, with higher expression in the high-risk group, suggesting that high-risk patients might benefit more from immunotherapy. We used TIDE scores to evaluate the potential clinical efficacy of immunotherapy in different risk groups. Typically, higher TIDE scores indicate a higher likelihood of tumor immune escape and poorer efficacy of immune checkpoint inhibitors. The results suggested (Figure 8D) that high-risk patients had lower TIDE scores than low-risk patients, implying that high-risk patients might achieve better efficacy from immunotherapy. And the IMvigor 210 cohort revealed that (Figure 8E), patients who responded to anti-PD-L1 immunotherapy had higher risk scores in our risk model compared to those who did not benefit from the treatment. Additionally, we evaluated the mutation frequency of HBV-HCC patients between different risk groups (Figures 8F, G), showing that TP53 had the highest mutation frequency, with a significant difference in mutation frequency between high-risk and low-risk groups. However, we found that there was no significant difference in TMB between high-risk group and low-risk group (Figure 8H).

Figure 8. (A) correlation between the signature and immunecells. (B) comparison of ESTIMATE scores, Stromal scores, and Immune scores between the high-risk and low-risk groups. (C) Differential expression of immune checkpoints between the high-risk and low-risk groups. (D) Differences in TIDE score between high-risk and low-risk groups. (E) The ability of model to predict patient responses to immunotherapy in IMvigor 210. (F) The frequency of mutations in the high-risk and (G) low-risk groups. (H) Differences in TMB between high-risk and low-risk groups.

We assessed the predictive capability of our gene signature for drug treatment in HBV-HCC patients. High-risk groups showed greater sensitivity to drugs such as Bortezomib, Camptothecin, Cisplatin, Cytarabine, Epothilone B, Etoposide, Gemcitabine, Paclitaxel, Sorafenib, and Tipifarnib. Low-risk groups were more sensitive to Axitinib, Erlotinib, Metformin, PD 0332991 (Palbociclib), Roscovitine, and Temsirolimus (Figure 9).

Figure 9. Drug sensitivity analysis between the high-risk and low-risk groups. (A) Bortezomib, (B) Camptothecin, (C) Cisplatin, (D) Cytarabine, (E) Epothilone B, (F) Etoposide, (G) Gemcitabine, (H) Paclitaxel, (I) Sorafenib, (J) Tipifarnib, (K) Axitinib, (L) Erlotinib, (M) Metformin, (N) PD 0332991 (Palbociclib), (O) Roscovitine, (P) Temsirolimus.

We conducted RT-qPCR analysis on 10 pairs of clinical samples from HBV-HCC patients. The results showed (Figure 10A) that the mRNA expression of DLAT in tumor samples was significantly higher than that in normal samples. Additionally, we conducted immunohistochemical analysis on another 10 pairs of clinical samples from HBV-HCC. Our immunohistochemical analysis revealed that DLAT protein expression in HBV-HCC tumor tissues was markedly elevated compared to adjacent normal liver tissues (Figure 10B).

Figure 10. (A) Verify the mRNA expression levels of DLAT. (B) IHC staining was performed on tumor tissues and adjacent normal tissues to verify the expression of DLAT.

In recent years, the incidence of liver cancer has been on the rise (Siegel et al., 2022), with hepatocellular carcinoma (HCC) constituting the majority of liver cancer diagnoses and deaths (Llovet et al., 2021). A history of HBV infection has become one of the most common risk factors for HCC. Recently, Xiaoguang Liu and colleagues proposed a new regulatory form of cell death, disulfidptosis, which could be utilized in cancer treatment (Liu et al., 2023). Research has indicated the significant role of disulfidptosis in various tumors, including lung cancer (Wang et al., 2024), renal cell carcinoma (Xu et al., 2023). However, in the study of hepatocellular carcinoma (HCC), research related to disulfidptosis remains limited. Zhao et al. developed (Zhao et al., 2023) a tumor prognosis signature based on disulfidptosis to assist in clinical decision-making. However, they simultaneously noted that the impact of viral infection was not considered in their study. Given that the etiology of potential liver disease may influence prognosis, they suggested that the prognostic assessment could be biased. Therefore, the prognostic role of disulfidptosis and its value in guiding treatment strategies for HBV-HCC patients require further investigation. We conducted additional research on HBV-HCC based on the differential expression of BHDRGs and identified two subgroups within the HBV-HCC patient population that exhibited different clinical characteristics and survival outcomes. Additionally, we constructed a prognostic model for HBV-HCC patients to guide clinical treatment and facilitate risk stratification.

This study performed two clustering analyses based on BHDRGs, leading to the identification of two subtypes with distinct clinical characteristics and significantly different prognoses. The first clustering analysis provided a general classification of this population, while the second clustering analysis allowed for a more detailed subdivision of these subtypes, revealing finer differences among them. Notably, in the context of high-heterogeneity tumor research, the second clustering analysis aids in uncovering the pronounced heterogeneity of tumors and offers deeper insights into the roles and mechanisms of BHDRGs in the prognostic differences among the various subtypes. By employing the second clustering analysis, we can more accurately identify the different subtypes within tumors and design more targeted treatment strategies.

In the first clustering analysis, patients were categorized into two distinct subtypes, each exhibiting different clinical characteristics and prognostic outcomes. K-M survival analysis revealed a poorer prognosis for group A compared to group B in the BHDRGs subtypes. Immune infiltration analysis suggested higher immune cell abundance in group A samples, indicating a more active tumor immune microenvironment. To further explore the possible reasons for the survival performance differences between subtypes, GSVA pathway analysis was conducted. This analysis revealed the enrichment of pathways related to small cell carcinoma, prostate cancer, colorectal cancer, and bacterial infection in group A samples. In contrast, group B was mainly enriched in various metabolic pathways, such as those involving amino acids, arachidonic acid, and linoleic acid. Additionally, multiple ligand-receptor signaling pathways were enriched in group B. Mossmann et al. has indicated (Mossmann et al., 2023) that while the mechanisms and targets of metabolic reprogramming remain unclear, it is certain that metabolic reprogramming plays a significant role in the occurrence and progression of tumors. Furthermore, several studies (Kudo et al., 2020; Lin et al., 2024; Liu et al., 2018) suggest that metabolic reprogramming could serve as a new direction for cancer treatment.

To further analyze the connection between BHDRGs and HBV-HCC, GO functional enrichment analysis and KEGG pathway enrichment analysis were performed on interGene. The research findings indicate that interGenes are significantly enriched in various functions and structures, including extracellular structures, the extracellular matrix, and the cytoskeleton, and are involved in multiple metabolic and signaling pathways. Notably, Wu et al. indicated (Wu et al., 2023) the association of the extracellular matrix with the progression of various tumors, including hepatocellular carcinoma, pancreatic ductal adenocarcinoma, and breast cancer. Extracellular matrix stiffening promotes cell proliferation, epithelial-mesenchymal transition, cancer cell metastasis, and drug resistance (Pickup et al., 2014; Walker et al., 2018). The functional and pathway differences between these two subtypes are closely related to the mechanism of disulfidptosis, which may be one of the contributing factors to the differences in patient prognosis.

The five key genes in our risk model are DLAT, STC2, POF1B, S100A9, and CPS1. Previous studies have found that all of these genes are involved, to varying degrees, in the development and progression of various tumors. STC2 is a glycoprotein that is expressed in multiple tumor tissues, and several studies have indicated (Cheng et al., 2018; Qie and Sang, 2022) that overexpression of STC2 promotes cell proliferation and migration, and is associated with tumor growth, invasion, metastasis, and poor prognosis. Inhibiting STC2 overexpression may be a promising candidate for targeted liver cancer therapy. Lacombe et al. indicated (Lacombe et al., 2006) that POF1B is located in a critical region for normal ovarian function, encoding a protein that binds to non-muscle actin filaments and plays an important role in the occurrence and development of premature ovarian failure. Crespi et al. suggest that POF1B may function in regulating cell adhesiveness (Crespi et al., 2015). S100A9 is an important immune-related protein associated with inflammation or cancer, and is believed to promote the occurrence, development, and metastasis of tumors. Zhan et al. reported (Zhan et al., 2023) that HBV-induced activation of S100A9 triggers the RAGE/TLR4-ROS signaling pathway, leading to the formation of numerous neutrophil extracellular traps, thereby promoting the growth and metastasis of HBV-HCC cells. Inhibiting S100A9 significantly inhibits the growth and metastatic ability of HCC. CPS1 encodes the rate-limiting enzyme of the urea cycle, which can more effectively clear ammonia from the body. Aberrant overexpression of the CPS1 is associated with the rapid proliferation of tumor cells, increasing the activity of pyrimidine biosynthesis through metabolic reprogramming (Zhang et al., 2023). The overexpression of CPS1 promotes the proliferation of tumor cells by increasing the de novo biosynthesis of pyrimidine nucleotides, making CPS1 a novel target for anti-tumor drugs (Owusu-Ansah et al., 2023). Interestingly, Zhang et al. (2023) found that the urea cycle key enzyme CPS1 is often low in HCC, and is positively correlated with patient prognosis. The contradictory findings suggest that the impact of CPS1 on the prognosis of HBV-HCC patients may involve more complex mechanisms, which require further investigation. As one of the core components of PDC, dihydrolipoyl transacetylase (DLAT) is essential for glucose metabolism and the tricarboxylic acid (TCA) cycle (Park et al., 2018). Previous studies have indicated that DLAT plays a positive role in the occurrence, development, and metastasis of various tumors (Chen et al., 2022; Li et al., 2023). Increasing evidence suggests that abnormal PDC activity is correlated with the malignant progression of cancer and poor clinical prognosis (Anwar et al., 2021; Jiang et al., 2024). Furthermore, studies have indicated (Tsvetkov et al., 2022) a strong correlation between DLAT and cuproptosis. Recently, Li et al. provided the first evidence that (Li et al., 2023) high expression of DLAT accelerates the progression of hepatocellular carcinoma, offering a promising therapeutic target for HCC treatment. Additionally, we noted that DLAT exhibited the highest coefficient weight in our risk model. Encouraged by these results, we decided to conduct further research on DLAT. Analysis of clinical HBV-HCC samples indicates that DLAT expression levels are significantly higher in cancerous tissues compared to adjacent normal tissues.

We believe that the differential expression of these five signature genes may play a crucial role in the prognostic differences between genetic subtypes. Therefore, we scored HBV-HCC samples according to the model, obtaining corresponding risk scores. Patients were divided into high-risk and low-risk groups based on the median risk score. Kaplan-Meier survival analysis revealed a significant difference in prognosis between the high-risk and low-risk groups. Additionally, we found significant differences in risk scores between both BHDRGs and genetic subtypes, with higher risk scores observed in subtypes with poorer prognosis. We validated the performance of the signature using both internal and external independent datasets. ROC analysis also demonstrated its strong predictive ability. Compared with previous studies, our BHDRGs signature showed relatively good performance and holds potential as a prognostic tool for predicting outcomes in HBV-HCC patients.

To illustrate the cli