Patient demographics

The clinical characteristics of our cohort are described in Table 1 and Supplementary Table S1. For our NDMM patients, the mean follow-up time was 35.1 months, with a range of 1.23 to 88.4 months. Pharmacotherapy used by patients at the time of diagnosis, as well as the chemotherapies used over the course of follow-up for NDMM patients are reported in Supplementary Tables S2 and S3. Patients received various regimens that included alkylating agents, glucocorticoids, proteasome inhibitors, immunomodulators, anti-CD38 monoclonal antibodies, and autologous stem cell transplantation (ASCT). The most common regimens over the follow-up period were CyBorD (Cyclophosphamide, Bortezomib, and Dexamethasone, N = 41); RD (Lenalidomide and Dexamethasone, N = 30); Lenalidomide alone (N = 23); and DRD (Daratumumab, Lenalidomide, and Dexamethasone, N = 16) (Supplementary Table S3). The most common first-line regimens were CyBorD (N = 38), RD (N = 16), DRD (N = 6), VRD (Bortezomib, Lenalidomide and Dexamethasone, N = 9) (Table 1). Thirty patients received an ASCT during the follow-up period (Table 1).

Table 1 Summary of clinical characteristicsMinimal batch effect

Experimental samples demonstrate mixing between plates in PCA space (Supplementary Figure S1A). As well, control samples tightly clustered by type rather than plate in PCA space (Supplementary Figure S1B). Technical replicates of control samples across plates demonstrated consistency in PG counts with a coefficient of variation (CV) of 3.34% to 15.6%, and PG intensity with a CV of 2.99% to 3.91% (Supplementary Figure S1C,D). PVCA revealed that duration of sample storage, plate on which samples were processed, and the two-way interaction of these variables accounted for approximately 11% of observed variance, with ~ 89% of the variance remaining residual after accounting for these covariates (Supplementary Figure S1E).

PG count similar among plasma cell dyscrasias

We identified a total of 11,128 PGs across our patient cohort, with an average of 8976 PGs per sample (Supplementary Figure S2A; Supplementary Table S4). Out of these, 4194 PGs were common to every patient, and 7849 PGs were common to at least 75% of patients (Supplementary Figure S2B). No association was found between diagnosis and the number of PGs identified (Supplementary Figure S2A).

Plasma cell-associated proteins in BMIF show greater intensity in BMIF of NDMM compared to SMM and MGUS patients

In our comparison of NDMM to MGUS patients, we identified 28 PGs with differential intensity (Supplementary Table S5). Of these 28 PGs, seven (SH2B2, IRF4, SDC1, TNFRSF17, ISG20, CXCR4, PRKCZ) were significantly more intense in NDMM than in MGUS, while 21 were significantly more intense in MGUS. (Fig. 1A). When comparing NDMM to SMM patients, 6 PGs showed differential intensity, with only one (SDC1) being more intense in NDMM (Fig. 1B; Supplementary Table S6). Many of the PGs that were more intense in NDMM compared to MGUS or SMM were known plasma cell markers, such as IRF4, SDC1, and TNFRSF17. The intensities of these proteins were significantly correlated with bone marrow plasma cell burden (Supplementary Figure S3).

Fig. 1

Plasma cell associated proteins are more intense in MM BMIF than MGUS or SMM BMIF: (A and B) Differential intensity analysis of PGs between MM and MGUS (A), and MM and SMM (B). The analysis was conducted on proteins with less than 25% missingness, PG intensities were median normalized within samples, and t-tests were used to compare mean PG intensity between groups. Correction for multiple testing was performed using FDR

Risk associated BMIF PGs are enriched for coagulation

To identify prognostic signals in our proteomic data, we analyzed the association between overall survival (OS) and the intensity of each PG in NDMM samples. Initially, considering only PGs common to all NDMM samples, we identified 56 PGs associated with increased OS and 34 PGs associated with decreased OS (Fig. 2A; Supplementary Table S7). When we expanded our analysis to include PGs with up to 25% missingness, using BPCA to estimate PG intensity, we identified 85 PGs associated with increased survival and 109 PGs associated with reduced OS (Fig. 2B; Supplementary Table S8). These PGs spanned the intensity dynamic range (Supplementary Figure S4). Clustering patients based on the Z-scored intensity of the OS-associated PGs from the complete dataset revealed three groups with significantly different OS and progression-free survival (PFS) (Fig. 2C-E). Differential intensity analysis comprising the 15 longest- and 15 shortest-surviving NDMM patients identified 346 PGs associated with longer OS and 128 PGs associated with shorter OS (Supplementary Figure S5, Supplemental Table S9).

Fig. 2

BMIF PGs define high-, medium-, and low-risk MM groups: (A and B) Cox hazard ratios for OS were calculated on each PG without missingness (A) and on each PG with less than 25% missingness (B). Correction for multiple testing was by FDR. C) Heatmap of PGs without missingness that were significantly associated with OS. Intensities were Z-scored, and clustering was performed on Euclidian distance using the Ward’s method. (D and E) Kaplan Meier for both OS (D) and PFS (E). Trace colours correspond to the clusters in (C)

Among the PGs with significantly low hazard (HR < 0.5, p < 0.05) from both the complete and imputed datasets, as well as the PGs associated with long OS from the differential intensity analysis, coagulation (GO:0007596, WP558) was the most significantly enriched Gene Ontology term/Wikipathway (Fig. 3A-B; Supplementary Figure S6A-D; Supplemental Tables S10-21). Considering the PGs with significantly high hazard (HR > 2, p < 0.05), a mix of terms were identified, with less consistency between the analyses that included “chromosome segregation”, “ribosome biogenesis”, and “mRNA processing”, among others (Supplemental Tables S10-21). The proteins within the "Blood Coagulation" Gene Ontology term that were significantly associated with OS included coagulation factors 2, 7, 9, and 12, as well as KLKB1, SERPINA10, SERPIND1, SERPINC1, PROC, PROZ, and HGFAC. Given the prominence of the coagulation cascade in these results, we clustered our NDMM patients based on proteins within this cascade (F2, F3, F5, F7, F8, F9, F10, F11, F12, F13A1, F13B, VWF, KLKB1, KLK1, FGA, FGB, FGG) as indicated on WikiPathways (WP272). This identified three distinct groups characterized by globally low (group 1), medium (group 2), and high (group 3) intensity of coagulation factors, with significantly different OS (log-rank p = 0.00078) and PFS (log-rank p = 0.017) (Fig. 4A-C). The group with the worst outcome (group 1, N = 26) was characterized by low intensity of all coagulation factors and had a median survival of 16.2 months. Conversely, the group with the best outcome (group 3, N = 24) had a high intensity of all factors and a median survival of 67.2 months.

Fig. 3

Overall survival-associated PGs in BMIF are enriched for coagulation: Gene set enrichment (A) and Wikipathway analysis (B) on the PGs with hazard ratios < 0.5 and p < 0.05 after correction for multiple testing

Fig. 4

Low-intensity of coagulation cascade PGs in BMIF identified high-risk MM: (A) Heatmap of PGs within the coagulation cascade. Intensities were Z-scored, and clustering was performed on Euclidian distance using the Ward’s method. (B and C) Kaplan Meier for both OS (B) and PFS (C). Trace colors correspond to the clusters in (A)

In PCA space, considering all PGs, coagulation group assignments clustered together, and a comparator group consisting of MGUS (N = 8), MGRS (N = 1), and POEMS syndrome (N = 1), largely projected amongst coagulation group 3 (Supplementary Figure S7A). This comparator group was selected due to its low plasma cell burden and relatively normal bone marrow histology. We performed differential intensity analysis comparing each of coagulation groups 1 and 3 to the comparator group. Comparing coagulation group 1 to this group identified 190 PGs significantly more intense in MM coagulation group 1 and 295 PGs significantly more intense in the comparator group (Supplementary Figure S7B; Supplementary Table S22). The PGs that were significantly more intense in MM coagulation group 1 included all PGs found to be more intense in NDMM compared to MGUS (Fig. 1A, Supplementary Figure S7B). Comparing MM coagulation group 3 to this group identified 102 PGs significantly more intense in MM coagulation group 3 and 359 PGs significantly more intense in the comparator group (Supplementary Figure S7C; Supplementary Table S23). Of the PGs significantly more intense in the whole NDMM cohort compared to MGUS, all but CXCR4 were significantly more intense in coagulation group 3 MM (Supplementary Figure S7C). Overall, 481 PGs were differentially intense between coagulation group 1 MM and the comparator group, while 465 PGs were differentially intense between coagulation group 3 MM and the comparator group. Notably, coagulation factors 2, 5, 7, 8, 9, 10, 11, 12, 13B, FGA, FGB, FGG, and KLKB1 were significantly less intense in coagulation group 1 compared to the comparator group (Supplementary Figure S7D). Conversely, the only coagulation cascade PG significantly differentially intense in coagulation group 3 MM compared to the comparator group was von Willebrand factor (VWF), which was more intense in coagulation group 3 (Log2FC = 0.99) (Supplementary Figure S7D). Additionally, several anticoagulant proteins were also significantly more intense in MM coagulation group 3 compared to MM coagulation group 1 or the comparator group, including PROC, PROS1, PROZ, SERPINC1, SERPINA6, SERPINCA10, SERPINI1, and SERPIND1 (Supplementary Figure S8). Finally, coagulation-regulating extracellular matrix proteins (ECM) FN1 and TNC were also significantly less intense in MM coagulation group 1 compared to MM coagulation group 3 (Supplementary Figure S9).

Coagulation biomarkers in BMIF are independent of prognostic indicators

To evaluate whether the prognostic value of coagulation clusters is independent of current clinical evaluations, we first analyzed the hazard associated with coagulation cluster assignment using a multivariate Cox model. This model accounted for R-ISS stage, age, bone marrow plasma cell burden, monoclonal (M)-protein quantity, and platelet count. The coagulation cluster assignment remained significant (HR = 0.39, p = 0.002) and outperformed the R-ISS staging (Fig. 5A). Next, we assessed for associations between coagulation clusters and various clinical and biochemical parameters, including age at diagnosis, plasma levels of albumin, beta-2 microglobulin, creatinine and eGFR, INR, PTT, M-protein, the ratio of involved to uninvolved serum free light chains, bone marrow plasma cell burden, and platelet count; no significant associations were found (p > 0.05) (Fig. 5B). Finally, we applied the coagulation clustering to patients staged by R-ISS and identified high- and low-risk groups within R-ISS stage 2 patients (n = 54, log-rank p = 0.021). There were too few patients in RISS stages 1 (n = 4) and 3 (n = 9) for meaningful analysis (Fig. 5C).

Fig. 5

Coagulation groups in BMIF prognostically independent of current clinical workup: (A) Forrest plot of the hazard ratio, confidence interval, and significance of coagulation clusters assessed in a multivariate Cox hazard model. (B) Box plots of continuous clinical covariates showing their distribution in each coagulation group. Wilcoxon was used to compare means between groups, and p values were corrected using FDR. (C) Kaplan Meier for OS with coagulation clustering applied to each RISS group

We also investigated potential associations between coagulation cluster assignment and the use of anticoagulants, antiplatelets, antiglycemics, antihypertensives, or antilipidemics at the time of diagnosis. No significant associations were observed (Supplementary Figure S10A-D). Finally, we assessed for associations between coagulation group and the use of chemotherapeutic classes and ASCTs over the follow-up period; no associations were found (Figure S10E).