Efficient profiling of total RNA in single cells with STORM-seq

Abstract

Despite significant advances, current single-cell RNA sequencing (scRNA-seq) technologies often struggle with accurately detecting non-coding transcripts, achieving full-length RNA coverage, and/or resolving transcript-level complexity. Many are also difficult to implement or inaccessible without specialized liquid handlers, further limiting their utility. We present Single-cell TOtal RNA-seq Miniaturized (STORM-seq), a random-hexamer primed, ribo-reduced single-cell total RNA sequencing (sc-total-RNA-seq) protocol using standard laboratory equipment. Adapted as a kit, STORM-seq constructs sequence-ready libraries in one working day, producing the highest complexity scRNA-seq libraries to-date, robustly measuring transcript isoforms and clinically relevant gene fusions in single cells. STORM-seq faithfully reconstructs expression profiles of locus-level transposable elements (TEs), and provides high-resolution profiling of transient, low-abundance enhancer RNAs (eRNAs), offering a powerful tool to dissect single-cell gene regulatory networks in unprecedented detail. Applied to human fallopian tube epithelium, the improved transcriptional resolution reveals a putative progenitor-like population and intermediate cell states, shaped by TEs and non-coding RNAs.

Figure 1.. STORM-seq efficiently profiles total RNA in single-cells.

a) Overview of the STORM-seq library preparation protocol. b) Fragment structure and annotations of final STORM-seq libraries. c) STORM-seq is the fastest single-cell total RNA-seq method, similar in total time to Smart-seq3xpress. Hands on time: active handling/pipetting of cells/samples; Hands off time: samples/cells not being actively handled (e.g. PCR steps). d) STORM-seq has the highest mapping rates in single-cells compared to VASA-seq and Smart-seq-total. e) STORM-seq has full gene-body length coverage, similar to other plate-based approaches. Protein coding genes only shown for comparison purposes. f) STORM-seq identifies thousands of genes and transcripts per cell across read depths. g) Gene detection rate comparison demonstrates STORM-seq measures the most genes/cell compared to other methods (UMI count minimum of 1). h) Background genomic alignment percentages in single cells (100k reads/cell) demonstrates that STORM-seq aligns similar proportions of reads to known coding, non-coding, and intergenic TE space, similar to bulk total RNA-seq. i) Example genomic background alignment coverage across technologies shows regions with coverage are flanked by poly-A sequences. j) Sequence logo plots of the most abundant UMI per cell across technologies. Expected results are even representation for random UMIs for STORM (NNNNNN - 8 bp UMI), VASA-seq (NNNNNN - 6 bp UMI/UFI), Smart-seq-total (16xN - 16 bp UMI), and Smart-seq3xpress (NNNNNNNNWW - 8bp random + 2 bp W (A/T)). All technology comparisons and results shown are in HEK293T cells at subsampled sequencing depths shown.

Figure 2.. STORM-seq reconstructs cell-type specific regulatory elements and clinically relevant gene fusions in single-cells.

a) Single-cell transposable element (TE) expression representation (obs/exp) across single-cell technologies and bulk RNA-seq. Bulk total RNA-seq is considered the “gold standard”. STORM-seq reconstructs TE profiles in single cells similar to bulk total RNA-seq. b) STORM-seq quantifies locus-level TE expression similar to bulk total RNA-seq. c) STORM-seq reveals single cell TE-derived transcript expression heterogeneity in previously reported bulk-level TE-derived transcript candidates in K-562 (Shah et al., Nat. Genetics 2023), demonstrating that bulk-level discovery of TE-derived transcripts does not necessarily indicate ubiquitous expression across single cells. d) Discovery and characterization of additional TE-derived transcript expression heterogeneity found in bulk total RNA-seq. LTR1A2-PURPL is being shown as an example of consistent expression across replicates in bulk total RNA-seq, but heterogeneous expression across K-562 single cells. Transcript names beginning with “TU” are unannotated genes/transcripts found in Ensembl 101 annotations. TE-derived transcripts above and below LTR1A2-PURPL have been removed for aesthetic purposes but the full list can be found in Supplemental Table 3. e) STORM-seq reads spanning annotated PURPL transcripts and LTR1A2-PURPL TE-derived transcripts are shown in pseudobulk (combined and assigned annotated PURPL and LTR1A2-PURPL tracks) and single cells with associated splice junction spanning reads as a percent of the maximum spanning read depth found in the combined track. The LTR1A2-PURPL TE-derived transcript comprises the majority of PURPL expression as shown by the percentage of total junction spanning reads in the assigned LTR1A2-PURPL track. f) Stranded coverage profiles of pseudobulk STORM-seq and VASA-seq across intergenic, distal K-562 enhancers, with transient transcriptome sequencing (TT-seq) and PRO-cap coverage. g) Detection rates of single-cell distal eRNAs when subsampled to 150k reads/cell shows STORM-seq identifies approximately twice as many eRNAs per cell compared to VASA-seq. h) Proportion of single cells with a detected, known gene fusion in K-562 from CCLE. Cell type proportions are weighted by respective bulk RNA-seq technology detection sensitivity to allow total RNA and mRNA single-cell protocols to be comparable. Single cells were downsampled to 150k reads/cell. i) Circos plot showing genomic alterations of known gene fusions (CCLE) in K-562 recovered by STORM-seq.

Figure 3.. STORM-seq identifies a continuous differentiation trajectory from an unclassified progenitor/”dual-feature” cell population in primary human fallopian tube epithelium.

a) RNA velocity supports a continuous differentiation process from an unclassified fallopian tube progenitor (UCFP) cells (cluster 3) to an intermediate branch point (cluster 5) to terminally differentiated secretory and ciliated cells (clusters 7–8 and 4, respectively). Thickness and direction of arrows indicate velocity. b) Latent time inference of cellular differentiation and c) CytoTRACE pseudotime supports a similar differentiation trajectory from early UCFP/”dual-feature” cells towards late, secretory and ciliated cell types. d) PCA space embedding of patient 1 with fitted principal curves and pseudotime show similar lineage trajectories as shown above. e) PCA space embedding of patient 2 with fitted principal curves and pseudotime show consistency with patient 1 lineage fate from UCFP to differentiatied secretory and ciliated cell types. RNA velocity, latent time, CytoTRACE pseudotime were calculated using scvelo. densMAP embeddings shown for patient 1 (a-c). Principal curves and pseudotime for patients 1 and 2 were calculated using slingshot (d-e). f) Representative cell-type specific intergenic transposable element (TE) expression from patient 1 demonstrates lineage restricted TE expression, with example locus-level LINE expression within the ciliated cell lineage. g) Independent validation of the presence of the “dual-feature” UCFP cells using cyclic immunofluorescence across 3 additional patients. Cell type classification and quantification of “dual-feature” cells as a proprotion of total cells. h) Separate regions of interest (ROI) with zoomed insets to show representative simultaneous expression of epithelial (PanCK) and non-epithelial (CD31) cells. i-k) FACS analysis of matched patient single cell suspensions using the same antibodies used for STORM-seq patients 1 and 2 demonstrate similar proportions of UCF-P/”dual-feature” cell populations (EpCAM+/CD31+).