RNA sequencing data and corresponding clinical information were obtained from TCGA-LIHC (n = 315). scRNA-seq data for [20] samples were obtained from GSE189903 and GSE242889. Samples lacking survival information or with a survival duration of $\le$30 days were excluded. Those exhibiting $>$30% missing values were also removed to ensure data quality. For RNA-seq data, transcripts per million values were extracted and normalized to log₂(TPM + 1) transformation. Replicate gene measurements were averaged. Batch effects between the TCGA and GEO datasets were corrected using principal component analysis (PCA). Only samples possessing both RNA-seq data and clinical annotations were retained for subsequent analyses.

Disulfidptosis-related gene set compiled from the literatures were integrated and defined as disulfidptosis feature genes (DFGs) and part of the gene set participating in cellular uptake and regulation of the actin cytoskeleton. The RNA-seq expression matrix was subsequently divided into mRNA and lncRNA expression matrices. LncRNAs co-expressed with DFGs were selected based on Pearson correlation coefficients ($|$R$|$ $>$ 0.4, adj. p 0.05) [16, 17].

Unsupervised clustering of TCGA-LIHC samples was performed with the non-negative matrix factorization (NMF) algorithm through the DFG-associated lncRNA expression matrix. NMF is a clustering method widely used for cancer molecular subtyping using gene expression data. This study uses the NMF algorithm to distinctly categorize disulfide death subtypes and analyzes their biological traits. Some categories prominently display disulfide death characteristics, while others are compared to the core subgroup from various perspectives. Kaplan-Meier survival curves were generated to assess the prognostic significance of identified subtypes.

A disulfidptosis score (Dis-score) was constructed via PCA utilizing DFGs. To verify the stability of NMF clustering results under data sampling variation, Bootstrap resampling combined with a multi-dimensional stability index was used for quantitative analysis. Based on the pre-processed lncRNA expression matrix (no missing values, all positive, filtered low expression features) and the optimal clustering rank (rank = 3, corresponding to C1, C2, and C3 subtypes), 1000-fold re-sampling was performed, with a subset of samples consistent with the original sample size extracted each time. For each resampling subset, the Brunet algorithm was used to perform NMF clustering (each clustering was repeated three times to avoid a local optimal solution), and the sample clustering label and consensus matrix were recorded. The core stability index was further calculated: Adjusted Rand Index (ARI) between iterations was used to quantify the similarity of any two resampled clustering structures (range –1~1, ARI $\ge$0.7 indicates structural stability). Finally, the ARI distribution density map is visualized.

Based on transcriptome data from TCGA-LIHC, we collected clinical information and multiple immune index data from different samples, including: (1) immune cell infiltration level (B cell, M1/M2 macrophages, neutrophils, and NK cell abundance as estimated by the QUANTISEQ method); (2) expression of immune checkpoint molecules (e.g., ICB-CD274, ICB-CTLA4); (3) immune-escape related score (TIDE, Dysfunction, Exclusion). Missing values in continuous variables were filled by medians, and samples with missing categorical variables were excluded. Spearman’s rank correlation analysis was used to evaluate the correlation strength between Dis-score and immune indexes, and the correlation coefficient (r) and significance (p-value) were calculated.

Three datasets of scRNA-seq raw data obtained from the GEO database were processed using the Seurat package (v4.3.0). The details for data preprocessing, including quality control (QC), normalization, selection of highly variable genes, and annotation are as follows:

(1) Raw Data Quality Control: To exclude low-quality cells, empty droplets, and cell aggregates, cells were filtered based on the number of detected genes per cell between 200 and 6000, and mitochondrial gene proportion 20%.

(2) Normalization, Data integration, and Dimensionality Reduction: SCTransform was applied for variance stabilization transformation, and simultaneous batch effect removal. To further eliminate technical variation across multiple samples, data integration was performed using the Harmony algorithm (v0.1.1). Following preprocessing, PCA was conducted based on highly variable genes.

(3) Cell Clustering: The first 16 principal components were used to construct a shared nearest neighbor (SNN) graph. Cell clustering was performed using the Louvain algorithm with a resolution parameter set to 0.4. Non-linear dimensionality reduction was achieved via Uniform Manifold Approximation and Projection (UMAP) for the visualization of clustering results.

(4) Cell Annotation: Cell type annotation was performed according to the expression patterns of canonical marker genes: epithelial cells (EPCAM, KRT19) [21], T cells (CD3D, CD3E), B cells (CD79A, MS4A1), myeloid cells (CD68, CD163), endothelial cells (PECAM1, VWF), and fibroblasts (COL1A1, ACTA2) [21, 22, 23].

To quantify the potential disulfidptosis status of tumor cells in HCC at the single-cell level, we constructed a scoring model based on gene characterization. First, we selected 15 genes explicitly associated with the core process of disulfidptosis based on published literature. These covered several key aspects, such as cystine uptake, cytoskeletal rearrangement, and metabolic adaptation (Table 1). The standard procedure for the R package “Seurat” was used to extract all annotated tumor cells (2238 in total) from the integrated single-cell data for subsequent analysis. For each tumor cell, the average Z-score for the expression of these 15-disulfidptosis-related genes was calculated at the cellular level. This was defined as the raw disulfidptosis death score of the cell. Prior to the calculation, all genes were normalized for overall expression levels within the tumor cell population.

Using the disulfidptosis phenotype of HCC tumor cells as the primary grouping criterion, we systematically characterized the metabolic reprogramming features of high/low-disulfidptosis score tumor cells across energy metabolism, biosynthesis, and antioxidant regulation. We constructed nine key metabolic pathway gene sets spanning three functional modules: energy metabolism (glycolysis, oxidative phosphorylation, tricarboxylic acid cycle), biosynthesis (pentose phosphate pathway, fatty acid metabolism, amino acid metabolism), and redox balance (glutamine metabolism, antioxidant system, disulfidptosis-related metabolism). The average expression of each metabolic pathway gene set was calculated and normalized, and Z-scores were applied to derive normalized metabolic pathway scores. Student’s t-test was used to compare pathway scores and compute fold changes between high- and low-risk groups, while the Wilcoxon rank-sum test was used to determine statistical significance. Heatmap visualization, cluster analysis, and functional grouping were integrated to comprehensively examine the significantly altered metabolic pathways.

Macrophage subsets were extracted from the total cell population for secondary clustering analysis. Three core subsets were defined based on the expression of characteristic markers: M1 macrophages (high expression of NOS2, CD80, CD86), M2 macrophages (high expression of CD163, MRC1, MS4A4A), and monocytes (high expression of CD14, FCGR3A).

Single-cell differentiation trajectories were constructed using the Monocle3 package (v1.3.4) with the following workflow: Seurat objects were converted to CellDataSet objects, followed by preprocessing and dimensionality reduction based on highly variable genes; UMAP was then used for non-linear dimensionality reduction, and cell states were identified via Louvain clustering; finally, monocytes were selected as the trajectory starting point, and pseudotime values were calculated using the order_cells function. Pearson correlation analysis was performed to evaluate the dynamic expression of disulfidptosis-related genes along pseudotime, and the statistical significance of correlation coefficients was tested.

To decipher cellular crosstalk between “High-Disulf tumor cells and M2 macrophages” in the HCC TME, we employed the CellChat package integrated with the human CellChatDB database to reconstruct communication networks. We specifically focused on interaction networks between tumor cells with high disulfidptosis risk and M2 macrophages. To systematically characterize directional interactions, five core functional communication axes and their corresponding ligand–receptor pair databases were predefined, encompassing recruitment and polarization, survival support, stress regulation, metabolic crosstalk, and promotion of invasion. Using communication probability (prob) as the core metric, interaction pairs within each axis were ranked in descending order of probability. The average probability, maximum probability, and number of interactions per axis were computed. One-way ANOVA was applied to assess differences in communication strength across axes, with a significance threshold of p 0.05. Finally, network diagrams and bubble plots were used to visualize the specificity of these crosstalk events.

A ferroptosis score was constructed based on a functionally validated 10-gene signature derived from previous studies. This encompassed core antioxidant regulators (GPX4, SLC7A11, FSP1, NRF2), iron metabolism-related genes (TFRC, FTH1), a lipid metabolism regulator (ACSL4), an inflammation-associated gene (PTGS2), and signaling mediators (STAT3, CD44). All genes were stably detected in the current dataset. The ssGSEA algorithm and parameters identical to those used for disulfidptosis scoring (implemented via the GSVA package) were applied to compute ferroptosis scores at the cell level, with normalized enrichment score (NES) values serving as the final quantitative metric, ensuring comparability between the two scoring systems.

Pearson correlation analysis was conducted to assess the linear relationship between disulfidptosis and ferroptosis scores in the tumor cell population. Within the M2 macrophage subpopulation—a core functional subset—the disulfidptosis score was further decomposed into three phases: Stage1 (DiSulfide_Stress), Stage2 (Core_Regulation), and Stage3 (Downstream_Execution). Pearson correlation analysis and correlation strength difference tests (Z-test) were subsequently performed between each phase-specific score and the ferroptosis score. One-way ANOVA was used to compare the overall distribution of ferroptosis scores across four groups: tumor cells, M1 macrophages, M2 macrophages, and monocytes.

To elucidate the significant positive correlation observed between the two scores in M2 macrophages, we introduced key intersection genes involved in both ferroptosis and disulfidptosis: RAC1, GPX4, NRF2, and SLC7A11. Among these, NRF2 and SLC7A11 were considered potential upstream regulators in a ternary correlation analysis. Pairwise correlation coefficients were computed for all gene pairs, and partial correlation analysis was performed using the ppcor package (v1.1).

For experimental verification of the co-expression axis, samples of tissue sections from 49 patients with primary HCC were collected at the First Affiliated Hospital of Anhui Medical University from 2023 to 2025. Multiplex immunofluorescence staining was performed to validate enhanced GPX4-RAC1 co-localization in M2 macrophages (CD68⁺, CD163⁺) versus other subsets, supporting their co-regulatory function.

Samples from the TCGA database with complete survival data (overall survival time and status) and GPX4/RAC1 expression values were included in the analysis. GPX4 and RAC1 expression levels were quantified using normalized transcriptomic data. The co-expression score was calculated using PCA. The first principal component (PC1) derived from GPX4 and RAC1 expression values was used as an alternative co-expression metric. Kaplan-Meier survival curves were generated to compare overall survival (OS) between high and low co-expression groups, with the log-rank test employed to assess statistical significance. Univariate and multivariate Cox proportional hazards models were constructed to evaluate the prognostic value of GPX4–RAC1 co-expression. The multivariate model was adjusted for established clinical prognostic factors, including pTNM stage and tumor grade. Hazard ratios (HRs) with 95% confidence intervals (CIs) were reported. SCT-standardized GPX4 and RAC1 gene expression data at the single cell level were extracted for the identified M2 macrophage subsets. Pearson and Spearman correlation analyses were performed to assess the correlation strength and significance of gene expression. Scatter plots and linear fitting were used to visualize the correlation results and verify the robustness of the correlation.

All statistical analyses were performed using R software (v4.4.2). Wilcoxon rank-sum test, or analysis of variance was selected for intergroup comparisons based on data distribution characteristics and the homogeneity of variance test results. Pearson or Spearman methods were used for correlation analysis. False discovery rate (FDR) control was applied to correct for multiple comparisons. Data visualization was performed using the ggplot2 package (v3.4.2), and network visualization was achieved using the igraph package (v1.4.2). Statistical significance was defined as p 0.05, with special cases indicated separately.

Through the integration of TCGA-LIHC transcriptomic data with an established set of DFGs (Table 1), 12 significantly correlated lncRNAs were identified using Pearson correlation coefficients (Supplementary Fig. 1). NMF clustering stratified patients into three subgroups: Subtype 1 (n = 99), Subtype 2 (n = 132), and Subtype 3 (n = 84) (Fig. 1A and Supplementary Fig. 1B). Subtype 1 exhibited significantly superior OS compared to Subtype 3 (p = 0.04), while Subtype 2 (intermediate) demonstrated intermediate prognosis (Fig. 1B). Analysis of the immune infiltration score revealed significant differences in the infiltration of myeloid cells (macrophages and neutrophils), B cells and T cells (CD8⁺ T cells and Tregs) between the three subtypes (Fig. 2E). With 1000 Bootstrap re-samplings performed, subset consistent with the original sample size were extracted each time. The ARI index shows the distribution characteristics of the similarity of cluster structures between any two resampling iterations, with an ARI range of –1 to 1. The closer to 1, the more consistent the clustering structure is indicated (the red dotted line in the figure represents the stability threshold of 0.7). The stability results of our classification show that the average ARI between iterations is 0.631 (Fig. 1E).

Fig. 1.

Identification of HCC molecular subtypes based on DR-LncRNAs and the construction of Dis-score. (A) Identification of three HCC subtypes based on expression profiles of 12 DR-LncRNAs. Cophenetic correlation coefficient for k = 2–7 for the cohort. (B) Survival differences (OS) of three subtypes from TCGA cohort. (C,D) Dis-score was constructed based on PCA algorithm. Dis-score in the box plot was compared using the Wilcoxon test. And the area under curve (AUC) can evaluate the model discrimination degree. AUC $>$0.6: May make sence; (E) shows the density map of the similarity of cluster structure (ARI) between iterations. The red dotted line in the figure represents the stability threshold of 0.7. The ARI index is close to the stable threshold, suggesting that the overall structure of NMF clustering has moderate consistency under data sampling variation. (F) Volcano plot of differential expression of genes between subtype 1 and 3. (G) Immune score of the samples in different subtypes was compared using the Wilcoxon test and visualized by the box plot. (H,I) The correlation analysis results of Dis-score with multiple immune cell infiltration scores, immune checkpoint expression levels, and immunotherapy response scores in different samples. (H) heatmap, (I) histogram. (J) The box plot showed the difference in the degree of immune cell infiltration and the expression of immune checkpoints between the high and low score groups. -, no significance, *p 0.05, **p 0.01. ***p 0.001, ****p 0.0001. NMF, non-negative matrix factorization; PCA, principal component analysis; HCC, hepatocellular carcinoma; LncRNAs, long non-coding RNAs; Dis-score, disulfidptosis score; OS, overall survival; ARI, Adjusted Rand Index.

Fig. 2.

Differences in biological characteristics between disulfide-dead subtype and disulfide-stable subtype. (A) Differentially expressed biological pathways in disulfidptosis subtypes (KEGG database). (B–D) WGCNA analysis in the disulfidptosis subtype. WGCNA co-expression network identified key modules based on lncRNAs and disulfidptosis related genes. The clustered modules of WGCNA, related to Dis-score, OS and macrophage M2. (B) Cluster trees corresponding to Pearson-based WGCNA. (C) The hierarch clustering on the clinical features. (D) Cell infiltration and all WGCNA modules. (E) Based on the Tip algorithm, the box plot shows the expression differences of different immunotherapy response targets between the three subtypes. *p 0.05, **p 0.01. ***p 0.001, ****p 0.0001. WGCNA, weighted gene co-expression network analysis.

A composite scoring system (Dis-score) was derived via PCA of DFGs, with principal component 1 (PC1) accounting for 68% of variance. The Dis-score was significantly higher in Subtype 1 versus Subtype 3 (p-value = 9.6 $\times$ 10^-7; Fig. 1C,D). Receiver operating characteristic (ROC) analysis showed high discriminative power of the Dis-score for disulfidptosis subtypes (AUC = 0.71), supporting clinical utility. However, it showed significant positive correlations with multiple immune markers, especially immune checkpoint genes (HAVCR2, CTLA4, and PDCD1), immune cell infiltration (B cells, M1/M2 macrophages, and MDSC), and the TIDE score. The volcano plot highlights pathways enriched by differentially expressed genes between Subtype 1 and Subtype 3, with upregulated pathways like actin regulation and ECM-receptor interaction (Fig. 1F). Additionally, a box plot using the Wilcoxon test shows that Subtype 1 has the highest M2 macrophage infiltration level (Fig. 1G). The infiltration scores of B cells, M1/M2 macrophages, TIDE scores, and most immune checkpoints were significantly higher in the high Dis-score group compared to the low Dis-score group (Fig. 1H–J).

Metabolomic validation revealed that Subtype 1 was enriched in glucose (log₂FC = 3.1), glutamine (log₂FC = 2.8), and Fe³⁺ (log₂FC = 4.5), whereas Subtype 3 accumulated pyruvate (log₂FC = –2.3) and acetate (log₂FC = –1.9). This pattern reflects the reliance of Subtype 1 on the glycolysis-glutaminase axis to support NADPH regeneration (Supplementary Fig. 1E). KEGG pathway enrichment analysis revealed significant differences in pathway activities between Subtype 1 (Disulfidptosis subtype) and Subtype 3 (Stabilization subtype) (Fig. 2A). Subtype 1 showed marked enrichment in pathways related to extracellular matrix (ECM)-receptor interaction, focal adhesion, the p53 and mTOR signaling pathways, and regulation of the actin cytoskeleton. In contrast, Subtype 3 was characterized by significant enrichment in oxidative phosphorylation and peroxisome pathways. This indicates a shift toward efficient energy metabolism and redox homeostasis, which may contribute to stability and cell survival. The above findings highlight distinct molecular mechanisms underlying the two subtypes, with Subtype 1 favoring disulfidptosis through ECM and stress signaling activation, while Subtype 3 promotes cellular stability via metabolic conservation and antioxidant responses.

We next analyzed the co-expression relationship between DFGs and lncRNAs in HCC, as revealed by WGCNA (Fig. 2B–D). The turquoise module was found to contain most of the genes related to disulfidptosis, and showed a high correlation with the Dis-score (Fig. 2B). The midnight blue module contained 64 lncRNAs and showed significant positive correlations with the Dis-score, M2 macrophage infiltration, B cell infiltration, Treg cell infiltration, and survival status. Furthermore, there was little intersection between the midnight blue and turquoise modules (Fig. 2E).

scRNA-seq of tumor and adjacent non-tumor tissues from 21 HCC patients identified 58,842 cells (Fig. 3A–C) clustered into 9 major cell types (Fig. 3D,E). Standardized analysis was conducted using the Seurat package. Low-quality cells were eliminated through QC to ensure data reliability (Supplementary Fig. 2), and samples were successfully de-batched and merged based on harmony. Based on the results of PCA (Fig. 3A–C and Supplementary Fig. 3), the first 15 principal components were selected to construct a shared nearest neighbor (SNN) map between cells. The Louvain algorithm was used for cell clustering, with the resolution parameter set to 0.4 (Fig. 3D, Supplementary Fig. 4A,B). Finally, 9 cell types with a clear boundary and uniform internal structure were obtained by clustering and annotation (Fig. 3E, Supplementary Fig. 4C,D). Markers used in the annotation process showed significant expression distribution (Supplementary Fig. 5).

Fig. 3.

Identification of cell clusters and analysis of cell population characteristics. (A–C) The PCA algorithm is used to reduce the dimension of the data, and then the harmony algorithm is used to remove the batch effect. (D,E) UMAP dimensionality reduction clustering was performed on 21 amples from GEO database, and 9 kinds of cell subsets were obtained. (F) The UMAP showed the intensity of disulfidptosis scores in cell levels in HCC tissues. (G) The violin plot showed the difference in the distribution of disulfidptosis scores after tumor cell grouping. (H) The bubble diagram representeed the difference in the expression of core genes between different tumor cell populations in tumor tissues. (I) The enrichment analysis result of differential genes obtained between high and low score groups. (J) Subpopulation and monocle analysis of MDCSs. (K) Monocle analysis of the top 5 significantly differentially expressed genes in pseudo-time. (L) The distribution of marker expression between different cell subsets. PCA, principal component analysis; UMAP, Uniform Manifold Approximation and Projection.

Using tumor-specific markers (high EPCAM expression and exclusion of normal epithelial traits), 2238 pure tumor cells were analyzed for disulfide death correlation. Of 15 related genes, 14 were detectable, except CHOP (Fig. 3H and Supplementary Fig. 6A). Disulfide death scores were calculated using ssGSEA, with the median (–0.03) dividing cells into high (1119 cells) and low (1119 cells) score groups (Fig. 3G). Enrichment analysis of the differentially expressed genes revealed changes in metabolism and transmembrane transport (Fig. 3I). UMAP visualization showed these groups were spatially scattered and interwoven, indicating that disulfide death heterogeneity is not directly linked to cell clustering but may reflect tumor cell adaptation to the microenvironment (Fig. 3F). Differential expression analysis was conducted on these groups (FDR 0.05, $|$FC$|$ $>$1.5) (Supplementary Fig. 6C). In the High_Disulf_group, 303 genes showed significant expression changes: 238 were up-regulated and 65 down-regulated, indicating notable gene expression activation in tumor cells with a high disulfidptosis score.

To further analyze the differentiation potential and functional heterogeneity of macrophage subsets in myeloid cells, a focused analysis was performed based on the bone marrow-derived suppressor cell (MDSC) population. In total, 5806 related cells were included. Subgroup segmentation and in-depth characterization were performed by combining marker gene integrity with the scoring model. With the core marker gene set expression and function score, the macrophage-related cells in the MDSC population were divided into three categories: M1 macrophages (n = 1980, 34.1%), M2 macrophages (n = 2237, 38.5%), and monocytes (n = 1589, 27.4%) (Fig. 3J–L). In the early stage of pseudo-time (Supplementary Fig. 7A–D), the cells showed high expression of the monocyte marker CD14. The middle stage showed enrichment with M1 macrophages (high expression of CD86) and M2 macrophages (high expression of CD163) (Fig. 3K and Supplementary Fig. 7H). Timing analysis of core genes showed that the key ferroptosis gene GPX4 was synergistically up-regulated with the disulfide death regulatory gene RAC1 (pseudo time 0.00 $\to$ 0.75). No significant fluctuation was observed in the expression of SLC7A11 (Supplementary Fig. 7E–G).

From the complete data set (58,842 cells), the core cell population was subjected to directional filtering. This resulted in 1119 high-disulfide tumor cells and 2237 M2 macrophages, giving rise to a total of 3356 cells for communication analysis (Fig. 4A,B). Three groups of secretory signal communication pairs with relative intensity differences were detected in the M2$\to$tumor direction, corresponding to the GALECTIN pathway (LGALS9-PTPRC, LGALS9-P4HB) and CypA pathway (PPIA-BSG). The average communication intensity in the tumor$\to$M2 direction was higher than in the M2$\to$tumor direction (Fig. 4D). A total of 11 communication pairs were detected in the tumor$\to$M2 direction. These covered a variety of signal modes and core pathways, with the majority being secretory signals. Core interaction pairs were migration inhibitory factor (MIF) pathway-related MIF-(CD74 + CXCR4), MIF-(CD74 + CD44), SPP1-CD44, MDK-LRP1, MDK-NCL, C3-C3AR1, ANXA1-FPR3, PPIA-BSG, and VTN-PLAUR. The core pairs covered secretory signals, ECM-receptor signals, and other modes (Fig. 4C). Additionally, in the tumor$\to$M2 direction, the top five pairs of interaction signals according to communication intensity were: MIF-(CD74 + CXCR4) $>$ MIF-(CD74 + CD44) $>$ PPIA-BSG $>$ MDK-NCL $>$ MDK-LRP1. The average communication probability of the MIF pathway was higher than that of other pathways. Moreover, the two groups of communication pairs in this pathway represent core interactions that have key regulatory significance. At the same time, the CypA pathway (PPIA-BSG) and GALECTIN pathway in the M2$\to$tumor direction were both statistically verified and included in the core network, although their intensity was low (Fig. 4D,E).

Fig. 4.

Multi-dimensional characteristics analysis between M2 macrophages and tumor cell subsets. (A) Tumor cells and M2 macrophages were included in the analysis. (B) The intensity of the bidirectional interaction between M2 macrophages and tumor cells was compared. (C) The heat map illustrated the interaction intensity between the two cell types in various directions. (D) The histogram indicated the top eight signaling axes in terms of the total interaction intensity between the two cell types. (E) The expression levels of key receptor and ligand genes across different cell populations were examined. (F) Differentially expressed metabolic-related biological processes were identified in tumor cell populations characterized by disulfide-induced cell death. (G) The top four differential metabolic processes were presented. NS, no significance.

We next analyzed the 2238 pure tumor cells (1119 High-Disulf group cells and 1119 Low-Disulf group cells) that had previously been rigorously filtered. Focusing on the core dimensions of “energy metabolism-biosynthesis-antioxidant regulation”, feature gene sets were constructed for 9 classical metabolic pathways (Table 2 and Fig. 4F). Each of these pathways showed significant differences between High- and Low-Disulf tumor cells (all pathways p 2 $\times$ 10^-16, except for glutamine metabolism, p = 1.5 $\times$ 10^-6) (Fig. 4G and Supplementary Files). The overall pattern was characterized by “specific upregulation of disulfidptosis-associated metabolism and general downregulation of traditional energy metabolism”.

Specifically, all five core genes in the Disulfidptosis_Metabolism pathway (SLC7A11, SLC3A2, GLUD1, G6PD, RAC1) were highly expressed in the High-Disulf group, with SLC7A11, G6PD, and RAC1 showing the most pronounced upregulation. The remaining 8 pathways were downregulated to varying degrees, with glutamine metabolism showing the most significant downregulation. Although the Antioxidant System pathway was downregulated overall (FC = 0.89), core antioxidant genes GPX4 and TXNRD1 were highly expressed in the High-Disulf group (Fig. 4G and Supplementary Files).

The ferroptosis score exhibited a significant heterogeneous distribution across different cell populations (Fig. 5G). The median ferroptosis score in the High-Disulf group cells was significantly higher than in the Low-Disulf group cells. Within other compartments, myeloid cells had the highest ferroptosis score, followed by M1 macrophages (Fig. 5A–E). The disulfidptosis score was divided into three stages and subsequently correlated with the ferroptosis score (Fig. 5F). All three stages showed significant positive correlations (all p 0.001), with a clear hierarchy in correlation strength: Stage2_Core_Regulation showed the strongest correlation with the ferroptosis score (r = 0.65). Expression of core genes in this stage (e.g., RAC1, NCKAP1) showed synchronized upregulation (Supplementary Fig. 8C–E and G) with the ferroptosis core gene GPX4. Stage1_DiSulfide_Stress ranked next, while Stage3_Downstream_Execution showed a relatively weaker correlation (r = 0.522). Cytoskeleton-related genes (e.g., ACTB, CFL1) in this stage showed strong associations with RAC1 and GPX4 (Supplementary Fig. 8F,H).

Fig. 5.

The results of correlation analysis between disulfide death and ferroptosis in multiple stages of HCC. (A–D) Correlation analysis of disulfide death score and ferroptosis score in different cells. (E) The distribution of ferroptosis scores in tumor cells with high and low disulfide death scores was different. (F) The correlation between the corresponding scores of the three stages of disulfide death and the ferroptosis score was demonstrated. (G) The UMAP map showed the distribution trend of the ferroptosis score in all cell populations. (H,I) Correlation analysis and strength comparison between SLC7A11, RAC1 and GPX4. (J) The correlation result was drawn from SCT-standardized data. (K) Immunofluorescence results in macrophages. Scale bar = 10 µm. CD68 red immunofluorescence, CD163 green immunofluorescence, GPX4 yellow immunofluorescence, RAC1 blue immunofluorescence. (L) Correlation analysis of GPX4-RAC1 expression based on fluorescence scores of 46 patients. (M,N) KM curve analysis of Dis-score in clinical features. (M) divided by the median core, (N) divided by the best cutoff values.

Based on the strong correlation (r = 0.653) between the core regulatory stage (Stage2) of disulfidptosis and ferroptosis in M2 macrophages identified previously (Fig. 5F), we further analyzed the core effector molecules of each cell death mode. Direct association analysis of the key ferroptosis suppression gene GPX4 and the core disulfidptosis regulatory gene RAC1 was performed to clarify inter-molecular regulatory relationships. Pearson correlation analysis of SCT-normalized expression data from 2237 M2 macrophages revealed a significant, moderate positive correlation between GPX4 and RAC1 expression (r = 0.444, p 2.2 $\times$ 10^-16, Fig. 5J). The analysis included 2237 valid cells (no NA values), ensuring sufficient statistical power. SLC7A11 is a candidate molecule involved in both ferroptosis inhibition (cystine transport) and disulfidptosis regulation (redox homeostasis). Association analysis showed only weak correlations between SLC7A11 and GPX4 (r = 0.08, p = 0.003), and between SLC7A11 and RAC1 (r = 0.029, p = 0.12), indicating weak or non-significant associations (Fig. 5H and Fig. 5I). Furthermore, bulk RNA-seq data and clinical information from the TCGA database were used to construct the co-expression score. Kaplan-Meier analysis revealed that co-expression of GPX4-RAC1was a risk factor for OS (Fig. 5M,N).

Immunofluorescence staining of M2 macrophages using CD163 and CD68 as markers revealed significant co-expression of GPX4 and RAC1 within this cell population, further validating the results of pseudo-time trajectory analysis (Fig. 5K,L).