Bioinformatics analysis of intrinsic drivers of immune dysregulation in multiple myeloma to elucidate immune phenotypes and discover prognostic gene signatures

Abstract

Multiple myeloma (MM) progression is driven by immune dysregulation within the tumor microenvironment (TME). However, myeloma-intrinsic mechanisms underlying immune dysfunction remain poorly defined, and current immunotherapies show limited efficacy. Using RNA-seq data from 859 MM patients (MMRF-CoMMpass), we integrated xCELL, CIBERSORT, and ESTIMATE algorithms to deconvolute immune-stromal dynamics. Consensus clustering identified immune subtypes, followed by differential gene analysis and LASSO-Cox regression to construct a prognostic model validated in an independent cohort (GSE19784, N = 328). Immune Subtype Classification: Two subgroups emerged: Multiple myeloma-associated immune-related cluster 1 (N = 482): Immune-dysfunctional TME with Th2 cell enrichment, preadipocyte accumulation, and CXCL family suppression, linked to poor survival (P < 0.001). Multiple myeloma-associated immune-related cluster 2 (N = 377): Immune-active TME with cytotoxic CD8 + T/NK cell infiltration and favorable outcomes. Prognostic Gene Signature: Ten immune-related genes (UBE2T, E2F2, EXO1, SH2D2A, DRP2, WNT9A, SHROOM3, TMC8, CDCA7, and GPR132) predicted survival (The One-year AUC = 0.682 and The Over 5-years AUC = 0.714). We define a myeloma-intrinsic immune classification system and a 10-gene prognostic index, offering a framework for risk-stratified immunotherapy. Integration with flow cytometry could optimize precision treatment in MM.

Keywords: Immune dysfunction; Immunosuppress; Multiple myeloma; Prognostic signature; Tumor microenvironment.

Fig. 1

Identification of two immune subtypes (MAIC1/MAIC2) in multiple myeloma training cohort. (A) Consensus clustering matrix for K = 2 subgroups. (B) Kaplan–Meier survival curves showing significantly worse outcomes in MAIC1 (brown, N = 482) versus MAIC2 (blue, N = 377) (P = 0.0066, log-rank test). (C) Integrated heatmap of clinical features (MM stage, gender, race) and immune cell composition. Black boxes highlight cell types with significant inter-subtype differences (***P < 0.001, log-rank test).

Fig. 2

Immune cell infiltration patterns in MM training cohort. (A) Stacked bar plot of 22 immune cell proportions across 859 patients (CIBERSORT). (B) Correlation matrix of immune cells (ssGSEA). Red/blue gradients indicate positive/negative correlations (range: − 1 to 1). (C) Differential immune cell abundance between MAIC1 and MAIC2 (Wilcoxon test). Significance levels: *P < 0.05, **P < 0.01, ***P < 0.001.

Fig. 3

Immune microenvironment characterization. (A) Violin plots comparing immune/stromal scores (ESTIMATE algorithm) between subgroups. (B) Heatmap of 28 immune infiltration features (ssGSEA). Red/blue indicates high/low activity (z-score range: − 2 to 2).

Fig. 4

Differentially expressed immune-related genes (DE-IRGs). (A) Volcano plot demonstrating the differentially expressed immune-related genes (DE-IRGs) between the two MAIC subgroups, identifying genes with significant upregulation (log2 fold change > 1 and p-value < 0.05). The plot represents the magnitude and significance of gene expression changes, with the size of the dots indicating the magnitude of differential expression. (B) KEGG pathway enrichment analysis of MAIC1-upregulated genes, identifying significant biological pathways associated with the tumor microenvironment and immune response. The analysis highlights functional categories enriched for genes regulated in MAIC1 subgroup. (C) KEGG pathway enrichment analysis of MAIC2-upregulated genes, with bubble size indicating the number of genes contributing to each pathway. (D) Gene Ontology (GO) enrichment analysis for MAIC1-upregulated genes, focusing on biological processes, molecular functions, and cellular components. (E) Detailed GO enrichment results for MAIC1-upregulated genes, showing specific pathways and their relevance to multiple myeloma progression. (F) Gene Set Enrichment Analysis (GSEA) results for MAIC2 group, identifying predicted biological pathways enriched in immune-active state. (G) GSEA results for MAIC1 group, identifying biological pathways enriched in immune-dysfunctional state.

Fig. 5

Prognostic model based on 10 DE-IRGs. (A) GO terms of prognostic DE-IRGs. (B) LASSO coefficient selection (λ = 0.021). (C) Univariable Cox regression results. (D) Co-expression network of 10 genes. (E and F) Survival and ROC curves in training/validation cohorts.

Fig. 6

Expression levels of 10 prognostic genes. Box plots comparing MAIC1 (red) vs. MAIC2 (blue) (P < 0.001 for all genes).