Gene expression profiles were downloaded from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/gds/). GSE24080, GSE4581 and GSE2658 were obtained using the GPL570 platform (Hanamura et al., 2006; Shi et al., 2010). GSE6477 was obtained using the GPL96 platform (Chng et al., 2007). GS136337 was obtained using the GPL27143 platform (Danziger et al., 2020). Multiple Myeloma Research Foundation (MMRF) CoMMpass study offers the expression profiles of MM patients and clinical information (including survival statistics). The normalized was conducted by the “limma” R package.

Thirteen cuproptosis-related genes were obtained from previous studies (Supplementary Table S1) (Tsvetkov et al., 2022) (Oliveri, 2022; Tang et al., 2022). 1,399 oxidative stress-related genes were extracted from GeneCards (https://www.genecards.org) with a relevance score ≥ 7 (Supplementary Table S2) (Wu et al., 2021).

The flowchart for this study was shown in Figure 1.

Using Pearson correlation analysis, cuproptosis -related oxidative stress genes were identified and a co-expression network was created based on the cutoff point (Pearson |R| > 0.3 and p < 0.05). A total of 419 cuproptosis-related oxidative stress genes were identified.

With p < 0.05 and |log2 Fold Change|≥ 0.6 as cut-off points, differentially expressed genes (DEGs) between MM patients and normal donors were identified using the “limma” package in R language. To illustrate DEGs, heat maps and volcano were created using the R package “ggplot2” and “pheatmap”. We then took the intersection of DEGs and cuproptosis-related oxidative stress genes to obtain differentially expressed cuproptosis -related oxidative stress genes (DECROGs).

For the DECROGs, the gene ontology (GO) annotation and Reactome pathway analysis were applied to investigated their functions and pathway enrichment utilizing the R package clusterProfiler (Yu et al., 2012). The cutoff criterion was the false-discovery rate (FDR). Terms with FDR <0.05 represented significant enrichment. GO enrichment consists of three components: biological process (BP), molecular function (MF) and cellular component (CC).

We identified prognosis-related genes by the univariate Cox regression analysis (p < 0.05). These genes were entered into the Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis to screen key genes (p < 0.05) through the R package glmnet (Friedman et al., 2010). Then, multivariate Cox regression was used to analyzed these genes and establish the risk score model. Risk scores of all samples were calculated according to the following formula: risk score = ∑Coef* Exp. The “Coef” refers to the coefficient of each mRNA from the multivariate Cox analysis and “Exp” indicates the expression level of each mRNA. We validated the risk model using the GSE2658, GSE136337, GSE4581 and MMRF datasets. The median risk score was used as the cut-off value to divide the high-risk and low-risk groups, The difference in overall survival between the high- and low-risk groups was evaluated using the Kaplan-Meier survival analysis and log-rank test. The R package “survivalROC” was utilised to produce the time-dependent Receiver Operating Characteristic (ROC) curves.

We constructed a nomogram prognostic prediction model based on risk scores and clinical pathological feature using the “rms” R package. Concordance index (C-index), calibration curve, and Decision Curve Analysis (DCA) were used to assess the prognostic performance of the established nomogram.

Immune infiltration analysis was analyzed by single-sample gene set enrichment analysis (ssGSEA) method. This involved the enrichment scores of 16 immune cell and 13 immune function categories. The “GSVA” package was utilized for the analysis (Hänzelmann et al., 2013). The correlation between the immune cell scores and the signature was investigated by Pearson correlation analysis, and the “corrplot” package was used for the analysis.

The data on drug sensitivity was obtained from the Genomics of Drug Sensitivity in Cancer (GDSC) database (Yang et al., 2013). “OncoPredict” package was used to calculate the 50% inhibiting concentration (IC50) values in the high- and low-risk groups (Maeser et al., 2021).

R software (version 4.2.0) and corresponding packages were carried out for statistical analyses. The Kaplan-Meier method and log-rank test were used to assess prognosis. Pearson method was applied for correlation analysis. The Wilcoxon test was used to compare data between two groups, while the Kruskal-Wallis H test was used for data comparison among three groups. All tests with p-value <0.05 indicated statistical significance.

The correlation between cuproptosis-related mRNAs and oxidative stress genes was analyzed by the Pearson correlation coefficient method (Pearson |R| > 0.3 and p < 0.05). A total of 419 cuproptosis-related oxidative stress genes were identified. The sankey plot demonstrated the correlation between cuproptosis-related genes (CRGs) and oxidative stress genes (OSGs) (Figure 2A). Then, we set a p < 0.05, and [log2FoldChange (log2FC)] > 0.6 as the cutoff point. Based on this criteria, we identified 2,298 DEGs in the GSE6477 dataset comparing multiple myeloma with healthy donors. The volcano plot of DEGs was illustrated in Figure 2B. Finally, 2,298 DEGs were intersected with 419 cuproptosis-related oxidative stress genes and thus we obtained 76 differentially expressed cuproptosis-related oxidative stress genes (DECROGs). The Venn diagram was shown in Figure 2C and the heatmap of the 76 DECROGs are shown in Figure 2D.

GO enrichment and Reactome analyses were conducted for exploring the potential function of 76 DECROGs. The result was demonstrated in Figures 2E, F. The results of GO enrichment analysis were categorized into the following three parts. For biological process (BP), DECROSGs were enriched in response to oxidative stress, and lipid localization. For cellular component (CC), the genes were associated with lumen and membrane. For molecular function (MF), the genes were associated with antioxidant activity and lipid-protein binding. Reactome pathway analysis demonstrated that the DECROGs were mainly associated with cell respiratiory electron transport, interleukin-12 signaling and neutrophil degranulation, defective intrinsic pathway for apoptosis.

Then, we investigated the prognostic significance of the 76 genes using univariate Cox regression. Consequently, a total of 26 DECROGs associated with prognosis were chosen (p < 0.05) (Table 1). Subsequently, the 26 DECROGs described above were integrated into the LASSO regression model (Figures 3A, B). Then stepwise multivariate Cox proportional hazard regression analysis was used and eight genes were identified to construct a prognostic model for MM patients (Figure 2C). We obtained eight genes, including RNASE3, APOE, CCNB1, MIF, FOXO1, KIT, PLA2G4A, ECG1 to build the risk model for MM patients. Among them CCNB1, MIF, PLA2G4A may be regarded as oncogenes, whereas RNASE3, APOE, FOXO1, KIT, and EGR1 may be tumor suppressor genes. The formula for calculating out each patient’s risk score is as follows: R i s k s c o r e = EXP R N A S E 3 * − 0.077 + EXP A P O E * − 0.214 + EXP C C N B 1 * 0.270 + EXP M I F * 0.290 + EXP F O X O 1 * − 0.238 + EXP K I T * − 0.069 + EXP P L A 2 G 4 A * 0.192 + EXP E G R 1 * − 0.104

According to the median risk scores of 558 samples, high-risk subgroups (n = 279) and low-risk subgroups (n = 279) were stratified. The heat map showed the expression levels of the eight DECROGs in the high- and low-risk groups (Figure 4A). Figure 4B shows the distribution of risk scores and survival status between the high-risk and low-risk groups. The mortality rate of patients increased with the increasing risk score.

Kaplan-Meier analysis showed that overall survival was worse in high-risk patients than in low-risk patients (p < 0.0001, Figure 4C). We evaluated the efficacy of the prognostic model by time-dependent ROC curves, and the area under the curves (AUC) was 0.653 for 1-year survival, 0.678 for 3-year survival, and 0.643 for 5-year survival (Figure 4D), indicating a moderate performance of this model.

The prognostic significances of this 8-gene signature were further vilified in 3 external independent GEO datasets and MMRF datasets with over 2,000 MM patients (Figure 5). The results showed that the overall survival of patients in the high-risk groups was significantly worse than that of patients in the low-risk group. The AUC of the risk model was >0.6, proving the performance of this model. In summary, we established an eight-gene risk model that exhibited acceptable performance in these datasets.

We evaluated the correlation between the risk score and other clinical characteristics using univariate and multivariate Cox regression analysis with the GSE24080 dataset to determine whether the risk score was an independent prognostic factor. The univariate analysis revealed that AGE, β2M, ALB, CRP, LDH, HGB, and risk score were substantially associated with the overall survival in MM patients (Figure 6A). Then, multivariate analysis confirmed that AGE, β2M, ALB, LDH, and risk score were significant independent prognostic factors (Figure 6B). We established a nomogram using the four independent factors to predict 1-,3- and 5- year overall survival rate in the GSE24080 dataset (Figure 6C). We evaluated the nomogram in terms of calibration, discrimination, and DCA curve. The calibration curve demonstrated that the predictive performance of the nomogram was generally consistent with the actual survival time in terms of the 1-, 3- and 5-year overall survival rate (Figure 6D). The C-index (concordance index) is mainly used to evaluate the discrimination between predictive outcomes of the nomogram and actual outcomes in survival analysis. The nomogram’s C-index was at 0.727 (95% CI 0.706–0.749), showing that the nomogram has good discrimination (Figure 6E). The five-year DCA curves showed that the nomogram was a good predictor of survival in MM patients (Figure 6F).

We further assessed the connection between the risk score and clinical features. The international staging system (ISS) is one of the earliest validated risk stratification for MM patients (Greipp et al., 2005). We first compared the predictive performance of this eight-gene risk model with ISS. In both GSE24080 and MMRF datasets, the time-dependent AUC of our model were higher than those of ISS, indicating that the predictive ability of our risk model was superior to ISS (Figures 7A, B). We then hypothesize that this risk model can further optimize the predictive performance of ISS for patient outcomes. We evaluated differences in overall survival between the high- and low-risk groups of three stages stratified by ISS. For the ISS stage I, the difference between the two groups was insignificant in both GSE24080 and MMRF datasets (Supplementary Figures S1A, S1B). For the ISS stage II, patients in the high-risk group had a worse prognosis than patients in the low-risk group in GSE24080. However, the difference between the two groups was not apparent in MMRF (Supplementary Figure S1C, S1D). As shown in Figures 7C, D, stage III patients were clearly split into two groups with varying survival rates, and the high-risk group had a worse prognosis. As shown in Figures 7E, F, the risk scores of GSE24080 and MMRF also increased significantly with increasing tumor stage (Kruskal–Wallis test p < 0.05, Figures 7E, F). When stage III patients were compared with stage I and II patients, respectively, there were significant differences in risk scores (Wilcoxon test p < 0.05). However, the difference in risk scores among stage I and stage II patients was not significant (Wilcoxon test p = 0.49, Wilcoxon test p = 0.089, respectively, Figures 7E, F). In conclusion, the eight-gene risk score model can be further used for stage III patients to predict the prognostic outcome more accurately.

Multiple myeloma (MM) is an exceptionally complicated and heterogeneous disease, which is characterized by various genetic alterations (Bustoros et al., 2022). Del (17p), del (13q) and translocation t (Palumbo and Anderson, 2011; Lipchick et al., 2016) were considered high-risk Chromosomal abnormalities by the IMWG (Palumbo et al., 2015). Amplification of 1q (amp1q) are associated with worse prognosis (D'Agostino et al., 2020). The prognostic significance of t (Schieber and Chandel, 2014; Lipchick et al., 2016) MM remains debatable (Bal et al., 2022). To confirm the predictive capacity of the 8-gene signature in patients with and without these genetic alterations, we selected the GSE136337 dataset containing these alterations for analysis. For genetic variation of amp1q, we only analyzed patients without amp1q, as there were only two patients with amp1q. We compared the survival of MM patients without these genetic indicators and found that patients in high-risk group had a worse prognosis (Figure 8). For patients with genetic indicators, including del (17p), del (13q), t (Palumbo and Anderson, 2011; Lipchick et al., 2016), and t (Schieber and Chandel, 2014; Lipchick et al., 2016), there was no difference between the high-risk group and low-risk group (Supplementary Figure S2). As a result, the MM patients without genetic abnormalities could be efficiently stratified by high- and low-risk using our eight-gene risk model.

To further investigate the relationship between risk score and immune infiltration, we used ssGSEA to evaluate the enrichment scores of different immune cells and immunological functions in the high-risk and low-risk groups. Compared to the high-risk group, most immune cell components were higher in the low-risk group (such as aDCs, macrophages, pDCs, Th1 cells, DCs, B cells, mast cells, Th2 cells, T helper cells, TIL, Treg), except for a lower proportion of CD8⁺ T cells. Moreover, except for MHC class I, the immune functional scores, such as APC co-inhibition, T cell co-inhibition, neutrophils, Type I IFN response, CCR, cytolytic activity, APC co-stimulation, HLA, T cell co-stimulation, inflammation-promoting, checkpoint, and parainflammation were lower in the high-risk group than in the low-risk group (Figures 9A, B). We further analyzed the correlation between the enrichment scores for immune cell and risk scores. The results showed a negative correlation between the risk score and the enrichment scores of aDCs, macrophages, pDCs, Th1 cells, DCs, B cells, mast cells, Th2 cells, T helper cells, TIL, and Treg, while the risk score and the enrichment score of CD8⁺ T cells were positively correlated (Figure 9C). These findings may suggest that the immune system of low-risk patients is more active.

To identify suitable drugs for high-risk patients, we compared the chemosensitivity of high- and low-risk groups based on IC50 values using the OncoPredict algorithm. A lower IC50 value indicates higher drug sensitivity. We found that patients in the high-risk group exhibited greater resistance to eight drugs, including SB216763 (a GSK3 inhibitor), Doramapimod (a p38 MAPK inhibitor), PF-4708671 (an S6K1 inhibitor), BMS-754807 (an IGF-1R/IR inhibitor), Selumetinib (a MEK inhibitor), NU7441 (a DNAPK inhibitor), Ribociclib (a CDK4/CDK6 inhibitor), and JQ1 (a BET inhibitor). On the other hand, high-risk patients showed greater sensitivity to eight drugs, including MIRA-1 (a TP53 inhibitor), GDC0810 (an ESR1/ESR2 inhibitor), Dihydrorotenone (a mitochondrial inhibitor), I-BRD9 (a BRD9 inhibitor), WIKI4 (a TNKS1/TNKS2 inhibitor), Fulvestrant (an ESR inhibitor), Linsitinib (an IGF1R inhibitor), and BI-2536 (a PLK1/PLK2/PLK3 inhibitor) (Figure 10). We believe that these drugs may be beneficial for the subsequent treatment of high-risk patients.