Scalable single-cell RNA sequencing from full transcripts with Smart-seq3xpress

Abstract

Current single-cell RNA sequencing (scRNA-seq) methods with high cellular throughputs sacrifice full-transcript coverage and often sensitivity. Here we describe Smart-seq3xpress, which miniaturizes and streamlines the Smart-seq3 protocol to substantially reduce reagent use and increase cellular throughput. Smart-seq3xpress analysis of peripheral blood mononuclear cells resulted in a granular atlas complete with common and rare cell types. Compared with droplet-based single-cell RNA sequencing that sequences RNA ends, the additional full-transcript coverage revealed cell-type-associated isoform variation.

Fig. 1. Scalable full-transcript coverage scRNA-seq with Smart-seq3xpress.

a, Schematic of nanoliter cDNA synthesis reactions performed in wells of 384-well PCR plates with 3 µl of hydrophobic overlay. b, Illustration of reduced-volume experiments with the lysis, RT and PCR volumes used. c, The number of genes detected per HEK293TF cell at each reaction volume, when sampling 100,000 sequencing reads (n = 100, 19, 32 and 28 cells, respectively). P value represents a two-sided t-test between the 10-µl and 1-µl conditions. d, Influence of hydrophobic overlays on miniaturized cDNA synthesis (1 µl total volume). For each compound, boxes depict the number of genes detected per HEK293FT cell (n = 17, 34, 39, 28, 25, 24, 28, 38 and 70, respectively), subsampled at 200,000 sequencing reads per cell. e, Replacement of the bead-based cDNA cleanup by dilution in single HEK293FT (n = 58 and 52, respectively) cells. Box plots show the number of genes detected per cell and condition (at 100,000 reads) with P value for a two-sided t-test across conditions. f, Tagmentation complexity using 0.1 µl of ATM Tn5 enzyme per HEK293FT cell in relation to input cDNA. The median number of detected genes as a function of raw sequencing reads (n = 51, 53, 54, 53, 53 and 52 cells for 25, 50, 75, 100, 200 and 500 pg, respectively). g, Tagmentation complexity for varying amounts of cDNA input. Complexity was summarized as unique aligned and gene-assigned UMI-containing read pairs per 400,000 raw reads and HEK293FT cell (n = 49, 51, 51, 50, 51 and 44). h, Schematic outline of the Smart-seq3 and Smartseq3xpress workflows. i, The number of genes detected with Smart-seq3xpress after variable amounts of pre-amplification PCR cycles. Median number of genes is reported as a function of raw sequencing reads in HEK293FT cells (n = 93, 98, 108, 113, 102, 114 and 118 cells for 10, 12, 13, 14, 15, 16 or 20 cycles, respectively). j, Fraction of UMI-containing reads to internal reads for HEK293FT cells prepared with Smartseq3xpress (KAPA HiFi; 12 PCR cycles), at a variable range of TDE1 Tn5 amounts (n = 64 cells each). k, Fraction of UMI-containing reads to internal reads for HEK293FT cells prepared with Smartseq3xpress (SeqAmp; 12 PCR cycles), at a variable range of TDE1 Tn5 amounts (n = 60 cells each). l,m, Optimization of RT and PCR conditions across 376 experimental conditions on HEK293FT cells. Colors indicate particular experimental conditions: Smart-seq3xpress with Smart-seq3 TSO (purple; n = 912), 52 °C RT/alternate TSO implementation (yellow; n = 74), fixed spacer TSO variant (blue; n = 45), FLASH-seq TSO variant (green; n = 55), Smart-seq3xpress with improved TSO (pink; n = 63) and all other conditions (gray; n = 21,707). Scatter plots denote the level of artifactual TSO-UMI reads and RNA counting errors (l) as well as a percentage of ribosomal RNA (rRNA) mapped reads and number of detected genes in 100,000 reads after removal of strand invasion reads (m). n, Benchmarking of Smart-seq3 variants. Box plots show the number of genes detected per HEK293FT cell in full-volume Smart-seq3 (ref. ), low-volume Smart-seq3 and Smart-seq3xpress implementations, at the indicated read depths (n = 109–110, 18–27, 9–170, 20–55 and 9–63 cells, depending on the cells available at the given sequencing depths). The box plots (in c, d, e, j, k and n) show the median and first and third quartiles as a box, and the whiskers indicate the most extreme data points within 1.5 lengths of the box. cSt, centistoke.

Fig. 2. Application of Smart-seq3xpress to hPBMCs.

a, Dimensional reduction (UMAP) of 26,260 hPBMC transcriptomes produced with Smartseq3xpress (KAPA, four donors; SeqAmp with improved TSO, three donors) colored and annotated by cell type. EM, effector memory; CM, central memory; NK, natural killer; ILC, innate lymphoid cell; HSPC, hematopoietic stem and progenitor cell; MAIT, mucosal-associated invariant T cell. b, Smartseq3xpress-based TCR reconstruction (TRaCeR) overlayed onto UMAP. c, QC of TCR reconstructions obtained with Scirpy, enumerating the number of T cells with certain types of TCR reconstructions. d, Benchmarking of Smart-seq2, Smart-seq3 and Smart-seq3xpress (SeqAmp, improved TSO) in primary hPBMCs. Each cell was downsampled to 100,000 reads, and the number of detected genes from exon-mapping reads is shown for representative cell types: B cells (n = 73, 366 and 859, respectively), CD4⁺ T cells (n = 261, 1,270 and 1,847, respectively), CD8⁺ T cells (n = 76, 272 and 913, respectively) and NK cells (n = 73, 352 and 601, respectively). P values indicate results of two-sided t-tests between the Smart-seq3 and Smart-seq3xpress. e, Differential gene expression analysis (Wilcoxon test, P_adj < 0.01) between naive CD4⁻ T cell cluster (n = 2,476) and clonal CD4⁻ T cell cluster (n = 682). Indicated are the top five TCR genes driving the clonal CD4⁻ T cell cluster separation. f, Dot plot showing expression of selected marker genes for MAIT, gamma-delta and clonal CD4⁺/CD8⁺ T cells in all annotated clusters, with size of the dot denoting the detection of a gene within the cells of the cluster and color scale denoting the average expression level. g, Analysis of captured transcribed genetic variation in donor-matched Smart-seq3xpress and 10x Genomics 3′ version 3.1 data. For each cell passing QC (n = 2,938 and 9,846, respectively), the number of SNPs with alternate allele coverage per cell are indicated (left) as well as the average SNP coverage normalized by the sequencing depth (right). h, Frequency of RNA-velocity-informative fully spliced reads in donor-matched Smart-seq3xpress and 10x Genomics 3′ version 3.1 data. For each cell in representative cell types—B cells (n = 404 and 642, respectively), CD4⁺ T cells (n = 1,320 and 1,317, respectively), CD8⁺ T cells (n = 417 and 1,181, respectively) and NK cells (n = 441 and 498, respectively)—we summarized the percentage of reads spanning exon–exon junctions, with nominal P values for a two-sided t-test. i, Differential splicing isoform analysis using BRIE2. Shown is a volcano plot of tested skipped-exon events color-coded by significant variation when testing for any cell type (ELBO gain >20). ELBO gain is a surrogate for the Bayes factor (that is, the likelihood ratio of two hypotheses). The x axis denotes the effect size on the distribution of PSI values in a cell type. j, Overlays of color-coded PSI values inferred for each sequenced cell (n = 26,260) in genes with significant cell type splicing variation (PTPRC, GUSBP11 and ISG20). The box plots (in d, g and h) show the median and first and third quartiles as a box, and the whiskers show the most extreme data points within 1.5 lengths of the box. HSPC, hematopoietic stem and progenitor cell.

Extended Data Fig. 1. Optimization of low-volume cDNA synthesis.

(a) Reduction of reaction volumes in single K562 cells. Shown are the number of genes detected per cell at each reaction volume when sampling 100,000 sequencing reads (n = 63, 39, 55, 53 cells, respectively). P-value represents a t test between the 10 µL and 1 µL conditions. (b) Coefficient of variance of cells at scaled volumes in both HEK293FT (n = 100, 19, 32, 28 cells, respectively) and K562 cell lines (n = 63, 39, 55, 53 cells, respectively). (c-d) For both (c) HEK293FT and (d) K562 cells the influence of creating standard volume Smart-seq3 libraries without (HEK293FT n = 37; K562 n = 26) and with an overlay (VaporLock; HEK293FT n = 15; K562 n = 31) was compared. Boxplots show genes detected and p-values show result of a two-sided t-test. (e) Replacement of the bead-based cDNA cleanup by dilution in single K562 (n = 57, 38, respectively) cells. Shown are the number of genes detected per cell for each condition at 100,000 reads with a p-value for a two-sided t-test within cell types. (f) Coefficient of variance between cells having received cDNA clean-up or dilution in HEK293FT (n = 58, 52, respectively) and K562 (n = 57, 38, respectively). (g) For indicated ratios of PCR volume to RT volume, boxplots show genes detected, genes with UMIs detected and UMIs captured, downsampled by sequenced reads. For PCR ratios 0.1x (n = 43), 0.2x (n = 48), 0.3x (n = 48), 0.4x (n = 47), 0.5x (n = 47), 0.6x (n = 45) HEK293FT cells were analyzed. (h) Boxplots for each PCR extension time (2, 3, 4, 6 min) using KAPA HiFi Hot Start polymerase are shown as detected genes binned by their transcript length at 250,000 reads per cell. (n = 29, 24, 29, 10; respectively) (i) Number of genes detected binned by transcript length for each extension time 3 min (n = 64), 4 min (n = 64), 6 min (n = 60) using SeqAmp polymerase at 250,000 reads per cell. The boxplots (in a-i) show the median, first and third quartiles as a box, and the whiskers show the most extreme data points within 1.5 lengths of the box.

Extended Data Fig. 2. Investigating tagmentation complexity.

(a) Systematic investigation of tagmentation complexity was performed by varying cDNA input with constant Tn5 amount, varying input and Tn5 amounts, varying Tn5 amounts with constant cDNA input, scaling of reaction volumes and Tn5 amounts with constant cDNA input, scaling of reaction volumes with constant cDNA and Tn5 amounts. For each boxplot, shown is the library complexity in terms of unique gene-assigned UMI-anchored read pairs (unique per-molecule cut-sites) from 400,000 raw sequencing reads. Each condition contains between 22 and 73 HEK293FT cells as annotated above each box. (b) Concordance of gene expression levels between HEK293FT cells tagmented using 0.05 µL ATM Tn5 (n = 23 cells) and 0.1 µL ATM Tn5 (n = 11 cells) (mean UMI counts over 15,915 genes). (c) Tagmentation complexity using 100 pg cDNA per single HEK293FT cell in relation to the enzyme amount (ATM Tn5). For each dot, the median number of detected genes is calculated from the indicated number of raw sequencing reads and plotted from n = 55, 54, 56, 58, 52, 51 cells (0.025 µL, 0.05 µL, 0.075 µL, 0.1 µL, 0.2 µL, 0.5 µL). (d) For varying amounts of Tn5 enzyme (see (c)), tagmentation complexity was summarized as unique aligned and gene-assigned UMI-containing read-pairs per 400,000 raw reads per HEK293FT cell (n = 53, 46, 56, 57, 52, 50, respectively). (e) The use of in-house Tn5 relative to performance of ATM Tn5 was compared in HEK293FT cells. Each boxplot shows the number of genes detected at 500,000 raw reads (n = 182, 52, 24, 226, 41, 53, 50, respectively). (f) Influence of post PCR clean-up for Smartseq3xpress. Genes detected after treating preamplified cDNA libraries to reduce “contaminants”, such as dNTPs and oligonucleotides, before going into tagmentation. All conditions (dilution alone (n = 152, 129), ExoSAP IT-express (“ExoSAP”; n = 254, 178), or ExoI + Fast-AP (n = 259, 167)) were either diluted in 9 or 19 µL of water. (g) Influence of post PCR clean-up for UMI captures. Pairwise comparisons of mean expression estimates (mean UMIs per gene) for evaluated clean-up conditions. The boxplots shown in a, d, e and f indicate the median, first and third quartiles as a box, and the whiskers show the most extreme data point within 1.5 lengths of the box.

Extended Data Fig. 3. Performance of pre-amplification polymerases for Smartseq3xpress.

(a) Mapping statistics of the six different tested commercial polymerases: KAPA HiFi Hot Start (n = 384), NEBNext Q5 (“NEBNext”, n = 384), NEBNext Ultra II Q5 (“Q5UltraII”, n = 384), Platinum II (n = 384), SuperFi II (n = 384), SeqAmp (n = 360). Dotplots show percentage of mapped reads to exons, introns, intergenic regions and unmapped with the sequenced read depth. (b) Number of Genes and UMIs detected per cell relative to sequencing depth for all six polymerases tested. (c) Fraction of UMI-containing reads versus internal reads for each of the tested polymerases. (d) Lines show median number of genes detected and median number of unique tagmentation sites after downsampling of sequenced reads using 0.5 µM, 1.0 µM and 2.0 µM TSO concentration in RT, combined with either 0.5 µM or 1 µM of each forward and reverse PCR primer. The replicate numbers for each of the panels are: TSO = 0.5 µM, PCR = 0.5 µM: n = 58, 57, 62, 63, 52, 55 cells; TSO = 0.5 µM, PCR = 1.0 µM: n = 51, 60, 61, 45, 51, 55 cells; TSO = 1.0 µM, PCR = 0.5 µM: n = 57, 58, 61, 52, 51, 57 cells; TSO = 1.0 µM, PCR = 1.0 µM: n = 51, 56, 63, 49, 52, 55 cells; TSO = 2.0 µM, PCR = 0.5 µM: n = 62, 59, 61, 53, 54, 51 cells; TSO = 2.0 µM, PCR = 1.0 µM: n = 49, 53, 59, 50, 55, 49 cells (KAPA, PlatinumII, SeqAmp, NEBNext, Q5UltraII, SuperFiII, respectively). (e) Molecular spike-ins were used to assess the accuracy of each polymerase in mRNA molecule counting, based on the counting difference between internal molecular spikes counts and Smart-seq3xpress UMIs at indicated amounts of TSO and PCR primer concentrations. Colored lines indicate the mean counting difference between the unique spike identifiers and quantified UMIs when sampling the spike at the indicated mean expression levels for each of the polymerases, shaded by the standard deviation. The counting differences are expressed in absolute deviance or relative to the mean molecule number. (f) Rate of base conversions in aligned reads relative to the reference genome. For every polymerase, we compute the average fraction of transitions and transversions, shown as boxplots over all cells, n = 328, 367, 315, 343, 302, 322 (for KAPA, PlatinumII, SeqAmp, NEBNext, Q5UltraII, SuperFiII). Boxplots indicate the median, first and third quartiles as a box, and the whiskers show the most extreme data point within 1.5 lengths of the box.

Extended Data Fig. 4. TSO strand-invasion artifact and improved TSO designs.

(a) Schematic representation of the search procedure for artifactual TSO-UMI reads. For every aligned 5’ read, a 20 bp window of reference genome sequence upstream of the alignment start position is considered. Within this sequence, we search for the UMI sequence allowing up to 1 mismatch. (b) Strand invasion leads to shortened captured molecules. We grouped reads by the presence or absence of TSO-UMI match in the 20 bp upstream window and retrieved the closest annotated transcription start site (TSS) of the assigned gene. Shown is the distance to TSS as the cumulative percentage of reads analyzed. (c) Influence of TSO design (x-axis), oligo-dT primer amount (columns) and forward/reverse PCR primer amounts (rows) on the occurrence of strand-invasion artifacts (y-axis). TSO concentration is indicated in color (red: 0.75 uM, teal: 1 µM). Replication was performed over many cells per condition as annotated above each respective box. (d) Shown are all evaluated UMI sequences incorporated into the Smart-seq3 TSO with their base composition in terms of random or stable bases. For each TSO, we show sequencing metrics in HEK293FT cells (numbers of cells per condition annotated in the Figure), in terms of the frequency of artificial TSO priming, Number of genes detected after discarding TSO primed molecules (100,000 raw reads) and the accuracy of the UMI counting as assayed by molecular spikes. Every box is colored by the coding capacity of the associated random bases in the TSO. Boxplots in c and d indicate the median, first and third quartiles as a box, and the whiskers show the most extreme data point within 1.5 lengths of the box.

Extended Data Fig. 5. Evaluation of candidate reaction conditions in human PBMCs and HEK293FT cells.

(a) Candidate TSO sequences and previous TSO (“Smartseq3 Original”) were evaluated in human PBMC samples with a change to 0.125 µM oligo-dT primer. At indicated sequencing depths, we show the median number of detected genes (left), median number of detected genes after discarding TSO-priming artifact reads (middle), and median number of detected UMIs after discarding TSO-priming artifact reads (right). (b) Investigation of the optimal amount of oligo-dT primer (colors) in the context of new TSO (“Smartseq3xpress improved”; left; n = 74, 71, 74, 83 cells for 0.0625 µM, 0.125 µM, 0.25 µM, 0.5 µM, respectively) and previous TSO (“Smartseq3 Original”; right; n = 293, 286, 313, 319 cells for 0.0625 µM, 0.125 µM, 0.25 µM, 0.5 µM, respectively) in PBMCs. Shown are the median number of detected genes per cell after discarding TSO-priming artifact reads. (c) Frequency of TSO artifact reads in PBMC cells for new TSO (“Smartseq3xpress Improved”; left; n = 64, 63, 53, 58 cells, respectively) and previous TSO (“Smartseq3 Original”; right; n = 266, 256, 262, 251 cells, respectively), with colors denoting the amount of oligo-dT primer used. (d) Benchmarking of new Smart-seq3xpress Improved TSO in HEK293FT cells. At the indicated sequencing depth, we show the number of genes in internal+UMI reads (left) and TSO-artifact filtered UMI reads (middle). (e) The reduction in the occurrence of TSO primed strand-invasion artifacts is shown as a boxplot for n = 94, 330 HEK293FT cells (improved TSO, original TSO). (f) Benchmarking of Smartseq3 variants and Smart-seq3xpress iterations. Shown are the number of UMIs captured in HEK293FT cells in the full-volume Smart-seq3 (n = 110 cells), low-volume Smart-seq3 (n = 27 cells) and Smart-seq3xpress iterations (n = 170, 55, 63 cells for KAPA, SeqAmp and SeqAmp improved TSO, respectively) at the indicated read depths. (g) Reproducibility over cells visualized by cell-to-cell correlation for Smartseq3 and Smartseq3xpress (n = 107, 62, respectively). Two-sided t-test p-value < 2.2e-16. (h) Representative correlation in captured molecules between a Smartseq3 cells and a Smartseq3xpress cell. Boxplots in c, e, f and g show the median, first and third quartiles as a box, and the whiskers show the most extreme data point within 1.5 lengths of the box.

Extended Data Fig. 6. Experimental improvements to the Smart-seq3xpress workflow.

(a) Introduction of a 3D-printed adapter to facilitate pooling of final libraries by centrifugation. Pictures showing reservoir and 3D printed frame/holder and assembly to facilitate pooling of 384 well plates quickly via gentle centrifugation. (b) For Smart-seq3 and Smart-seq3xpress (42 °C and 52 °C reverse transcription implementations), we show the workflow in the library preparation process with the associated timings. Any steps that require hands-on work are shaded with pattern.

Extended Data Fig. 7. Analysis of PBMC samples using Smart-seq3xpress.

(a) Distribution of cell type abundances as a percentage of all cells from each of the 7 donors. (b) FACS-based (index sorting) data overlayed onto UMAP embeddings. “Top Left” shows annotated UMAP based on gene expression. “Top middle”, shows classification of all sorted live cells based on their expression of CD4 and CD8. (DnT = double negative CD4- and CD8-, DpT = Double positive CD4 + and CD8 + , NA = cells without staining). “Top Left” all recorded, and index sorted live cells are categorized based on their expression of CD45RA and CCR7 overlayed over UMAP embeddings (CM = CD45RA- CCR7 + , EM = CD45RA- CCR7-, TEMRA = CD45RA + CCR7-, TN = CD45RA + CCR7 + , NA = cells without stainings). “Lower panels” show protein levels for antibody-staining against CD45RA, CD45RO and CCR7 overlayed over UMAP embeddings of sequenced transcriptomes as log-scaled mean fluorescence intensity (mFI).

Extended Data Fig. 8. Downsampling of Smart-seq3xpress PBMC data.

Sequencing data generated for 26,260 human PBMCs was downsampled from a median depth of 258,000 read pairs per cell to steps of 75% (~193,000 reads), 50% (~129,000 reads), 25% (~64,000 reads), 10% (~26,000) of reads while retaining the relative abundances of per-cell coverage. (a) At each of the downsampling depths, we repeated the Seurat workflow of normalization, clustering and dimensionality reduction using UMAP. (b) We tracked the assigned cluster identity of every cell in the dataset over the downsampling depths and visualized the flow of clustered cells by connecting cluster labels (nodes) with line widths scaled to the number of cells.

Extended Data Fig. 9. Direct comparison of Smart-seq3xpress and 10x Genomics 3’ v3.1 data.

Human PBMCs from the same donor were processed and sequenced using both protocols. We sampled the same number of sequenced reads (approx. median 140,000 per cell) and cells (3,000) and applied an identical data analysis workflow using Seurat (see Methods). (a) Shown are the UMAP embeddings for Smartseq3xpress identifying the number of clusters (left), annotated cluster names (middle) and cluster identities by reference-based identification via Azimuth cell type predictions (level 2 annotations) (b) 10xv3.1 UMAP embeddings showing clusters found (left), annotated cluster names (middle) and reference based predicted annotations by Azimuth (prediction level 2). (c-d) Shows the full the donor matched data without downsampling or cell subsetting after QC for both Smartseq3xpress (n = 3187) and 10xv3.1 (n = 6483). (left) clusters detected, (middle) annotated cell types based on clusters, (right) reference-based azimuth cell type predictions (prediction level 2). (e-f) Reference-based annotation of all human PBMC generated in this study using Smartseq3xpress by Azimuth (Hao et al., 2020) using the human PBMC reference dataset available at https://azimuth.hubmapconsortium.org. (e) Smartseq3xpress UMAP is colored and annotated by Azimuth cell type annotation level 2. (f) Smartseq3xpress UMAP is colored and annotated by Azimuth cell type annotation level 3.

Extended Data Fig. 10. PBMC splicing analysis using Smart-seq3xpress.

(a) Percent-spliced-in (PSI) for skipped exon events in the CD27-AS1, ETHE1 and POP5 genes are overlayed as color scale on the UMAP embedding of sequenced PBMC cells (n = 26,260). (b) For genes with significant cell-type specific splicing patterns (shown in Figures. 3j, Supplementary Fig. 18a), scaled expression levels (log-normalized read counts) are annotated as colors.