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

We used the “edgeR” R package to identify differentially expressed genes (DEGs) between HBV-HCC samples and normal tissues with a history of HBV in the TCGA database under the criteria |log2FC|>0.585 and Padj<0.05. By intersecting the genes included in GSE45114 with DRGs and DEGs, we define the resulting gene set as BHDRGs refers to the disulfidptosis-related genes in hepatitis B virus-associated hepatocellular carcinoma. We merged the TCGA HBV-HCC transcriptome data with the GSE45114 HBV-HCC transcriptome data using the “sva” and “limma” R packages and averaged gene expression levels for duplicate samples. Subsequently, we reassessed the differential expression of BHDRGs using the Wilcoxon rank-sum test. We evaluated the prognostic significance of BHDRGs using Kaplan-Meier (K-M) survival curves and univariate Cox analysis. Finally, we visualized the results.

To reveal the genomic alterations of BHDRGs, we processed CNV data of HBV-HCC patients and assessed the copy number variation status of BHDRGs in HBV-HCC patients.

Using the expression levels of BHDRGs and survival data, we performed an initial consensus clustering analysis of the samples using the “ConsensusClusterPlus” R package with the PAM algorithm. The optimal number of clusters (K) was determined using clustering heatmaps and the PAC method, and samples were divided into different subtypes. We evaluated the distinction between subtypes using principal component analysis (PCA). Kaplan-Meier survival analysis was performed on different subtypes using the “Survival” and “Survminer” R packages, and heatmaps were drawn combining clinical characteristics. To evaluate the differences in enriched metabolic pathways between subtypes, we performed GSVA (Gene Set Variation Analysis) using the “GSVA” and “GSEABase” R packages with files “c2.cp.kegg.symbols.gmt” and “c5.go.symbols.gmt” from MSigDB. We then assessed immune infiltration levels between subtypes using ssGSEA analysis. To explore potential reasons for the existence of different subtypes in HBV-HCC patients, we identified differential genes (interGenes) between groups using the “limma” R package under the conditions |log2FC| > 1 and P < 0.05.

To analyze the relevant biological functions and structures of interGenes and to identify the corresponding enriched pathways, we utilized the R packages “ClusterProfiler,” “org.Hs.eg.db,” and “enrichplot” to conduct GO and KEGG analyses of the upregulated and downregulated genes of interGenes. The results were subsequently visualized using R.

We conducted a univariate Cox regression analysis on interGenes to identify genes significantly associated with prognosis (uniSigGenes). To investigate more thoroughly the interactions and consistency among BHDRGs and to reveal the finer differences among these subtypes, we performed a more refined clustering analysis of the patients using uniSigGenes.Using the “limma”, “ConsensusClusterPlus”, “survminer”, and “survival” R packages, we performed clustering analysis on the samples based on the expression levels of uniSigGenes using the PAM algorithm, determined the optimal K value and divided the samples into different gene subtypes (genecluster) according to the PAC method and clustering analysis results. Kaplan-Meier survival analysis was performed on gene subtypes combined with survival data.

Next, we randomly divided the samples into validation and test groups in a 1:1 ratio. Using the “glmnet” R package, we performed 10-fold cross-validation and 1000-cycle Least Absolute Shrinkage and selection Operator (LASSO) analysis in the test group, followed by multivariate COX regression analysis to determine the signature genes and establish the prognostic model. We calculated risk scores for samples based on the expression levels and related coefficients of the signature genes. The risk score (RS) formula is as follows: RS = β 1 × expression of  gene 1 + β 2 × expression of  gene 2 + . . . + β n × expression of  gene n .

We performed Kaplan-Meier survival analysis and time-dependent receiver operating characteristic (ROC) curve analysis for high and low-risk groups, and calculated the area under the ROC curve (AUC) to evaluate the model’s performance. The risk score formula was applied to the validation set and the total dataset. Patients were grouped based on risk scores, and heatmaps of signature gene expression, risk score distribution, and survival curves were drawn to validate the prognostic value of the risk scores. We selected the GSE14520 dataset for external validation of our model. The Kaplan-Meier survival curve and the ROC curve were used to evaluate the model’s performance within this external dataset. Additionally, we collected risk models from previous relevant literature to compare with our proposed BHDRGs model. The efficacy of each model was assessed using Kaplan-Meier survival analysis and ROC curve analysis.

Using the “rms” and “regplot” R packages, we constructed an HBV-HCC prognostic nomogram with the help of riskScores, age, gender, stage (American Joint Committee on Cancer-Tumor Node Metastasis staging), AFP (Alpha-Fetoprotein), grade (tumor pathological grading). We plotted calibration curves for 1, 3, and 5 years to verify the accuracy of the nomogram, ROC curves to validate the discrimination of the nomogram, and decision curve analysis (DCA) curves to assess the clinical applicability of the nomogram.

We conducted an analysis of immune cell infiltration based on risk scores, tumor microenvironment (TME) analysis, immunotherapy-related analysis, and tumor mutation burden (TMB) analysis.

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.

We obtained drug half-maximal inhibitory concentration (IC50) data from the Genomics of Drug Sensitivity in Cancer (GDSC) database. Subsequently, we performed drug sensitivity analysis using the “pRRophetic” R package to predict differences in drug therapy based on tumor gene expression levels between the high-risk and low-risk groups. We also conducted a t-test to compare the IC50 values between these two groups. The analysis results were visualized using the “ggplot2” R package.

We collected tumor and adjacent non-tumor tissues from 20 patients with HBV-HCC at the First Affiliated Hospital of Wenzhou Medical University. The study protocol was approved by the Ethics Committee of the First Affiliated Hospital of Wenzhou Medical University (reference number: 2020-074). Written informed consent was obtained from all participants, and the study adhered to the standards outlined in the Declaration of Helsinki.

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.

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.

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

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.

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

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: RS = expDLAT * 0.43798077   5786097 + expSTC 2 * 0.24041938   9109836 + expPOF 1 B * 0.17373046   5002983 + expS 100 A 9 * 0.19007896   8663933 − expCPS 1 * 0.11323419   4208144

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.

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

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

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