A single-cell atlas characterizes dysregulation of the bone marrow immune microenvironment associated with outcomes in multiple myeloma

Abstract

Multiple myeloma (MM) remains incurable despite advances in treatment options. Although tumor subtypes and specific DNA abnormalities are linked to worse prognosis, the impact of immune dysfunction on disease emergence and/or treatment sensitivity remains unclear. We developed an Immune Atlas of MM by generating profiles of 1,397,272 single cells from the bone marrow (BM) of 337 newly diagnosed participants and characterized immune and hematopoietic cell populations. Cytogenetic risk-based analysis revealed heterogeneous associations with T cells of BM, with 17p13 deletion showing distinct enrichment of a type 1 interferon signature. The disease progression-based analysis revealed the presence of a proinflammatory immune senescence-associated secretory phenotype in rapidly progressing participants. Furthermore, signaling analyses suggested active intercellular communication involving a proliferation-inducing ligand and B cell maturation antigen, potentially promoting tumor growth and survival. Lastly, using independent discovery and validation cohorts, we demonstrated that integrating immune cell signatures with known tumor cytogenetics and individual characteristics significantly improves stratification for the prediction of survival.

Fig. 1. Overview of the Immune Atlas design, workflow and participant characteristics.

a, Overview of the Immune Atlas study design, discovery (n_participants = 263) and validation (n_participants = 74) participant cohorts, sample processing and analysis workflow. b, Clinical characteristics of participants (n_participants = 263) in the discovery cohort. The forest plot illustrates the effect of various clinical features on PFS. Error bars display the 95% confidence interval (CI). c, Dot plot depicting the cross-section of samples based on ASCT and frontline treatment, where the dot size indicates the number of participants and dot color indicates the type of treatment regimen. d, Bar chart showing the total number of participants with each of the six genetic events used for risk classification. e, UpSet plot showing the distribution of the major cytogenetic abnormalities comprising the Davies-based HR myeloma definition. f, UpSet plot showing the intersection of participants categorized as SR or HR and NP or RP at baseline. g, Overview of progression group categorization and study design for the discovery cohort. The participants with a progression event within the first 18 months following therapy were classified as RPs (n_participants = 67). Participants with durable remission or no observed progression for at least 4 years were classified as NPs (n_participants = 83). Participants with a progression event between 18 months and 4 years were classified as RPs (n_participants = 71). The participants who exited the study before 4 years of disease diagnosis without experiencing a progression event were classified as Inc (n = 42). h, Kaplan–Meier curves display the survival analysis for participants categorized on the basis of risk stratification (HR versus SR), transplant as a frontline treatment, treatment type and ISS staging. Participants lacking ISS stage information at baseline or WGS information for cytogenetic risk stratification were omitted from the relevant figure panels. The P values were estimated using a log-rank test. a and g were created with BioRender.com. Source data

Fig. 2. Single-cell Immune Atlas of samples from participants with MM.

a, UMAP embedding of 1,149,344 CD138⁻ BMME cells collected from participants with MM. A total of 106 clusters were observed, spanning five major compartments defined by canonical markers: T and NK cells, B cells and erythroblasts, myeloid cells, erythrocytes and plasma cells. Populations identified as doublets are colored gray. b, Feature plots displaying the normalized gene expression for a selection of lineage-specific markers. UMAPs and per-aliquot cluster compositions to depict the effects of batch correction for major lineages are shown in Extended Data Fig. 2. c, A stacked bar chart displaying the average per-participant cell type composition at baseline. Clusters are colored by their lineage and shaded by subtype. Doublet populations are omitted. d, UMAP of the T lymphocyte and NK compartment. Cells are colored by individual cell type, with shaded boundaries representing regions predominantly containing CD4⁺ (purple), CD8⁺ (orange) or NK (green) cells. The color for specific cell types is included in the corresponding dot plots (f–h). An extended dot plot for precise annotation of different T and NK cell subtypes is shown in Extended Data Fig. 3a. UMAPs for only the CD8⁺ and CD4⁺ T cells are also shown in Extended Data Fig. 3b,c. e, Feature plots displaying the normalized gene expression per cell for markers to distinguish CD4⁺, CD8⁺ and NK cells. f–h, Dot plots displaying the average scaled expression of select marker genes used for precise cluster annotation. Expression is visualized on a red–blue color scale, with the size of each dot corresponding to the percentage expression. Dot plots are split by lineage into NK cells (f), CD8⁺ T cells (g) and CD4⁺ T cells (h). The colored triangle next to the cluster name matches the cluster color in the corresponding UMAP (d). Percent.mt refers to the percentage of mitochondrial transcripts. i, UMAP of the myeloid compartment. Cells are shaded by their subtype. Doublet populations are colored gray. j, Dot plot displaying the average scaled expression of select marker genes for precise cluster annotation in the myeloid compartment. Expression is visualized on a red–blue color scale, with the size of each dot corresponding to the percentage expression. The triangle next to the cluster name matches the cluster color in the corresponding UMAP. k, UMAP of the B cell and erythroblast compartment. Cells are colored by their lineage (B cells, cyan; erythrocytes, red; others, dark blue), shaded by subtype. Doublet populations are colored gray. l, Dot plot displaying the average scaled expression of select marker genes used for precise cluster annotation in the B cell and erythroblast compartment. Expression is visualized on a red–blue color scale, with the size of each dot corresponding to the percentage expression. Source data

Fig. 3. BMME alterations associated with cytogenetic abnormalities.

a, BM immune cell type and subtype abundances comparing HR versus SR participants or stratifying by individual HR abnormalities. Proportions were normalized to the total number of nonmalignant cells per participant. Colors indicate the coefficient of a linear model fitted to logit-transformed proportion adjusting for processing site, with orange and blue indicating cell populations with higher and lower abundances, respectively. Shapes denote two-sided P values (circles, not significant; diamonds, P < 0.05; squares, P < 0.01). b, Expression of marker genes representing CD4⁺ and CD8⁺ T cell states in participants stratified by composite risk (that is, HR versus SR) or individual cytogenetic risk abnormalities. The colors indicate z score normalization, with positive values indicating higher expression levels (red) in participants with the risk event compared to the dataset average, whereas negative values (blue) indicate lower expression levels. c, Box plots show per-participant proportions for CD8⁺ Teff HLA cells as a fraction of total CD8⁺ cells, stratified by combined Davies risk (HR, n_participants = 123; SR, n_participants = 108) or the presence of 1q21 gain (n_participants = 72), NSD2 t(4;14) (n_participants = 12), their combination (n_participants = 18) or neither (n_participants = 127). Two-sided P values were computed using a linear model on logit-transformed proportions adjusting for site. d, Pseudotime trajectory of the CD8⁺ T cells, with arrows indicating the paths along the trajectory originating at CD8⁺ Tn cells. The cell types with high and lower proportions in HR as compared to SR are shown with shades of orange and blue colors, respectively. Proportions are shown as log odds ratios relative to total CD8⁺ T cells. e, Putative dysfunctional T cell signature (Supplementary Table 4) enrichment in CD8 Teff HLA⁺ cells from participants with HR NDMM (n_participants = 123) as compared to SR (n_participants = 108). The significance of the difference in signature levels was determined using the Wilcoxon rank-sum test, two-sided. The dashed red line indicates the median dysfunctional signature score for standard-risk patients. f, IFN-I response signature levels across major cell compartments, grouped by composite or individual HR abnormalities. Per-participant IFN-I response scores across each compartment are in Extended Data Fig. 4g. g, Correlation between IFN-I response signature scores between plasma cells of CD138⁺ (bulk RNA-seq, GSVA) and CD138⁻ (pseudobulk scRNA-seq) compartments. Participants with greater than 50 plasma cells in scRNA-seq and matching bulk data were included (n_participants = 108). Significance of the correlation was calculated using a linear model, adjusting for processing site using a two-sided test. h,i, IFN-I response signature scores derived from bulk RNA-seq, comparing participants with and without 1q21 copy-number alterations (h; normal, n_participants = 408; gain, n_participants= 213; Amp, n_participants = 39) and participants with and without TP53 loss-of-function mutations (i; none, n_participants = 571; partial, n_participants = 55; complete, n_participants = 24). Participants included in the CoMMpass registry with available risk information and CD138⁺ bulk RNA-seq data were analyzed. The significance of differences in enrichment was evaluated using pairwise two-sided Student’s t-tests between groups. j, Dot plot summarizing differential abundance results across cell populations, including plasma cells. The marker shape indicates two-sided P values and color represents the log odds ratios from a linear model fitted to logit-transformed proportions, with positive values denoting enrichment correlated with plasma cell levels. Proportions were computed as a fraction of all cell populations, excluding doublets. Differential abundance was assessed using the plasma cell percentage as both a continuous covariate and a categorical covariate (≥20% versus <20%). k–m, Scatter plots illustrating the relationship between BM plasma cell percentages, as estimated by flow cytometry before CD138 isolation (x axis) and the abundance of CD8⁺ Teff memory cells (k), BM-resident NK cells (l) and fibroblasts (m). n, Comparison of exhaustion-related markers TIGIT and TOX expression across CD8⁺ T cells between participants with less than (blue; n = 189) and greater than (orange; n_participants = 74) 20% plasma cells. In the box plots, bounds of the box represent the 25th and 75th percentiles, with the center displaying the median. Whiskers extend to 1.5× the interquartile range (IQR) beyond the bounds of the box. Source data

Fig. 4. Single-cell level alterations in the BMME of MM RPs.

a, Stacked bar chart of mean per-participant cell type proportions at baseline in RP versus NP participants. Clusters are colored by their major cell type and shaded by individual clusters. The average proportion of major cell types is shown on the graph, normalized as a fraction of all cells, excluding doublets. b, Proportions of T, B and myeloid cells per participant by progression groups (RP, n_participants = 67; NP, n_participants = 83), calculated from total immune cells. Two-sided P values were calculated using a linear model adjusting for processing sites. Nonsignificant P values > 0.05 are not shown. c, CD14⁺ monocyte trajectory with differential abundance results. Arrows indicate lineage paths from immature Neutrophil_RPS/RPL. Circles represent clusters, with labels adjacent to each center and size representing the number of cells within a cluster. Colors correspond to the log odds ratio for RP and NP participants, computed as a fraction of myeloid cells. d, A volcano plot displaying the differentially expressed genes between NP and RP participants in CD14⁺ monocytes. The x axis displays the natural log fold change and the y axis displays the −log₁₀ two-sided BH-adjusted P value. Significantly differentially expressed genes associated with inflammation and IFN-I response pathways are shown in red and blue, respectively. e, A bar plot displaying GSEA for differentially expressed genes shown in d. The x axis shows the NES, with positive values indicating association with NP-enriched markers and negative values indicating association with RP-enriched markers. Pathways are colored by the sign of the NES and shaded by −log₁₀FDR; those with BH-adjusted P values < 0.1 were considered significant. f, Box plots displaying the distribution of CD8⁺ cells (left) and selected significantly enriched T cell subtypes (right) as a fraction of CD3⁺ T cells. Open circles represent individual participants. The difference in proportions between RPs (n_participants = 67) and NPs (n_participants = 83) was assessed using a linear model (two-sided P value). g, A volcano plot displaying the differentially expressed genes between NP and RP participants in CD3⁺ T cells. The x axis displays the natural log fold change and the y axis displays the −log₁₀ two-sided BH-adjusted P value. Select genes are highlighted and colored on the basis of their associated function. h, Bar plot displaying GSEA results for the differentially expressed genes shown in g. The x axis shows the NES, with positive and negative values indicating association with the NP-enriched and RP-enriched markers, respectively. The pathways are colored by the sign of the NES and shaded by the −log₁₀FDR; those with BH-adjusted P values < 0.1 were considered significant. i, The trajectory of CD8⁺ cells along with differential abundance results. Arrows indicate lineage paths from CD8⁺ Tn cells. Circles represent clusters, with labels shown adjacent to each center and size representing the number of cells within a cluster. The lineage highlighted in red corresponds to the trajectory associated with cytotoxic cells. Colors correspond to the log odds ratio for RP and NP participants, computed as a fraction of all CD8⁺ T cells. j, A density plot showing the distribution of cells by pseudotime along the cytotoxicity lineage, from the original cluster (CD8_Tn, low pseudotime) to differentiated cytotoxic clusters (high pseudotime). k, Smoothed normalized expression along the cytotoxicity lineage’s pseudotime for five Tn-associated genes (blue) and five cytotoxicity-associated genes (red), weighted by Slingshot’s lineage assignment weight. l–o, Survival plots displaying the relationship between PFS and the participant’s average cytotoxicity signature (l), Tn signature (m), exhaustion signature (n) or putative T cell dysfunction signature (o) scores across CD3⁺ T cells. Significance values were determined using two-sided P values from a Cox proportional hazards (PH) model. For the survival curves, participants were binned into groups with ‘high’ (brown) or ‘low’ (blue) expression, with the cutoff determined using maximally selected rank statistics. In the box plots, the middle bar represents the median, lower and upper hinges correspond to first and third quartiles and upper whiskers extend to the largest value no further than 1.5× the IQR. Source data

Fig. 5. Pathway and systems biology analysis to decipher mechanisms of poor outcomes in MM.

a, Top: differential cell population abundances stratified by composite cytogenetic risk (HR versus SR), progression (RP versus NP) and progression within the subsets of SR participants (SR: RP versus NP) and HR participants (HR: RP versus NP). The color indicates the linear model coefficient fitted to logit-transformed proportions, with orange indicating higher abundance in the first group of each comparison and blue indicating lower abundance. Shapes indicate the two-sided P value for the coefficient, with circles representing no significant difference, diamonds representing P < 0.05 and squares representing P < 0.01. Proportions were computed as a fraction of all immune cells excluding plasma cells, erythroid cells and doublets. Bottom: average normalized scores for select immune signatures (Supplementary Table 4) across the various cell populations. b, Bar graph of differentially enriched markers within the CD34⁺ HSC population. The log₂ fold change values are relative to RPs, with overexpressed genes in orange and downregulated genes in blue. c, Heatmap of intercellular communication depicting key patterns of outgoing (top) and incoming (bottom) signaling between cell types. All cell types, including plasma cells, were included, although some populations were combined to simplify the interpretability of the cell communication analysis (Supplementary Table 6). Each row corresponds to a ligand–receptor pair. The heatmaps show relative strength of outgoing signals (top; ligand expression) and the corresponding incoming signals (bottom; receptor expression) by each cell type. d, Chord diagram indicating the IFNγ signaling network in all cells. Chords are colored by the ‘sender’ cell type (ligand) and point toward the ‘receiver’ cell type (receptor). e–h, Average expression of IFNγ in T cells (e,f) and IFNγR2 in CD14⁺ monocytes (g,h) and their associations with outcome in SR participants. Box and violin plots (e,g) compare the per-participant average expression between SR-NP and SR-RP participants, with each dot representing a participant. P values were calculated using a two-sided Wilcoxon rank-sum test. In the box plots, the middle bar represents the median, lower and upper hinges correspond to first and third quartiles and upper whiskers extend to the largest value no further than 1.5× the IQR. Kaplan–Meier curves (f,h) display the association between expression level and PFS, stratified by median expression (high, above the median; low, below the median). Hazard ratios and two-sided P values were estimated using CoxPH models. i, Heatmap showing normalized average AUC scores for transcriptional regulons on selected myeloid populations. Additional columns display hazard ratios and two-sided P values, from CoxPH models fitted on average participant AUC scores categorized into high and low activity using a cutpoint approach. j, Survival plots display survival associations between regulon activity in the myeloid compartment and participant outcomes, where high E2F8 regulon expression (bottom) and low IRF7 regulon expression (top) are associated with poor outcomes. The two-sided P values from Cox models are shown. k, Feature plots of the per-cell AUC values for IRF7 (left) and E2F8 (right) TF regulons across key myeloid populations. Blue indicates low activity (or AUC) and red indicates high activity. The color bar represents the regulon enrichment score. Source data

Fig. 6. Terminally differentiated and senescent T cells predict poor outcomes in an independent validation cohort of 74 participants.

a, Characteristics of the validation cohort. Left: scRNA-seq analysis on CD138⁻ BM of 74 participants with NDMM yielded 247,928 high-quality cells. Middle: PFS curves comparing the discovery (gray) and validation (red) cohorts. Dashed lines indicate the median survival time for each cohort. The adjacent box plot indicates the number of RPs and NPs in the validation cohort. Right: number of participants from the validation cohort with different HR abnormalities per the Davies risk definition. prob., probability. b, UMAP embedding of 247,928 CD138⁻ BMME cells from the validation cohort (Fig. 2a). Major cell types are shown in consistent colors with the discovery cohort, with shades representing different cell states and subtypes. c, Correlation of cellular abundances between discovery and validation cohorts. Points represent individual cell types and subtypes, with colors corresponding to the Kendall correlation coefficient. The shaded region represents the 95% CI. d, Top: differential abundance analysis of RPs vs NPs in the validation cohort. Colors indicate the log odds ratio derived from a linear model on logit-transformed proportions, adjusting for the study site, with orange indicating higher abundance in RP and blue indicating lower abundance. Shapes indicate the two-sided P value for the coefficient, with circles representing no significant difference, diamonds representing P < 0.05 and squares representing P < 0.01. Proportions were estimated as a fraction of all immune cells. Bottom: average normalized signature scores for select immune signatures (Supplementary Table 4) across immune populations. e,f, Survival analysis of CD8⁺ Teff HLA⁺ cell abundance, as a fraction of all CD8⁺ cells, for OS (e; P = 0.07, log-rank test) and PFS (f; P = 0.011, log-rank test). g, Volcano plot of the differentially expressed genes across CD3⁺ T cells in the validation cohort between RPs (right) and NPs (left). The x axis shows the batch-corrected log₂ fold change, with positive values corresponding to higher expression in RP participants and negative values corresponding to higher expression in NP participants. The y axis shows the −log₁₀ BH-adjusted P value based on a two-sided test using a linear model fit to log-normalized expression. Vertical dashed lines mark the log₂ fold change ± 0.1 and the horizontal dashed line marks adjusted P = 1 × 10⁻⁵⁰. Certain genes are highlighted on the basis of their functional role (red, cytotoxic or cytolytic; blue, IFN-I; green, Tn cell; yellow, stress). h, Pseudotime trajectory of CD8⁺ T cells, with arrows indicating the paths along the trajectory originating at CD8⁺ Tn cells. Circles represent clusters and colors indicate the log odds ratio of proportion as a fraction of CD8⁺ T cells between RP and NP participants, with orange showing higher abundance in RP and blue showing higher abundance in NP. i,j, Survival plots of PFS associated with the participant’s average cytotoxicity signature score (i) or average dysfunction signature score (j) across all CD3⁺ T cells. CoxPH models were fit on continuous signature scores, with hazard ratios and two-sided P values reported. For the survival curves, participants were binned into groups with high or low expression, with the cutoff determined using maximally selected rank statistics. Source data

Fig. 7. Prediction of MM progression by integration of cytogenetic risk along with immune signatures.

a, Schematic of variables tested (immune signatures, cytogenetics and clinical variables (covariates)) and the three regression strategies used (elastic net, logistic regression and Cox), followed by bootstrap validation used for model selection. b, Receiver operating characteristic curves for progression prediction models based on single clusters, clinical variables and cytogenetics or Immune Atlas variables alone and in combination. Curves are colored by model. The labels indicate subclusters (SubC) and covariates (CoV), including age, batch, site, ISS and cytogenetic. c–e, Kaplan–Meier curves showing the separation of participants by predicted PFS based on age, ISS stage and batch (c), cytogenetics, age, ISS stage and batch (d) and Immune Atlas signatures, age, ISS stage and batch (e). f,g, Kaplan–Meier curves showing the separation of participants when cytogenetic risk scores are combined with the best 11 predictive Immune Atlas subclusters (f) or with all 83 subclusters (g). h, Importance of immune subclusters for predicting the progression. The clusters with better and poor MM outcomes are shown with blue and red colors, respectively. The red dashed line marks the P value threshold of 0.1 from the ANOVA Wald chi-squared test. i, PFS predictive model with 11 predictive immune clusters, excluding ASCT, in the discovery cohort stratifying participants by high versus low risk (AUC = 0.80). j, Validation of the PFS predictive model based on 11 immune clusters and clinical covariates (excluding ASCT) on an independent validation cohort of 74 participants with NDMM. The model demonstrates excellent performance in stratifying participants at higher risk of progression from a low-risk category, achieving an AUC of 0.94. All survival curves display the two-sided P value from a log-rank test. k, Summary of the key cellular subtypes and signaling pathways comprising the MM BMME and their association with participant outcomes. Within the aging BM, a state of chronic inflammation, known as ‘inflammaging’, results in altered lymphoid and myeloid cell populations, enabling immune escape of malignant plasma cells. Within the T cell compartment, participants with MM showing poor outcomes exhibit a shift toward immunosenescent and late-activated CD8⁺ T cells, producing type 2 interferon (IFN-II) that drives senescence-associated and immunosuppressive phenotypes in myeloid compartment. In contrast, participants with MM showing better outcomes display highly proliferative Tn and Tcm CD8⁺ subsets, in addition to enriched Th populations driven by increased MHC-II antigen presentation among myeloid cells. T cell and myeloid populations in these participants appear to be stimulated by IFN-I, in contrast to participants with poor outcomes exhibiting enrichment of IFN-II signaling. This difference in IFN stimulation appears to be linked to participant outcomes, in part, through the differential expression of BAFF by IFN-I-stimulated monocytes and APRIL by IFN-II-stimulated monocytes. Notably, BAFF preferentially binds to mature B cells to promote survival, potentially enhancing B cell-mediated immunity and leading to improved outcomes. Conversely, APRIL has been shown to inhibit memory B cell function and promote malignant plasma cell survival. This dysregulation is further highlighted in the shift from B cell development toward increased myelopoiesis in participants with poor outcomes. k was created with BioRender.com. Source data

Extended Data Fig. 1. Overview of clinical characteristics of the validation cohort of 74 NDMM patients.

(a) Summary of clinical characteristics. The forest plot illustrates the effect of various clinical features on progression-free survival (PFS). Error bars display a 95% confidence interval. (b) Bar chart showing the total number of patients with each of the six cytogenetic events used for risk stratification. (c) UpSet plot showing the distribution and overlap of the major cytogenetic abnormalities comprising the Davies-based high-risk myeloma definition between patients. (d) UpSet plot showing the intersection of patients categorized as standard-risk (SR) or high-risk (HR) and non-progressor (NP) or rapid progressor (RP). (e) Kaplan-Meier curves depict survival outcomes for patients categorized based on risk stratification (HR vs SR), transplant as a frontline treatment, treatment type, and ISS staging. Two-sided p-values from a log-rank test are displayed. Patients lacking ISS stage information at baseline or who lack WGS information for cytogenetic risk stratification were omitted from the survival analysis. Source data

Extended Data Fig. 2. Evaluation of Batch Effects in the scRNA-seq data.

UMAP projections and per-sample cluster compositions across immune and tumor cell compartments, illustrating the effects of batch correction using Harmony approach. Each panel (a-f) corresponds to distinct cellular compartments, including (a) all cells, (b) T and natural killer cells, (c) myeloid cells, (d) B cells, erythroblasts, and progenitors, (e) erythrocytes, and (f) plasma cells, further subdivided into 6 subpanels (x.1-x.6). In panel (a), each cluster is colored by their lineage group, with different shades distinguishing the clusters (see Fig. 2a-c). In panels (b-f), clusters are colored by their sub-compartment colors (see Fig. 2d-i, Extended Data 3). Subpanels (x.1-x.4) display UMAPs before (x.1-x.2) and after (x.3-x.4) Harmony batch correction, colored by cell type (x.1, x.3) with doublets marked in gray or by processing site (x.2, x.4) with individual sample aliquots distinguished by lighter or darker shades of the respective color (Emory: blue, Mayo: red, MSSM: green, WUSTL: purple). (x.5) Stacked bar charts highlighting the sample aliquots contributing to each (a.5) lineage group or (b.5-f.5) cluster, where the size of each segment is proportional to the number of cells coming from a given sample in the respective population and colored as described in panels x.2 and x.4. (x.6) Stacked bar charts displaying the cellular composition within individual aliquots, where the size of each segment is proportional to the cells associated with that cluster as a fraction of all cells in the corresponding compartment, colored by cluster as described in panels x.1 and x.3. Source data

Extended Data Fig. 3. Extended Cell Type Annotation information.

(a) Dot plot displaying additional markers for CD4⁺ T, CD8⁺ T, and natural killer cell annotations. Scaled expression values for each gene are visualized on a red-blue color scale, with the size of each dot representing the percent expression of the corresponding gene. The colored triangle next to each cell name corresponds to the T and NK cell clusters (see Fig. 2d). (b-c) UMAPs displaying subclusters of CD4⁺ T cells (b) and CD8⁺ T cells (c). (d-e) UMAP (d) and dot plot displaying the top differentially expressed markers (e) for mature erythrocyte populations. (f) Feature plot displaying plasma cell markers. UMAP embeddings correspond to those displayed in Fig. 2a. (g-h) UMAP (g) and dot plot displaying the top differentially expressed markers (h) for the plasma cell compartment. (i) Scatter plot showing the relationship between logit-transformed plasma cell proportions in CD138^neg scRNA-seq data and plasma cell fractions estimated via flow cytometry on unsorted aspirates. Each dot represents an individual sample, colored by processing site. The black line with p and R² values represents the line of best fit average across processing sites, where dashed color lines represent fits for individual processing sites. (j) Box plot depicting the plasma cell proportion in the scRNA-seq data for each patient (n_patient=263). Patients are binned based on whether their plasma cell fraction is ≥20% (n_patient=74) or <20% (n_patient=189), as estimated by flow cytometry on unsorted samples. Two-sided p-values comparing the groups is estimated via a linear model. In the box plots, bounds of the box represent the 25^th and 75^th percentile, with the center displaying the median. Whiskers extend to 1.5*IQR beyond the bounds of the box. (k-m) Analyses to assess for malignancy of the plasma cells in the CD138^neg scRNA-seq data. (k) UMAP highlighting cells with driver overall and individual gene mutations in red and inferred copy number in purple. UMAP embeddings correspond to those displayed in Fig. 2a (l) UMAP of all cell types of RP (left) or NP (right) cohort samples showing various driver mutations. (m) CCND1 expression across cells in the RP vs NP patient cohorts (top). CCND1 expression in cells from patients with mutations (mut) and/or translocations (Tx) determined based on analysis of WES and WGS data (bottom). Unadjusted p-values from a wilcoxon rank-sum test between RP and NP samples is displayed if significant. Source data

Extended Data Fig. 4. Summary of Cytogenetic Risk Associated Immune Alterations in Multiple Myeloma.

(a) Survival based on CD8 Teff HLA⁺ cell abundance as a fraction of CD8⁺ T cells (p = 0.036, log-rank test). Cell abundances were batch-corrected using regression residuals. Cut-off was determined by maximally selected rank statistics at the 82% quantile (191 low, 39 high). (b) Progression-free survival analysis from regression of the above described CD8 Teff HLA⁺ metrics (p = 0.062, log-rank test) with cutoff set at the 89% quantile (205 low, 25 high). (c-d) Scaled expression of ‘dysfunctional’ T cell signature genes, including exhaustion and senescence markers (Supplementary Table 4). Dot plot (c) showing percent of cells expressing each gene with dot size and average expression with color (red = high, blue = low) and tile plot (d) showing average expression across cluster. Rows are clusters, and columns are genes, grouped by signature category. (e-f) Trajectory plots depict the predicted differentiation of CD8⁺ T cell subtypes from CD8⁺ naïve T cells. Arrows indicate directionality and dots represent trajectory clusters colored by batch-adjusted log-odds abundance and sized by cell count. Comparisons shown for patients with both t(4;14)[NSD2] and 1q21 gain (e) and TP53 complete loss (f) compared to those without them (orange = high, blue = low). (g) Per patient heatmap of overall IFN-I response signature scores. Red corresponds to high signature score, blue corresponds to a lower signature score. Signature scores are normalized within each cell lineage. Patients with no cells of a specific type will have a grey bar for their IFN-I signature score. Patient tumor cytogenetics are displayed in a title map to the left of the signature score plot. (h) Dot plots of differential cellular abundance analysis for patients with 1q21 gain in combination with other high-risk abnormalities. Each row corresponds to a comparison between 1q21 and other cytogenetic events. “HR_No1q” = high risk without 1q21 gain; “HR+1q21” = high risk with 1q21 gain. Colors indicate log-odds ratios, and shapes indicate two-sided p-values comparing cluster proportions from a linear model (circle = ns, diamond = p < 0.05, square = p < 0.01). (i) Dot plot of differential cellular abundance analysis for patients with partial or complete loss of TP53 via mutation or copy number loss. Partial loss is defined as either monoallelic loss of 17p13 or one non-synonymous mutation of TP53. Complete loss is defined as biallelic loss of 17p13 or monoallelic loss of 17p13 with mutation. (j) Dot plot summarizing differential cellular abundance analysis for patients with Chromothripsis or APOBEC events. (k-m) Box plots illustrating the relationship between bone marrow plasma percentages (>=20%, n_patient = 74; <20%, n_patient = 189), as estimated via flow cytometry before CD138 isolation (x-axis), and the abundance of CD8⁺ T effector memory cells (n_patient=263) (k), BM-resident NK cells (l), and fibroblasts (m). Two-sided p-values for each comparison were computed using a linear model. In the box plots, bounds of the box represent the 25^th and 75^th percentile, with the center displaying the median. Whiskers extend to 1.5*IQR beyond the bounds of the box. Source data

Extended Data Fig. 5. Differential abundance analyses for immune subpopulations comparing rapid progressors versus non-progressors.

(a-d) Summary differential abundance volcano plots and individual cluster box plots, computed as either a fraction of (a) all immune cells, (b) B cells, (c) myeloid cells, or (d) CD3⁺ cells, comparing RP and NP patients. Each panel consists of a volcano plot in the top left, displaying the log-odds ratio (x-axis) derived from a site-adjusted linear model associating progression group with logit-transformed cell proportions, and the y-axis displays the -log10 of the two-sided p-value from the linear model. Points are colored according to cluster identity in Fig. 2. Box plots depict the per-patient proportion for each cluster within each compartment, with individual patients shown as open circles (Rapid Progressors (RP, n_patient = 67)=Orange, Non-Progressors (NP, n_patient = 83)= Blue). Two-sided p-values are displayed above the box plot if the populations are significantly different between RP and NP (p < 0.05) as assessed using a linear model. In the box plots, bounds of the box represent the 25^th and 75^th percentile, with the center displaying the median. Whiskers extend to 1.5*IQR beyond the bounds of the box. (e) Comparison of differential abundance profiles between triplet therapy patients and the overall cohort. The x and y axes display log-odds ratio change in proportion between RP (positive) and NP (negative) patients across triplet therapy patients (x-axis) or all patients (y-axis) as a fraction of immune cells. Source data

Extended Data Fig. 6. Selected gene regulators and their survival associations in myeloid populations.

Kaplan-Meier curves (top) and feature plots displaying the per-cell AUC scores (bottom) in selected myeloid populations for STAT1 (a), BRCA1 (b), E2F1 (c), E2F8 (d), IRF1 (e), IRF2 (f), IRF7 (g), IRF9 (h), and KLF9 (i) regulons. Gene regulator analysis and regulon identification were performed using pySCENIC on selected myeloid clusters. Kaplan-Meier curves are based on the average AUC value across the selected myeloid populations. The cut-point approach was implemented to stratify patients into either high or low regulon expression groups. Two-sided p-values indicate the significance of a CoxPH model fitted to the regulon AUC score. Source data

Extended Data Fig. 7. UMAP and dot plot comparisons of validation and discovery cohorts across cell compartments.

UMAP projections and stacked bar charts comparing cellular composition between discovery and validation cohorts across major immune cell compartments. Each panel (a-d) corresponds to a distinct cellular compartment including (a) all cells, (b) NK/T cells, (c) myeloid cells and (d) B cells, erythroblasts and progenitor cells, further subdivided into 5 subpanels (x.1-x.5). Subpanels (x.1-x.2) display UMAPS of both discovery and validation cells colored by cell type (x.1) or sample (x.2). Subpanels (x.3-x.4) show the same UMAPS separated by discovery (x.3) and validation (x.4) cohorts, colored by cell type. Subpanels (x.5) contain stacked bar charts representing the distribution of cells within each cluster across individual samples and separated by cohort. Samples from the discovery cohort are colored different shades of blue, while samples from the validation cohort are colored different shades of red. (e-i) Validation cohort dot plot showing the previously described cluster markers for (e) CD8⁺ cells, (f) CD4⁺ cells, (g) natural killer cells, (h) myeloid cells, and (i) B cells, erythroblasts, and progenitor cells. Source data

Extended Data Fig. 8. Additional Model Diagnostics and Predictive Modeling of Overall Survival and Progression in MM Using Immune Signatures, Cytogenetics, and Clinical Features.

Panels (a-c) show model AUCs in the discovery cohort across immune compartments: (a) AUCs by immune cell compartments and variable combinations, (b) Boxplots showing the distribution of AUCs for models integrating a single cell type with various covariates, along with the AUC for the top models integrating all clinical covariates and either 7, 11, 34, or all immune populations, and (c) Boxplots showing the AUCs for models derived from individual immune populations, grouped by cellular compartment. In the box plots, bounds of the box represent the 25^th and 75^th percentile, with the center displaying the median. Whiskers extend to 1.5*IQR beyond the bounds of the box. Whiskers extend to 1.5*IQR beyond the bounds of the box. (d-i) Receiver operating characteristic (ROC) and Kaplan-Meier (KM) analysis for overall survival (OS) prediction in the discovery cohort. (d) ROC curves for models with single immune subclusters (SubC), clinical covariates (CoV), cytogenetics, and combinations. Covariates include age, batch, site, ISS stage and cytogenetics. KM curves depict predicted OS based on (e) clinical covariates (f) cytogenetics + clinical covariates, (g) Immune Atlas Signature + clinical covariates, and (h) the top 20 predictive immune subclusters + clinical covariates. (i) The importance of immune subclusters for predicting the OS colored by favorable (blue) or poor (red) OS association. (j-o) ASCT’s contribution to survival prediction in discovery and validation cohorts. KM curves showing ASCT association with OS in discovery (j) and validation (k) cohorts. Predictive models including ASCT, cytogenetics, clinical variables, and top immune features in discovery (l) and validation (m). Equivalent models excluding ASCT as a variable, still stratifying high- versus low-risk patients in discovery (n) and validation (o) cohorts. (p-t) OS prediction in the validation cohort. (p) A forest plot based on a multivariate CoxPH illustrating the bias of ASCT and ISS for OS. Two-sided p-values from the CoxPH model are displayed. (q) Box plot of bootstrapped AUCs for models using various immune compartments (n_patient=71), with the integrative model as a superior option. However, the AUC for OS remained below 0.75 in general. (r) Box plot of bootstrapping applied to an integrative model of feature selection based on AUC to identify the optimal model for OS prediction using immune populations, clinical information, cytogenetics, and ASCT. In the box plots, bounds of the box represent the 25^th and 75^th percentile, with the center displaying the median. Whiskers extend to 1.5*IQR beyond the bounds of the box. (s) KM curves show the effect of ASCT on the prediction of OS. (t) Integrative model using the top 20 immune signatures combined with clinical and cytogenetic information in the prediction of OS. Source data

Extended Data Fig. 9. Prediction of MM progression by integration of IMWG cytogenetics risk along with immune signatures.

(a) Diagram illustrating the type of variables that were tested (immune signatures, IMWG cytogenetics, and clinical covariates) followed by the three regression strategies used (elastic net, logistic regression, and Cox) with bootstrap validation for model selection. (b) Receiver operating characteristic (ROC) curves for progression prediction models based on single clusters, clinical variables, and cytogenetics or immune atlas variables alone and in combination are shown and colored based on the specific group of models. The labels include subclusters (SubC) and covariates (CoV), which include age, batch, site, ISS, and cytogenetics. Kaplan-Meier curves showing the separation of patients with high or low scores for prediction of PFS are shown for (c) demographics-based, (d) IMWG 2024 high-risk criteria, and (e) Immune Atlas signatures. Kaplan-Meier curves show the separation of patients when cytogenetic risk scores are combined with the (f) best 10 predictive immune atlas subclusters or (g) 57 subclusters. (h) Volcano plot displaying the importance of immune subclusters for predicting the progression based on the coefficients from the modeling. The clusters with better and poor MM outcomes are shown with blue and red colors, respectively. Receiver operating characteristic (ROC) curves for using 10 subclusters and excluding ASCT for discovery (i) and validation (j) cohorts. (k) Dot plot showing the average scaled expression of marker genes for the “MatNeut2” phenotype reported in Jong et al., 2024. Expression is visualized on a red-blue color scale, with the size of each dot corresponding to the percent expression of marker genes. Expression is normalized relative to the average cluster expression across all CD14⁺ monocyte clusters. Source data