Comprehensive analysis of prognosis of cuproptosis-related oxidative stress genes in multiple myeloma

1 Introduction

Multiple myeloma, the second most common hematologic malignancy, accounts for 1.8% of all malignancies. It is a heterogeneous malignant plasma cell disorder characterized by aberrant proliferation of mature B cells (Mitsiades et al., 2007; Rajkumar et al., 2014). The first-line treatment for multiple myeloma (MM) mainly consists of the combination of bortezomib, and dexamethasone with either cyclophosphamide or doxorubicin. With the introduction of more effective proteasome inhibitors, such as carfilzomib and more potent immunomodulatory drugs such as pomalidomide, as well as the development of a new class of monoclonal antibody therapies, the clinical outcome of relapsed or refractory MM patients has improved significantly (Joshua et al., 2019). However, MM remains incurable with a high recurrence rate. The prognosis of MM patients is still poor (Palumbo and Anderson, 2011; Sonneveld et al., 2016). Therefore, performing risk stratification for patients and finding a new prognostic biomarker is crucial to improve prognostic accuracy and direct therapeutic treatment.

Copper ion plays a crucial role in numerous biological processes, including mitochondrial respiration, redox signaling, kinase signaling, cell wall remodeling, oxidative stress responses and other processes (Rae et al., 1999; Kim et al., 2008; Yruela, 2009). Dysregulation of copper plays a key role in cancer and mitochondria in many diseases. Cuproptosis is a recently identified form of programmed cell death (PCD) distinguished from known death mechanisms like apoptosis necroptosis, pyroptosis, and ferroptosis. Cuproptosis occurs through the binding of the intracelluar copper to lipoylated components of the tricarboxylic acid (TCA) cycle in mitochondria (Tsvetkov et al., 2022). This leads to aggregation of lipoylated protein and subsequent loss of Fe-S cluster-containing proteins, which results in acute proteotoxic stress and ultimately cell death. Although cuproptosis has received increased attention, the mechanism and role of cuprotosis in multiple myeloma remain unclear.

Oxidative stress occurs due to the excessive production of reactive oxygen species (ROS) or the failure of antioxidants to eliminate ROS inadequately. It is an essential factor in driving tumorigenesis and cancer progression (Schieber and Chandel, 2014; Lipchick et al., 2016). The production and degradation of ROS levels are strictly controlled in normal cells. Mitochondria are thought to be primary sources of ROS (Holmström and Finkel, 2014). The dysregulation of ROS induces mitochondrial DNA damage (Cadenas and Davies, 2000). Additionally, accumulating DNA damage eventually results various genomic alterations and initiates tumorigenesis. Only a few studies describe genetic events in multiple myeloma cells that primarily affect intracellular redox status during its progression (Lipchick et al., 2016). ACA11, an orphan box H/ACA class small nucleolar RNA, was shown to inhibit oxidative stress in multiple myeloma cell, which was upregulated in MM patient with t (Palumbo and Anderson, 2011; Lipchick et al., 2016) chromosomal translocation (Chu et al., 2012). The relationships between oxidative stress with survival in MM patients needs further investigation.

The hypothesis of this study is that combining cuproptosis and oxidative stress can provide a more accurate prognostic value for patients with multiple myeloma (MM). To test this hypothesis, we developed a novel 8-cuproptisis-associated-oxidative stress signature to predict the survival outcomes of MM patients by analyzing data obtained from the Gene Expression Omnibus database. This signature was helpful for risk stratification and prognosis. Furthermore we established a prognostic nomogram that could accurately predict overall survival of MM patients. Our findings may shed light on the underlying mechanisms of MM progression and suggest that the signature may serve as a promising prognostic marker for multiple myeloma patients.

2 Materials and methods2.1 Data acquisition

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.

FIGURE 1. Flow chart of this study.

2.2 Cuproptosis -related oxidative stress genes

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.

2.3 Differential gene expression analysis

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).

2.4 Functional enrichment analysis of 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).

2.5 Construction and validation of a prognostic risk model

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.

2.6 Construction and evaluation of a predictive nomogram

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.

2.7 Immune microenviroment

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.

2.8 Drug sensitivity 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).

2.9 Statistical analyses

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.

3 Results3.1 DECROGs identification and functional enrichment analysis

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.

FIGURE 2. Identification of candidate DECROGs. (A) Sankey diagram showed the correlations between CRGs and OSGs (B)The volcano plot between healthy donors and MM patients. (C) Venn diagram identifying DEGs correlated with cuproptosis-related oxidative stress genes. (C) Heatmap of the expressions of 76 DECROGs. (E) Top 10 classes of GO enrichment terms in biological process (BP), cellular component (CC), and molecular function (MF). (F) Top 10 classes of Reactome pathway enrichment terms.

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.

3.2 Exploration of the prognostic DECROGs in MM

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:

Riskscore=EXPRNASE3*−0.077+EXPAPOE*−0.214+EXPCCNB1*0.270+EXPMIF*0.290+EXPFOXO1*−0.238+EXPKIT*−0.069+EXPPLA2G4A*0.192+EXPEGR1*−0.104

TABLE 1. 26 potential genes based on univariate Cox regression analysis.

FIGURE 3. LASSO regression analysis of 26 genes. (A) LASSO regression analysis of 26 genes. (B) Cross-validation of results of LASSO regression analysis. (C) Multivariable Cox regression analysis of genes determined by lasso regression analysis.

3.3 Establishment and validation of risk model for predicting overall survival

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.

FIGURE 4. Risk score model based on eight-gene signature in GSE24080. (A) Gene expression heat map for eight prognostic genes in the high-risk and low-risk groups. (B) MM patients were divided into high-risk and low-risk groups based on the median risk score. Scatter plot of survival time and status in the high-risk and low-risk groups. (C) Kaplan–Meier curves of overall survival. (D) Assessment of the predictive ability of the model by time-dependent ROC analysis.

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.

FIGURE 5. Validation of the eight-gene risk score model in the testing datasets. (A) The Kaplan–Meier OS curves for patients in the high- and low-risk groups in the MMRF cohort (p = 0.00026, top). ROC curves showed the prognostic performance of the prognostic model in the MMRF cohort. (bottom) (B) OS in GSE4581 (p < 0.0001, top), ROC curve in GSE4581 (bottom). (C) OS in GSE2658 (p < 0.0001, top), ROC curve in GSE4581 (bottom). (D) OS in GSE136337 (p = 0.0072, top), ROC curve in GSE136337 (bottom).

3.4 Construction and validation of the nomogram prognostic model

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).

FIGURE 6. Nomogram construction based on the eight-gene signature and prognostic value of genes. Forest plot of univariate (A) and multivariate (B) Cox regression analysis of the clinical features in GSE24080. (C) Nomogram of risk score and other clinical factors for predicting MM1-, 3-, and 5-year overall survival in GSE24080. (D) The calibration plot (D), C-index (E), and DCA curves (F) of the nomogram in GSE24080.

3.5 Correlation of the risk model with clinical characteristics

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.

FIGURE 7. Validation of the eight-gene risk score model in ISS (A) The C-index of our risk model and ISS in GSE24080. (B) The C-index of our risk model and ISS in MMRF. (C) Kaplan–Meier curves of MM patients in stage III of ISS in GSE24080 (p < 0.001). (D) Kaplan–Meier curves of MM patients in stage III of ISS in MMRF cohort (p = 0.043). (E) Boxplot showed the difference in riskscore for different stages in GSE24080 (p < 0.001). (F) Boxplot showed the difference in risk scores for different stages in MMRF (p < 0001).

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.

FIGURE 8. Validation of the eight-gene risk score model in patients without genetic indicators by Kaplan–Meier curves. (A) MM patients without t (Palumbo and Anderson, 2011; Lipchick et al., 2016) (p = 0.002). (B) MM patients without amp1q (p = 0.007). (C) MM patients without del (17p) (p = 0.009). (D) MM patients without del (13q) (p = 0.016) (E) MM patients without t (Schieber and Chandel, 2014; Lipchick et al., 2016) (p = 0.014). MM patients were divided into high-risk and low-risk groups by the median risk score.

3.6 Immune microenvironment of high- and low-risk groups

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.

FIGURE 9. Characteristics of immune cell and immune function in different risk groups. (A) Comparison of the 16 immune cell scores calculated by ssGSEA between different risk groups.(B) Comparison of the 13 immune-related functions by ssGSEA between different risk groups. (C) Relations between immune cells and the risk score (*p < 0.05; **p < 0.01; ***p < 0.001)

3.7 Chemotherapy sensitivity prediction

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.

FIGURE 10. Chemotherapeutic sensitivity prediction. High-risk group showed lower sensitivity to SB216763 (A), Selumetinib (B), NU7441 (C), PF-4708671 (D), BMS-754807 (E), Doramapimod (F), JQ1 (G) and Ribociclib (H), and higher sensitivity to WIKI4 (I), I-BRD9 (J), BI-2536 (K), MIRA-1 (L), Dihydrorotenone (M), Linsitinib (N), Fulvestrant (O) and GDC0810 (P).

4 Discussion

Multiple myeloma is a highly heterogeneous hematologic malignancy, and the clinical outcomes of MM patients vary widely. What they have in common is that most cases are incurable (Siegel et al., 2019). Cell death is essential for maintaining of organismal homeostasis, and preventing excessive proliferation. Cuproptosis is a newly discovered form of death, which is closely associated with mitochondrial metabolic activity (Tsvetkov et al., 2022). Oxidative stress plays a critical role in various stages of tumorigenesis and cancer progression. Excessive ROS causes mitochondrial DNA damage, which associated with initiates tumorigenesis (Cadenas and Davies, 2000). Both cuprotosis and oxidative stress are closely related to mitochondrial homeostasis. Oxidative stress also plays an indispensable role in cell death, which can result in ferroptosis and autophagy (Chen et al., 2022; Ma et al., 2022). We combined the two in our analysis.

In our study, we first identified cuproptosis-related oxidative stress genes and then obtained DEGs between healthy donors and MM patients, after which the two were taken to intersect to obtain 76 DECROGs. A novel 8 gene signature as a prognostic marker for MM was developed by univariate, LASSO, and multivariable Cox analysis. Patients were divided into the high- and low-risk groups according to the median cut-point of risk score. There are significant differences in the Kaplan–Meier survival curve between the high-risk and low-risk groups in GSE24080. Similar results were observed in other three GEO and MMRF cohorts.

The eight genes in the risk score model: CCNB1 (Cyclin B1) is associated with the cell cycle and mitosis. In previous study, comparing to replased MM patient, the level of CCNB1 was downregulated in newly diagnosed MM patients (Dementyeva et al., 2013). MIF (Macrophage Migration Inhibitory Factor) has been shown to accelerates the development of the disease in multiple myeloma. MM patients with high MIF had a worse prognosis (Wang et al., 2020; Xu et al., 2021). A number of studies have shown that PLA2G4A (placental phospholipase A24A is involved in the biological procedure of stress response in various diseases (Brien et al., 2017; Zhang et al., 2018; Mishra et al., 2022). PLA2G4A is overexpressed in different solid cancers such as hepatocellular carcinoma (Fu et al., 2017), prostate cancer (Patel et al., 2008), breast cancer (Chen et al., 2017), cervical cancer (Xu et al., 2019), as well as malignant hematologic diseases, such as AML (acute myeloid leukemia) (Runarsson et al., 2007; Bai et al., 2020; Zhang et al., 2022) In multiple myeloma patients, expression of the PLA2G4A was higher compared to healthy individuals and high gene expression of PLA2G4A is associated with adverse outcomes (Mahammad et al., 2021). A previous study showed that FOXO1 (Forkhead Box O1) activation could inhibit the tumor growth and induce cell autophagy and cell death. In Chronic myeloid leukemia, FOXO1 is upregulated in drug-resistant cells with BCR-ABL1 kinase mutation (Wagle et al., 2016). In MM, FOXO1 mediates cell apoptosis and inhibits tumor cell growth (Liu et al., 2016). RNASE3 (Ribonuclease A Family Member 3), also known as ECP(Eosinophil cationic protein), is mainly involved in human immune function (Ostendorf et al., 2020). A previous study described the protective role of ECP under oxidative stress. It inhibits ROS-induced apoptosis in cardiomyocytes via PI3K-Akt pathway (Ishii et al., 2015). In hematological cancers, it was significantly upregulated in CML patients’ PBMCs compared with healthy controls (Yao et al., 2022). APOE (apolipoprotein E) is a known immune suppressant, which is significantly regulated in many human tumors, but its role in multiple myeloma has not been defined (Wu et al., 2022). KIT (KIT proto-oncogene) encodes a receptor tyrosine kinase. Its role in MM needs to be further investigated. EGR1 (early growth response protein 1) is characterized as a tumor suppressor in multiple myeloma (Chen et al., 2010).

A nomogram integrates multiple independent prognostic factors, and scores of each factor are calculated based on their contribution to the outcome. Then a total score for an individual patient can be calculated. Finally, the outcome for a given patient can be predictive. Therefore, we can specify a more precise treatment strategy based on the results of the nomogram (Iasonos et al., 2008). In this study, we evaluated the relations between DECROGs and patients’ clinical characteristics, and multivariate Cox regression analysis showed that risk score was one of the independent prognostic factors indicating that the DECROGs could serve as a reliable predictor of survival. Then we constructed a nomogram to predict the outcome of MM patients. The C-index, calibration curve and DCA curve demonstrated that the nomogram’s predicted value was in close agreement with the actual outcome.

The International Staging System (ISS) is the standard for staging of myeloma. Therefore, we wanted to evaluate whether the model could improve prognosis prediction in MM patients by combining with the ISS. The outcomes demonstrated that the model could optimize ISS to some extent by further differentiating patients with ISS stage III. Chromosomal change are present in the plasma cell of almost 100% of MM patients (Avet-Loiseau et al., 2009), which resulting the internal heterogeneity of MM. We evaluate the role of our model in patients with or without genetic alterations. In patients with genetic alterations, including del (17q), t (Palumbo and Anderson, 2011; Lipchick et al., 2016), del (13q), amp1q, and t (Chu et al., 2012; Schieber and Chandel, 2014), the difference was not significant, whereas in patients without these alterations, our model allows further stratification of patients by prognosis. The high-risk group demonstrated shorter overall survival than the low-risk group. Above all, the model can increase the prediction accuracy significantly combining ISS and genetic alterations.

As immune cells play a crucial role in the development of multiple myeloma, we further analyzed the relationship between the model and immune infiltration (García-Ortiz et al., 2021). In multiple myeloma, cell-mediated immunity is suppressed (Kyle and Rajkumar, 2004). Our study found that compared to the low-risk group, most immune cells were lacking in the high-risk group, indicating an overall decrease in immune activity. The reduction of these immune cells may lead to increased susceptibility to infection and the development of other diseases, and may also affect the patient’s response to multiple myeloma treatment, which is unfavorable for prognosis (Allegra et al., 2021).

Drug sensitivity research is beneficial for discovering potential therapeutic drugs. We used the Oncopredict tool to evaluate the IC50 of different drugs in the high-risk and low-risk groups by analyzing cancer drug sensitivity genomics (GDSC). High-risk patients were found to be more sensitive to MIRA-1, GDC0810, Dihydrorotenone, I-BRD9, WIKI4, Fulvestrant.1, Linsitinib, and BI-2536, providing candidate drugs for the treatment of multiple myeloma.

There are inevitably some limitations in the research work, which we hope to solve in the future.

First, this is a retrospective study and needs further prospective studies to confirm our results. Second, to increase the sample size, we employed numerous datasets. Biases between platforms are unavoidable, which may cause differences in the results. Moreover, the role of the eight-gene signature in the pathogenesis of MM remains to be addressed based on further experimental research.

In conclusion, we established an eight-gene risk model based on oxidative stress genes associated with cuproptosis. The risk model and prognosis are significantly correlated. Combining the risk model with clinical features, we established a novel nomogram that may better predict the survival in MM patients more accurately. This study offers a novel perspective on understanding multiple myeloma and provides potential targets for the diagnosis and treatment of multiple myeloma.

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found here: https://www.ncbi.nlm.nih.gov/geo/ https://research.themmrf.org/.

Author contributions

TL and QW designed and managed the study. TL, LY and YH drafted the manuscript. All authors contributed to the article and approved the submitted version.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2023.1100170/full#supplementary-material

References

Allegra, A., Tonacci, A., Musolino, C., Pioggia, G., and Gangemi, S. (2021). Secondary immunodeficiency in hematological malignancies: Focus on multiple myeloma and chronic lymphocytic leukemia. Front. Immunol. 12, 738915. doi:10.3389/fimmu.2021.738915

PubMed Abstract | CrossRef Full Text | Google Scholar

Avet-Loiseau, H., Li, C., Magrangeas, F., Gouraud, W., Charbonnel, C., Harousseau, J-L., et al. (2009). Prognostic significance of copy-number alterations in multiple myeloma. J. Clin. Oncol. 27 (27), 4585–4590. doi:10.1200/JCO.2008.20.6136

PubMed Abstract | CrossRef Full Text | Google Scholar

Bai, H., Zhou, M., Zeng, M., and Han, L. (2020). PLA2G4A is a potential biomarker predicting shorter overall survival in patients with non-M3/NPM1 wildtype acute myeloid leukemia. DNA Cell Biol. 39 (4), 700–708. doi:10.1089/dna.2019.5187

PubMed Abstract | CrossRef Full Text | Google Scholar

Bal, S., Kumar, S. K., Fonseca, R., Gay, F., Hungria, V. T., Dogan, A., et al. (2022). Multiple myeloma with t(11;14): Unique biology and evolving landscape. Am. J. cancer Res. 12 (7), 2950–2965.

PubMed Abstract | Google Scholar

Brien, M., Larose, J., Greffard, K., Julien, P., and Bilodeau, J. F. (2017). Increased placental phospholipase A2 gene expression and free F2-isoprostane levels in response to oxidative stress in preeclampsia. Placenta 55, 54–62. doi:10.1016/j.placenta.2017.05.004