Abstract Background Disulfidptosis is a newly discovered type of cell death. We aim to identify hub genes associated with both disulfidptosis and immune infiltration in hepatocellular carcinoma (HCC) patients, and to develop an individualized risk prediction model. Methods The TCGA-LIHC cohort was utilized as the training set to identify molecular subtypes associated with disulfidptosis and to perform immune infiltration analysis. WGCNA, univariate Cox, and LASSO algorithm were employed to select hub genes for constructing the prognostic model. ICGC-LIRI cohort was utilized as an independent testing set. Validation of the expression of hub genes was performed in vitro using qRT-PCR and Western blot. Results Cluster 1 was identified as the disulfidptosis associated molecular subtype, characterized by higher expression of disulfidptosis related genes (DRGs) and immune infiltration levels. ANXA2, MSC, and ST6GALNAC4 were identified as hub genes for calculating the risk score. The high-risk group were more likely to benefit from immunotherapy, targeted therapy and chemotherapy. A prognostic model was developed combining clinicopathological factors with satisfactory predictive accuracy. The hub genes were found upregulated in HCC cell line. Conclusions Our findings provide valuable theoretical support for prognostic prediction and the evaluation of therapeutic outcomes in relation to disulfidptosis and immune infiltration in HCC, highlighting the importance of conducting in-depth research on disulfidptosis-related mechanisms. Keywords: Disulfidptosis, Immune infiltration, Immunotherapy, Hepatocellular carcinoma, Prognosis prediction model Graphical abstract [27]Image 1 [28]Open in a new tab 1. Background Cell death is a continuous event that takes place throughout the lifespan of an individual. Over the years, several modes of cell death have been identified, including apoptosis, necroptosis, autophagy, ferroptosis, and cuproptosis. These modes of cell death play significant roles in various biological processes that contribute to the occurrence and progression of HCC. In HCC, down-regulating the expression of anti-apoptotic proteins cFLIPL and cIAP-1 can be achieved by knocking out RIPK1/TRAF2, ultimately promoting apoptosis [[29]1]. Lan et al. found that inhibiting NUPR1 can promote apoptosis and necroptosis of HCC [[30]2]. Lin et al. discovered that under hypoxic conditions, the downregulation of METTL3 can activate the autophagy-related pathway, promote HCC resistance to sorafenib, and increase the expression of angiogenesis-related genes [[31]3]. The LIFR/NF-κB/LCN2 axis has been found to induce HCC and affect HCC sensitivity to ferroptosis inducers [[32]4]. The aforementioned cell death pathways have significant implications in HCC, with various genes and signaling pathways being identified as potential therapeutic targets. Despite these advancements, the prognosis for HCC remains unsatisfactory. Disulfidptosis is a recently identified cell death pathway that exhibits close association with the actin cytoskeleton. In cells characterized by high levels of SLC7A11, an abnormal disulfide bonding of actin cytoskeleton proteins occurs under conditions of glucose starvation. Consequently, there is a collapse of f-actin mediated by SLC7A11, ultimately resulting in cell death. Given its close connection to the actin cytoskeleton, this newly discovered death pathway holds promise as a potential target for tumor therapy [[33]5]. However, the role of disulfidptosis in the development and progression of HCC is still unclear, and further research is needed to improve the poor prognosis of HCC. The infiltration of immune cells within the tumor microenvironment (TME) of HCC plays a vital role in assessing the prognosis of HCC patients and predicting the effectiveness of immunotherapy [[34]6]. In recent years, the emergence of immunotherapy, including immune checkpoint inhibitors, has proved to be a powerful tool in improving the prognosis of HCC patients. However, the variable response rates among patients pose a challenge. Therefore, the identification of HCC patients who are more likely to benefit from immunotherapy can aid in the development of personalized treatment strategies [[35]7]. In this study, we utilized publicly available databases such as The Cancer Genome Atlas (TCGA) and International Cancer Genome Consortium (ICGC) to conduct an expression analysis of disulfidptosis-related genes (DRGs) in both HCC tissue and normal tissue. Subsequently, we employed bioinformatics techniques to identify hub genes that exhibit significant associations with the prognosis of HCC. Additionally, we calculated risk scores to predict immune infiltration and drug sensitivity. By integrating clinicopathological factors, we constructed a prognostic model that demonstrated favorable predictive accuracy. Overall, our study provides a theoretical foundation and highlights a new avenue of research regarding the role of disulfidptosis in HCC. 2. Results 2.1. Identification of disulfidptosis-related molecular subtypes In the TCGA-LIHC cohort, we conducted differential analysis of gene expression be-tween HCC and normal tissues, and found that most DRGs were upregulated in HCC ([36]Fig. 1A and B). Subsequently, we performed consensus clustering analysis based on DEGs and identified two clusters ([37]Fig. 1C–E). Cluster 1 had higher expression of DRGs and could be considered a molecular subtype related to disulfidptosis ([38]Fig. 1F). In addition, we conducted survival analysis and observed that patients in Cluster 1 had significantly worse overall survival rates than the patients in the other subtype group ([39]Fig. 1G). By using ESTIMATE algorithm, we conducted immune infiltration analysis and found a significant difference in immune score between clusters, with Cluster 1 having significantly higher levels of immune score than Cluster 2 ([40]Fig. 1H). Fig. 1. [41]Fig. 1 [42]Open in a new tab Identification of disulfidptosis-related subtype. (A–B) Heatmap of 24 DRGs and volcano plot of DEGs. (C–E) Unsupervised clustering based on DEGs. (F) Expression of 24 DRGs between clusters. (G) Survival analysis. (H) Comparison of immune score using ESTIMATE. ∗p < 0.05, ∗∗p < 0.01, and ∗∗∗p < 0.001. ns indicates p > 0.05. 2.2. Identification of hub genes related to disulfidptosis and immunity Considering the notable variations in DRGs and immune infiltration between the two clusters, we conducted WGCNA to identify potential hub genes among these DEGs that are significantly associated with disulfidptosis and immune infiltration ([43]Fig. 2A and B). The results showed a significant correlation between the black and salmon modules and both clustering and immune infiltration ([44]Fig. 2C–E). Consequently, we obtained a total of 240 genes (151 from black module and 89 from salmon module) for further analysis. Fig. 2. [45]Fig. 2 [46]Open in a new tab WGCNA was performed to identify hub modules. (A) Optimal power selected using scale-free topology fit index and network connectivity. (B) Dendrogram of DEGs in HCC. (C) Correlation between modules and both clustering and ESTIMATE scores. (D) Correlation between GS for Cluster 1 and MM for salmon module. (E) Correlation between GS for Cluster 1 and MM for black module. 2.3. Calculation of risk score based on DIRGs We performed univariate Cox regression analysis on the genes within these two modules and identified 26 genes that exhibited significant associations with prognosis ([47]Fig. 3A). Notably, all of these genes were found to be prognostic risk factors. Subsequently, we conducted LASSO-Cox regression analysis using these 26 prognosis-related genes and selected three DIRGs to calculate a risk score ([48]Fig. 3B and C). By using the median risk score as a threshold, patients were categorized into high-risk and low-risk groups ([49]Fig. 3D and E). We observed notable disparities in the expression of 24 DRGs between these two groups, with the majority of DRGs demonstrating significantly heightened expression in the high-risk group ([50]Fig. 3F–H). Consistent outcomes were observed in the testing set. Considering previous research indicating that elevated expression of SLC7A11 is involved in disulfidptosis, we performed correlation analysis. As anticipated, all three DIRGs displayed substantial positive correlations with SLC7A11 ([51]Fig. 3I–N). Fig. 3. [52]Fig. 3 [53]Open in a new tab Identification of DIRGs. (A) Univariate Cox regression analysis of hub genes in black and salmon modules. (B–C) LASSO analysis to identify DIRGs. (D–E) Survival information, risk score and DIRGs expression between high- and low-risk groups. (H) Correlation analysis performed on DRGs, DIRGs, and risk scores. (I–K) Correlation analysis performed on SLC7A11 and DIRGs in TCGA cohort. (L–N) Correlation analysis performed on SLC7A11 and DIRGs in ICGC cohort. ∗p < 0.05, ∗∗p < 0.01, and ∗∗∗p < 0.001. ns indicates p > 0.05. 2.4. Clinicopathological features and biological functional analysis between groups We performed survival analysis on the high- and low-risk groups, which uncovered a significantly worse prognosis in the high-risk group ([54]Fig. 4A and B). Additionally, Time-ROC analysis demonstrated satisfactory predictive accuracy at various time points for HCC ([55]Fig. 4C and D). Furthermore, we examined the disparities in clinicopathological factors between the two groups and identified a correlation between the risk score and T stage, suggesting that patients with higher risk scores tend to have more advanced T stages ([56]Fig. 4E and F). Fig. 4. [57]Fig. 4 [58]Open in a new tab Clinicopathological analysis and biological functional analysis. (A–B) Survival analysis between two groups. (C–D) Time-ROC analysis of TCGA and ICGC cohorts. (E–F) Comparison of risk score among different T stages. (G–H) GO and KEGG analysis of risk score in TCGA cohort. (I–K) GO, and KEGG analysis of DEGs in high-risk group in TCGA and ICGC cohorts. (K) GSEA analysis of DEGs in high- and low-risk groups. ∗p < 0.05, ∗∗p < 0.01, and ∗∗∗p < 0.001. - indicates p > 0.05. By performing GO analysis, we found that risk scores are associated with multiple biological processes, including immunity, cytoskeleton, cell adhesion, and others ([59]Fig. 4G, [60]Fig. S2). KEGG analysis revealed enrichment in signaling pathways related to immunity and tumors ([61]Fig. 4H, [62]Fig. S2). Next, we conducted differential analysis between the test and training sets. GO analysis of the DEGs indicated their close association with processes such as immune cell migration, activation, and differentiation, as well as additional processes ([63]Fig. 4I). KEGG analysis showed significant enrichment in signaling pathways related to immunity and tumors ([64]Fig. 4J). GSEA demonstrated that the upregulated DEGs in the high-risk group exhibited significant enrichment in adhesion- and immune-related signaling pathways, such as positive regulation of leukocyte cell-cell adhesion and regulation of T cell receptor signaling pathway. Conversely, the downregulated DEGs in the high-risk group presented notable enrichment in metabolism-related signaling pathways such as fatty acid catabolic process ([65]Fig. 4K). 2.5. Tumor mutation burden analysis We examined the nucleotide mutation landscapes and copy number variations using the R package maftool and the online tool Gene Set Cancer Analysis (GSCA). Different degrees of copy number variations have been observed in 24 DRGs and 3 DIRGs, with copy number heterozygous amplification/deletion being more frequently observed ([66]Fig. 5A). 19 out of the 24 DRGs exhibited varying degrees of single nucleotide variations. Missense mutation was the most prevalent type of deleterious mutation, primarily characterized by the C > T base substitution ([67]Fig. 5B). Top 5 DRGs with the highest frequency of deleterious mutations were FLNB, TLN1, FLNA, LRPPRC, and IQGAP1. Nevertheless, there was no notable disparity in TMB between the high- and low-risk groups ([68]Fig. 5C, D, [69]Fig. S3). Subsequently, we conducted a survival analysis based on TMB by selecting the optimum cutoff value that produced the most significant survival difference. Patients in high TMB group demonstrated considerably worse prognosis ([70]Fig. 5E). Moreover, survival analysis incorporating both TMB and risk scores unveiled patients exhibiting high levels of both factors experienced the most unfavorable prognosis, whereas those with low scores and TMB tended to have the most favorable prognosis ([71]Fig. 5F and G). Fig. 5. [72]Fig. 5 [73]Open in a new tab TMB analysis of DRGs and LIHC patients. (A) CNV analysis of DRGs and three DIRGs. (B) SNV analysis of DRGs and three DIRGs. (C–D) TMB analysis of high- and low-risk groups. (E–G) Survival analysis based on TMB and risk score. ∗p < 0.05, ∗∗p < 0.01, and ∗∗∗p < 0.001. ns indicates p > 0.05. 2.6. Immune infiltration analysis In order to elucidate the distinctions in immune microenvironments between the groups, we performed immune infiltration analysis. The results revealed that the high-risk group exhibited elevated stromal score, immune score, and ESTIMATE score, along with reduced tumor purity ([74]Fig. 6A and B). The immune cell landscapes in training sets and testing sets were estimated by CIBERSORT ([75]Fig. S4). The ssGSEA analysis demonstrated substantial upregulation of 24 immune-infiltrating cell types in the high-risk group, such as CD4, CD8, dendritic cells, macrophages, and monocytes ([76]Fig. 6C and D). Based on the correlation analysis, we discovered significant associations between DIRGs, the risk score, and the majority of immune cells, with ANXA2 displaying the most notable correlation ([77]Fig. 6E). CNP0000650 scRNA seq revealed considerable expression of ANXA2 in diverse constituents of TME, encompassing immune cells, epithelial cells, and mast cells ([78]Fig. 7A). Analysis of [79]GSE125449 using online tool TISCH displayed remarkably expression of ANXA2 in HCC cells and endothelial cells ([80]Fig. 7B). The expression of MSC and ST6GALNAC4 can be detected in a subset of cells of TME. The online tool scAtlasLC further confirmed these findings ([81]Fig. 7C). It is worth noting that, analysis using TISCH revealed remarkably elevated ANXA2 expression patterns across various components of TME, including monocytes/macrophages, tumor cells, endothelial cells, and fibroblasts, among patients receiving PD1/CTLA4 immunotherapy ([82]Fig. 7D). There is no significant difference in the expression of MSC and ST6GALNAC4 in most cells before and after immunotherapy. Fig. 6. [83]Fig. 6 [84]Open in a new tab Immune infiltration analysis. (A–B) ESTIMATE scores between groups. (C–D) ssGSEA analysis between groups. (E) Correlation analysis of DIRGs, risk score and immune infiltrating cells. ∗p < 0.05, ∗∗p < 0.01, and ∗∗∗p < 0.001. ns indicates p > 0.05. Fig. 7. [85]Fig. 7 [86]Open in a new tab ScRNA-seq analysis. (A) ScRNA-seq of DIRGs in CNP0000650. (B) ScRNA-seq of DIRGs using TISCH2. (C) ScRNA-seq of DIRGs using scAtlasLC. (D) Expression of DIRGs between PD-L1/CTLA-4 group and control group. 2.7. Therapeutic efficiency analysis We conducted the correlation analysis of the risk score and immunomodulators, including immunostimulators, immunoinhibitors, and major histocompatibility complex (MHC) genes. A significant positive correlation was observed between the risk score and most of the immunomodulators. Top 10 immunomodulators from each of the three categories were selected for visualization ([87]Fig. 8A–C). Considering that the primary approach in immunotherapy is immune checkpoint inhibitors, we also examined the expression differences of immunoinhibitors between the two groups. The majority of immunoinhibitors showed significantly higher expression in the high-risk group ([88]Fig. 8E and F). Using the Genomics of Drug Sensitivity in Cancer database (GDSC) database and the pRRophetic package, we compared the drug sensitivities of the high- and low-risk groups by estimated IC50. Out of the 138 drugs, 89 had lower IC50 values in the high-risk group, while 18 drugs had lower IC50 values in the low-risk group ([89]Fig. 8D). Many commonly used drugs in systemic therapy for HCC, such as cisplatin, fluorouracil, gemcitabine, and doxorubicin, exhibited lower IC50 values in the high-risk group ([90]Fig. 8G). Additionally, through the publicly available dataset [91]GSE109211, we found that patients with high-risk score exhibit a higher response rate to sorafenib. ([92]Fig. 8H). The Sankey diagram visualized the survival information, clustering results, risk scores, and TMB status of TCGA-LIHC patients ([93]Fig. 8I). Fig. 8. [94]Fig. 8 [95]Open in a new tab Therapeutic effect analysis. (A) Top 10 immunostimulators exhibiting the highest correlation with risk score. (B) Top 10 immunoinhibitors exhibiting the highest correlation with risk score. (C) Top 10 MHC genes exhibiting the highest correlation with risk score. (D) Estimated IC50 of 138 drugs. (E–F) Expression of immunoinhibitors between high- and low-risk groups. (G) Estimated IC50 of commonly used drugs in HCC. (H) Comparison of risk score between responder and non-responder group in [96]GSE109211. (I) Sankey diagram of LIHC patients. ∗p < 0.05, ∗∗p < 0.01, and ∗∗∗p < 0.001. ns indicates p > 0.05. 2.8. Establishment of individualized risk model We performed a multivariate analysis of risk score and several clinicopathological factors, such as T stage, gender, and age. The findings demonstrated that T stage and risk score were significant independent prognostic factors for HCC prognosis ([97]Fig. 9A). By utilizing the rms package, we constructed a personalized risk prognosis model for HCC patients by integrating age, T stage, and risk score ([98]Fig. 9B). The calibration curves for both the training and testing datasets indicated a high level of predictive accuracy ([99]Fig. 9C and D). Fig. 9. [100]Fig. 9 [101]Open in a new tab Establishment of risk model. (A) Multivariate analysis of risk score and clinicopathological factors. (B) Nomogram built for predicting overall survival of HCC patients. (C–D) Calibration curve plots comparing actual and predicted OS. (E) qRT-PCR results of three DIRGs. (F) WB results of three DIRGs. ∗p < 0.05, ∗∗p < 0.01, and ∗∗∗p < 0.001. ns indicates p > 0.05. We proceeded to validate the expression of three DIRGs in vitro. As anticipated, ANXA2, MSC, and ST6GALNAC4 all displayed significantly elevated expression levels at both the mRNA and protein levels in HCC cell line SNU-387 compared to normal liver cell line HL-7702 ([102]Fig. 9E and F). 3. Discussion Abnormalities in the actin cytoskeleton of cells contribute to processes, such as invasion and metastasis, in HCC [[103]8,[104]9]. Liu et al.'s study unveiled a novel finding that the actin cytoskeleton's vulnerability to disulfide bond stress can induce disulfidptosis [[105]5]. The treatment based on disulfidptosis is still in the research stage. According to literature reports, glucose depletion inducers can trigger disulfidptosis, while disulfide reducing agents can inhibit disulfidptosis [[106]5,[107]10]. The development of therapies capable of inducing this recently identified type of cell death could potentially enhance the unfavorable prognosis of HCC.In recent times, immunotherapy, exemplified by immune checkpoint inhibitors, has emerged as one of the foremost systemic treatment approaches for HCC. Nevertheless, the response rate remains approximately 20 % [[108]11]. This study was conducted to address the need for a risk scoring system that predicts overall survival and response rate to systemic therapy. The implementation of such a system will facilitate the creation of individualized treatment plans in clinical settings, ultimately leading to improved patient outcomes and treatment efficacy. Bioinformatic analysis was conducted using two independent databases, TCGA and ICGC, to improve the study's reliability. A remarkably upregulation of DRGs was observed in HCC tissues, suggesting a potential involvement of disulfidptosis in the pathogenesis and progression of HCC. Subsequently, through unsupervised clustering analysis of DEGs in HCC, Cluster 1 characterized by higher expression of DRGs, including SLC7A11, was identified as the molecular subtype associated with disulfidptosis. Furthermore, Cluster 1 exhibited worse overall survival and higher immune infiltration. These findings support our hypothesis that DRGs and relevant pathways may be implicated in immune infiltration and overall survival in HCC patients. By conducting WGCNA, we identified two hub modules significantly associated with disulfidptosis and immune infiltration. Through univariate analysis and LASSO-Cox regression analysis, we ultimately obtained 3 DIRGs closely related to HCC prognosis, which were used to construct a risk scoring model. The majority of DRGs were highly expressed in the high-risk group, which had worse overall survival and higher immune infiltration. Cluster 1 had a higher risk score, validating it as a subtype related to disulfidptosis ([109]Fig. S5). Three DIRGs were identified as hub genes that determine molecular subtypes related to disulfidptosis. They exhibit a significant positive correlation with SLC7A11, suggesting their potential crucial role in the induction of disulfidptosis. By detecting these three genes, we can accurately predict the status of disulfidptosis and determine the level of immune infiltration in individual patients with HCC. The risk score was associated with the T staging of tumors, with patients at higher T stages tending to have higher risk scores, indicating a potential correlation between the risk score and malignancy of HCC cells. Further functional analysis revealed significant correlations between our risk score and biological processes such as immunity, cellular cytoskeleton, and cell adhesion. Additionally, enrichment was observed in cancer related signaling pathways. Therefore, we speculate that molecular pathological changes happened in biological processes such as cellular cytoskeleton and cell adhesion in HCC patients with high-risk scores, leading to the upregulation of disulfidptosis in key regulatory cells in the TME, thereby affecting the invasion, metastasis, drug resistance and other processes of HCC. TMB provides a quantitative estimate of the total number of mutations in the coding regions of the tumor genome. These mutations are processed into neoantigens and presented to T cells, triggering an anti-tumor immune response. Therefore, tumor cells with higher TMB levels may be more easily recognized by the immune system, leading to a stronger immune response to ICIs [[110]12]. TMB can serve as a biomarker, particularly in predicting patient responses to ICIs [[111]13]. In this study, we found that patients who had both high TMB and high risk scores had the worst clinical prognosis, indicating that TMB and risk score serve as prognostic factors for HCC patients. Furthermore, immune infiltrating cells, which play a crucial role in tumor development within TME [[112]14,[113]15], were found to be universally upregulated in the high-risk group. The advent of immunotherapy has significantly transformed the treatment approach for HCC patients, particularly those in the intermediate to advanced stages. The combination of immunotherapy and targeted therapy has now become a first-line option for intermediate to advanced stage HCC, including combinations such as atezolizumab + bevacizumab [[114]11,[115]16]. The ESTIMATE algorithm estimates the degree of immune infiltration in tumors based on tumor RNA sequencing data, but it imposes a significant economic burden on patients. In this study, we developed a risk scoring model using only 3 DRGs, which accurately predicts levels of immune infiltration and expression of immunomodulators. The cost of testing the reagents for these three DIRGs is relatively inexpensive. The high-risk group exhibits elevated levels of immune infiltration, suggesting that patients with higher risk scores are more likely to benefit from immunotherapy. In recent years, a study assessing the efficacy of sorafenib in combination with intraarterial FOLFOX (SoraHAIC) has shown positive outcomes as a first-line treatment for patients with portal vein thrombosis [[116]17]. The selection of effective targeted drugs and perfusion chemotherapy drugs will undoubtedly contribute to improving treatment efficacy. By comparing estimated IC50 between groups and analyzing data from [117]GSE109211, we found that the patients with high-risk score display increased sensitivity to commonly used chemotherapy drugs and targeted drugs for HCC, such as fluorouracil, cisplatin, gemcitabine, sorafenib, and doxorubicin. The aforementioned results indicate that our risk scoring system plays a crucial role in guiding the selection of systemic therapeutic drugs. Previous research has demonstrated ANXA2's oncogenic role in HCC. Elevated ANXA2 expression is indicative of advanced HCC stage, unfavorable prognosis, and is correlated with the progression, metastasis, and development of chemotherapy resistance in HCC [[118][18], [119][19], [120][20]]. Our study yielded intriguing findings. Analysis using TISCH revealed remarkably elevated ANXA2 expression patterns across various components of TME among patients receiving PD1/CTLA4 immunotherapy. Integrating our findings, it is plausible to propose that the combination of immunotherapy with disulfidptosis inducers may enhance the response rate. ANXA2 will be a potential molecular biomarker for evaluating the feasibility and efficacy of this combined therapy. Nevertheless, this study has some limitations. Further mechanistic research to investigate the involvement of disulfidptosis in the development of HCC is in urgent need, especifically the role of DIRGs in intercellular communication within TME and the signaling pathways in which DIRGs are involved in regulating immune cells in HCC progression. Due to the limited number of HCC samples, more samples are needed to confirm the generalizability of this model. Moreover, as a retrospective study, it lacks several clinicopathological features that could contribute to the construction of a prognostic risk model. This could potentially underestimate the clinical applicability of the study. Samples with more comprehensive clinicopathological factors can contribute to the construction of a more effective prognostic model. In the subsequent validation of the DIRG model, stratified analyses will be conducted based on backgrounds such as HBV and HCV, with the aim of developing a more clinically meaningful and personalized model. 4. Conclusions We identified disulfidptosis-related molecular subtype in HCC. We developed a three-gene risk scoring system that demonstrated high accuracy in predicting prognosis, effectively reflecting the disulfidptosis status, immune infiltration level, and identifying sensitive drugs. Our research can serve as a valuable theoretical foundation for the utilization of disulfidptosis inducers and the implementation of immunotherapy. Thorough mechanistic research on disulfidptosis in HCC holds the potential for advancing personalized diagnostic and therapeutic approaches, leading to enhanced prognosis. 5. Methods 5.1. Acquisition and processing of data DRGs were obtained from previous studies, totaling 24 genes. RNA sequencing data and clinical data were obtained from the TCGA-LIHC cohort. After excluding samples with incomplete survival information, we obtained a total of 354 HCC samples and 49 normal samples with complete RNA sequencing data and clinical survival information, which were used as the training set for subsequent analysis. Somatic mutation data were obtained from the UCSC database, which allowed us to obtain a total of 341 samples with somatic mutation data for further analysis. The data from the ICGC-LIRI cohort were used as the testing set and processed in the same way. [121]Fig. S1 showed graphical abstract of this research. 5.2. Differential analysis and unsupervised clustering We performed differential analysis on TCGA-LIHC patients using the DESeq2 package. We set the threshold for differentially expressed genes (DEGs) at an absolute logFC >1 and FDR (False Discovery Rate) adjusted p < 0.05. We identified disulfidptosis-related molecular subtypes in LIHC patients by performing unsupervised clustering using the R package consensusclusterplus on the expression profiles of these DEGs. We compared the expression of DRGs between the clusters. Additionally, we compared the differences in clinicopathological factors and immune infiltration levels between the subtypes. 5.3. Immune infiltration and biological functional analysis We used the Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) algorithm to analyze immune infiltration in patients from both TCGA and ICGC cohorts. The degree of immune cell infiltration was evaluated by stromal score, immune score, ESTIMATE score, and tumor purity. We used the Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) algorithm to evaluate the composition of 22 immune cells, and the R package GSVA and single-sample gene set enrichment analysis (ssGSEA) methods to evaluate the expression levels of 28 immune cells [[122]21,[123]22]. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were used for biological function analysis and signaling pathway enrichment analysis. The Gene Set Enrichment Analysis (GSEA) software was used for functional enrichment analysis, with a threshold set at p < 0.05, false discovery rate (FDR) < 0.25, and normalized enrichment score (NES) > 1. 5.4. Weighted gene co-expression network analysis (WGCNA) We used the R package WGCNA to perform WGCNA on differentially expressed genes related to disulfidptosis. We set the scale-free topology fit index to 0.9 and combined network connectivity to obtain an optimal power value of 3 [[124]23]. Subsequently, we established a hierarchical clustering dendrogram to identify module membership and associated each module with ESTIMATE results and molecular subtypes to obtain key modules related to both immune and disulfidptosis. Then gene significance (GS) and module membership (MM) were calculated for each module. 5.5. Least absolute shrinkage and selection operator (LASSO) cox analysis After obtaining key modules related to immune infiltration and disulfidptosis, we conducted univariate Cox regression analysis on the genes in these modules to obtain immune infiltration- and disulfidptosis-related genes (DIRGs) related to HCC prognosis. We then used the glmnet R package to perform LASSO algorithm to identify hub genes and calculate risk scores. The risk score is calculated according to the following formula: [MATH: Riskscore=∑δ=13(Expδ∗Coeδ) :MATH] , where Expδ represents the expression level of each DIRG, and Coeδ represents the corresponding coefficient. Subsequently, patients were divided into high- and low-risk groups based on the median risk score. This formula was used to process data of testing set and obtain corresponding risk groups. 5.6. Acquisition and processing of single-cell RNA sequencing (scRNA-seq) data We selected twelve scRNA-seq samples (P08, P09, P10, P11, P12, P13, P14, P15, P16, P17, P18, P19) of primary liver cancer from the CNP0000650 dataset available on the China National GeneBank DataBase (CNGBdb) for analysis. The samples were filtered based on the criteria: min.cells = 3, min.features = 200 and percent.mt < 5. Data normalization was performed using the LogNormalize method. The FindVariableFeatures function was used to identify the top 2000 genes exhibiting the most significant differences. Principal component analysis (PCA) was carried out utilizing the RunPCA function. The appropriate data dimension, determined as 15 based on the ElbowPlot, was utilized to define cell groups through common cell annotations. Clustering and visualization were performed employing the UMAP method. To validate the expression of DIRGs, we utilized the online tool Single-cell Atlas in Liver Cancer (scAtlasLC). To explore the expression of DIRGs after immune therapy, we employed the online tool Tumor Immune Single-cell Hub 2 (TISCH2). 5.7. Cell culture and quantitative real time polymerase chain reaction (qRT-PCR) HL-7702 (human liver cell line) and SNU-387 (human liver cancer cell line) were obtained from Nanjing Cobioer Biotechnology Corporation and were cultivated in Dulbecco's minimum essential media (DMEM) containing 10 % fetal bovine serum (FBS, WISENT, Canada). The mRNA expression of three DIRGs was detected using qRT-PCR. Total RNA was extracted using TRIzol (T9108, Takara, Dalian, China) and reverse transcribed into complementary DNA using PrimeScript™ RT Master Mix (ES Science, Shanghai, China). SYBR Green Master Mix (ES Science, Shanghai, China) was used for qRT-PCR amplification of the three DIRGs. Finally, the mRNA expression levels were normalized to the ACTB gene,and calculated by 2^−ΔΔCt method. The primers used for PCR amplification can be found in [125]Supplementary Material Table S1. 5.8. Western blot Total protein was extracted from cells using RIPA lysis buffer containing a mixture of protease inhibitors (Beyotime Biotech, Shanghai, China). The protein lysates were separated using SDS-PAGE and then transferred onto PVDF membranes, which were blocked with 5 % skim milk at room temperature for 1 h. The membranes were then incubated overnight at 4 °C with primary antibody anti-beta Actin (ab8227, Abcam, Shanghai, China). Afterwards, the membranes were incubated with recommended dilution of secondary antibodies (Anti-ST6GALNAC4, ab127016, Anti-ANXA2, ab41803, Anti-MSC, ab221097, Abcam, Shanghai, China) in the blocking buffer for 1 h at room temperature. The original, unedited image of the blots in [126]Fig. 9 can be found in Supplementary Material. 5.9. Statistical analysis All statistical analyses were performed using R software version 4.2.2 or GraphPad Prism 9. For normally or approximately normally distributed data, comparisons between two groups were conducted using t-test, while for skewed data, the Wilcoxon test was used. Kaplan-Meier curves were generated for survival analysis visualization. Spearman correlation coefficient was used for correlation analysis. P-value <0.05 was considered statistically significant. CRediT authorship contribution statement Zhe Xu: Writing – original draft, Validation, Software, Methodology, Formal analysis, Conceptualization. Chong Pang: Writing – original draft, Visualization, Validation, Conceptualization. Xundi Xu: Writing – review & editing, Supervision, Methodology, Data curation, Conceptualization. Ethical statement Not applicable. Data availability statement The datasets analysed during the current study are available from the corresponding author on reasonable request. Consent for publication Not applicable. Funding This work was supported by the National Natural Science Foundation of China (81670111) and Hunan Provincial Key Research and Development Program (2019SK2242). Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments