MULTI-seq: sample multiplexing for single-cell RNA sequencing using lipid-tagged indices

Abstract

Sample multiplexing facilitates scRNA-seq by reducing costs and identifying artifacts such as cell doublets. However, universal and scalable sample barcoding strategies have not been described. We therefore developed MULTI-seq: multiplexing using lipid-tagged indices for single-cell and single-nucleus RNA sequencing. MULTI-seq reagents can barcode any cell type or nucleus from any species with an accessible plasma membrane. The method involves minimal sample processing, thereby preserving cell viability and endogenous gene expression patterns. When cells are classified into sample groups using MULTI-seq barcode abundances, data quality is improved through doublet identification and recovery of cells with low RNA content that would otherwise be discarded by standard quality-control workflows. We use MULTI-seq to track the dynamics of T-cell activation, perform a 96-plex perturbation experiment with primary human mammary epithelial cells and multiplex cryopreserved tumors and metastatic sites isolated from a patient-derived xenograft mouse model of triple-negative breast cancer.

Figure 1:. MULTI-seq demultiplexes cell types, culture conditions, and time points for single-cell and single-nucleus RNA sequencing.

(A) Diagram of the anchor/co-anchor LMO and CMO scaffolds (black) with hybridized sample barcode oligonucleotide (red). LMOs and CMOs are distinguished by their unique lipophilic moieties (e.g., lignoceric acid, palmitic acid, or cholesterol). (B) Schematic overview of a proof-of-concept single-cell RNA sequencing experiment using MULTI-seq. Three samples (HEKs and HMECs with and without TGF-β stimulation) were barcoded with either LMOs or CMOs and sequenced alongside unlabeled controls. Cells were pooled together prior to scRNA-seq. Next-generation sequencing produces two UMI count matrices corresponding to gene expression and barcode abundances. (C) Cell type annotations for LMO-labeled cells demonstrate separation between HEKs (pink), MEPs (cyan), and LEPs (dark teal) in gene expression space (see Fig. S2A). Ambiguous cells positive for multiple marker genes are displayed in grey. n = 6,186 MULTI-seq barcoded cells. (D) MULTI-seq sample classifications for HEKs (dark red), unstimulated HMECs (green), and TGF-β-stimulated HMECs (blue) match cell state annotations. Cells classified as doublets (black) predominantly overlap with ambiguously-annotated cells. n = 6,186 MULTI-seq barcoded cells. (E) TGF-β-stimulated HMECs (blue) exhibited elevated TGFBI expression relative to unstimulated HMECs (green). *** = Wilcoxon rank sum test (two-sided), p <= 10⁻¹⁶. n = 1,950 MULTI-seq barcoded HMECs. Data are represented as mean ± SEM. (F) Single nucleus MULTI-seq sample classification proportions for each cell type identified by clustering in gene expression space (see Fig. S2E–G). n = 5,894 MULTI-seq barcoded nuclei. (G) MULTI-seq sample classifications illuminate temporal gene expression patterns in Jurkat cells following activation with ionomycin and PMA for varying amounts of time. Time-point centroids in gene expression space are denoted with larger circles. n = 3,709 Jurkat nuclei. (H) Violin plots of gene expression marking different stages of Jurkat cell activation. n = 3,709 Jurkat nuclei.

Figure 2:. MULTI-seq barcoding of multiplexed HMEC culture conditions

(A) Barcode UMI abundances (left) and doublet classifications (right) mapped onto barcode space. MULTI-seq barcode #3 is used as a representative example. Doublets localize to the peripheries of sample groups in large-scale sample multiplexing experiments. n = 25,166 cells. (B) Cell state annotations demonstrate separation between MEPs (cyan) and LEPs (dark teal) in gene expression space (left, see Fig. S5A). Ambiguous cells positive for multiple marker genes are displayed in grey. MULTI-seq classifications grouped by culture composition (right) — e.g., LEP-alone (blue), MEP-alone (green), and both cell types together (dark red) — match cell state annotations. Discordant region where annotated MEPs are classified as doublets by MULTI-seq is indicated with arrows. n = 25,166 cells. (C) MULTI-seq doublet classifications (left) and computational predictions produced by DoubletFinder (right) largely overlap in gene expression space. Discordant region where DoubletFinder-defined doublets that are classified as singlets by MULTI-seq indicated with arrows. n = 25,166 cells. (D) MEP co-culture induces LEP proliferation and TGF-β signaling. Clusters corresponding to resting (black) and proliferative (green) LEPs are identifiable in gene expression space (Fig. S5B). Projecting sample classification densities onto gene expression space for co-cultured LEPs (dark red, top left) and LEPs cultured alone (blue, top right) illustrates that co-cultured LEPs are enriched in the proliferative state (table, bottom left). Co-cultured LEPs also express more TGFBI than LEPs cultured alone. Each point represents an average of LEPs grouped according growth factor condition. *** = Wilcoxon rank sum test (two-sided), p = 3.1×10⁻⁶. n = 32 signaling molecule condition groups. Data are represented as mean ± SEM. (E) Hierarchical clustering and heat map analysis of resting LEPs grouped by treatment. Emphasized genes are known EGFR signaling targets. RNA UMI abundances are scaled from 0–1 for each gene. Values correspond to the average expression within each signaling molecule treatment group. Dendrogram labels: E = EGF, W = WNT4, A = AREG, I = IGF-1, R = RANKL, C = Control.

Figure 3:. PDX sample multiplexing demonstrates low-RNA cell detection, reveals immune cell proportional shifts and classical monocyte heterogeneity in the progressively metastatic lung.

(A) Schematic overview of PDX experiment. (B) MULTI-seq sample classifications (WT, early, mid, late tumor progression) mapped onto barcode space. Replicate tissues are denoted as ‘A’ or ‘B’. n = 10,427 cells. (C) MULTI-seq classifications facilitate low-RNA and low-quality cell deconvolution. CellRanger discards cells barcodes with low RNA UMI counts (red dotted line). Gene expression profiles for classified low-RNA cells reflect established immune cell types (top right, see Fig. S6F). Unclassified low-RNA cells resemble low-quality single-cell transcriptomes (bottom right, see Table S4). n = 2,580 (classified), 583 (unclassified) cells. (D) Cell state annotations (top) and tumor stages (bottom) for lung immune cells in gene expression space. Mono. = monocyte, C = classical, NC = non-classical, Mac. = macrophage, DC = dendritic cell, pDC = plasmacytoid DC. Cells with undeterminable annotations displayed in grey. n = 5,965 cells. (E) Statistically-significant shifts in lung immune cell type proportions for each tumor stage relative to WT. Two-proportion z-test with Bonferroni multiple comparisons adjustment, * = 0.05 > p > 10⁻¹⁰; ** = 10⁻¹⁰ > p > 10⁻²⁰; *** = p < 10⁻²⁰. n = 44 tumor-stage/cell type groups. Statistically-insignificant proportional shifts omitted. (F) Subsetted classical monocyte gene expression space overlaid with sample classification densities corresponding to tumor stage. Inset illustrates heterogeneity within late-stage classical monocytes characterized by differential expression of Thbs1 and Cd14. n = 2,496 (all), 1,087 (inset) cells.