A total of 1,653 AML samples and 337 normal samples were included in this study. AML samples included The Cancer Genome Atlas-Acute Myeloid Leukemia (TCGA-LAML) cohort and seven GEO cohorts (GSE10358-GPL570, GSE12417-GPL96, GSE12417-GPL570, GSE37642-GPL96, GSE37642-GPL570, GSE71014-GPL10558), GSE14468-GPL570) (Supplementary Table S1), and normal samples were Genomic tissue expression (GTEx)-whole blood cohort. For GEO cohort data from the affymetrix platform, raw “CEL” files were downloaded and normalized with the use of the robust multiarray averaging (RMA) method, whereas microarray data from the other platforms were directly downloaded with a normalized matrix file. TCGA-LAML and GTEx RNA-seq data (RSEM TPM) were downloaded from the UCSC XENA database (https://xenabrowser.net/datapages/). The “human.gtf” file was adopted to raw matrix annotation. All data were analyzed using R x64 and the associated R Bioconductor software package, and the data information is shown in Supplementary Table S1. Ten DRGs were retrieved from the study by Gan et al., of which, SLC7A11, SLC3A2, RPN1, and NCKAP1 are drivers, and GYS1, NDUFS1, OXSM, LRPPRC, NDUFA11, and NUBPL are suppressors.

We used the ssGSEA algorithm to calculate enrichment scores for disulfidptosis drivers and suppressors, defined as disulfidptosis positive score and negative score, respectively, and subtracted the negative score from the positive score to obtain disulfidptosis activity score (Subramanian et al., 2005).

Based on the expression profiles of 10 DRGs, the “ConsensusCluster” package was used to perform unsupervised clustering of the TCGA-LAML dataset by consensus clustering method (Wilkerson and Hayes, 2010), and two cluster subtypes with significant differences were obtained. We performed 1,000 replicates to ensure stable and reliable clustering.

WGCNA can assess patterns of gene expression correlations and perform methods for visualization of co-expression networks. We used the “WGCNA” software package to identify genes associated with disulfidptosis scores in the TCGA-LAML cohort (Langfelder and Horvath, 2008). Pearson correlation analysis was used for adjacency matrix formation for all paired genes, and a scale-free topology of the adjacency matrix was implemented based on the optimal soft threshold power. Then, we further transform the adjacency matrix into the topological overlap matrix (TOM). Based on the TOM difference measure, a minimum module size of 30 and a cut height of 0.2 were set to partition genes with similar expression patterns into the same modules by average linkage hierarchical clustering. Then, the correlation between module eigengenes (MEs) and disulfidptosis score was evaluated, and the modules that met the study purpose were determined according to the degree of correlation.

For the target module genes identified by WGCNA, Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was used to identify gene functions. The Gene Set Variation Analysis (GSVA) algorithm was used to calculate the activity scores of individual gene sets in different samples, and the Gene Set Enrichment Analysis (GSEA) algorithm was used to evaluate the difference in pathway activity between patients with high- and low-risk score groups to analyze the biological differences among patients with different risk scores. We used q value < 0.05 as a threshold for significant enrichment.

We used the ESTIMATE algorithm (Yoshihara et al., 2013) to evaluate the immune and stromal scores for each AML sample and applied the CIBERSORT algorithm (Newman et al., 2015) to determine the proportion of immune cell subsets in each sample.

We performed univariate Cox regression analysis of disulfidptosis activity score related genes based on p < 0.01 was used to identify the DASRGs significantly related to AML prognosis. In order to limit the influence of multicollinearity between variables, the least absolute shrinkage and selection operator (LASSO) regression analysis was used to further reduce the dimension and screen out the optimal variables to prevent instability caused by model estimation distortion, so as to construct an accurate prognostic risk score model. The risk score for each sample was obtained by multiplying the expression value of each model gene with its corresponding coefficient and adding it. Then, the risk scores of all patients were ranked, and AML patients were divided into high-risk score group and low-risk score group based on the optimal cut-off value, and the differences in clinicopathological factors and biological characteristics between the two groups were further analyzed.

We download the somatic mutation data from the TCGA database (https://portal.gdc.cancer.gov/), and compared the mutation differences in high- and low-risk score groups. The “pRRophetic” package was used to predict the half maximal inhibitory concentration (IC50) of AML samples to commonly used therapeutic drugs (Geeleher et al., 2014). A smaller IC50 value indicates a better treatment effect. We further used the SubMap (https://cloud.genepattern.org/gp) algorithm to predict the response of different risk score groups to anti-PD-1 and anti-CTLA4 immune checkpoint inhibitors.

We downloaded AML single-cell sequencing data containing 21 cell types (GSE116256) from the GEO database, as well as another group of AML Single-cell sequencing data in the context of PD-1 blocking (GSE198052). We referred to previously published literature related to single cells (Jiang et al., 2022; Zheng et al., 2023). The 10 × scRNA-seq data were processed by R software according to a standardized procedure. The original gene expression matrix was introduced into the “Seurat” package for processing, only genes expressed in at least three single cells, and cells with unique molecular identifiers (UMI) counts <200 were removed. Moreover, only cells expressing more than 1,500 genes and less than 6,000 genes were included. The percentage of mitochondrial or ribosomal genes was calculated for each cell, and cells with more than 20% mitochondrial gene expression were considered low quality cells and also not subjected to downstream analysis. Then, normalized counts were obtained by using the library size normalization of the original matrix, and the top 2,000 genes with a large coefficient of variation were obtained by using the “FindVariableFeatures” function. After z-score processing, principal component analysis (PCA) was performed based on high-variable genes. The uniform manifold approximation and projection (UMAP) algorithm was used to realize the visualization of clustering. Cell types were referred to the annotation file provided by van Galen et al. (2019).

Clinical samples were collected in accordance with the Declaration of Helsinki and institutional guidelines, and informed consent was obtained from each patient and healthy volunteer. Ethical approval was obtained from the Ethics Committee of the Second Affiliated Hospital of Nanchang University [No. review. (2018) No. (092)], and all experimental protocols and methods were performed in accordance with relevant protocols and regulations. Five samples from patients with newly diagnosed chronic myeloid leukemia without any previous treatment, five samples from patients in blast crisis, and five normal samples from healthy volunteers were collected according to the World Health Organization classification of tumors of hematopoietic and lymphoid tissues. Detailed details of sample collection, next-generation sequencing, and processing procedures are available in our previous publications (Li et al., 2020).

Wilcoxon and Kruskal–Wallis tests were used for between-two and multiple-group comparisons, respectively. The “survminer” package was used to determine the optimal cut-off value. The number of patients in a single risk group was set to be no less than 30% of the total population. Kaplan-Meier survival analysis was performed using the log-rank test. The receiver operating characteristic (ROC) curve was used to evaluate the specificity and sensitivity of the risk score, and the area under the curve (AUC) was determined. Bilateral p < 0.05 was considered statistically significant.

In comparison with normal samples, DRGs exhibited upregulated expression in AML samples, suggesting potential crosstalk between DRGs and AML (Figure 1A). Expression correlation analysis indicated positive correlations among most DRGs, with suppressor NDUFA11 showing a negative correlation with drivers SLC7A11 and NCKAP1. Additionally, NCKAP1 and GYS1 exhibited a negative correlation (Figure 1B), suggesting the presence of an antagonistic regulatory mechanism in disulfidptosis. Cox regression analysis identified DRGs as risk factors, except for NUBPL (Figure 1C). K-M curve analysis demonstrated that high expression groups of SLC7A11, SLC3A2, OXSM, NDUFA11, and NDUFS1 had significantly worse prognoses than low expression groups, while the high expression group of NUBPL had a significantly better prognosis (Supplementary Figures S1A, B). PPI network analysis revealed NDUFS1 as the hub gene among the suppressors (Figure 1C). Furthermore, a set of TFs with potential regulatory roles with DRGs was identified (Figure 1D). Using the GSVA algorithm, disulfidptosis positive, negative, and activity scores were computed. Compared to normal samples, AML samples exhibited higher disulfidptosis positive and activity scores and lower negative scores, indicating AML cell resistance to disulfidptosis (Figure 1E). K-M curve analysis showed that patients in both high positive and activity score groups had worse prognoses than those in low score groups, with the opposite observed for the negative score (Figures 1F–H). The occurrence of disulfidptosis was found to be influenced by glucose starvation, with glucose metabolic pathways such as citrate cycle (TCA cycle), glycolysis gluconeogenesis, pentose and glucuronate interposition, and pentose phosphate pathway enrichment scores being significantly positively correlated with disulfidptosis negative score. Conversely, the TCA cycle was significantly negatively correlated with the positive score (Figure 1I), implying that increased glucose metabolism activity is unfavorable for disulfidptosis occurrence.

The GSE116256 cohort included bone marrow samples from 16 patients with AML and 5 healthy participants, encompassing 21 cell types (Figure 2A). Among these, six types of AML malignant cells were identified as HSC-like, progenitor-like, granulocyte-monocyte-progenitor (GMP)-like, promonocyte-like, monocyte-like, and conventional dendritic cell-like. Expression profile analysis revealed that disulfidptosis suppressor genes like LRPPRC and NDUFA11 exhibited higher detectable expression rates in single cells, particularly in GMP-like cells (Figure 2B). Furthermore, the calculation of disulfidptosis scores in various cells showed that AML malignant cells, especially GMP-like cells, had higher negative scores (Figures 2C–F). Additionally, malignant cells demonstrated higher positive scores, resulting in no significant difference in activity scores between them and normal cells (Figure 2G).

Correlation analysis of immune infiltration showed that a high disulfidptosis activity score was associated with reduced infiltration of M2 macrophages and memory B cells, and increased infiltration of neutrophils. Conversely, high disulfidptosis negative scores were associated with increased infiltration of memory B cells, M1 macrophages, mast cells, and resting CD4⁺ T cells, along with reduced infiltration of monocytes and activated CD4⁺ T cells (Figure 2H). Analysis of TME characteristics demonstrated that higher disulfidptosis activity scores were associated with higher immune and stromal scores and lower tumor purity, whereas the opposite trend was observed for negative scores (Figure 2I). This suggests the involvement of disulfidptosis in antitumor responses within the AML TME. Moreover, the disulfidptosis activity score was significantly positively correlated with the expression of most immune checkpoints, with the opposite seen for the negative score (Figure 2J). Pathway analysis revealed that the activity score was positively correlated with the enrichment score of most tumor marker pathways (Figure 2K). Among these pathways, proliferation-related signaling pathways like MYC targets V1/V2 and DNA repair were negatively correlated with the activity score, indicating a potential inhibitory effect of disulfidptosis on cell viability. Correlation analysis with clinicopathological factors showed that white blood cell (WBC) count exhibited a negative correlation with both the positive score and activity score, while age showed a positive correlation with the activity score and a negative correlation with the negative score (Figure 3A). In the comparison of categorical variables, the activity score and positive score increased with an increase in cytogenetic risk, with the highest values observed in the French-American-British (FAB) classification of M6 and M7 (Figures 3B, C). Moreover, disulfidptosis scores did not show significant differences among patients (Figure 3D).

The weighted correlation network analysis (WGCNA) analysis was performed on The Cancer Genome Atlas—Acute Myeloid Leukemia (TCGA-LAML) dataset to identify additional potential DRGs. The cluster dendrogram displayed an increase in color depth corresponding to the magnitude of the disulfidptosis score (Figure 3E). Figures 3F, G show the scale-free fit exponent and average connectivity analysis for various soft threshold powers. The blue and purple module feature genes were combined with a cut height of 0.2 (Figures 3H, I). A soft threshold power of β = 5 (unscaled R2 = 0.9) was selected to categorize the top 5,000 genes, sorted by standard deviation, into 14 independent co-expression modules (Figure 3J). The correlation plot of the module-trait relationship indicated that the magenta gene module, comprising 152 genes, exhibited the highest correlation with the disulfidptosis activity score (Figures 3J, K; Supplementary Table S2). KEGG analysis highlighted that these genes were mainly enriched in signaling pathways such as metabolic pathways, hematopoietic cell lineage, platelet (PLT) activation, and proteoglycans in cancer (Figure 3L).

We performed a consensus cluster analysis based on the expression of DRGs retrieved in the study of Gan et al. According to the distribution characteristics of the cluster plots and considering the small size of samples in the TCGA-LAML cohort (Supplementary Figure S2). We chose cluster number 2 as optimal. Two disulfidptosis-related molecular subtypes, Cluster A and Cluster B, were thus identified (Figure 4A). Cluster B demonstrated a significantly worse prognosis than Cluster A (Figure 4B). The expression of RPN1, GYS1, SLC3A2, and NDUFA11 was significantly higher in Cluster B than in Cluster A (Figure 4C). Additionally, Cluster A exhibited a higher disulfidptosis negative score, while Cluster B displayed a higher activity score (Figure 4D). TME analysis indicated that Cluster B possessed higher stromal and immune scores (Figure 4E), primarily due to a greater proportion of monocyte infiltration. In contrast, Cluster A was enriched in more naive B cells, resting memory CD4⁺ T cells, and eosinophils (Figure 4F). Higher expression of PD-1, TNFRSF9, and CD86 was observed in Cluster B, suggesting the potential presence of immunosuppression in this subtype (Figure 4G). Difference analysis results showed that the activity of 50 tumor marker pathways in Cluster B exceeded that in Cluster A (Figure 4H).

Univariate Cox regression analysis on potential DRGs identified by WGCNA analysis revealed significant correlations between AML prognosis and GCLM, PLEKHH3, NEO1, CSF1, ST6GALNAC4, AK1, and SLC14A1 (Figure 5A). To reduce dimensionality and construct a prognostic risk score model, LASSO regression analysis was performed using six genes, excluding PLEKHH3 (Figures 5B, C) (Supplementary Table S3). Patients with AML were stratified into high- and low-risk score groups based on the optimal cut-off value (Figure 5D). Compared to the low-risk score group, the high-risk score group exhibited a higher number of deceased patients (Figure 5E), higher expressions of GCLM, NEO1, CSF1, ST6GALNAC4, and AK1, and lower expression of SLC14A1 (Figure 5F). Survival analysis demonstrated that the high-risk score group had a significantly worse prognosis than the low-risk score group (Figure 5G). ROC curve analysis revealed high area under the curve (AUC) values for the risk score at 1, 3, and 5 years (0.779, 0.714, and 0.778, respectively), indicating the robust prognostic value of the risk score (Figure 5H). Univariate and multivariate Cox regression analyses confirmed the risk score as an independent factor for predicting AML prognosis (p < 0.001) (Figures 5I, J).

The prognostic prediction ability of the risk score model was further validated across seven AML cohorts, with the high-risk score group consistently exhibiting a significantly worse prognosis than the low-risk score group (Figures 6A–G). Univariate Cox regression analysis confirmed the prognostic power of risk score model (p < 0.05) (Figure 6H). A clinical nomogram was constructed by combining clinicopathological factors significantly associated with AML prognosis, namely, age and cytogenetic risk (Figure 6I). Calibration curve analysis demonstrated the consistency between observed overall survival (OS) and predicted OS (Figure 6J). ROC curve analysis showed high AUC values for the nomogram at 1, 3, and 5 years (0.784, 0.769, and 0.871, respectively), confirming its strong prognostic value (Figure 6K).

Further exploration was conducted to explore potential reasons for the significant prognostic differences between the high- and low-risk score groups. The high-risk score group had a higher count of monocytes, while the low-risk score group had more memory B cells and resting and activated mast cells (Figure 7A). Notably, immune checkpoint expression, including PD-L1, CTLA1, LAG3, PD-1, PD-L2, CD80, CD86, and TNFRSF9, was significantly upregulated in the high-risk score group (Figure 7B). GSEA enrichment analysis revealed heightened activity in immune-related signaling pathways, such as antigen processing and presentation, B cell receptor signaling, chemokine signaling, cytokine-cytokine receptor interaction, and neutrophil extracellular trap formation in the high-risk score group (Figure 7C). Conversely, metabolism-related pathways, including ascorbate and aldarate metabolism, glycosaminoglycan biosynthesis-heparan sulfate/heparin, hedgehog signaling, lipoic acid metabolism, and pentose and glucuronate interposition, were enriched in the low-risk score group (Figure 7D). The high-risk score group displayed a higher gene mutation rate, with DNMT3A, NPM1, FLT3, TP53, and RUNX1 being the most frequently mutated genes (Figures 7E, F). Risk scores also significantly differed across various clinical characteristics, with patients of advanced age, male gender, and worse cytogenetic risk demonstrating higher risk scores (Figures 7G–I). In terms of the French-American-British (FAB) classification, patients with the M3 type had the lowest risk scores, while patients with the M5-M7 type exhibited higher risk scores (Figure 7J). Notably, risk scores did not differ significantly across white blood cells (WBC) and platelet (PLT) count groups (Figures 7K, L).

Sensitivity to commonly used AML drugs was predicted, with the low-risk score group exhibiting lower IC50 values for cytarabine and midostaurin, indicating greater sensitivity to these drugs. No significant difference in sensitivity to doxorubicin was observed (Figures 8A–C). Predictive results for immunotherapy showed that the high-risk score group was more responsive to anti-PD-1 treatment, with significantly upregulated PD-1 expression (Figure 8D). We further analyzed the differences in risk scores among different patients in a set of AML single-cell sequencing datasets after anti-PD-1 treatment. Figures 8E–G shows AML cell distribution, treatment response versus non-response population, and expression levels of risk scores in patients with AML, respectively. Analysis of the risk scores revealed that the single-cell risk score and average risk score were higher in patients with AML who did not respond to anti-PD-1 therapy compared to those who responded (Figures 8H, I). In the TCGA-LAML cohort, patients with both high PD-1 expression and high-risk scores exhibited the worst prognosis, while those with low levels of both had the best prognosis (Figure 8J). These results suggest that the risk score can predict the sensitivity of patients to chemotherapy and immunotherapy.

In a real-world clinical cohort, transcriptome sequencing was performed on samples from five normal individuals, five myeloid leukemia chronic-phase patients, and five myeloid leukemia acute-phase patients. Compared to the normal samples, the expression of CSF1 was upregulated in the chronic and acute phase samples, while AK1, NEO1, and SLC14A1 were downregulated. The expression of GCLM was upregulated in the chronic phase samples, while ST6GALNAC4 showed no significant change (Figure 9A). Risk scores increased with the progression of myeloid leukemia (Figure 9B). Correlation analysis with clinicopathological factors indicated positive correlations between the risk score and WBC, PLT, hemoglobin (HB), and age, while red blood cell (RBC) count showed a negative correlation (Figure 9C). Although the small sample size did not yield significant correlations (p < 0.05), disulfidptosis negative scores increased and activity scores decreased with the progression of myeloid leukemia (Figure 9D), confirming the presence of disulfidptosis resistance. Notably, there were no significant differences in disulfidptosis scores between gender groups (Figure 9E). Positive and active scores were positively correlated with WBC, RBC, PLT, HB, and age, while negative scores exhibited negative correlations with WBC, PLT, and HB (Figure 9F).