Overloading And unpacKing (OAK) - droplet-based combinatorial indexing for ultra-high throughput single-cell multiomic profiling

Abstract

Multiomic profiling of single cells by sequencing is a powerful technique for investigating cellular diversity. Existing droplet-based microfluidic methods produce many cell-free droplets, underutilizing bead barcodes and reagents. Combinatorial indexing on microplates is more efficient for barcoding but labor-intensive. Here we present Overloading And unpacKing (OAK), which uses a droplet-based barcoding system for initial compartmentalization followed by a second aliquoting round to achieve combinatorial indexing. We demonstrate OAK's versatility with single-cell RNA sequencing as well as paired single-nucleus RNA sequencing and accessible chromatin profiling. We further showcase OAK's performance on complex samples, including differentiated bronchial epithelial cells and primary retinal tissue. Finally, we examine transcriptomic responses of over 400,000 melanoma cells to a RAF inhibitor, belvarafenib, discovering a rare resistant cell population (0.12%). OAK's ultra-high throughput, broad compatibility, high sensitivity, and simplified procedures make it a powerful tool for large-scale molecular analysis, even for rare cells.

Fig. 1. Principle and performance of OAK in single cell profiling of multiple molecular modalities.

a Schematic of OAK’s scRNA-Seq workflow. mRNA hybridizes with poly-dT oligos within droplets in fixed cells or nuclei. Following reverse transcription and emulsion break, cells/nuclei are pooled and re-distributed into aliquots for secondary indexing via PCR. TSO: template switch oligo. pR1: primer binding sequence for TrueSeq Read 1. Schematic created in BioRender. Darmanis, S. (2023) BioRender.com/w63n572b Simulation of droplets containing zero (blue), one (magenta), and more than one (yellow) cell, at various cell numbers loaded per channel. Pink and green indicate cell number ranges in regular Chromium and OAK, respectively. c OAK results at different cells per channel. Image scale bars: 75 µm. Mean cells per droplet with 95% confidence margins of error are presented. d-e, Number of genes in K562 (d) and NIH/3T3 cells (e) vs. reads per cell. Green: 150,000 cells loaded; yellow: 450,000 cells loaded, same as in c. f-g, Number of genes in K562 (f) and NIH/3T3 cells (g). ~15000 reads per cell. Boxplots: center lines are medians; limits denote Q1 (lower) and Q3 (higher) quartiles; whiskers extend to 1.5 times the interquartile range (IQR) or last data points if within limits. K562: regular Chromium NextGEM 3’ RNA-Seq, n = 6022; OAK scRNA-Seq with 150,000 cells loaded, n = 3647; scifi-RNA-seq, n = 1617. NIH/3T3: OAK scRNA-Seq with 150,000 cells loaded, n = 691; sci-RNA-seq, SPLiT-seq, sci-CAR, and Paired-seq, n = 868 each. h Percentage of human bronchial epithelial cells assigned to each sample hashtag (n = 9) by standard Chromium and OAK. Each dot is a sample hashtag. i Percentage of fragments overlapping TSS. Chromium: Chromium’s standard Multiome ATAC + Gene Expression, n = 4484; OAK_FA: OAK’s multiome with formaldehyde (FA), n = 1835; OAK_MeOH: OAK’smultiome with methanol (MeOH), n = 2903. Boxplots: center lines are medians; limits denote Q1 and Q3; whiskers extend to 1.5 times IQR or last data points if within limits. j–k, Number of genes (j) and ATAC fragments (k) using OAK with FA or MeOH, or Chromium’s standard multiome. Source data are provided as a Source Data file.

Fig. 2. OAK paired snRNA-Seq and snATAC-Seq on the human peripheral retina.

a Uniform Manifold Approximation and Projection (UMAP) of snRNA-Seq data with table of number and percentage of each cell type. The color of the dot in the table indicates the position in the UMAP. b Heatmap displaying OCRs in each cell type. c Chromatin tracks in major cell types for the genomic region spanning the ARR3 gene with a Ridge plot (expression values are normalized and log transformed) indicating gene expression of ARR3 from snRNA-Seq data. d Chromatin tracks in bipolar cell types for the genomic region including the TSS of DOK5, with a Ridge plot indicating DOK5 expression level. e Significant transcription factors by weighted gene activity for each major cell type. Differential activity was determined using a one-sided t-test adjusted for multiple comparisons. Source data are provided as a Source Data file.

Fig. 3. OAK single-cell lineage tracing and transcriptome profiling for melanoma cells during belvarafenib treatment.

a Diagram of lineage tracing experiment. IPC-298 cells labeled with lineage barcodes were sampled for scRNA-Seq on Days 0, 10, 20, and 90. Belvarafenib treatment commenced following Day 0 subculture collection. Schematic created in BioRender. Darmanis, S. (2023) BioRender.com/l09z998b Fold change in cell count for each lineage at each time point. Cell counts from Day 0 served as the baseline. Enriched (yellow) includes lineages with over tenfold increase from Day 0 to Day 20. Resistant (enriched) refers to the lineage enriched on Day 20 and resistant on Day 90. Stable refers to lineages neither depleted nor enriched. c Volcano plot depicting differentially expressed genes on Day 20 between depleted and enriched lineages. P values are calculated by the Wilcoxon rank-sum method (two-sided) with the benjamini-hochberg correction method. Genes with adjusted p values lower than 1e-8 and log2 fold changes beyond ±0.5 are labeled. d Violin plots for FN1 expression level (normalized and log-transformed) in cells within depleted and enriched lineages. e Fold changes between the depleted and the enriched lineages, with specific genes labeled the same as in (d). Green dashed lines denote ±1.5-fold changes. f PROGENy pathway scores for each cell at Day0 (n = 42), Day 10 (n = 59), Day 20 (n = 275), and Day 90 (n = 4827) within the resistant lineage. P values are calculated using the Mann-Whitney-Wilcoxon test (two-sided) with Bonferroni adjustment. Boxplots’ center lines represent medians. Boxplots: center lines are medians; limits denote Q1 and Q3; whiskers extend to 1.5 times IQR or last data points if within limits. g De-differentiation and differentiation scores for cells within the resistant lineage. Source data are provided as a Source Data file.