The overall study workflow is presented in Figure 1. Gene expression profiles of MM patients were retrieved from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) for the datasets GSE2658 and GSE57317, while normal sample data were obtained from GSE118985. All three datasets were generated using the Affymetrix Human Genome U133 Plus 2.0 Array (GPL570) platform. Raw microarray data were merged and batch-corrected using the “Combat” algorithm in the SVA package, yielding a unified dataset comprising 559 MM samples (GSE2658), 55 MM samples (GSE57317), and 68 normal samples (GSE118985) (see Additional file 1). To further assess the prognostic value of the model, dataset GSE136337 from GEO was used as an external validation cohort. DRGs were extracted based on the gene list reported by the Li research group (6).

Interactions among DRGs were analyzed using the R packages “iGraph”, “Psych”, “Reshape2”, and “RColorBrewer”. Correlations among DRGs were examined with the “corrplot” package. Additionally, the “RCircos” package was used for localizing the position of genes on chromosomes. To explore protein-protein interactions (PPI), DRGs were submitted to the STRING database (version 12.0) (https://cn.string-db.org/), and the resulting network was analyzed (7). Somatic alteration frequencies and chromosomal mutation sites in MM were extracted from the cBioPortal for Cancer Genomics (http://www.cbioportal.org/). Furthermore, differential expression between MM and normal samples was assessed using the “limma” R package.

This study evaluated a total of 24 DRGs: ACTB, ACTN4, CAPZB, CD2AP, DSTN, FLNA, FLNB, GYS1, INF2, IQGAP1, LRPPRC, MYH10, MYH9, MYL6, NCKAP1, NDUFA11, NDUFS1, NUBPL, PDLIM1, RPN1, SLC3A2, SLC7A11, TLN1, and OXSM. The rationale for including each gene, encompassing its function, directionality, and relevance to disulfidptosis, is summarized in Additional File 2. Data were collected and analyzed using the “ConsensusClusterPlus” R package, which stratified the 614 tumor samples into two clusters (Cluster A and Cluster B). Survival differences between the clusters were assessed using the “survival” R package. Additionally, differences in clinical features between clusters were evaluated with the “pheatmap” R package.

The “GSVA” and “GSEABase” R packages were used to investigate the underlying biological processes and signaling pathways. The datasets “c5.go.v2024.1.Hs.symbols.gmt” and “c2.cp.kegg_legacy.v2024.1.Hs.symbols.gmt” were obtained from the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb/) as collections of annotated gene sets. A false discovery rate (FDR) < 0.05 was applied as the statistical threshold for clustering. Furthermore, single-sample gene set enrichment analysis (ssGSEA) was performed to quantify differences among the clustering classes, enabling a detailed comparison of their biological characteristics.

To quantitatively assess the relationship between MM and disulfidptosis, a disulfidptosis-related prognostic model was developed. The “limma” R package was employed to identify two distinct clusters, with filtering criteria set at FDR < 0.05 and absolute log₂ fold change > 1. The complete set of MM samples (n = 614) was randomly divided into a training set (n = 307) and a validation set (n = 307) at a 1:1 ratio. DEGs between these subtypes were subjected to univariate Cox proportional hazards regression, with p < 0.05 considered statistically significant for prognostic relevance. LASSO-Cox regression was subsequently performed using the “glmnet” R package to construct a prognostic risk model associated with cell death (8). The optimal penalty parameter (λ) was determined via 10-fold cross-validation; λ.min, corresponding to the minimum cross-validation error, was selected for model construction, while λ.1se was also recorded as a reference for model robustness. Risk scores for disulfidptosis were calculated based on the expression profiles of key genes using the derived model equation. The training set was used to establish the prognostic model, and the validation set was employed to assess its predictive performance. The risk score was calculated according to the following formula:

In the risk score formula, exp represents gene expression, and coef denotes gene risk coefficients. Samples were stratified into high-risk and low-risk groups based on the median risk score, following standard risk classification criteria. Kaplan-Meier (K-M) survival analysis was performed to evaluate the prognostic significance of the risk model, and receiver operating characteristic (ROC) curves were generated to assess its predictive accuracy. These analyses were conducted using the R packages “survival”, “survminer”, “pheatmap”, and “timeROC”. Principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) were also performed using the “Rtsne” and “ggplot2” packages to visualize differences between the two risk groups. Missing values were addressed via complete-case analysis, and potential outliers were examined using standardized residuals and leverage statistics, with exclusions applied only when justified. The proportional hazards assumption of the Cox models was tested using Schoenfeld residuals, and no significant violations were observed. Additionally, pathway enrichment scores for key biological processes were calculated using GSVA, and their correlations with the disulfidptosis risk score were assessed by Spearman correlation.

A nomogram was constructed based on MM patient data to predict disease risk using the calculated risk scores. The results were visualized with the R packages “rms” and “rmda”. To further evaluate the clinical utility and robustness of the proposed model, its prognostic performance was compared with that of previously published MM risk models using the “survival” and “survcomp” packages. Decision curve analysis (DCA) was performed with the “rmda” package, and clinical impact curves (CICs) were generated to depict the number of patients classified as high risk and the corresponding true positives across different threshold probabilities.

To analyze differences in the TME between the high-risk and low-risk groups associated with cell death, the “estimate” R package was used to calculate the immune score, stromal score, and estimate score. Gene set enrichment analysis (GSEA) was performed using the “clusterprofiler” R package to explore the underlying biological pathways and mechanisms. The CIBERSORT algorithm (https://cibersortx.stanford.edu/) was implemented in R to evaluate 22 immune cell types per TME sample (9), alongside an analysis of their relationships in immune cell infiltration. MCP-counter and xCell were additionally employed to provide orthogonal assessments of immune and stromal cell populations, thereby enhancing robustness and capturing bone marrow-specific features. Spearman correlation analysis was conducted to assess the relationship between the risk model gene signatures and immune infiltration levels. To investigate potential factors contributing to tumor immune evasion, the tumor immune dysfunction and exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu/) was applied. A heatmap was generated to visualize correlations between immune checkpoint genes and risk scores. Immune checkpoint genes significantly associated with risk scores (p < 0.05) were further subjected to FDR correction and multivariable Cox regression to evaluate their independent prognostic significance, with hazard ratios and 95% confidence intervals (CIs) reported and visualized in forest plots. K-M survival analysis was also performed using the “survival” and “survminer” R packages to assess their prognostic relevance.

Expression and sensitivity data for targeted drugs were obtained from the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/). To identify statistically significant differences in drug response, a threshold of P < 0.001 was applied. The “oncopredict” and “parallel” R packages were used to assess differences in drug sensitivity between the two risk groups (10). This analysis allowed prediction of individualized treatment responses and highlighted potential therapeutic agents for risk-stratified MM patients.

The MM cell lines RPMI8226 and U266 were obtained from Baidi Biotechnology Co., Ltd. (Zhejiang, China). Cells were cultured in RPMI 1640 medium (GIBCO, USA) supplemented with 20% fetal bovine serum (FBS) and 1% penicillin/streptomycin at 37°C in a 5% CO₂ atmosphere. Bone marrow aspirates were collected from 13 MM patients and 4 healthy controls at Weifang People’s Hospital. Normal plasma cells and patient-derived myeloma plasma cells were isolated from bone marrow mononuclear cells using CD138⁺ magnetic bead selection. The clinical characteristics of the MM patients are summarized in Additional File 3. This study was approved by the Ethics Committee of Shandong Second Medical University, and informed consent was obtained from all participants. Total RNA was extracted using TRIzol reagent (Sparkjade, China) following the manufacturer’s protocol. Gene expression was quantified by qRT-PCR using the SYBR Green method on an Applied Biosystems 7500 Fast Real-Time PCR system. Relative expression levels were calculated using the 2^-△△Ct method, with GAPDH serving as the housekeeping gene for normalization. The primers used for qRT-PCR are listed in Table 1.

WB analysis was subsequently performed. Equal amounts of cell lysates were loaded onto SDS-PAGE gels for electrophoresis, and proteins were then transferred onto polyvinylidene fluoride (PVDF) membranes. Membranes were blocked with 5% skim milk in Tris-buffered saline containing 0.1% Tween-20 (TBST) and incubated overnight with primary antibodies. Following this, membranes were incubated with HRP-conjugated goat anti-rabbit secondary antibody for 1 h. Specific protein bands were visualized using an enhanced chemiluminescence reagent and an Amersham Imager 600 (GE Healthcare Biosciences, Pittsburgh, PA, USA). Band intensities were quantified using ImageJ software (v.1.50a) (https://imagej.net/software/imagej/).

All data analyses and visualizations were performed via R (version 4.4.1) and GraphPad Prism 9.5. For quantitative variables, Student’s t test was used to evaluate differences between groups with normally distributed data, whereas the Wilcoxon test was used for skewed data, with a P value of less than 0.05 considered statistically significant. In this context, the levels of significance are represented as follows: * P < 0.05, ** P < 0.01, *** P < 0.001, and **** P < 0.0001. All genome-wide and multi-gene P values were adjusted by the FDR method, with FDR < 0.05 considered significant.

Chromosomal mapping of the 24 DRGs revealed their genomic distribution (Figure 2A). Regulatory network analysis demonstrated their interactions, suggesting potential prognostic relevance in MM (Figure 2B). Correlation analysis identified strong positive associations (e.g., ACTB-MYH9, r = 0.41) and negative associations (e.g., ACTN4-CD2AP, r = -0.26) among DRGs (Figures 2C, D). Genomic profiling detected missense and splice-site mutations in MM samples (Figure 2E). Differential expression analysis of GEO datasets revealed significant downregulation of ACTN4, CAPZB, CD2AP, DSTN, GYS1, LRPPRC, MYL6, PDLIM1, RPN1, and SLC3A2 in MM tissues compared with controls, whereas ACTB, FLNA, FLNB, INF2, IQGAP1, MYH9, NCKAP1, NUBPL, and SLC7A11 were upregulated (Figure 2F), implicating DRGs in MM pathogenesis.

Consensus clustering of 614 MM samples (GSE2658 and GSE57317) identified two robust subtypes (k = 2), which were validated by cumulative distribution function (CDF) curves (Figures 3A–C). Cluster A (n = 306) and Cluster B (n = 308) exhibited clear separation in principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) plots (Figures 3D, E). Survival analysis revealed significantly longer overall survival (OS) in Cluster A compared with Cluster B (P < 0.001; Figure 3F). Cluster A was characterized by upregulation of CD2AP, DSTN, FLNB, NCKAP1, NDUFS1, and NUBPL, whereas Cluster B overexpressed MYL6, NDUFA11, PDLIM1, SLC3A2, and TLN1 (Figures 3G, H).

ssGSEA indicated higher immune infiltration in Cluster B, with CD56dim NK cells and monocytes significantly enriched (P < 0.001; Figure 4A). Immune cell correlation networks highlighted interactions with cell death-related genes (Figure 4B). Gene set variation analysis (GSVA) revealed that Cluster B was enriched in pathways related to prion diseases, viral myocarditis, DNA repair, and replication, whereas Cluster A exhibited enrichment in guanylyltransferase activity and histone ubiquitination pathways (Figures 4C, D).

From 690 initial DEGs, univariate Cox regression identified 121 prognostic genes. DEGs were determined exclusively within the training set, and a Lasso-Cox model with 10-fold cross-validation was applied to select the optimal penalty parameter (λ.min ≈ 0.051), narrowing the candidates to 29 genes. Subsequent multivariate regression yielded a 9-gene prognostic signature: Risk score = (0.4500 × TPST2) + (0.1514 × HIF1A) + (0.1571 × KIF21B) + (-0.6397 × MCPH1) + (-0.5123 × MAST4) + (0.3354 × ANXA2) + (-0.4341 × ALG14) + (-0.2931 × PQLC3) + (0.4405 × RANGAP1). The results of the multivariate regression analysis are shown in Supplementary Figure 1 and Additional File 4. High-risk patients exhibited significantly worse OS in the combined cohort (P < 0.001), training set (P < 0.001), and validation set (P < 0.05; Figures 5C–H). ROC analysis confirmed the prognostic accuracy of the model, with areas under the curve (AUCs) of 0.674, 0.786, and 0.727 at 1, 3, and 5 years, respectively (Figures 5I–K). Validation in the external dataset GSE136337 yielded AUCs of 0.646, 0.705, and 0.735 at 1, 3, and 5 years, respectively (Supplementary Figure 2), demonstrating high predictive performance. Furthermore, PCA and t-SNE analyses revealed clear separation between the high- and low-risk groups, confirming distinct risk stratification (Figures 5L–O).

To elucidate the biological implications of the disulfidptosis-related signature, we examined its correlation with key pathways involved in MM pathogenesis and disulfidptosis mechanisms. Spearman correlation analysis was performed between the risk score and pathway enrichment scores representing proteasome activity, endoplasmic reticulum (ER) stress, oxidative phosphorylation, and actin cytoskeleton organization. Notably, the risk score was significantly positively correlated with proteasome activity (r = 0.34, P < 0.001) and oxidative phosphorylation (r = 0.28, P < 0.001) (Supplementary Figure 3). Conversely, a significant negative correlation was observed with the ER stress response (r = -0.16, P < 0.001), and a weak but statistically significant inverse correlation was found with actin filament organization (r = -0.09, P < 0.05). These findings suggest that high-risk patients, as defined by our model, exhibit enhanced proteasomal function and metabolic activity, accompanied by a dampened ER stress response and subtle alterations in actin cytoskeletal dynamics.

To improve prognostic precision in MM patients, we constructed a nomogram incorporating the disulfidptosis-related risk score (Figure 6A). Bootstrap validation with 1,000 resamples demonstrated excellent concordance between predicted and observed outcomes, as indicated by the calibration curve (average absolute error = 0.022), confirming the model’s robustness (Figure 6B). Comparative analysis revealed significantly lower risk scores in Cluster B compared with Cluster A (P < 0.001), highlighting the model’s discriminative ability (Figure 6C). Differential gene expression patterns between the high- and low-risk groups further validated the biological relevance of the cell death-related model (Figure 6D). Decision curve analysis (DCA) indicated that the nomogram provided superior net benefit across a range of threshold probabilities, and clinical impact curves (CICs) confirmed that predicted high-risk cases closely matched actual outcomes, supporting its potential clinical utility (Supplementary Figure 4).

The Sankey diagram illustrated consistency in patient stratification across classification methods, with Cluster A predominantly comprising high-risk patients and Cluster B enriched for low-risk cases (Figure 6E). In benchmark comparisons with four established MM prognostic models (11–14), our disulfidptosis-based model demonstrated superior performance, exhibiting significant survival discrimination across all comparative models (P < 0.05; Supplementary Figures 5A–D) and higher predictive accuracy (C-index: 0.709 vs. Sun: 0.681, Zhu: 0.647, Zhang: 0.631, Wang: 0.585; P < 0.05; Figures 6F, G; Supplementary Figures 5E–G). Collectively, these results establish the disulfidptosis signature as a clinically actionable prognostic tool with enhanced stratification power relative to existing models.

We applied the Estimation of Stromal and Immune cells in Malignant Tumor Tissues Using Expression Data (ESTIMATE) algorithm to assess TME characteristics, including stromal, immune, and combined estimate scores. Comparative analysis revealed significantly higher immune scores in high-risk patients (P < 0.001), suggesting increased immune cell infiltration (Figure 7A).

To further investigate functional differences between high- and low-risk groups, GSEA was performed to identify pathways associated with risk stratification. High-risk samples were significantly enriched in pathways related to the cell cycle, DNA replication, progesterone-mediated oocyte maturation, proteasome activity, and pyrimidine metabolism. In contrast, low-risk samples exhibited enrichment in pathways associated with ascorbate and aldarate metabolism, lysosome function, neuroactive ligand-receptor interaction, retinol metabolism, and vascular smooth muscle contraction (P < 0.05; Figures 7B, C).

To characterize the immune landscape of the TME, we applied the CIBERSORT algorithm to deconvolute 22 immune cell subsets across all samples. This analysis consistently revealed enrichment of memory B cells and plasma cells throughout the cohort, with high-risk patients exhibiting significantly higher infiltration of memory B cells, CD8⁺ T cells, and activated dendritic cells compared with low-risk patients (P < 0.001; Figure 7D; Supplementary Figure 6). Correlation analysis identified the strongest negative associations between risk scores and both memory B cells and activated dendritic cells, whereas naive B cells and plasma cells showed the most pronounced positive associations; Spearman correlation further confirmed a significant relationship between the gene signature and infiltrating immune cells (Figure 7E). Sensitivity analyses using MCP-counter and xCell revealed consistent patterns: high-risk samples exhibited increased tumor/plasma-cell signals alongside reduced B-lineage cells, CD8⁺ T cells, cytotoxic lymphocytes, and NK cells. MCP-counter additionally highlighted alterations in neutrophils, monocytic lineages, endothelial cells, and fibroblasts. Notably, plasma-cell signatures displayed an inverse relationship with normal B-cell signals, consistent with the expectation that malignant plasma-cell expansion suppresses normal B-cell states (Supplementary Figure 7). Collectively, these findings show that the tumor microenvironment in high-risk patients is characterized by an abundance of malignant plasma cells and a scarcity of key effector immune cells, thereby fostering an immunosuppressive milieu associated with a poor prognosis.

TIDE analysis revealed significantly higher scores in high-risk patients (P < 0.001), suggesting increased potential for immune evasion and reduced responsiveness to immunotherapy (Figure 7F). Hotspot analysis further demonstrated correlations between the disulfidptosis-related gene signature and immune checkpoint genes (Figure 7G). K-M survival analysis indicated that higher expression of checkpoint genes such as CD70 and LAIR1 was associated with improved survival within the high-risk group, whereas elevated expression of CD40, CD200R1, CD274, PDCD1LG2, TNFRSF4, and TNFRSF14 was correlated with significantly better outcomes in intermediate- and high-risk patients compared with low-risk individuals (P < 0.05; Figure 7H). Multivariable Cox regression analysis of all immune checkpoint genes is shown in Supplementary Figure 8. Notably, CD274 exhibited a hazard ratio (HR) of 0.73 (95% CI, 0.50-1.06), indicating a non-significant trend toward reduced mortality with higher expression, consistent with the K-M results but not reaching statistical significance. This observation may reflect biological heterogeneity, as CD274 can mark an inflamed, T-cell-rich microenvironment associated with better prognosis, measurement variability, or residual confounding factors such as tumor burden or prior treatment, and thus should be interpreted cautiously.

To evaluate the clinical applicability of our risk model, we assessed correlations between disulfidptosis-related risk scores and therapeutic responses using the oncoPredict algorithm for 198 FDA-approved and experimental antitumor agents. To ensure reliable drug sensitivity predictions, expression data were subjected to strict preprocessing and quality control. PCA revealed no obvious batch effects, and expression distributions were highly consistent across samples (Supplementary Figure 9), confirming dataset comparability and suitability for downstream oncoPredict modeling. Significant differences in predicted drug sensitivity were observed between risk groups (P < 0.05). Specifically, low-risk patients exhibited increased sensitivity to WIKI4, WEHI-539, AZD6738, Linsitinib, AZD7762, and I-BRD9, whereas high-risk patients were more responsive to Doramapimod (BIRB 796), SB216763, NU7441, SCH772984, JQ1, and IAP_5620 (Figure 8). These findings highlight the potential of the risk model for precision oncology applications, with risk stratification predicting distinct therapeutic vulnerabilities. The observed differential sensitivity profiles suggest mechanistic differences in drug response pathways between disulfidptosis-defined subgroups.

To validate the clinical relevance of our prognostic signature, we performed systematic molecular characterization of the nine identified key genes. Previous studies have reported upregulation of ANXA2 and HIF1A in MM, whereas MAST4 exhibits tumor-suppressive downregulation. We therefore focused experimental validation on the remaining six genes (MCPH1, PQLC3, TPST2, ALG14, RANGAP1, and KIF21B), whose roles in MM pathogenesis remain uncharacterized. qRT-PCR analysis of MM cell lines (RPMI8226 and U266) compared with normal plasma cells revealed significant downregulation of ALG14, MCPH1, and PQLC3, and marked upregulation of TPST2, RANGAP1, and KIF21B (Figure 9A). These expression changes were confirmed in paired primary MM samples (n = 13) by paired t-test, with all comparisons reaching statistical significance (P < 0.05; Figure 9B). Considering MCPH1’s established tumor-suppressor function in other malignancies, we performed WB analysis, which demonstrated reduced MCPH1 expression in MM cell lines relative to controls (Figure 9C; uncropped images are shown in Supplementary Figure 10), suggesting that loss of MCPH1 may contribute to MM pathogenesis. This multi-level validation confirms the biological relevance of our prognostic signature and identifies MCPH1 as a potential novel tumor suppressor in MM.