Single cell characterization of myeloma and its precursor conditions reveals transcriptional signatures of early tumorigenesis

Abstract

Multiple myeloma is a plasma cell malignancy almost always preceded by precursor conditions, but low tumor burden of these early stages has hindered the study of their molecular programs through bulk sequencing technologies. Here, we generate and analyze single cell RNA-sequencing of plasma cells from 26 patients at varying disease stages and 9 healthy donors. In silico dissection and comparison of normal and transformed plasma cells from the same bone marrow biopsy enables discovery of patient-specific transcriptional changes. Using Non-Negative Matrix Factorization, we discover 15 gene expression signatures which represent transcriptional modules relevant to myeloma biology, and identify a signature that is uniformly lost in abnormal cells across disease stages. Finally, we demonstrate that tumors contain heterogeneous subpopulations expressing distinct transcriptional patterns. Our findings characterize transcriptomic alterations present at the earliest stages of myeloma, providing insight into the molecular underpinnings of disease initiation.

Fig. 1. The landscape of normal and abnormal plasma cells at single cell resolution.

a Overview of cohort and experimental setup, including the number of samples per disease stage, sex, sample preparation batch, and whether the sample was fresh or stored frozen prior to 10x sequencing. CD138+ bone marrow cell fractions were isolated and are analyzed in this study. b, c UMAP representation of plasma cells colored by disease stage (b) and sample ID (c). Cells similar in expression profile are placed nearby in this embedding. d Results of Leiden clustering of all cells. Seven clusters were merged to define a single cluster of healthy plasma cells. e Sample composition of Leiden clusters, by disease stage and sample ID (colors match the legends given in (b) and (c), respectively). The majority of clusters each consist of cells from a single sample. f Violin plots showing distribution of expression of genes commonly upregulated in patients with translocations (y-axis), along with annotations of the cytogenetic alterations detected in samples by clinical iFISH assay (top).

Fig. 2. In silico dissection of transcriptional differences in normal and abnormal plasma cells within patient samples.

a The number of cells (top) and estimated purity of each sample with 95% confidence intervals (bottom). Sample purity was estimated using two orthogonal methods: clustering of individual samples (blue; the fraction of cells labeled abnormal per sample is plotted) and our Bayesian hierarchical purity model (orange; the mode of posterior sample purity is plotted). Source data are provided as a Source Data File. b UMAP localization of individual cells labeled normal or abnormal. c Cartoon schematic of our differential expression analysis. We run two DE analyses: First, we compare all abnormal (purple) vs. all normal (yellow) cells using limma-voom. Next, we compare patients’ abnormal cells to their own normal cells, controlling for inter-patient variability. Samples with 100% normal or abnormal cells were excluded from the within-patient analysis. d Volcano plot of limma-voom DE results for abnormal vs. normal cell populations. Orange denotes genes with q-value < 0.1. The 4 most significantly up- and downregulated genes and other selected genes are annotated. e Pseudobulk expression of DEGs detected between abnormal and normal pseudosamples using limma-voom (z-scored per gene). Each column represents the normal or abnormal cells from a given sample. Color annotations denote disease stage (top), normal or abnormal (second), paired columns coming from the same sample (third; matching colors denote that columns correspond to the same sample; black denotes that there was no paired sample), and whether IgH translocation or hyperdiploidy was detected in that sample by iFISH (bottom). f Quantification of DEGs uniquely discovered using within-patient DE. The venn diagram represents the overlap of DEGs found using limma-voom and our within-patient DE approach. The bar plot describes the number of DEGs found per sample using within-patient DE (right side) and the number of abnormal and normal cells per sample (left side). g Volcano plot of 1760 DEGs uniquely discovered using our within-patient DE approach. The y-axis represents the maximum -log₁₀(q-value) of the gene across samples included in the within-patient analysis, and x-axis represents the maximum log₂(fold change). The color and size of a dot denote the number of samples for which that DEG was detected, with blue dots representing DEGs detected in just one sample.

Fig. 3. Bayesian non-negative matrix factorization uncovers gene signatures which capture myeloma cell biology across disease stages.

a Top genes for nine representative gene signatures. The importance score, plotted on the x-axis, is based on both the strength of the gene’s contribution to the signature and its specificity to the signature (see Methods). b A signature with top contribution from CCND1 is discovered and is most active in samples with t(11;14), as expected. c, d We discover a ‘normal plasma cell signature’ that is active in normal plasma cells across disease stages and downregulated in abnormal cells from MM and precursor conditions. We visualize this signature’s activity by showing its mean activity ± s.e.m. for the normal and abnormal populations within each sample (c) and on a UMAP plot (log scale) (d). Mean activities were compared between groups, with *** denoting q < 0.001 for group differences (abnormal cells from SMM (n = 12) and MM (n = 8), respectively, significantly differed from NBM (n = 9)). e Validation on external dataset: our NMF algorithm run on external CD138+ single cell data from MGUS, SMM, MM and healthy donors independently discovers a gene signature similar to our normal plasma cell signature, with shared top genes CD27, CD79A, and JSRP1. f After labeling cells in that dataset as normal or abnormal, we discover that this signature follows the same pattern as in our data, with high activity in normal cells and a significant decrease in activity in abnormal cells across disease stages. Mean activities ± s.e.m. across cells in normal and abnormal portions of samples are shown, with *** denoting q < 0.001 for group differences (abnormal cells from SMM (n = 5) and MM (n = 13), respectively, significantly differed from NBM (n = 11)). Source data are provided as a Source Data File.

Fig. 4. IFN-inducible signature is correlated between CD138+ and microenvironment cells, and gene signatures exhibit intratumor heterogeneity.

a Mean activity ± s.e.m. of CD138+ IFN-inducible signature across normal and abnormal plasma cell populations. Both normal and abnormal plasma cells exhibit significantly increased activity of the interferon-inducible signature in MM vs. NBM (q = 5.2 × 10⁻³ and q = 3.2 × 10⁻⁴, resp.). Source data are provided as a Source Data File. b Mean activity per sample of IFN-inducible signature discovered in CD138+ cells (top), T cells (middle) and CD14+ monocytes (bottom). Mean expression levels for the ten genes with the highest values in the W matrix for each signature are also shown. Expression of additional interferon-inducible genes IFI27 and IFI6 is shown for CD138+ samples (see Supplementary Note 2). CD138+ samples from patients were limited to abnormal cells before calculating means. NMF signature results and expression data for T cells and monocytes were taken from Zavidij et al.. c Subpopulations within patient tumors heterogeneously express gene signatures. Cells from a given MM sample were projected onto a UMAP plot based on expression of highly variable genes, and colored by the activity level of NMF signatures determined to be heterogeneously expressed in that sample (see Methods).