Effective ribosomal RNA depletion for single-cell total RNA-seq by scDASH

Cell culture and sorting

The human colon cancer cell line HCT116 (CCL-247™) was purchased from the American Tissue-type Culture Collection (ATCC). HCT116 cells were maintained in Dulbecco’s Modified Eagle Medium (DMEM) (Thermo Fisher Scientific, Cat# 11965118) supplemented with 10% Fetal Bovine Serum (Thermo Fisher Scientific, Cat# 16000044) and 1% Penicillin/Streptomycin (Thermo Fisher Scientific, Cat# 15070063). Cells were cultured to 70–80% confluence, trypsinized and collected. Single cells were sorted using BD FACSAria™ IIIu cell sorter (BD Biosciences) into 200 µl 8-tube strip containing lysis buffer prepared according to MATQ-seq protocol (Sheng et al., 2017). 100 cells were sorted for each mini-bulk sample.

Library construction and sequencing

Single-cell or mini-bulk RNA-seq libraries were constructed as instructed by the MATQ-seq protocol. cDNA libraries were quantified using a Qubit 3 Fluorometer with dsDNA HS Assay Kit (Thermo Fisher Scientific, Cat# Q32851). Quality check on cDNA library size (Fig. S1) was performed using Fragment Analyzer (Advanced Analytical) with HS NGS Fragment Kit (1–6000 bp) (Agilent, Cat# DNF-474-1000). Individual libraries were pooled using Nextseq500/550 High Output Kit v2.5 (Illumina, Cat# 20024906). As an optional step, pooled libraries were size-selected using Agencourt AMPure XP magnetic beads (Beckman Coulter, Cat# A63882). Pooled libraries were then subjected to scDASH treatment as illustrated in Fig. 1, followed by enrichment amplification with P5 and P7 primers. The detailed scDASH protocol can be found in Supplemental File (scDASH Protocol). Illumina libraries were sequenced using a Nextseq500/550 (Illumina) to obtain approximately 1.5 million paired sequencing reads for each library.

In silico design and generation of sgRNA sequences

To ensure the utility of the sgRNA library in targeting inter-chromosomal loci, we obtained the reference sequences of non-clustered 45S rRNA genes (18S: NR_003286.2; 5.8S: NR_003285.2; 28S: NR_003287.2) and used this in parallel with a recently assembled sequence of ribosomal 45S RNA cluster (GenBank accession KY962518.1) (Kim et al., 2018) to generate sgRNA candidates in silico. We used Benchling (Benchling, 2020) to design sgRNA on both positive and negative strands of the target sequences, with specificity score (off-target score) and efficiency score (on-target score) calculated for each sgRNA. The models used for scoring guide sequences in terms of specificity and efficiency were proposed by Hsu et al. (2013) and Doench et al. (2016) respectively. Only the common sgRNA candidates generated from both non-clustered reference sequences and KY962518.1 were considered. We selected only those commonly identified sgRNAs with both off-target and on-target scores 0.5 standard deviation higher than the mean. We selected 58 sgRNAs against the 7,096 bp-long sequence of 18S, 5.8S and 28S rRNA genes. The sgRNA target sites are approximately spaced out by 100 bp intervals. A complete list of the sgRNA sequences can be found in Table S1.

In vitro sgRNA synthesis

Transcription templates of each sgRNA were split into halves and purchased as DNA oligos. The first half (54–55 nt) contains a T7 RNA polymerase promoter, a 20-nt target-specific spacer region and a 14-nt overlapping scaffold sequence. G was added upstream if spacer sequence does not begin with a G, in order to satisfy the sequence requirements of the T7 RNA polymerase promoter. The entire sequence is as follows: 5′  TTCTAATACGACTCACTATA(G)(crRNA-20)GTTTTAGAGCTAGA 3′. The second half (80 nt) constitutes a constant region required for the association of sgRNA with Cas9 protein: 5′  AAAAGCACCGACTCGGTGCCACTTTTTCAAGTTGATAACGGACTA GCCTTATTTTAACTTGCTATTTCTAGCTCTAAAAC 3′. The two halves of the sgRNA transcription templates were assembled and amplified with a PCR reaction using KAPA HiFi™ HotStart DNA Polymerase (Roche, Cat# 7958935001). The resulting 58 amplicons (120-121 bp) were pooled and purified using NucleoSpin^® Gel and PCR Clean-up Kit (Macherey Nagel, Cat# 740609.250) and their sizes were verified by agarose gel electrophoresis.

In vitro transcription was performed using the HiScribe™ T7 Quick High Yield RNA Synthesis Kit (NEB, Cat# E2050S) in accordance with the protocol provided. We performed overnight incubation to maximize yield. Samples were then cleaned up using RNA Clean and Concentrator-5 (Zymo Research, Cat# R1013), and typically 60–90 µg of RNA were obtained after the clean-up. We then verified transcript size using Fragment Analyzer Automated CE System (DNF 472).

In vitro Cas9 digestion

CRISPR/Cas9 treatment was performed as described by Gu et al. (2016) with certain adaptations. Formation of ribonucleoprotein (RNP) complex was achieved by mixing 90 ng of SpCas9 nuclease (NEB, Cat# M0386S) and 2700 ng of sgRNAs in Cas9 buffer. This corresponds to 100-fold excess of Cas9 protein and 3000-fold excess of sgRNAs relative to the rRNA substrate (i.e., 0.9 ng of rRNA substrate if using 1 ng of cDNA as input). Incubation duration for Cas9 digestion of sample cDNA was extended to 2 h and 30 min. Digested samples were purified using Agencourt AMPure XP magnetic beads (Beckman Coulter, Cat# A63882) and either directly subject to PCR enrichment or stored at −20 °C.

Detailed protocols for in vitro sgRNA synthesis and Cas9 digestion are supplied in the Supplemental File.

Quantitative PCR

Following PCR enrichment with P5 and P7 primers and subsequent clean-up, qPCR analysis was performed on the cDNA library using Quantagene q225 (Kubo Technology). A set of primer sequences specific to 18S, 5.8S and 28S was designed such that the minimum distance between the primer recognition regions and Cas9 target site is 20 nt (Table S2). Levels of rRNA measured in the cDNA library were normalised to GAPDH levels. Relative quantification of rRNA cDNA levels was performed using the delta delta CT method (Pfaffl, 2001).

Data processing and analysis

All sequencing raw data have been deposited on Gene Expression Omnibus: GSE158880. Demultiplexing and adapter trimming of raw reads were performed using BaseSpace (Illumina). Read quality was verified using FastQC (Andrews, 2010) where the mean Phred quality score across each base position in the read is > 30. Quality-checked reads were quantified with kallisto (Bray et al., 2016) using a concatenated reference sequence comprising Gencode v.29 human reference transcriptome and KY962518.1 rRNA sequence. Per-nucleotide coverage of 45S rRNA sequence was calculated from .bam files using the depth command in SAMtools. The calculated depth at each position was then normalised to the total depth. Gene-level read count was obtained by importing kallisto-generated transcript abundance to tximport (Soneson, Love & Robinson, 2015). Library size difference was normalised by random downsampling of gene-level read count using custom R scripts. This random downsampling was repeated three times for each sample, and either the mean or median read count was taken for different analyses (see Results). All data processing and analysis code are available at GitHub repository: https://github.com/lscdanson/scDASH.

Graphic design and visualisation tools

Figure 1 was created with BioRender.com (, 2020). All figures in the results section were generated using ggplot2 (Wickham, 2016).

scDASH can effectively deplete rRNA genes in scRNA-seq libraries

The overall scDASH approach is depicted in Fig. 1. The main differences between the previously described DASH and scDASH are: (1) the workflow and reaction conditions to generate the cDNA library follows that of the single cell transcriptome amplification protocol; (2) modified reaction conditions to achieve optimal Cas9-mediated depletion efficiency and specificity, in particular we found that the optimal concentration of Cas9 and sgRNA required in scDASH is higher than previously described for DASH (see Methods); (3) the specific library of sgRNA that we have designed to deplete rRNA. This rRNA-targeting sgRNA library design is non-trivial. Of the human cytoplasmic rRNA species, 5S rRNAs, ribosomal 45S RNA clusters, and non-clustered ribosomal 45S RNA genes, the 45S rRNA are the most abundant and targeting this group of rRNA sequences for depletion with the DASH approach would be the most effective. The 45S rDNA that encode for the ribosomal 45S RNA clusters are made up of 30-40 tandem repeats that are transcribed as a single unit to produce the precursor rRNA that is then later processed to generate the 18S, 5.8S, and 28S rRNA (Agrawal & Ganley, 2018). Due to its highly repetitive nature, a full reference sequence for the 45S rRNA was not available until recently, with the use of TAR cloning and long-read sequencing technology generated the first 44,838 bp-long reference of the human chromosome 21 rRNA genes (GenBank accession KY962518.1) (Kim et al., 2018). We used this rRNA reference to design a library of sgRNAs in silico that should enable targeted rRNA cleavage when applied in our scDASH protocol. However, due to the same technical challenges, the 45S rRNA clusters from other chromosomes have not yet been sequenced, and it is unclear what degree of sequence homology there is between the 45S rRNA clusters from different chromosomes within and between individuals. It has been reported that human rDNA copy number could be as high as 1,590 copies per individual (Parks et al., 2018), therefore if the sequence variation between the chromosomal rDNA clusters is large, it is possible that sgRNA designed using the 45S reference from chr21 may not effectively target those rDNA sequences from other chromosomes. On the other hand, other literature has reported that intra-individual variations with high allelic frequency are rare in human rDNA loci (Parks et al., 2018), which suggests that sequence homology between chromosomal 45S clusters may be sufficiently high, but this assumption has not been tested experimentally before this study. The resulting 58 sgRNAs chosen for this study are spaced out approximately in 100 bp intervals and were specifically selected against the 7,096 bp-long sequence of 18S, 5.8S and 28S rRNA genes in the KY962518.1 reference (Figs. 2B, 2C). Specific details of the sgRNA design rationale can be found in Materials and Methods.

We first examined the efficacy of scDASH in depleting cytoplasmic rRNA from scRNA-seq libraries. We designed primers against the sgRNA target sites on each of the 18S, 5.8S and 28S rRNA genes (Fig. S2), then performed qPCR using these primers to assess the quantity of sequence fragments from these target site in both untreated libraries and scDASHed libraries. Comparisons of the quantity of rRNA fragments in scRNA-seq libraries before and after scDASH, as measured by qPCR, indicates that scDASH treatment results in > 99% depletion in all three sgRNA-targeted rRNA gene loci (Fig. 2A).

Next, we comprehensively evaluated the performance of scDASH by comparing the sequencing results from 8 single-cell and 4 mini-bulk transcriptomes, half of which were subjected to scDASH and half were untreated. Reads that mapped to individual nucleotides within the relevant window from position 3,000 to 14,000 on the KY962518.1 reference sequence were quantified, as illustrated in Figs. 2B and 2C respectively for mini-bulk and single-cell samples. In untreated samples, there are multiple regions showing peaks with 8000x coverage, largely localized to within the 18S and 28S rRNA coding regions. Such a high coverage is presumably an aggregate total contributed by an average of 315 rRNA loci. Depletion with scDASH resulted in roughly 70% reduction in total coverage in this rRNA region for the mini-bulk libraries, and 55% reduction for the single-cell libraries. This is consistent with our qPCR measurement results. While most of the peaks were dramatically reduced in the DASHed samples, a few loci consistently appear to be less efficiently depleted. The reason for this is unclear; steric hinderance of Cas9 cleavage due to formation of secondary structures in-situ, or the sequence variations between those loci from different chromosomes are two possible speculative reasons that can be further investigated in a future study.

scDASH is highly specific in depleting cytoplasmic rRNA

Having established the applicability of scDASH in depleting cytoplasmic rRNAs in scRNA-seq libraries, we wondered whether the depletion is specific to the targeted 18S, 5.8S and 28S rRNA genes; in other words, we wanted to assess the off-target effects of this approach. To do so, we assessed the correlation in gene expression levels between scDASH treated and untreated sample groups, modelling the relationship between depleted rRNA and undisrupted gene expressions using linear regression (Fig. 3). By plotting the untreated sample on the x-axis, and the scDASH treated sample on the y-axis, genes that are unaffected by the treatment should fall close to the unity line as their expression should be highly correlated between treated and untreated samples; whereas genes that are depleted by scDASH will appear further right of the unity line, as their expression in the treated sample will be much lower than in the untreated sample. In scDASH treated single-cell RNA-seq libraries, we observed that rRNA genes fell to the right of the unity line as expected (Fig. 3A), while the overall agreement in gene expressions remains high with a coefficient of determination (r²) of 0.93. This suggests that the gene expression in the DASHed scRNA-seq libraries are largely similar to that of untreated samples (Fig. 3A), suggesting a negligible off-target depletion. It was previously reported that low expression genes in single-cell data can be highly variable among replicates due to dropout events, and that they can only be reliably detected with bulk or single-cell ensemble datasets (Wu et al., 2014). Therefore, to further validate our results, we performed the same off-target assessment experiment in mini-bulk samples of 100 cells per experiment (Fig. 3B). As expected, a higher (r²) value (0.97) was observed in mini-bulk sample pairs while the three targeted rRNA genes remain deviated from the regression line (Fig. 3B). This indicates that scDASH depletion is generally highly specific, and that the lower (r²) value seen in single-cell measurements is largely due to the greater level of technical variation that is expected in single-cell RNA-seq experiments.

Although we did not observe any coding genes or major classes of non-coding transcripts to be severely affected by off-target depletion, there were a few other ‘non-informative’ RNA species, that were depleted after scDASH treatment. First, the mitochondrially encoded rRNA genes MT-RNR1 and MT-RNR2 were not specifically targeted, but were found to be significantly depleted when compared to untreated samples. This is because over half of the sgRNAs we designed share a 100% sequence identity to different regions of MT-RNR1 and MT-RNR2, presumably a result of their sequence similarity with cytoplasmic rRNAs. In experiments where these mitochondrial rRNA are of-interest, those sgRNA that target regions with homology to the mitochondrial rRNA can be removed from the scDASH sgRNA pool, thus preserving these fragments in the final sequencing library while the trade-off may be a larger portion of rRNA fragments being sequenced as well. Another set of rRNA that were untargeted but depleted are three rRNA genes that do not originate from the 45S rRNA cluster: FP671120.3, FP236383.2 and FP236383.3. As shown in Fig. 3, these genes stand out as outliers in close proximity to the targeted 45S rRNA genes. To provide a quantitative measure of the depletion effect, we computed the z-score of these three untargeted but depleted rRNA genes. More specifically, we first calculated the difference in normalised read counts for every single gene between DASHed and untreated replicates (ΔDU), as well as between untreated replicates (ΔUU). Then, we computed a two-sample t distribution of the mean difference between ΔDU and ΔUU of all genes. The three off-target genes present a z-score ranging from 63–65, while 28S rRNA gene gives a z-score of 224. This shows that even though the untargeted genes were depleted, the depletion effect is still limited in extent when compared to the targeted rRNA genes. FP671120.3 and FP236383.2 belong to a class of non-coding RNA known as NoRC associated RNA, which function as antisense against rRNA genes when ribosome synthesis is not desired. These RNA species were depleted since they share a 100% sequence identity with multiple rRNA-specific regions on KY962518.1. FP236383.3 is classified as a novel subclass of long non-coding RNA (lncRNA) called long intervening non-coding RNA (lincRNA). Different regions on FP236383.3 exhibit a 99.6–100% sequence identity to 18S, 5.8S and 28S rRNA genes. Although the role of lincRNA remains elusive, it has been speculated that FP236383.3 could be a by-product of the concerted evolution of rRNA genes (Robicheau et al., 2017).

scDASH enables more informative whole transcriptome profiling at equal sequencing depth

We further evaluated the effect of cytoplasmic rRNA depletion on the detectability of other RNA types. Specifically, we calculated the number of reads detected for each RNA type at gene-level before and after scDASH treatment, while accounting for sequencing library depth. In the untreated libraries, consistent with past analyses (Eun, 1996), rRNA genes accounted for the largest proportion of reads (78%) in single-cell transcriptomes (Table 1; Fig. 4A). Following scDASH depletion, rRNA-aligned reads were brought down to 34% compared to without scDASH (Fig. 4B). Correspondingly, reads that align to coding genes expanded by three-fold from 20% to 61% after scDASH rRNA depletion (Figs. 4A, 4C). The detectable gene expression level of other RNA types including pseudogenes, lncRNAs, and microRNAs were all significantly elevated upon treatment with scDASH (Figs. 4D, 4E, 4F). The numbers of genes detected also appeared to increase with scDASH treatment, although the increase was not statistically significant (Fig. S3). Similar increases in detectability of non-rRNA species were also seen in the mini-bulk RNA-seq datasets (Fig. S4). Collectively, this illustrates the effectiveness of scDASH in freeing up sequencing space for non-rRNA species-of-interest in scRNA-seq applications.