A Monocyte-Driven Prognostic Model for Multiple Myeloma: Multi-Omics and Machine Learning Insights

Introduction

Multiple myeloma (MM) is an incurable malignancy characterized by the proliferation of malignant plasma cells (PCs) within the bone marrow (BM).^1,2 Despite advancements in treatment, MM exhibits high recurrence and drug resistance rates, with a median survival time of only 5 to 6 years.³ Furthermore, its prevalence continues to rise.⁴ In recent years, immunotherapy, including immunomodulatory drugs (IMiDs), proteasome inhibitors (PIs), and immune checkpoint inhibitors (ICIs), has emerged as a promising therapeutic approach for MM.^5,6 Unfortunately, ICIs have yielded disappointing clinical outcomes in MM due to the complex tumor niche composition and highly suppressive immune microenvironment.⁷ Elucidating the composition of the tumor microenvironment (TME) in MM is essential for optimizing treatment strategies and enhancing the efficacy of ICIs.

The pathophysiology of MM is influenced by BM infiltration of monoclonal PCs, which leads to complex interactions between tumor cells and immune cells within the BM-TME.^8,9 Monocytes serve as a bridge between innate and adaptive immunity, influencing the tumor microenvironment through cytokine production, antigen presentation, and differentiation into macrophages or dendritic cells. Their dysfunction has been implicated in immune evasion and treatment resistance in various malignancies, including MM.^10,11 Alterations in the BM monocytes are evident even at the precursor stage of MM, including a decline in monocyte function and impaired chemotaxis.^12,13 Moreover, peripheral blood monocyte count serves as a dynamic prognostic biomarker in MM,¹⁴ and the proportion of monocytes in the BM has been found to increase in direct correlation with tumor cell load.¹⁵ These features manifest monocytes are associated with patient prognosis and response to immunotherapy.

Here, we first employed Mendelian randomization (MR) analysis to identify immune cells with significant causal associations to MM. Flow cytometry further confirmed that monocytes were significantly increased in MM. Subsequently, differentially expressed genes related to monocytes were identified using single-cell RNA sequencing (scRNA-seq). Based on 482 prognostically consistent monocyte-related genes, we developed a monocyte-related gene prognostic signature (MGPS) by integrating 101 machine learning algorithms and transcriptome sequencing data. In both the training and validation cohorts, MGPS demonstrated consistent and superior predictive performance in forecasting overall survival (OS) and response to immunotherapy. Overall, our study offers a promising tool for guiding the clinical management and personalized treatment of MM.

Materials and Methods

Clinical Samples and Flow Cytometry Analysis

This study was approved by the ethical committee of the Institute of the Third Xiangya Hospital of Central South University. Written informed consents were obtained from patients and healthy donors before sample collection. BM mononuclear cells (BMNCs) were obtained from 11 healthy controls (HCs) and 21 MM patients (Supplement Table 1) simultaneously and isolated by Ficoll density-gradient centrifugation. BMNCs were subsequently stained and analyzed by flow cytometry. We used CD3, CD4, and CD8-specific markers to gate T cells; CD19-specific markers to gate B cells; CD3, CD16, and CD56-specific markers to gate NK cells; and CD33, CD14, and HLA-DR-specific markers to gate monocytes. All samples were analyzed using a Navios flow cytometer (Beckman Coulter) and FlowJo 10.8.1 software (see Figure S1 for gating strategies). Detailed antibody information is provided in Supplementary Table 2.

Mendelian Randomization Analysis

Exposure and Outcome Data Acquisition

The 731 immunophenotypes GWAS summary statistics for each immune trait can be accessed from the GWAS Catalog, with accession numbers ranging from GCST90001391 to GCST90002121.¹⁶ The original GWAS on 731 immunophenotypes utilized data from 3,757 European individuals. Detailed information on the 731 immunocyte phenotypes is provided in Supplement Table 3. The GWAS summary statistics for MM were obtained from the FinnGen study posted in 2022 (https://r9.risteys.finngen.fi/endpoints/O15_PRE_OR_ECLAMPSIA), which included 1,249 cases and 299,952 controls.

Selection of IVs

There are three core assumptions for MR design: (1) genetic variants directly affect exposures; (2) genetic variants are not associated with potential confounders; (3) genetic variants affect outcomes only via the effects on exposures. If these assumptions are met, MR analyses can infer causality without the bias from unmeasured confounders.¹⁷ To identify robust and independent single nucleotide polymorphisms (SNPs) associated with exposures, we employed a multi-step approach.

First, SNPs were selected from the GWAS data with a significance threshold of p < 5×10⁻⁶. Next, SNPs were retained based on linkage disequilibrium, using a threshold of r² > 0.001 and a distance of 10,000kb between SNPs. Subsequently, the PhenoScanner database was used to check each SNP for associations with potential confounders and other outcome-related traits.¹⁸ Palindromic SNPs with intermediate allele frequencies were removed. The F statistic for each instrument was above the threshold of 10, indicating that all SNPs are robust.¹⁹ Figure 1 illustrates the overall study design.

Figure 1 Flow chart of the study.

Analysis

The primary analysis for the MR study was conducted using the inverse variance weighted (IVW) method. To assess potential violations of model assumptions, we performed comprehensive sensitivity analyses, including MR-Egger, weighted median, simple mode, and weighted mode approaches. Subsequently, the results were subjected to sensitivity analyses such as the heterogeneity test and horizontal pleiotropy test.²⁰ To control for the impact of multiple comparisons, we applied a multiple test adjustment using the false discovery rate (FDR) method. Causality was concluded if the IVW method yielded a P_FDR < 0.05, the direction of the causal effect was consistent across the four methods in sensitivity analyses, and there was no evidence of heterogeneity or pleiotropy.

Single-Cell RNA-Seq Data Collection and Analysis

Single-cell RNA-seq datasets (GSE124310)⁸ of BMNC including healthy individuals (normal bone marrow, NBM = 9) and 23 patients (smoldering multiple myeloma, SMM = 11; monoclonal gammopathy of undetermined significance, MGUS = 5; MM = 7) were analyzed according to the workflow (https://github.com/hbctraining/scRNA-seq)²¹ in R (version 4.3.0). For data quality control (QC), we retained cells with less than 10% mitochondrial gene content and genes expressed in at least three cells within the 500–7000 expression range. Expression data were normalized using the “Log-Normalize” function in the “Seurat” package. Principal component analysis (PCA) was performed on the top 3000 variable genes. Uniform manifold approximation and projection (UMAP) was employed for dimensionality reduction and cluster identification. Clusters were identified with a resolution of 0.2 and annotated according to marker genes corresponding to different cell types. We performed pseudobulks processing on the singlecell data and then used the “DESeq2” package for analysis. The count matrix was extracted and the differential expression genes (DEGs) between MM and NBM monocytes were identified using the “DESeq” package. Monocyte-DEGs (moDEGs) with a p-value of less than 0.05 were selected for further analysis.

Bulk RNA-Seq Data Collection and Analysis

The MMRF‐CoMMpass dataset (NCT01454297), comprising bulk RNA‐seq data for BM samples from MM patients, including those with newly diagnosed multiple myeloma (NDMM) and relapsed/refractory multiple myeloma (RRMM), was downloaded from the UCSC Xena (https://xenabrowser.net/datapages/). Samples with incomplete survival information were excluded, resulting in a final cohort of 844 MMRF-MM patients. The MMRF cohort was randomly divided into an internal discovery set and an internal validation set in a 7:3 ratio, with the distribution of survival characteristics balanced between the two sets. Additionally, two datasets from the GEO database, GSE136337 (n=436 NDMM samples)²² and GSE24080 (n=556 NDMM samples),²³ were included as independent sets for external validation (Supplement Table 4).

Construction of Prognostic Signature by Integrative Machine Learning Approaches

To develop a reliable prognostic model with high predictive accuracy, we integrated 10 classical algorithms: least absolute shrinkage and selection operator (LASSO), Stepwise Cox, random forest (RSF), CoxBoost, elastic network (Enet), gradient boosting machine (GBM), survival support vector machine (Survival-SVM), supervised principal components (SuperPC), ridge regression and partial least squares regression for Cox (plsRcox). The signature generation procedure was as follows:^24,25 (a) Univariate Cox regression identified the prognostic moDEGs in the MMRF-discovery dataset; (b) A total of 101 algorithm combinations were performed to match prediction models based on the leave-one-out cross-validation (LOOCV) framework in the MMRF-discovery dataset; (c) All models were evaluated in the MMRF-validation dataset, GSE136337 and GSE24080 cohorts; (d) The predictive performance of each model was assessed by calculating the Harrell’s concordance index (C-index) of in all validation cohorts, and the model with the highest average C-index was considered optimal. Consequently, we established a final model, termed the Monocyte-related Gene Prognostic Signature (MGPS). In addition, we collected model indices from previous researchers and compared the MGPS with the previous models.

Extraction of Total RNA from Bone Marrow Samples and RT-qPCR

To confirm the expression pattern of the MGPS signature genes in MM more precisely, BMNC samples of MM patients and healthy controls were used for this study. 1 mL of each sample was taken separately, and all intracellular RNA was extracted by Trizol reagent, and the quality of the extracted RNA were detected by nanodrop (Thermo scientific). The extracted RNA was reverse transcribed to cDNA according to manufacture instructions to detect the following genes expression. The BlazeTaq™ SYBR^® Green qPCR Mix2.0 kit (Genecopoeia) and the following reaction system were utilized to perform RT-qPCR reactions next. Primer sequences are shown in Supplement Table 5. The CT values of each gene were counted, and the relative expression of characteristic genes was analyzed according to the 2-ΔCt method using GAPDH as the internal reference gene.

Statistical Analysis

All data cleaning, analyses and result visualization were performed with R (version 4.3.0). The chi-squared test was applied to compare categorical variables, and the Wilcoxon rank-sum test or t-test was employed to compare continuous variables. The “survminer” package was used to determine the optimal cut-off value. Cox regression and Kaplan–Meier analyses were performed using the “survival” package. The “timeROC” package was used to conduct the Receiver Operating Characteristic (ROC) curve analysis. The nomogram was generated by the “rms” package, and the accuracy of the nomogram was evaluated via ROC curves and calibration curves, and its net clinical benefit was evaluated through decision curve analysis (DCA). The Sankey plot was visualized using the SankeyMATIC online tool (https://sankeymatic.com/build/). Gene Ontology (GO) analysis was conducted using the “clusterProfiler” package. Gene Set Enrichment Analysis (GSEA) was executed using the “enrichr” and “GSEABase” package. Gene Set Variation Analysis (GSVA) was performed using the “GSVA” package. In the Mendelian randomization analysis, FDR-corrected p values were used to evaluate confirmatory positive results. A two-sided p-value of < 0.05 was considered statistically significant (*P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001, NS: not significant).

Results

Complex Causal Association Between 731 Immune Phenotypes and Multiple Myeloma

To elucidate the causal relationship between immune phenotypes and MM at the genetic level, we conducted a bidirectional MR analysis. The workflow of our study is depicted in Figure 1. The causal effects of immune cells on MM are presented in Figure 2A, showing that the trait of the following seven immune cells traits are positively correlated with the development of MM (OR>1, P_FDR<0.05), including B cell panel: BAFF-R on CD20⁻; Monocyte panel: CD14 on CD14⁺CD16⁻ monocyte, HLA-DR on CD14⁺ monocyte; Myeloid cell panel: CD33^bright HLA-DR⁺CD14⁻%CD33^bright HLA-DR⁺, CD33 on CD14⁺ monocyte, CD33 on CD33^brightHLA-DR⁺CD14^dim, CD33 on CD33^dimHLA-DR⁺CD11b⁺, CD33 on CD33^dimHLA-DR⁺CD11b⁻, CD33 on CD66b⁺⁺ myeloid cell, CD33 on Mo MDSC, CD33 on Im MDSC, CD33 on CD33^bright HLA-DR⁺, CD33 on CD33^brightHLA-DR⁺CD14⁻, CD11b on Mo MDSC, CD11b on CD33^brightHLA-DR⁺CD14^dim; TBNK panel: CD8^dim AC, HLA-DR⁺CD4⁺ AC, CD3 on NKT, FSC-A on HLA-DR⁺NK, SSC-A on NK; Treg panel: CD127 on CD45RA⁺CD4⁺. The remaining 14 traits reduces the incidence of MM (OR<1, P_FDR<0.05), including B cell panel: IgD⁺CD38⁻ AC, IgD⁻CD38^bright%B cell, CD25 on IgD⁺CD38^bright, CD25 on IgD⁻CD38^bright, CD38 on transitional; cDC panel: Plasmacytoid DC %DC, CD62L on CD62L⁺ DC; Maturation stages of T cell panel: CD3 on TD CD8^bright; Myeloid cell: CD66b on CD66b⁺⁺ myeloid cell; TBNK panel: CD4⁺%T cell, B cell % CD3⁻ lymphocyte; Treg panel: Secreting Treg AC, CD25^high AC, CD28 on CD39⁺CD8^bright. For details of 731 characterization immunophenotypes, refer to Supplement Table 3.

Figure 2 Mendelian randomization analysis results of forest plots and flow cytometry analysis results. (A) The causal association between immune cell traits on multiple myeloma. (B) The causal association between multiple myeloma and immune cell traits. (C) Flow cytometry violin plot shows the proportion of B cells in CD3− lymphocytes from HC and MM patients (HC: n=10; MM: n=10). Flow cytometry analysis shows the MFI expression of CD14 (D) and HLA-DR (E) on CD14+monocyte from HC (n=11) and MM patients with different tumor burdens (Low burden: n=11; High burden: n=10). NS: not significant, *p < 0.05, **p < 0.01.

The causal effects of MM on immunocytes are shown in Figure 2B. Specifically, MM was identified as a risk factor for three immunophenotypes (OR>1, P_FDR<0.05): CD66b on CD66b⁺⁺ myeloid cell, CD45 on CD66b⁺⁺ myeloid cell, CD11b on CD66b⁺⁺ myeloid cell (Myeloid cell panel). Conversely, MM acted as a protective factor for three immunophenotypes (OR<1, P_FDR<0.05): CD27 on T cell (B cell panel), CD62L on CD62L⁺ DC (cDC panel), CD28 on secreting Treg (Treg panel).

Among the immune cell-related risk factors for MM, HLA-DR and CD14 were the most frequently observed (12/21, Figure 2A). To further confirm immune cell variations in MM, we conducted a flow cytometry analysis. Our findings indicated a decreased proportion of B cells within CD3⁻ lymphocytes in MM patients (Figure 2C), which aligns with their protective role in the onset of MM according to MR results. In MM patients with lower tumor infiltration (MM cell < 10%), we consistently observed a significant increase in the MFI of CD14 on CD14⁺ monocyte in the BM (Figure 2D). Similarly, in these patients, we also found a significant increase in the MFI of HLA-DR on CD14⁺ monocytes in the BM (Figure 2E). No significant variations were detected in other immunophenotypes between the two groups (Figure S2). These findings validate that the increased MFI of CD14 and HLA-DR on CD14⁺ monocytes is causally associated with MM, particularly in the early stages of MM when the tumor burden is lower. Both membrane proteins are indicative of the immunosuppressive state of monocytes in cancer.^26,27 Therefore, we further analyzed the monocyte gene characteristics in the MM BM-TME through scRNA-seq.

Identification of Monocyte Differential Expression Genes by scRNA-Seq

We obtained the scRNA-seq data (GSE124310) for analysis after a series of cleaning and quality control procedures, which included 23,191 cells. Based on marker genes for distinct cell types (Supplement Table 6), we annotated the cells into 13 major categories, including CD4 cells, CD14 monocytes, CD8 cells, NK cells, plasma cells, B cells, hematopoietic stem cells, CD16 monocytes, pDCs, mDCs, pre-B cells, hematopoietic progenitor cells, and plasmablasts (Figure 3A). The doughnut chart illustrates the variation in cell proportions across different stages of the disease (Figure 3A). Additionally, we isolated the monocyte subsets and defined four subsets after further dimensionality reduction,²⁸ including CD14⁺CD16⁻, OxPhos, HLA-DR^high and CD14^lowCD16⁺ monocytes. Notably, our results showed that HLA-DR^high monocytes increased significantly with the progression of MM disease, consistent with the findings from our MR and flow cytometry results (Figure 3B). Furthermore, we identified 1,447 differentially expressed genes (moDEGs) (p < 0.05) between the MM and NBM groups for subsequent analysis (Figure 3B and Supplement Table 7).

Figure 3 Identification of Monocyte-related candidate genes and construction of MGPS based on machine learning. (A) UMAP visualization of the immune cell distribution across different stages (NBM, MGUS, SMM, MM). The proportion of various immune cell types is shown for each stage. (B) Volcano plot identifying differentially expressed genes in monocytes (moDEGs) between MM and NBM. UMAP visualization illustrates the distribution of monocyte subsets across different groups. (C) Venn diagram displaying 1,218 overlapping genes among the CoMMpass-MMRF dataset, GSE136337, GSE24080, and moDEGs. (D) C-indexes of 101 machine-learning algorithm combinations across three cohorts. (E) Coefficients of 18 model genes derived from the MGPS using the CoxBoost and RSF algorithms. (F) Overall survival (OS) analysis in MM patients based on the expression levels of each signature gene in the MMRF cohort.

Construction of a Prognostic Signature by Machine Learning Approaches

To ensure that the moDEGs are consistently expressed across multiple datasets, we first identified the intersection of the moDEGs with the MMRF, GSE136337 and GSE24080 cohorts, and yielded 1,218 overlapping genes (Figure 3C). Subsequently, to develop prognostic signatures based on these overlapping moDEGs, we use the univariate Cox regression analysis in the MMRF-discovery dataset, which identified 482 prognostic moDEGs associated with overall survival (OS) (Supplement Table 8). Next, we integrated 101 prediction models using 10 machine learning algorithms, including LASSO, Stepwise Cox, RSF, CoxBoost, Enet, GBM, Survival-SVM, SuperPC, Ridge and plsRcox. To assess the robustness of these models and identify the most effective prognostic signature, we employed a tenfold cross-validation approach on the MMRF-validation cohort (n=253) and two external validation cohorts (GSE136337, n=436; GSE24080, n=556). Finally, among 101 models generated, CoxBoost plus RSF was chosen based on its superior concordance index (C-index) across validation datasets, ensuring stability and predictive accuracy (Figure 3D). These combined algorithms identified 18 key genes, which were used to develop a highly reliable prognostic model termed the Monocyte-related Gene Prognostic Signature (MGPS) (Figure 3E).

Of these 18 genes, 11 were positively correlated with a higher risk of poor prognosis in MM, while the remaining 7 were negatively correlated (Figure 3F) (Supplement Table 9). Therefore, to determine the association between the MGPS gene expression levels and clinical outcomes in MM patients, we conducted survival analyses for each model gene using the MMRF dataset. As expected, for the 11 positively correlated genes (coef > 0), our results showed that patients with higher expression had significantly shorter OS (Figure 3F). On the contrary, for the 7 negatively correlated genes (coef <0), higher expression levels of five genes, except for the ENY2 and SNRPC with minimal coefficients (−0.01 and −0.009, respectively), were associated with better OS (Figure 3F). This indicates that the level of model genes is indeed closely related to the survival outcomes of MM patients. Moreover, in the GSE136337 cohorts, most genes exhibited similar trends (Figures S3 and S4A).

MGPS Exhibits Robust and Stable OS Predictive Performances

To comprehensively assess the robustness of MGPS, we calculated MGPS scores for each patient and categorized them into either low-risk or high-risk groups based on the median MGPS score. The MGPS score was formulated based on the coefficients and categorical values of expression level as follow: MGPS score = (0.0582 x PFDN2) + (0.2206 x TUBA1B) + (0.1787 x ERH) + (0.0483 x FAM49B) + (0.1253 x HNRNPR) + (0.1223 x RHOC) + (0.2697 x ANAPC11) + (0.0369 x IL32) + (0.1522 x GLIPR2) + (0.0723 x ZNF90) + (0.1402 x YWHAZ) - (0.1041 x PECAM1) - (0.0101 x ENY2) - (0.2086 x PSAP) - (0.0086 x SNRPC) - (0.2099 x IL16) - (0.1319 x PCGF5) - (0.2628 x ERP29).

To interrogate the difference in survival between low and high-risk groups, we first conducted the Kaplan–Meier analysis based on the discovery MMRF cohort, which demonstrated that the low-risk group had significantly better OS than the high-risk group (p<0.0001, Figure 4A), and consistent results were observed in the validation MMRF cohort (p<0.0001, Figure 4B), GSE136337 cohort (p<0.0001, Figure 4C) and GSE24080 cohort (Figure S4B). Additionally, the number of patient deaths increased progressively with higher MGPS scores in all cohorts (Figure 4D–F). The AUCs of 1-, 2-,3- and 4-year OS were 0.793, 0.819, 0.808 and 0.926 in the discovery MMRF cohort; 0.718, 0.802, 0.836 and 0.798 in the validation MMRF cohort; and 0.675, 0.679, 0.688 and 0.670 in the GSE136337 cohort (Figure 4G–I). The calibration curves further confirmed the good predictive performance of MGPS (Figure S5A and S5B). Taken together, these data support the idea that MGPS could accurately and robustly predict the prognosis of MM patients, suggesting that MGPS may become an attractive tool for clinical practice.

Figure 4 Survival analysis and predictive performance evaluation of MGPS. Overall survival analysis in low- and high-risk MM patient groups across the Discovery MMRF cohort (A), Validation MMRF cohort (B), and GSE136337 cohort (C). Distribution of MGPS according to survival status and time in the Discovery MMRF(D), Validation MMRF cohort (E), and GSE136337 cohort (F). (G–I) Receiver operating characteristic (ROC) curves showing the prediction of 1-,2-,3-, and 4-year OS in low- and high-risk groups across the cohorts. (J) Multivariate Cox regression analysis shows that the MGPS score serves as an independent prognostic factor for MM patients. (K) Time-dependent ROC analysis predicting 1-, 2-, 3-, and 4-year OS for Cox models with or without MGPS scores in the GSE136337 cohort. (L) Nomogram constructed for predicting OS in MM patients and (M) Decision curve analysis (DCA) for 5-year survival prediction. (N) MGPS compared with R-ISS and other published molecular signatures. ***p < 0.001.

Subsequently, to assess the independent predictive value of MGPS, we performed a univariate Cox regression analysis of OS combined with the clinical features of the GSE136337 cohort (Supplement Table 10). Significant variables identified in the univariate analysis were further included in a multivariate Cox regression analysis (Figure 4J). Remarkably, the risk scores based on MGPS showed statistically significant in both the univariate (p=6.2e-10) and multivariate analyses (Supplement Table 10). These findings indicate that MGPS is an independent prognostic factor for MM patients (HR 2.72, 95% CI: 1.84–4.0, p<0.001).

To enhance the clinical applicability of MGPS, we developed an OS-nomogram incorporating MGPS and clinical characteristics (Figure 4L). The AUCs of the nomogram were 0.783, 0.765, 0.751 and 0.756 at 1-, 2-, 3- and 4-year intervals, respectively, outperforming models without MGPS scores (0.769, 0.742, 0.695, 0.695) (Figure 4K). Decision curve analysis (DCA) further demonstrated that the nomogram provided a greater net clinical benefit than models without MGPS scores (Figures 4M and S5C). These findings suggest that the MGPS-based nomogram exhibits robust and reliable predictive performances. To compare the prognostic efficacy of MGPS with R-ISS and existing MM molecular signatures, we integrated previous studies that used different biologically significant features, such as glycolysis,²⁹ cuproptosis,³⁰ ubiquitin,³¹ lactylation,³² ferroptosis,³³ hypoxia,³⁴ mitophagy,³⁵ inflammatory,³⁶ necroptosis³⁷ etc. Notably, MGPS exhibited better C-index performance than almost all models in the MMRF, GSE136337, and GSE24080 datasets (Figure 4N). In conclusion, these findings confirm the idea that MGPS is a more effective prognostic model for MM.

Evaluation of the Clinical Features of MGPS

Next, we investigate the distribution of 18 MGPS genes across different MGPS-based risk groups in the training and validation cohorts (Figure 5A–C and E). Five genes—PCGF5, PECAM1, PSAP, ERP29, and IL16—were consistently enriched in the low-risk group across all three cohorts, while the remaining genes were predominantly enriched in the high-risk group. Furthermore, we compared several classical clinical characteristics of MM across the two MGPS risk groups, finding that patients in the high-risk group exhibited higher levels of β2m, LDH, and plasma cell percentage, but lower albumin levels (Figure 5D and F).

Figure 5 The clinical features and the immune landscape associated with MGPS in MM. Distribution of clinical characteristics and the expression levels of model genes according to the MGPS risk score in the Discovery MMRF cohort (A), Validation MMRF cohort (B), and GSE136337 cohort (C). Differences in clinical characteristics between high-risk and low-risk patient groups (D and F). Sankey diagram illustrating the interrelationships between risk group, age, gender, ISS stage, and plasma cell percentage in MM patients within the GSE136337 cohort (E). Spearman correlation analysis between model genes and inhibitory immune checkpoints or their ligand genes across the GSE136337 (G), MMRF (H), and GSE24080 (I) cohorts. TIDE score and T cell exclusion analysis between different risk groups (J) (NS: not significant, *p < 0.05, **p < 0.01, ***p < 0.001).

Predictive Value of MGPS in Immune Escape and Immunotherapy Response

Based on our previous evidence demonstrating that high MGPS is associated with poor prognosis in MM, we sought to explore the potential of the MGPS in predicting the response to immunotherapy and immune escape. Given that the ICIs have been employed in preclinical and clinical trials as the primary strategy for immunotherapy, we compared the Spearman correlation of 18 MGPS-genes with several immune checkpoint receptors and ligands. The results show that low-risk enriched genes, like PCGF5, PECAM1 and PSAP, were significantly positively correlated with the expression of ligands genes (Figure 5G–I). This suggests that the low-risk group may be more responsive to immunotherapy. Subsequently, the tumor immune dysfunction and exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu/login/) validated that higher MGPS scores were associated with an increased probability of immune evasion (Figure 5J), indicating that immune checkpoint blockade therapy may be less effective in these patients. Collectively, these findings suggest that patients in the low-risk group identified by MGPS are more likely to benefit from immunotherapy.

The Correlation of the MGPS with Single-Cell Characteristics

We next investigated the expression of MGPS genes in immune cells within the MM BM-TME at the single-cell level, revealing that MGPS genes are broadly expressed across various immune cell types, including pre-B cell, DC cell, monocyte, hematopoietic progenitor cell (Figure 6A). Then, the monocyte subsets were extracted and the distribution of MGPS genes was visualized using UMAP (Figure 6B). Notably, our results revealed that 6 genes positively correlated with a higher risk of poor prognosis in MM, were increased with the progression of MM disease, including ANAPC11, GLIPR2, YWHAZ, TUBA1B, ERH, and FAM49B. Besides, PSAP was highly expressed across the entire monocyte population (Figure S6). Subsequently, we conducted trajectory analysis and calculated the contribution of MGPS genes during monocyte development. The results showed that ERP29 and ZNF90 are highly expressed during the initial phases of monocyte development, GLIPR2 and ENY2 are mainly expressed in the intermediate stage,ANAPC11 and ERH are expressed in both the early and late stages. Other signature genes are primarily expressed mainly in the late stages of monocyte development (Figure 6C).

Figure 6 The correlation of the MGPS with single-cell characteristics. (A) Violin plot displaying the expression levels of model genes across different immune cell types. (B) UMAP visualization of the expression distribution of model genes specifically within monocytes. (C) Cell trajectory and pseudotime analysis of the identified monocyte subsets. (D) Bar graph illustrating GO enrichment pathways associated with model genes. (E) Heatmap showing GSVA enrichment analysis pathways of differentially expressed genes between the two risk groups. (F) Gene enrichment in the high or low risk group was analyzed by GSEA (MSigDB: h.all.v2024.1.Hs.symbols.gmt). (G) RT-qPCR analysis of MGPS signature genes (coef< 0) relative expression levels between HC and MM patients. (H) RT-qPCR analysis of MGPS signature genes (coef> 0) relative expression levels between HC and MM patients. (NS: not significant, *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001).

To further investigate the molecular mechanisms of MGPS genes in MM, we performed the GO enrichment analysis. The results indicated that these genes are involved in several biological processes and cellular components, including protein homodimerization activity, extracellular exosome, melanosome, nucleus, extracellular space, and signal transduction (Figure 6D). Additionally, the GSVA analysis based on the KEGG pathways showed that DGEs between the high and low-risk groups exhibited significant differences (padj. < 0.05) in pathways related to ribosome (protein synthesis) and renin angiotensin system (Figure 6E). Moreover, GSEA demonstrated that the high-risk group was more enriched in the hypoxia signature, whereas the low-risk group exhibited enhanced DNA repair function and interferon response capability (Figure 6F). These findings suggest that a high MGPS score reflects a state of cellular functional impairment. To explore the role of MGPS genes in the development of MM, we conducted qPCR validation to assess the expression levels of the 18-genes signature across both HC and MM patients. The qPCR results showed a significant upregulation in the expressions of ERP29, IL16, PSAP, PCGF5, PECAM1, ENY2 and SNRPC in HC (Figure 6G). Higher expression of ERH, PFDN2, ANAPC11, RHOC, FAM49B, YWHAZ, ZNF90 and GLIPR2 accumulated in MM patients (Figure 6H). The evidence supports the notion that MGPS genes may participate in the development and progression of MM through the regulation of monocyte function.

Discussion

Over the past three decades, MM has shown a rising global incidence, with a ~1.36-fold increase from 1990 to 2019.³⁸ Despite significant advancements in understanding its pathophysiological mechanisms in recent years,^39–42 the causal relationship between MM and immunophenotypes remains unconfirmed. Our study is the first to investigate these causal relationships between MM and immune cell traits using extensive publicly available GWAS genetic data. We identified seven immune cell traits that exhibit significant causal effects on MM, while MM, in turn, demonstrates substantial causal effects on three immune cell traits.

Monocytes, which are involved in immune responses such as give rise to macrophages, the destruction of microbes, and tumor cells,¹⁰ have been identified as risk factors for MM in both Mendelian randomization and flow cytometry analyses. It has been reported that MM patients may have functional defects in monocytes,¹⁰ including the dysregulation of MHC-II molecules on CD14⁺ monocytes, resulting in T cell suppression in vitro.⁸ Furthermore, overexpression of IL21R in CD14⁺ monocytes has been shown to enhance osteoclast formation in MM.⁴³ Our study finds that the MFI of CD14 in CD14⁺CD16⁻ monocytes, along with the MFI of CD33 and HLA-DR in CD14⁺ monocytes, is positively associated with MM.

Using single-cell sequencing, we identified monocyte-related differential expression genes between MM patients and healthy individuals. These genes were further used to establish a consistent MGPS by applying 101 combination models derived from 10 machine-learning methods. Previous prognostic models in MM have relied on transcriptomic signatures or clinical staging, but often lack integration of immune cell-specific markers. Our approach uniquely combines multi-omics data with Mendelian randomization, ensuring a causality-driven selection of prognostic genes. Notably, our results displayed that the MGPS-based high-risk group is more likely to experience immune escape, potentially correlating with drug resistance,⁴⁴ while the low-risk group exhibits a better response to immunotherapy and is more likely to benefit from ICI therapy. Significantly, nonclassical monocytes have been shown to accumulate in the myeloma niche and contribute to CAR T-cell dysfunction, further linking monocyte activity to immune evasion.⁴⁵ These findings suggest that targeting monocyte-driven pathways could enhance therapeutic efficacy in MM by overcoming resistance and restoring immune function.

As anticipated, the TIDE algorithm and correlation analysis of immunosuppressive receptor genes suggest that MGPS may serve as a novel predictive biomarker closely associated with immunotherapy response. Moreover, we found that a high MGPS score reflects a state of cellular functional impairment, characterized by enriched hypoxia features, dysfunctional DNA repair mechanisms and a compromised interferon response. These findings provide biological evidence and insights into the adverse prognosis associated with high-risk groups. By identifying high-risk patients with an unfavorable immune landscape, MGPS may aid in refining treatment strategies, such as intensifying immunomodulatory therapies or personalizing immunotherapy regimens.

In the era of precision medicine, the Revised International Staging System (R-ISS),⁴⁶ with its limited parameters, fails to satisfy clinicians’ requirements for an ideal prognostic marker. Existing prediction models for MM often suffer from biases in algorithm selection or lack validation across multiple datasets,^47,48 resulting in suboptimal performance or overfitting. To address these limitations, we developed 101 models using data from multiple MM cohorts and ten commonly used machine learning algorithms. Our results showed that the RSF plus CoxBoost model was the optimal choice through rigorous evaluation. The derived MGPS is identified as an independent prognostic factor and demonstrates outstanding predictive accuracy across various datasets.

The MGPS genes were reported that notably associated with immune responses or tumor progression. PSAP was identified as highly expressed across the entire monocyte population and has previously been reported to be mainly involved in regulating inflammatory and immune responses.⁴⁹ Additionally, IL-32 is known to induce inflammatory cytokines such as TNF-alpha, interleukin-1beta, and interleukin-6 from monocytes/macrophages, and it synergizes with signals from pattern-recognition receptors.⁵⁰ Notably, PECAM1/CD31 functions as a checkpoint molecule modulated by FcγR-mediated signaling in monocytes, making it a potential target to enhance FcγR functions in antibody-mediated therapies.⁵¹ Besides, IL-16 has been found to stimulate the expression and production of pro-inflammatory cytokines in human monocytes.⁵² Although limited studies documented other genes in the differentiation and function of monocytes, they have been implicated in the pathogenesis of various tumors including MM.^53–57 Collectively, these findings underscore the significant impact of MGPS on monocyte function and MM immune status, highlighting its potential role in the clinical management and personalized treatment of MM patients.

While MGPS has been demonstrated to be a robust prognostic predictor for MM, several limitations must be acknowledged. First, training and test cohorts heterogeneity, such as differences in treatment protocols and sample collection, may introduce biases and slightly affect the predictive performance of the MGPS model. The Mendelian randomization data used in this study primarily derive from European populations, which may potentially limit the generalizability of the findings to other ethnic groups. Additionally, the biological functions of the MGPS genes in the context of MM pathogenesis and progression remain to be fully elucidated. Further functional validation in experimental models, including both in vitro genetic knockdown assays and in vivo animal studies, is essential to provide deeper insights into the mechanisms underlying the prognostic significance of MGPS and to confirm its potential as a clinically relevant biomarker. In future studies, we intend to explore MGPS performance in diverse patient populations and validate gene functions through mechanistic in vitro and in vivo experiments. Moreover, large-scale, multicenter prospective studies are also necessary to further confirm MGPS’s clinical utility.

Conclusion

Our findings reveal a complex interplay between immune cells and MM, providing immunological insights into the pathogenesis of MM. Furthermore, based on multi-omics analyses and machine learning algorithms, we established a robust and powerful monocyte-related prognostic signature for prognosis prediction and personalized immunotherapy in MM patients. By identifying high-risk patients with an unfavorable immune landscape, MGPS may help refine treatment strategies, such as intensifying immunomodulatory therapies or personalizing immunotherapy regimens. These results underscore the potential of MGPS as a clinical tool to improve risk stratification and guide therapeutic decisions in MM.