Integrative analysis of the prognostic value and immune microenvironment of mitophagy-related signature for multiple myeloma

Background Multiple myeloma (MM) is a fatal malignant tumor in hematology. Mitophagy plays vital roles in the pathogenesis and drug sensitivity of MM. Methods We acquired transcriptomic expression data and clinical index of MM patients from NCI public database, and 36 genes involved in mitophagy from the gene set enrichment analysis (GSEA) database. Least absolute shrinkage and selection operator (LASSO) Cox regression analysis was conducted to construct a risk score prognostic model. Kaplan–Meier survival analysis and receiver operation characteristic curves (ROC) were conducted to identify the efficiency of prognosis and diagnosis. ESTIMATE algorithm and immune-related single-sample gene set enrichment analysis (ssGSEA) was performed to uncover the level of immune infiltration. QRT-PCR was performed to verify gene expression in clinical samples of MM patients. The sensitivity to chemotherapy drugs was evaluated upon the database of the genomics of drug sensitivity in cancer (GDSC). Results Fifty mitophagy-related genes were differently expressed in two independent cohorts. Ten out of these genes were identified to be related to MM overall survival (OS) rate. A prognostic risk signature model was built upon on these genes: VDAC1, PINK1, VPS13C, ATG13, and HUWE1, which predicted the survival of MM accurately and stably both in training and validation cohorts. MM patients suffered more adverse prognosis showed more higher risk core. In addition, the risk score was considered as an independent prognostic element for OS of MM patients by multivariate cox regression analysis. Functional pathway enrichment analysis of differentially expressed genes (DEGs) based on risk score showed terms of cell cycle, immune response, mTOR pathway, and MYC targets were obviously enriched. Furthermore, MM patients with higher risk score were observed lower immune scores and lower immune infiltration levels. The results of qRT-PCR verified VDAC1, PINK1, and HUWE1 were dysregulated in new diagnosed MM patients. Finally, further analysis indicated MM patients showed more susceptive to bortezomib, lenalidomide and rapamycin in high-risk group. Conclusion Our research provided a neoteric prognostic model of MM based on mitophagy genes. The immune infiltration level based on risk score paved a better understanding of the participation of mitophagy in MM. Supplementary Information The online version contains supplementary material available at 10.1186/s12885-023-11371-7.


Introduction
Multiple myeloma (MM) is a fatal hematologic cancer featured with abnormal propagation of monoclonal plasma cells in bone marrow [1].Over the past decades, many effective treatments, such as bortezomib (injectable proteasome inhibitor), lenalidomide (oral immunomodulatory drug), chimeric antigen receptor-engineered T cells, and autologous stem cell transplantation [2][3][4][5], improve the outcome of MM patients.However, the medical need of MM remains unmet as a result of the significant heterogeneity [6].Various factors are involved in MM progress, mainly including genetic abnormalities [7], changes in bone marrow microenvironment (BM-ME) [8], and epigenetic alterations [9].Studies confirmed that immune cells in BM-ME and the dysregulation of genes involved in immune checkpoint are related to the immune infiltration level of BM-ME [10,11].The protection from BM-ME and high genetic instability guard MM cells against chemotherapies, or receptor-targeting drugs, which virtually leads to resistance and relapse [12].
Mitochondria is not only a crucial house of producing energy via oxidative phosphorylation, but also a center of producing cellular metabolites [13].Mitophagy is a process of lysosome-dependent mitochondrial autophagy, which protects cells against proapoptotic proteins, poisonous reactive oxygen species (ROS), and the unavailing hydrolytic action of adenosine triphosphate (ATP), induced by the depolarization of mitochondrial membrane or the mitochondrial DNA (mtDNA) changes [14,15].Accumulating evidences have recently confirmed that mitophagy is a double-edged sword in cancer development.On one hand, the decrease of mitophagy promotes the cancer progression [16].On the other hand, the increased mitophagy facilitates cancer cells proliferation and progression by defending cancer cells from apoptosis [17].Nonetheless, roles of mitochondrial dysfunction in the immune microenvironment and prediction of outcome in MM remain indistinct.
In our study, we performed the least absolute shrinkage and selection operator (LASSO) Cox regression analysis to build a five-mitophagy-related-gene prognostic risk signature model.The risk model revealed a great predicted value of actual survival probabilities.In addition, enrichment analysis identified the alteration of immune checkpoint and immune microenvironment based on the risk score.Finally, drug sensitivity analysis predicted latent drugs for treating MM.

Data acquisition
MM transcriptomic expression data and clinical features were acquired from the NCI Gene Expression Omnibus database (GEO).The raw data was normalized and log2 transformed.All detailed information of GEO dataset were shown in Supplementary Table 1.The overall design of this study was shown in Supplementary Fig. 1.Among them, GSE6477 and GSE13591 were applied to appraisal differentially expressed genes (DEGs).GSE9782 was used as the training set to establish the prognostic risk score model, and GSE24080 and GSE4204 were utilized as validation cohorts.GSE24080 was applied to perform univariate analysis and multivariate analysis for overall survival (OS) rate and construct the nomogram model.GSE6477 and GSE47552 were applied to evaluate the performance of risk score for MM diagnosis.

Establishment of the prognostic risk model
First, univariate cox regression analysis was applied to obtain the OS-related genes with p < 0.05, and ten genes (SLC25A4, VDAC1, RNF41, SLC25A5, PINK1, SQSTM1, VPS13C, ATG13, HUWE1, and OPTN) were significantly correlated with MM OS time.Next, we performed LASSO Cox regression to optimize the prognostic model by further compressing the genes and constructing the prognostic model by "glmnet" package (version 4.1-1), and five genes finally came into the risk score formula (VDAC1, PINK1, VPS13C, ATG13, and HUWE1).Moreover, MM patients were divided into two groups upon the optimal cutoff of the risk score with "Survminer" package (version 0.4.9).Receiver operating characteristic (ROC) curve was applied to estimate the prognostic value of risk score model in MM patients.

Conclusion
Our research provided a neoteric prognostic model of MM based on mitophagy genes.The immune infiltration level based on risk score paved a better understanding of the participation of mitophagy in MM.

Highlights
• The study reported a mitophagy-related genes signature in multiple myeloma.
• The mitophagy-related genes signature had ideal prognostic independence in multiple myeloma.
• A nomogram to predict the overall survival of multiple myeloma was built by combining the five-gene signature, LDH and ISS stage.• Immune infiltration was related to the mitophagy-related risk signature.

Construction of the nomogram
To assess the prediction value of risk score in MM, we performed univariate and multivariate cox regression analysis.Variables with p < 0.10 was supposed to the multivariate cox regression analysis, and p < 0.05 was considered as remarkable independent prognostic factors.Then, we used "rms" package (version 6.2-0) to construct the nomogram independent prognostic model.Finally, we assessed the predictive value of nomogram by ROC curve and calibration curve.

Estimation of immune infiltration
We computed the immune score, stromal score, tumor purity, and ESTIMATE score for each sample in GSE24080 by ESTIMATE algorithm [21].According to the previous study [22], single sample GSEA (ssGSEA) was conducted to evaluate the abundance of infiltrating cells by "GSVA" package (version 1.38.2).

Prediction of drug sensitivity
To achieve precise treatment upon mitophagy-related signature and identified potential drugs for MM, we utilized the Genomics of drugs sensitivity in cancer (GDSC, https://www.cancerrxgene.org/) to predict the chemotherapeutic response [23].R package "pRROphetic" was implemented to evaluate the half-maximal inhibitory concentration (IC 50 ) by ridge regression [24].The forecast precision was determined through 10-fold cross-validation using the GDSC training set.Additional methods are described in supplementary methods.

Statistical analysis
At least three dependent experiments were performed in qRT-PCR and values were presented as mean ± SD.We used unpaired Student's t-test and Mann-Whitney U test to judge the otherness in two groups.More than two groups, we used One-way ANOVA (for parametric data) and Kruskal-Wallis (for non-parametric data) test to contrast significance.ROC curve was performed to assess the diagnosis value of risk score and the prognostic value of nomogram model in MM.The Kaplan-Meier method with a two-sided log-rank test was performed to evaluate the OS of MM patients.SPSS 21 software (SPSS, Chicago, USA) and GraphPad Prism 8 were used for statistical analysis.P value lower than 0.05 was defined markedly different.

Identification of 15 mitophagy-related genes differentially expressed in MM
We acquired a total of 36 mitophagy-related genes from GSEA GOBP database (supplementary Table 2).The datasets of GSE13591 (including 5 normal plasma cells (NPC) and 133 MM patients) and GSE6477 (including 15 NPC and 73 MM patients) were used to analyze the genes involved in mitophagy which were differentially expressed in MM.We found 19 genes and 21 genes were differentially expressed in MM patients compared to human healthy donors in GSE6477 and GSE13591, respectively (Fig. 1A and B).We observed that SLC25A4, PHB2, CERS1, VPS13C, HUWE1, VDAC1, and SLC25A5 were significantly upregulated, while OGT, ATG13, PINK1, and OPTN were significantly downregulated in MM patients (supplementary Fig. 2A and 2B).Then, venn plot was used to identify 15 mitophagy-related genes overlapped in the two datasets (Fig. 1C).To ulteriorly probe the relationship among 15 mitophagy-related genes, we performed PPI network analysis (Fig. 1D).We found PINK1, SQSTM1 and VDAC1 were hub genes.Besides, we conducted correlation analysis to ascertain the relationship among these genes (Fig. 1E).The correlation between VDAC1 and SLC25A5 was significantly positive (r = 0.678), whereas the correlation between SQSTM1 and PHB2 was significantly negative (r = -0.416).Furthermore, we used cBioPortal, an online database, and found no frequent mutation of 15 mitophagyrelated genes in 211 MM samples (Fig. 1F) [25].

The mitophagy-related risk score was an independent prognostic element in MM
We applied GSE24080 database to performed univariate and multivariate cox regression analysis.We observed lactate dehydrogenase (LDH) (hazard ratio (HR), 2.134; 95% CI, 1.39-3.277;p = 0.001), international staging system (ISS) stage (HR, 1.83; 95% CI, 1.414-2.369;p < 0.0001), and risk score (HR, 1.601; 95% CI, 1.059-2.422;p = 0.026) were the independent prognostic elements for MM survival (Table 1).In order to visualize the prognosis prediction, we further constructed a nomogram model with the above independent prognostic features (Fig. 3A).Calibration curves are typically used to assess the agreement between predicted probabilities and observed event rates or frequencies in the real world, with the 45° line representing the best prediction scenario.Calibration ROC and curves were applied to estimate the precision and efficacy of the prognostic risk model.Decision curve analysis (DCA) is usually used to evaluate the clinical utility of different predictive models [26].We observed that the precision of the total points in nomogram for 1-year survival was 0.8240 (95%CI: 0.7570 to 0.8911, P < 0.0001), 2-year survival 0.7729 (95%CI: 0.6912 to 0.8545, P < 0.0001), and 3-year survival 0.6653 (95%CI: 0.5933 to 0.7372, P < 0.0001), respectively (Fig. 3B).The calibration curves and DCA revealed a great predicted value of actual survival probabilities for 1-year, 3-year, and 5-year, respectively base on the nomogram model (Fig. 3C and D).

The mitophagy-related risk signature was related to clinical characteristics of MM
Next, we used the database of GSE24080 to estimate the relation between the risk score and clinical characteristics (Table 2).The results showed that the value of risk score was related to LDH (p = 0.03), albumin (ALB) (p = 0.025), hemoglobin (HGB) (p = 0.02), beta-2 microglobulin (B2M) (p = 0.037), ISS stage (p = 0.022), the percentage of plasma cells in bone marrow biopsy (BMPC) (p = 0.018), and cytogenetic abnormalities (p < 0.0001).Cytogenetic abnormalities has been regarded as a poor prognostic factor for MM [27].Therefore, we choose GSE136337 as the validation dataset, which consists with several cytogenetic information.We observed that MM patients with cytogenetic abnormalities, including del13q  expression (supplementary Table 3).These results indicated mitophagy-related risk signature was associated with clinical characteristics of MM, especially cytogenetic abnormalities.
To further assess the applicability of the risk signature model for different MM conditions, we performed Kaplan-Meier survival analysis.The risk score showed stabilized predictive value in term of ISS stage (ISS I/ II, p = 0.0052, ISS III, p = 0.0079, Fig. 4A and B), gender (male, p = 0.0075; female, p = 0.0016, Fig. 4C and D), the condition of cytogenetic (no cytogenetic abnormalities, p = 0.022, with cytogenetic abnormalities, p = 0.027, Fig. 4E F).Among these, MM patients had the shorter OS in higher-risk score group.However, this phenomenon was found in young people (p = 0.00096, Fig. 4G), but not in older people (p = 0.06, Fig. 4H).In addition, higher risk score was related to more adverse OS in MM patients received various treatment, including bortezomib (p = 0.017, Fig. 4I), thalidomide (p = 0.0013, Fig. 4J), dexamethasone (DEX, p < 0.0001, Fig. 4K), and PS341 (proteasome inhibitor, p < 0.0001, Fig. 4L).Taken together, the prognostic model satisfied the personalized medical needs of MM.

Immune infiltration was related to the mitophagy-related risk signature
Due to the significantly enrichment of immune response in downregulated DEGs in MM patients, ESTIMATE was performed to assess the different immune infiltration levels of MM patients in different risk group.First, Kaplan-Meier analysis revealed that low stromal score, immune score and ESTIMATE score were associated with adverse OS, while high tumor purity was associated with poor OS (Fig. 6A).In addition, MM patients had lower levels of stromal score, immune score and ESTIMATE score in the high-risk group, and higher levels of tumor purity in both training cohort and validation cohort (Fig. 6B C).Furthermore, we conducted ssGSEA to probe the disparate immune cell subsets based on the risk score.The results found 19 of 28 immune cells were notably dysregulated in high-risk group (Fig. 6D).Activated B cell, central memory CD4 T cell, activated dendritic cell, central memory CD8 T cell, CD56 bright killer cell, effector memory CD8 T cell, immature B cell, nature killer cell, nature killer T cell, and type 17 T helper cell, participated in anti-tumor immune response, were all decreased in high-risk group.Besides, the levels of immature dendritic cell, regulatory T cell, MDSC, macrophage and mast cell, belonging to immunosuppressive cells, were also decreased in MM patients with high-risk score.Moreover, we explored the relation between risk score and immune checkpoint, which was regarded as potential therapy targets for MM [28].Intriguingly, the results showed that risk score was connected with immune checkpoint negatively, including IDO1, CD276, CD86, PD-L1, and PD-L2, while positively related to CD279 (Fig. 7A).Taken together, the prognostic risk score was connected with tumor immune infiltration levels and genes involved in immune checkpoint.Therefore, our research verified the immune checkpoint inhibitors, including PD-L1, IDO1, CD276, CD86 had potential clinical values for MM patients.

Validation for mitophagy-related genes in MM samples
To validate these results of bioinformatic analysis, we collected the mononuclear cells of NDMM patients and healthy donors in bone marrow.Then, the qRT-PCR was executed.We observed the expression levels of VDAC1 were markedly upregulated in MM patients than that in healthy donors (P = 0.0056, Fig. 7B), while the expression levels of PINK1 and HUWE1 were downregulated in MM patients (P = 0.001, P = 0.019, Fig. 7C and D).Although, the expression level of VPS13C was lower and ATG13 was higher in MM patients than these in controls, no evident differences were observed (P = 0.14, P = 0.838, Fig. 7E F), maybe for the reason of insufficient sample size.

Prediction of possible drugs for MM upon risk score signature
To further identified potential drugs for MM, we explored the estimated half-maximal inhibitory concentration (IC 50 ) between the two group on the bias of the Genomics of Drug Sensitivity in Cancer (GDSC).We found decreased IC 50 of 16 drugs in MM patients with high-risk score, including bortezomib and lenalidomide, which are widely used in MM therapy (Fig. 8).Besides, the estimated IC 50 of ABT.263 (bcl-2 inhibitor), ABT.888 (PARP inhibitor), AICAR (AMPK activator), ATRA (alltrans retinoic acid), dasatinib (tyrosine kinase inhibitor), AZD8055 (mTOR inhibitor), erlotinib (EGFR inhibitor), etoposide, MG.132 (proteasome inhibitor), parthenolide (HDAC inhibitor), rapamycin (mTORC1 complex inhibitor), and thapsigargin (ATPase inhibitor) were also lower in high-risk score group (Fig. 8), which might provide novel insights into MM treatment.To sum up, the risk score can be considered as an index to choose appreciate therapeutic targets for MM precisely.

Discussion
As the second most common hematological malignancy, MM remains incurable because of drug resistance and relapse [29].There is no excellent survival prediction for MM so far.Mitophagy is a vital cellular progress, which results in degradation of dysfunctional mitochondria.Several proteins have been reported involved in mitophagy, including PINK1, Parkin, OPTN, ATG13, and p62/ SQSTM1, and FUNDC1, etc. [30].Up to now, increasing researches have revealed that mitophagy acts essential part in the development and drug sensitivity of various cancers, especially in MM.It has been reported that Parkin/PARK2 carried mutation in glioma [31], lung cancer [32], and breast cancer [33].In addition, the expression of thioredoxin was notably increased in MM cells, which resistant to bortezomib.Increased thioredoxin increased resistance to bortezomib in MM via mitophagy inactivation [34].However, as far as we known, there is no research focus on the potential relationship between mitophagy-related genes and MM prognosis.
Herein, we constructed a five-mitophagy-related risk signature for MM.The prognostic model worked well in  In our study, the mitophagy-related risk signature consists of VDAC1, PINK1, VPS13C, ATG13, and HUWE1, which can accurately predict MM prognosis.MM patients in high-risk score group suffered more adverse OS.Our findings are consistent with previous studies.VDAC1(voltage-dependent anion channel 1), a main factor in the outer mitochondrial membrane, was upregulated in multiple cancers, including lung cancer [35], breast tumor [36], cervical tumor [37], and liver cancer [38].Overexpressed VDAC1 was a prognostic element and was related to immune infiltrates in breast cancer [36].Besides, VDAC1 involved in bromodomain inhibitor resistance in breast cancer [39].VDAC1 was upregulated in mesothelioma patients.The expression levels of VDAC1 were positive related to disease stage and high expression of VDAC1 had correlation with shorter survival rates.Loss of VDAC1 could inhibit cell propagation in mesothelioma cancer cells [40].PINK1 (PTEN-induced kinase 1)/Parkin was a classical axis to regulate mitophagy in mammalian cells [41].Previous researches uncovered that PINK1 and Parkin decreased in colorectal cancer [42].Loss of PINK1 restrained mitophagy, facilitated the Warburg effect, and promoted macrophages polarization in gastric cancer [43].Besides, another study found PINK1 was reduced in MM and regulated the MOB1B-mediated Hippo-YAP/ TAZ pathway leading to MM migration and homing [44].Matrine, a natural alkaloid, accelerated the apoptosis of liver cancer cells by intercepting the PINK1/Parkin pathways and restraining mitophagy [45].Consistent with our research, Li et al. identified an autophagy-related signature for MM prognosis, which contained PINK1, EIF2AK2, KIF5B, MYC, NRG2, and VEGFA [46].These finding revealed the important role of PINK1 in MM.Recent study revealed that VPS13C (vacuolar protein sorting 13 homolog C) was descended and associated with poor prognosis in skin cutaneous melanoma [47].ATG13 (autophagy related 13) is a member of autophagy initiation complex.ULK1 could active autophagy via phosphorylating ATG13 to inhibit the progress of breast cancer [48].Zhou et al. reported that estrogen receptor α (ERα) could bind the promoter of ATG13 leading to the increasing of mitophagy.Erα inhibitor oxabicycloheptene sulfonate could reduce the expression of ATG13, leading to restrain the viability of breast cancer cells [49].Niu et al. reveled that chemotherapy drug licochalcone A could active the upstream of autophagy ULK1/ATG13 complex, inducing hepatocellular carcinoma cells apoptosis by inducing autophagy [50].Previous studied demonstrated HUWE1 played disparate role in different tumors, as an oncogene in liver cancer, lung cancer, colon adenocarcinoma, and stomach adenocarcinoma, whereas an antioncogene in skin tumors, and thyroid cancer [51].The reason for the opposite function of HUWE1 in tumor may be the different tumor microenvironment and genetic backgrounds.The function and mechanism of these five mitophagy-related genes in various tumors had been overwhelmingly illustrated.However, the further study is urgently needed to be conducted in MM.
GSEA results indicated oncogene-related gene set, including MYC targets, E2F target, and mTORC1 signaling was significantly associated with high-risk group.Our findings were in line with the lower rapamycin IC 50 in the high-risk group by subsequently drug sensitivity analysis, revealing the potential role of mTOR inhibitors in MM therapy.In fact, there is a cross-talk between mTOR signaling pathway and ubiquitin proteasome system [52].The combination of these drugs offers multiple possibilities for treating MM.A clinical trial aims to RRMM (NCT00483262) uncovered the synergistic effect of bortezomib and mTOR inhibitor.In this clinical trial, 14 of 43 MM patients obtained partial response or better [53].
The results of ESTIMATE and ssGSEA disclosed that MM patients with high-risk score suffered lower levels of stromal score, immune score and ESTIMATE score, whereas higher levels of tumor purity.Besides, cell populations participated in anti-tumor immune response, such as, immature B cell, central memory CD4 T cell, type 17 T helper cell, central memory CD8 T cell, CD56 bright killer cell, effector memory CD8 T cell, nature killer cell, and nature killer T cell, were all markedly reduced in high-risk score group.Research has confirmed the essentials of mitochondrial dynamics in immune cells [54].Paul et al. revealed the elevation of mitophagy facilitated anti-tumor immunity in intestinal epithelial cells by activated CD8 T cell through cross-dressing of dendritic cells [55].Moreover, the activation of anti-tumor immune cells demands energy.The lessen antitumor activity of immune cells may be associated with mitophagy dysfunction [56].In addition, the levels of immature dendritic cell, regulatory T cell, MDSC, macrophage and mast cell, belong to immunosuppressive cells, were also decreased in high-risk score group.A study revealed that macrophages was as a harmful prognostic index in innate immunity [57].To sum up, the results indicated the tight correlation between mitophagy-related risk score and immunosuppression, which may explain the reason of poor prognosis in MM patients with high risk.New strategies for the treatment of tumors focus on immune checkpoint inhibitors.We ulteriorly probed the correlation between risk score and immune checkpoint.These results disclosed that risk score had negative association with immune checkpoint, including IDO1, CD276, CD86, PD-L1, and PD-L2, which indicated inhibitors targeting checkpoints, such as PD-L1 may be less valid in MM patients in high-risk score group.However, a clinical trial has been reported that RRMM patients with extramedullary disease obtained benefit from the combination of PD-L1 inhibitor avelumab with radiotherapy [58].Therefore, these findings suggested that mitophagy status and sensibility to immune checkpoints may be heterogeneous between primary and recurrent MM patients.
However, there are still some limitations in our study.First, the expression of five-mitophagy-genes were verified only in mRNA level, and the expression at protein level should be further clarified.Second, the nomogram model was not be applied in validation cohort, due to the lack of clinical information in the two independent datasets.Third, the function of the five-mitophagy-genes in MM has not yet to be cleared.Therefore, further experiments in vitro and in vivo need to be conducted.

Conclusions
In conclusion, we constructed a five-mitophagy-genes (VDAC1, PINK1, VPS13C, ATG13, and HUWE1) prognostic risk model, which as an independent element for MM OS, could estimate the survival of MM accurately and stably both in training and validation cohorts.The molecular landscape characteristics upon the risk score, including the regulatory pathways, immune level, and potential drug targets, improved our cognition for MM, which provide novel insights into MM treatment.

Fig. 2
Fig. 2 Construction of prognostic risk model based on mitophagy-related genes.(A).Forest plot of the univariate cox regression in GSE9782.(B).The LASSO Cox analysis identified mitophagy-related genes in the training cohort.(C).Partial likelihood deviance of different numbers of variables.Onethousand-fold cross-validation was applied for tuning penalty parameter selection.(D).Kaplan-Meier survival analysis in training cohort based on the risk score.(E).ROC curve for the prognostic risk signature.LASSO: least absolute shrinkage and selection operator; ROC: receiver operating characteristic

Fig. 3
Fig.3The mitophagy-related risk score was an independent prognostic factor in MM. (A).Nomogram predicting 1-, 3-, and 5-year survival for MM patients based on mitophagy-related genes risk score.To use this nomogram, the specific point for each variable of the patient lies on each variable axis.Draw a vertical line upward to determine the point at which each variable accepts; the sum of these points is located on the Total Points axis, and draw a vertical line down to the survival axis to determine the probability of 1-, 3-and 5-year overall survival.LDH: lactate dehydrogenase; ISS: international staging system.(B).ROC curve for the nomogram prognosis system.(C).DCA curve for the nomogram prognosis system.The x-axis represents the threshold probability, while the y-axis stands for the net benefits.(D).Calibration plot of the nomogram for 1-year, 3-year, and 5-year OS.The risk score in normal donors and MM patients with different stages in GSE6477 (E) and GSE47552 (F).(G).ROC curve for MM with different myeloma stages.DCA: decision curve analysis * p < 0.05, ** p < 0.01, *** p < 0.001

Fig. 4
Fig. 4 Survival analysis in different subgroup.Kaplan-Meier survival analysis of MM patients grouped by ISS stage I/II/III (A, B), male and female (C, D), no cytogenetic abnormalities and with cytogenetic abnormalities (E, F), age < 65 and ≥ 65 (G, H), bortezomib treatment and thalidomide treatment (I, J) in GSE24080.Kaplan-Meier survival analysis of MM patients stratified by dexamethasone treatment and PS341 treatment (K, L) in GSE9782.

Fig. 5
Fig. 5 Enrichment analysis of mitophagy-related risk signature.Gene ontology and KEGG pathway enrichment analysis of DEGs in GSE24080 dataset grouped by risk score.(A).The significantly enriched gene ontology biological process (GOBP) of downregulated DEGs.(B).The significantly enriched KEGG pathway of downregulated DEGs.(C).GSEA analysis results of KEGG gene set, HALLMARK gene set, and GOBP gene set.NES, normalized enrichment score; FDR, false discovery rate

Fig. 6
Fig. 6 Immune infiltration was related to the mitophagy-related risk signature.(A).Kaplan-Meier survival analysis upon stromal score, immune score, ESTIMATE score, and tumor purity.(B).The distribution of stromal score, immune score, ESTIMATE score, and tumor purity upon risk score in GSE24080.(C).The distribution of stromal score, immune score, ESTIMATE score, and tumor purity upon risk score in GSE9782.(D).The heatmap of the comparison in 28 immune-related gene sets upon risk score in GSE24080.* p < 0.05, ** p < 0.01, *** p < 0.001

Table 2
The correlation of prognostic risk signature and clinical characteristics

Table 3
Impact of risk-score on the most common cytogenetic abnormalities of multiple myeloma (GSE136337) Bold values indicate statistically significant p values less than 0.05