zUMIs - A fast and flexible pipeline to process RNA sequencing data with UMIs

zUMIs - A fast and flexible pipeline to process RNA sequencing data with UMIs Background: Single-cell RNA-sequencing (scRNA-seq) experiments typically analyze hundreds or thousands of cells after amplification of the cDNA. The high throughput is made possible by the early introduction of sample-specific bar codes (BCs), and the amplification bias is alleviated by unique molecular identifiers (UMIs). Thus, the ideal analysis pipeline for scRNA-seq data needs to efficiently tabulate reads according to both BC and UMI. Findings: zUMIs is a pipeline that can handle both known and random BCs and also efficiently collapse UMIs, either just for exon mapping reads or for both exon and intron mapping reads. If BC annotation is missing, zUMIs can accurately detect intact cells from the distribution of sequencing reads. Another unique feature of zUMIs is the adaptive downsampling function that facilitates dealing with hugely varying library sizes but also allows the user to evaluate whether the library has been sequenced to saturation. To illustrate the utility of zUMIs, we analyzed a single-nucleus RNA-seq dataset and show that more than 35% of all reads map to introns. Also, we show that these intronic reads are informative about expression levels, significantly increasing the number of detected genes and improving the cluster resolution. Conclusions: zUMIs flexibility makes if possible to accommodate data generated with any of the major scRNA-seq protocols that use BCs and UMIs and is the most feature-rich, fast, and user-friendly pipeline to process such scRNA-seq data. Keywords: single-cell RNA-sequencing; digital gene expression; unique molecular identifiers ; pipeline porate unique molecular identifiers (UMIs) to label individual Introduction cDNA molecules with a random nucleotide sequence before am- The recent development of increasingly sensitive protocols al- plification [ 7]. This enables the computational removal of am- lows for the generation of RNA-sequencing (RNA-seq) libraries plification noise and thus increases the power to detect expres- of single cells [1]. The throughput of such single-cell RNA-seq sion differences between cells [8, 9]. To increase the throughput, (scRNA-seq) protocols is rapidly increasing, enabling the pro- many protocols also incorporate sample-specific bar codes (BCs) filing of tens of thousands of cells [ 2, 3] and opening exciting to label all cDNA molecules of a single cell with a nucleotide se- possibilities to analyze cellular identities [4, 5]. As the required quence before library generation [10]. This allows for early pool- amplification from such small starting amounts introduces sub- ing, which further decreases amplification noise [ 6]. Addition- stantial amounts of noise [6], many scRNA-seq protocols incor- Received: 17 October 2017; Revised: 16 March 2018; Accepted: 15 May 2018 The Author(s) 2018. Published by Oxford University Press. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy059/5005022 by Ed 'DeepDyve' Gillespie user on 21 June 2018 2 zUMIs - RNA-seq with UMIs Table 1: Features of available UMI pipelines for the quantification of gene expression data. Compatible Open Quality UMI BC Down- UMI library Name Reference source filter collapsing Mapper detection Intron sampling protocols Cell Ranger [2] yes BC+UMI Hamming STAR A no yes [2] distance CEL-seq [15] yes BC+UMI Identity bowtie2 WL no no [15, 46] only dropEst [16] yes BC Frequency- TopHat2 or WL,top- yes no [2, 13, 19] based Kallisto n,EM Drop-seq- [13] no BC+UMI Hamming STAR WL,top-n no no [13, 15, 17] tools distance scPipe [47] yes BC+UMI Hamming subread WL,top-n no no [13, 17, 18, 46] distance umis [14] yes BC Frequency- Kallisto WL,top- no no [2, 13, 17–19, based n,EM 46, 48] UMI-tools [25] yes BC+UMI Network- BWA WL no no [17, 19] based zUMIs This work yes BC+UMI Hamming STAR A,WL,top-n yes yes [2, 3, 12, 13, distance 15, 17, 18, 21, 46, 48] We consider whether the pipeline is open source, has sequence quality filters for cell BCs and UMIs, mappers, UMI-collapsing options, options for BC de tection (A, automatically infer intact BCs; WL, extract only the given list of known BCs; top-n, order BCs according the number of reads and keep the top n BCs; EM, merge BCs with given edit distance), whether it can count intron mapping reads, whether it offers a utility to make varying library sizes more comparable via downsampling, and finally with which RNA-seq library preparation protocols is it compatible ally, for cell types such as primary neurons, it has been proven Implementation and Operation to be more feasible to isolate RNA from single nuclei rather than Filtering and mapping whole cells [11, 12]. This decreases mRNA amounts further so that it has been suggested to count intron mapping reads origi- The first step in our pipeline is to filter reads that have low- nating from nascent RNAs as part of single-cell expression pro- quality BCs according to a user-defined threshold (Fig. 1). This files [ 11]. However, the few bioinformatic tools that process RNA- step eliminates the majority of spurious BCs and thus greatly seq data with UMIs and BCs have limitations (Table 1). For ex- reduces the number of BCs that need to be considered for count- ample, the Drop-seq-tools is not an open source [13]. While Cell ing. Similarly, we also filter low-quality UMIs. Ranger is open, it is exceedingly difficult to adapt the code to The remaining reads are then mapped to the genome using new or unknown sample BCs and other library types. Other tools the splice-aware aligner STAR [22]. The user is free to customize are specifically designed to work with one mapping algorithm mapping by using the options of STAR. Furthermore, if the user and focus mainly on transcriptome references [14, 15]. Further- wishes to use a different mapper, it is also possible to provide more, the only other UMI-RNA-seq pipeline providing the utility zUMIs with an aligned bam file instead of the fastq file with the to also consider intron mapping reads, dropEst [16], is only appli- cDNA sequence, with the sole requirement that only one map- cable to droplet-based protocols. Here, we present zUMIs, a fast ping position per read is reported in the bam file. and flexible pipeline that overcomes these limitations. Transcript counting Findings Next, reads are assigned to genes. In order to distinguish exon zUMIs is a pipeline to process RNA-seq data that were mul- and intron counts, we generate two mutually exclusive an- tiplexed using cell BCs and also contain UMIs. Read-pairs are notation files from the provided gtf, one detailing exon posi- filtered to remove reads with low-quality BCs or UMIs based tions, the other introns. Based on those annotations, Rsubread on sequence and then mapped to a reference genome (Fig.1). featureCounts [23] is used to first assign reads to exons and af- Next, zUMIs generates UMI and read count tables for exon and terward to check whether the remaining reads fall into introns, exon+intron counting. We reason that very low input material in other words, if a read is overlapping with intronic and ex- such as from single nuclei sequencing might profit from in- onic sequences, it will be assigned to the exon only. The output cluding reads that potentially originate from nascent RNAs. An- is then read into R using data.table [24], generating count ta- other unique feature of zUMIs is that it allows for downsampling bles for UMIs and reads per gene per BC. We then collapse UMIs of reads before collapsing UMIs, thus enabling the user to as- that were mapped either to the exon or intron of the same gene. sess whether a library was sequenced to saturation or whether Note that only the processing of intron and exon reads together deeper sequencing is necessary to depict the full mRNA com- allows for properly collapse of UMIs that can be sampled from plexity. Furthermore, zUMIs is flexible with respect to the length the intronic as well as from the exonic part of the same nascent and sequences of the BCs and UMIs, supporting protocols that mRNA molecule. have both sequences in one read [2, 3, 12, 13, 15, 17, 18]aswell Per default, we only collapse UMIs by sequence identity. If as protocols that provide UMI and BC in separate reads [19–21]. there is a risk that a large proportion of UMIs remains under- This makes zUMIs the only tool that is easily compatible with all collapsed due to sequence errors, zUMIs provides the option to major UMI-based scRNA-seq protocols. collapse UMIs within a given Hamming distance. We compare Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy059/5005022 by Ed 'DeepDyve' Gillespie user on 21 June 2018 Parekh et al. 3 Figure 1: Schematic of the zUMIs pipeline. Each of the gray panels from left to right depicts a step of the zUMIs pipeline. First, fastq files are filtered according to user-defined bar code (BC) and unique molecular identifier (UMI) quality thresholds. Next, the remaining cDNA reads are mapped to the reference genome using STAR. Gene-wise read and UMI count tables are generated for exon, intron, and exon+intron overlapping reads. To obtain comparable library sizes, reads can be downsampled to a desired range during the counting step. In addition, zUMIs also generates data and plots for several quality measures, such as the number of detected genes/UMIs per BCe and distribution of reads into mapping feature categories. (A) (B) Drop−seq−tools vs zUMIs 1 < 17 UMI−tools vs zUMIs 1 < 17 UMI−tools vs zUMIs Hamming1 1.000 0.997 0.995 0.992 0.990 Drop−seq−tools UMI−tools zUMIs 1e+02 (C) (D) zUMIs − 1 < 17 zUMIs − Hamming1 UMI−tools − edit1 directional Drop−seq−tools − 1 < 17 1e+00 1e−02 Unfiltered zUMIs − 1 < 17 zUMIs − Hamming1 1e−04 Drop−seq−tools − 1 < 17 UMI−tools − edit1 directional 50000 100000 150000 200000 UMIs unfiltered Figure 2: Comparison of different UMI collapsing methods. We compared Drop-seq-tools and UMI-tools with zUMIs using our HEK dataset (227 mio reads). (A) Run time to count exonic UMIs using zUMIs (hamming distance = 0), UMI-tools (”unique” mode) and Drop-seq-tools (edit distance = 0). (B) Box plots of correlation coefficients of gene-wise UMI counts of the same cell generated with different methods. UMI counts generated using zUMIs (quality filter “1 base under phred 17” or ha mming distance = 1) were correlated to UMI counts generated using Drop-seq-tools (quality filter “1 base under phred 17” ) and UMI-tools (“directional adjacency” mode ). (C) Comparison of the total number of UMIs per cell derived from different counting methods to “unfiltered” counts. (D) Violin plots of gene-wise dispersion estimates with different quality filtering and UMI collapsing methods. Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy059/5005022 by Ed 'DeepDyve' Gillespie user on 21 June 2018 UMIs filtered Time taken (minutes) Dispersion Pearson's R 4 zUMIs - RNA-seq with UMIs the two zUMIs UMI-collapsing options to the recommended di- or UMIs enough to justify the extra cost. In our HEK-cell exam- rectional adjacency approach implemented in UMI-tools [25]us- ple dataset, the number of detected genes starts leveling off at ing our in-house example dataset (see Methods section). zUMIs 1 million reads. Sequencing double that amount would only in- identity collapsing yields nearly identical UMI counts per cell crease the number of detected genes from 9,000 to 10,600 when as UMI-tools, while Hamming distance yields increasingly fewer counting exon reads (Fig.3D). In line with previous findings [ 8, UMIs per cell with increasing sequencing depth (Fig.2C). Smith 14], the saturation curve of exon+intron counting runs parallel et al [25] suggest that edit distance collapsing without consid- to the one for exon counting, both indicating that a sequencing ering the relative frequencies of UMIs might indeed overreach depth of 1 million reads per cell is sufficient for these libraries. and overcollapse the UMIs. We suspect that this is indeed what happens in our example data, where we find that gene-wise dis- Output and statistics persion estimates appear suspiciously truncated as expected if several counts are unduly reduce to one, the minimal number zUMIs outputs three UMI and three read count tables: gene-wise after collapsing (Fig.2D). counts for traditional exon counting, one for intron and one for However, note that the above-described differences are mi- exon+intron counts. If a user chooses the downsampling option, nor. By and large, there is good agreement between UMI counts six additional count tables per target read count are provided. To obtained by UMI-tools [25], the Drop-seq pipeline [13], and zU- evaluate library quality, zUMIs summarizes the mapping statis- MIs. The correlation between gene-wise counts of the same cell tics of the reads. While exon and intron mapping reads likely is >0.99 for all comparisons (Fig. 2B). In light of this, we consider represent mRNA quantities, a high fraction of intergenic and un- the >3 times higher processing speed of zUMIs to be a decisive mapped reads indicates low-quality libraries. Another measure advantage (Fig.2A). of RNA-seq library quality is the complexity of the library, for which the number of detected genes and the number of iden- tified UMIs are good measures (Fig. 1). We processed 227 mil- Cell BC selection lion reads with zUMIs and quantified expression levels for exon In order to be compatible with well-based and droplet-based and intron counts on a Unix machine using up to 16 threads, scRNA-seq methods, zUMIs needs to be able to deal with known which took less than 3 hours. Increasing the number of reads as well as random BCs. As default behavior, zUMIs infers which increases the processing time approximately linearly, where fil- BCs mark good cells from the data (Fig.3A, 3B). To this end, we tering, mapping, and counting each take up roughly one third of fit a k-dimensional multivariate normal distribution using the the total time (Fig.3E). We also observed that the peak random R-package mclust [26, 27] for the number of reads/BC, where k access memory usage for processing datasets of 227, 500, and is empirically determined by mclust via the Bayesian informa- 1,000 million pairs was 42 Gb, 89 Gb, and 172 Gb, respectively. tion criterion. We reason that only the kth normal distribution Finally, zUMIs could process the largest scRNA-seq dataset re- with the largest mean contains BCs that identify reads originat- ported to date with around 1.3 million brain cells and 30 billion ing from intact cells. We exclude all BCs that fall in the lower read-pairs generated with 10xGenomics Chromium (see Meth- 1% tail of this kth normal distribution to exclude spurious BCs. ods section) on a 22-core processor in only 7 days. The HEK dataset used here contains 96 cells with known BCs and zUMIs identifies 99 BCs as intact, including all the 96 known Intron counting BCs. Also, for the single-nucleus RNA-seq from Habib et al. [12], zUMIs identified a reasonable number of cells; Habib et al. report Recently, it has been shown that intron mapping reads in RNA- 10,877 nuclei and zUMIs identified 11,013 intact nuclei. However, seq likely originate from nascent mRNAs and are useful for gene we recommend to always check the elbow plot generated by zU- expression estimates [31, 32]. Additionally, novel approaches MIs (Fig.3B) to confirm that the cutoff used by zUMIs is valid for a leverage the ratios of intron and exon mapping reads to infer given dataset. In cases where the number of BCs or BC sequences information on transcription dynamics and cell states [33]. To are known, it is preferable to use this information. If zUMIs is ei- address this new aspect of analysis, zUMIs also counts and col- ther given the number of expected BCs or is provided with a list lapses intron-only mapping reads as well as intron and exon of BC sequences, it will use this information and forgo automatic mapping reads from the same gene with the same UMI. To as- inference. sess the information gain from intronic reads to estimate gene expression levels, we analyzed a publicly available DroNc-seq dataset from mouse brain ([12]; see Methods section). For the Downsampling ∼11,000 single nuclei of this dataset, the fraction of intron map- scRNA-seq library sizes can vary by orders of magnitude, which ping reads of all reads goes up to 61%. Thus, if intronic reads complicates normalization [28, 29]. A straight-forward solution are considered, the mean number of detected genes per cell for this issue is to downsample overrepresented libraries [30]. increases from 1,041 for exon counts to 1,995 for exon+intron zUMIs has an built-in function for downsampling datasets to a counts. Next, we used the resulting UMI count tables to investi- user-specified number of reads or a range of reads. By default, gate whether exon+intron counting improves the identification zUMIs downsamples all selected BCs to be within three absolute of cell types, as suggested by Lake et al. [11]. The validity and ac- deviations from the median number of reads per BC (Fig.3C). Al- curacy of counting introns for single-nucleus sequencing meth- ternatively, the user can provide a target sequencing depth, and ods has recently been demonstrated [34]. Following the Seurat zUMIs will downsample to the specified read number or omit pipeline to cluster cells [35, 36], we find that using exon+intron the cell from the downsampled count table if fewer reads were counts discriminates 28 clusters, while we could only discrimi- present. Furthermore, zUMIs also allows the user to specify a nate 19 clusters using exon counts (Fig.4A, 4B). The larger num- multiple target read number at once for downsampling. This ber of clusters is not simply due to the increase in the counted feature is helpful if the user wishes to determine whether the UMIs and genes. When we permute the intron counts across RNA-seq library was sequenced to saturation or whether fur- cells and add them to the exon counts, the added noise actu- ther sequencing would increase the number of detected genes ally reduces the number of identifiable clusters (Fig. 4E). Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy059/5005022 by Ed 'DeepDyve' Gillespie user on 21 June 2018 Parekh et al. 5 1.5e+07 (A) (B) (C) 2.0e+08 0.9 1.0e+07 1.5e+08 0.6 1.0e+08 0.3 5.0e+06 5.0e+07 0.0 0.0e+00 0.0e+00 0 1000 2000 3000 4000 log10(Reads) Cell Index Cells (D) (E) Counting Filtering Mapping Summarizing total Exon Exon+Intron 0e+00 1e+06 2e+06 3e+06 200 400 600 800 1000 Sequencing Depth Million Read pairs Figure 3: Utilities of zUMIs. Each of the panels shows the utilities of zUMIs pipeline. The plots from A–D show the results from the example HEK dataset used here. (A) The plot shows a density distribution of reads per BC. Cell BCs with reads right of the blue line are selected. (B) The plot shows the cumulative read distribution in the example HEK dataset where the BCs in light blue are the selected cells. (C) The bar plot shows the number of reads per selected cell BC with the red lines showing upper and lower median absolute deviation (MAD) cutoffs for adaptive downsampling. Here, the cells below the lower MAD have very low coverage and are discarded in downsampled count tables. (D) Cells were downsampled to six depths from 100,000 to 3,000,000 reads. For each sequencing depth, the genes detected per cell are shown. (E) Runtime for three datasets with 227, 500, and 1,000 million read-pairs. The runtime is divided in the main steps of the zUMIs pipeline as follows: filtering, mapping, counting, and summarizing. Each dataset was processed using 16 threads (“-p 16”). We continue to further characterize the seven clusters that neurons [37], suggesting that the split of cluster 7 might be bio- were subdivided by the addition of intron counts (Fig.4D). First, logically meaningful (Fig.4E). we identify DE genes between the newly formed clusters. If we In order to evaluate the power gained by exon+intron count- count only exon reads, there appear to be, on average, only 10 ing in a more systematic way, we perform power simulations us- DE genes between the subgroups, while exon+intron counting ing empirical mean and dispersion distributions from the largest yields ∼10 times more DE genes, thus corroborating the signal and most uniform cluster (∼1,500 cells) [9]. For a fair comparison, found with clustering. The log2-fold changes of those additional we include all detected genes, which is equivalent to the num- DE genes estimated with either counting strategy are generally ber of genes detected with exon+intron counting. Also, since we in good agreement; especially large log2-fold changes are de- call a gene detected as soon as one count is associated, exon tected with both exon and exon+intron counting (Fig.4F). Genes counting is necessarily a subset of exon+intron. Thus, there are, that are detected as DE in only one of our counting strategies on average, 4 times more genes in the lowest expression quan- have small log2-fold changes, and there are more of these small tile for exon counting than for exon+intron counting (Fig.4H). For changes detected using exon+intron counting. those genes, expression is too spurious to be used for differen- Detecting more genes naturally increases the chance to also tial expression analysis; for exon+intron counting, we have, on detect more informative genes. Here, we cross-reference the average, 60% power to detect a DE gene in the first mean expres- gene list with marker genes for transcriptomic subtypes de- sion bin with a well-controlled false discovery rate (FDR) (Fig.4G). tected for major cell types of the mouse brain [37] and find that In summary, the increased power for exon+intron counting and ∼5% of the additional genes are also marker genes, which cor- probably also the larger number of clusters are due to better de- responds well to the general frequency of marker genes among tection of lowly expressed genes. Furthermore, we think that al- the detected genes (4%). In the same vein, we also detect propor- though potentially noisy, the large number of additionally de- tionally more DE genes with exon+intron counting compared to tected genes makes exon+intron counting worthwhile, espe- exon counting. Thus, including introns simply allows us to bet- cially for single-nuclei sequencing techniques that are enriched ter detect present transcripts, while leaving the proportions of for nuclear nascent RNA transcripts, such as DroNc-seq [12]. Ad- interest unaltered. Having a closer look at cluster 7, it was split ditionally, exon+intron counting may help in extracting as much into a bigger (7) and a smaller cluster (24) using exon+intron information as possible from low coverage data as generated counting (Fig.4A-C), we find one marker gene (Il1rapl2) to be in the context of high-throughput cell atlas efforts (e.g., 10,000– DE between the subclusters using exon+intron counting, while 20,000 reads/cell [38, 39]. Last, users should always exclude the Il1rapl2 had only spurious counts using exon counts. Il1rapl2 is possibility of intronic reads stemming from genomic DNA con- a marker for transcriptomic subtypes of GABAergic Pvalb-type tamination in the library preparation by confirming low inter- Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy059/5005022 by Ed 'DeepDyve' Gillespie user on 21 June 2018 Density Genes Detected Cumulative Reads Runtime (hours) Number of Reads 2 7 6 1 5 8 3 Exon+ Exon Intron 6 zUMIs - RNA-seq with UMIs (A) (B) (C) (D) Exon Exon+Intron 60 27 exPFC 7 GABA 7 CA 19 6 DG SMC 14 13 1 0 OPC 7777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777 7 ODC MG 9 ASC −25 −30 8 8 8 END 0.00 0.25 0.50 0.75 1.00 10 1000 −50 Cell Type Fraction −40 −20 0 20 40 −50 −25 0 25 50 Number of DE Genes (E) (F) Exon+fake−Intron reads Exon Exon+Intron (H) 7 2000 20 01234 Il1rapl2− normalized counts (I) 0.75 (G) 0.50 0.4 Both 0.25 Exon 0.0 0.00 Exon+Intron 0.15 −0.4 0.10 −0.8 −20 0.05 −1.2 0.00 −2 −1 0 Log2−Fold Change Exon only −20 0 20 strata of mean−expression Figure 4: Contribution of intron reads to biological insights. We analyzed published single-nucleus RNA-seq data from mouse prefrontal cortex (PFC) and hippocampus [12] to assess the utility of counting intron in addition to exon reads. We processed the raw data with zUMIs to obtain expression tables with exon reads as well as exon+intron reads and then used the R-package Seurat [35, 36] to cluster cells. With exon counts, we identified 19 clusters (A), and with exon+intron counts we identified 27 clusters (B). Clusters are represented as t-SNE plots and colored according to the most frequent cell-type assignment in the original article [12]: glutamatergic neurons from the prefrontal cortex (exPFC), GABAergic interneurons (GABA), pyramidal neurons from the hippocampal CA region (CA), granule neurons from the hippocampal dentate gyrus region (DG), astrocytes (ASC), microglia (MG), oligodendrocytes (ODC), oligodendrocyte precursor cells (OPC), neuronal stem cells (NSC), smooth muscle cells (SMC) and endothelial cells (END). Different shades of those clusters indicate that multiple clusters had the same major cell type assigned. If we randomly sample counts from the intron data and add them to the exon counting, the noise reduces the number of clusters and the Seurat pipeline can only identify 9–11 clusters (E). The composition of each cluster based on exon+intron is detailed in panel (C), and cells that were not assigned a cell type in [12] are displayed as empty. The boxes mark the clusters that were not split when using exon data only. For example, cluster 7 from exon counting, which mainly consists of GABAergic neurons, was split into clusters 7, 24 (506, 66 cells) when using exon+intron counting. In (D), we show the numbers of genes that were differentially expressed (DE) (limma p-adj <0.05) between the clusters found only with exon+intron counts. The panel numbers represent the exon counting cluster numbers and the y-axis the exon+intron counting cluster number. The log2-fold changes corresponding to these contrasts are also used in (G). Among the genes that were additionally detected to be DE by exon+intron −5 counting was the marker gene Il1rapl2 (limma p-adj = 10 ). In (F), we present a violin plot of the normalized counts for Il1rapl2 in cells of the GABAergic subclusters 7 and 24. Log2-fold changes calculated with exon+intron counts correlate well with exon counts (G). Note that for exon counting only, half as many genes could be evaluated as for exon+intron counting and thus only half of the exon+intron genes are depicted in (G). Large log fold changes (LFCs) are found to be significant with both counting strategies (purple points are close to the bisecting line). We conducted simulations based on mean and dispersion measured using exon cluster 0 (1,616 cells, ∼90% exPFC). In (I) we show the expected true positive rate and the false discovery rate for a scenario comparing 300 vs 300 cells. Results for exon and exon+intron counting were stratified into five quantiles according to the mean expression of genes, where stratum 1 contains lowly expressed genes and stratum 5 the most highly expressed genes. The numbers of genes falling into each of the bins using exon+intron and exon counting are depicted in (H). genic mapping fractions using the statistics output provided by Analyzed RNA-seq datasets zUMIs. HEK293T cells were cultured in DMEM high glucose with L- glutamine (Biowest) supplemented with 10% fetal bovine serum (Thermo Fisher) and 1% penicillin/streptomycin (Sigma-Aldrich) Conclusion in a 37 C incubator with 5% carbon dioxide. Cells were passaged zUMIs is a fast and flexible pipeline for processing raw reads to and split every 2 or 3 days. For single-cell RNA-seq, HEK293T cells were dissociated by incubation with 0.25% Trypsin (Sigma- obtain count tables for RNA-seq data using UMIs. To our knowl- edge, it is the only open source pipeline that has a BC and Aldrich) for 5 minutes at 37 C. The single-cell suspension was washed twice with phosphate-buffered saline, and dead cells UMI quality filter, allows intron counting, and has an integrated downsampling functionality. These features ensure that zUMIs were stained with Zombie Yellow (Biolegend) according to the is applicable to most experimental designs of RNA-seq data, manufacturer’s protocol. Single cells were sorted into DNA including single-nucleus sequencing techniques, droplet-based LoBind 96-well polymerase chain reaction (PCR) plates (Eppen- methods where the BC is unknown, as well as plate-based UMI- dorf) containing lysis buffer with a Sony SH-800 cell sorter in 3- methods with known BCs. Finally, zUMIs is computationally ef- drop purity mode using a 100-μmnozzle. Next, single-cell RNA- ficient, user-friendly, and easy to install. seq libraries were constructed from one 96-well plate using a slightly modified version of the mcSCRB-seq protocol. Reverse transcription was performed as described previously [40], with the only change being the use of KAPA HiFi HotStart enzyme for Methods Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy059/5005022 by Ed 'DeepDyve' Gillespie user on 21 June 2018 Log2−Fold Change Exon+Intron Exon+Intron cluster ID No. of FDR TPR Cluster ID (Exon+Intron) Genes Parekh et al. 7 PCR amplification of cDNA. Resulting libraries were sequenced ters. To illustrate that the additional clusters found by count- using an Illumina HiSeq1500 with 16 cycles in Read 1 to decode ing exon+intron reads are not spurious, we use intron-only UMI cell BCs (6 bases) and UMIs (10 bases) and 50 cycles in Read 2 to counts from the same data to add to the observed exon-only sequence into the cDNA fragment, obtaining ∼227 million reads. counts. More specifically, to each gene we add scran-size factor- Raw fastq files were processed using zUMIs, mapping to the hu- corrected intron counts from the same gene after permuting man genome (hg38) and Ensembl gene models (GRCh38.84). them across cells. We assessed the cluster numbers from 100 In addition, we analyzed data from 1.3 million mouse brain such permutations. cells generated on the 10xGenomics Chromium platform [2]. Se- quences were downloaded from the National Center for Biotech- Comparison of UMI collapsing strategies nology Information Sequence Read Archive under accession number SRP096558. The data consist of 30 billion read-pairs In order to validate zUMIs and compare different UMI collaps- from 133 individual samples. In these data, read 1 contains 16 bp ing methods, we used the HEK dataset described above. We ran for the cell BC and 10 bp for the UMI and read 2 contains 114 bp zUMIs (1) without quality filtering, (2) filtering for onebase un- of cDNA. zUMIs was run using default settings, and we allowed der Phred 17, and (3) collapsing similar UMI sequences within a 7 threads per job for a total of up to 42 threads on an Intel Xeon hamming distance of 1. To compare with other available tools, E5-2699 22-core processor. we ran the same dataset using the Drop-seq-tools version 1.13 Finally, we obtained mouse brain DroNc-seq read data [12] [13] and quality filter “1 base under Phred 17” without edit dis- from the Broad Institute Single Cell Portal [41]. This dataset con- tance collapsing. Last, the HEK dataset was used with UMI tools sists of ∼1,615 million read-pairs from ∼11,000 single nuclei. [25] in (1) “unique” and (2) “directional adjacency” mode with Read 1 contains a 12 bpcell BC and a 8 bpUMI and read 2 60 bpof edit distance set to 1. Also, we compared the output of zU- cDNA. MIs from the DroNc-seq dataset when using default parame- The two mouse datasets were mapped to genome version ters (“1 base under Phred 20”) to UMI-tools in (1) “unique,” (2) mm10 and applying Ensembl gene models (GRCm38.75). “directional adjacency,” and (3) “cluster” settings. For each set- ting and tool combination, we compared per-cell/per-nuclei UMI contents in a linear model fit. Power simulations and DE analysis We evaluated the power to detect differential expression with thehelpofthe powsimRpackage [9]. For the DroNc-seq dataset, Availability of source code and requirements we estimated the parameters of the negative binomial dis- Project name: zUMIs tribution from one of the identified clusters, namely, cluster Project home page: https://github.com/sdparekh/zUMIs 0, compromising 1,500 glutamatergic neuronal cells from the Operating system(s): UNIX prefrontal cortex (Fig.4D). Since we detect more genes with Programming language: shell, R, perl exon+intron counting (4,433 compared to 1,782), we included Other requirements: STAR >= 2.5.3a, R >= 3.4, Rsubread >= this phenomenon in our read count simulation by drawing mean 1.26.1, pigz >= 2.3 & samtools >= 1.1 expression values for a total of 4,433 genes. This means that License: GNU GPLv3.0 the table includes sparse counts for the exon counting. Log - Research Resource Identification Initiative ID: SCR 016139 fold changes were drawn from a gamma distribution with shape equal to 1 and scale equal to 2. In each of the 25 simulation it- erations, we draw an equal sample size of 300 cells per group Availability of supporting data and test for differential expression using limma-trend [42]on log counts per million (CPM) values with scran [43] library size All data that were generated for this project were submitted to correction. The true positive rate and FDR are stratified over the GEO under accession GSE99822. An archival copy of the source empirical mean expression quantile bins. code and test data are available via the GigaScience repository For the differential expression analysis between clusters, we GigaDB [45]. use the same DE estimation procedure as in the simulations: scran normalization followed by limma-trend DE-analysis (c.f. [44]). Abbreviations BC: barcode; DE: differentially expressed; FDR: false discovery Cluster identification rate; MAD: median absolute deviation; PCR: polymerase chain reaction; PFC: prefrontal cortex; scRNA-seq: single-cell RNA se- After processing the DroNc-seq data [12] with zUMIs as de- quencing; UMI: unique molecular identifier. scribed above, we cluster cells based on UMI counts derived from exons only and exons+introns reads using the Seurat pipeline [35, 36]. First, cells with fewer than 200 detected genes were fil- Competing interests tered out. The filtered data were normalized using the LogNor- malize function. We then scale the data by regressing out the The authors declare that they have no competing interests. effects of the number of transcripts and genes detected per cell using the ScaleData function. The normalized and scaled data are then used to identify the most variable genes by fitting a Funding relationship between mean expression (ExpMean) and disper- This work has been supported by the Deutsche Forschungsge- sion (LogVMR) using the FindVariableGenes function. The iden- tified variable genes are used for principle component analysis, meinschaft (DFG) through SFB1243 subprojects A14/A15. and the top 20 principle components are then used to find clus- ters using graph-based clustering as implemented in FindClus- Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy059/5005022 by Ed 'DeepDyve' Gillespie user on 21 June 2018 8 zUMIs - RNA-seq with UMIs ing and sequencing using droplet microfluidics. Nat Protoc Author contributions 2017;12(1):44–73. S.P. and C.Z. designed and implemented the pipeline. B.V. tested 21. Hochgerner H, Lonnerber ¨ g P, Hodge R, et al. STRT-seq-2i: the pipeline and helped in power simulations. All authors con- dual-index 5’ single cell and nucleus RNA-seq on an address- tributed to writing the manuscript. able microwell array. Sci Rep 2017;7(1):16327. 22. Dobin A, Davis CA, Schlesinger F, et al. STAR: ultrafast uni- versal RNA-seq aligner. Bioinformatics 2013;29(1):15–21. References 23. Liao Y, Smyth GK, Shi W. featureCounts: an efficient general 1. Sandberg R. Entering the era of single-cell transcriptomics in purpose program for assigning sequence reads to genomic biology and medicine. Nat Methods 2014;11(1):22–4. features. Bioinformatics 2014;30(7):923–30. 2. Zheng GXY, Terry JM, Belgrader P, et al. Massively parallel 24. Dowle M, Srinivasan A. data.table: Extension of ‘data.frame.” digital transcriptional profiling of single cells. Nat Commun 2017, https://CRAN.R-project.org/package=data.table,r package version 1.10.4. 2017;8:14049. 3. Rosenberg AB, Roco CM, Muscat RA, et al. Single-cell profiling 25. Smith TS, Heger A, Sudbery I. UMI-tools: modelling sequenc- of the developing mouse brain and spinal cord with split- ing errors in unique molecular identifiers to improve quan- pool barcoding. Science 2018;360(6385):176–82. tification accuracy. Genome Res. 2017. 4. Wagner A, Regev A, Yosef N. Revealing the vectors of cel- 26. Fraley C, Raftery AE. Model-based clustering, discrimi- lular identity with single-cell genomics. Nat Biotechnol nant analysis, and density estimation. J Am Stat Assoc 2016;34(11):1145–60. 2002;97(458):611–31. 5. Regev A, Teichmann SA, Lander ES, et al. The Human Cell 27. Fraley C, Raftery AE Enhanced Model-Based Clustering, Atlas. Elife, 2017; 6. Density Estimation and Discriminant Analysis Software: 6. Parekh S, Ziegenhain C, Vieth B, et al. The impact of ampli- MCLUST. , J. Classification , 2003, 20, 263–286. fication on differential expression analyses by RNA-seq. Sci 28. Vallejos CA, Risso D, Scialdone A, et al. Normalizing single- Rep 2016;6:25533. cell RNA sequencing data: challenges and opportunities. Nat 7. Kivioja T, Vah ¨ ar ¨ autio A, Karlsson K, et al. Counting absolute Methods 2017;14(6):565–71. 29. Evans C, Hardin J, Stoebel DM. Selecting between-sample numbers of molecules using unique molecular identifiers. Nat Methods 2012;9(1):72–4. RNA-seq normalization methods from the perspective of 8. Ziegenhain C, Vieth B, Parekh S, et al. Comparative anal- their assumptions. Brief Bioinform 2017. ysis of single-cell RNA sequencing methods. Mol Cell 30. Grun ¨ D, van Oudenaarden A. Design and analysis of single- 2017;65(4):631–43.e4. cell sequencing experiments. Cell 2015;163(4):799–810. 9. Vieth B, Ziegenhain C, Parekh S, et al. powsimR: power anal- 31. Hendriks GJ, Gaidatzis D, Aeschimann F, et al. Extensive os- ysis for bulk and single cell RNA-seq experiments. Bioinfor- cillatory gene expression during C. elegans larval develop- matics 2017;33(21):3486–3488. ment. Mol Cell 2014;53(3):380–92. 10. Ziegenhain C, Vieth B, Parekh, et al. Quantitative 32. Gaidatzis D, Burger L, Florescu M, et al. Analysis of intronic single-cell transcriptomics. Brief Funct Genomics 2018; and exonic reads in RNA-seq data characterizes transcrip- doi:10.1093/bfgp/ely009. tional and post-transcriptional regulation. Nat Biotechnol 11. Lake BB, Ai R, Kaeser GE, et al. Neuronal subtypes and diver- 2015;33(7):722–9. sity revealed by single-nucleus RNA sequencing of the hu- 33. La Manno G, Soldatov R, Hochgerner H, et al. RNA velocity in single cells. bioRxiv 2017;p. 206052. man brain. Science 2016;352(6293):1586–90. 12. Habib N, Avraham-Davidi I, Basu A, et al. Massively paral- 34. Lake BB, Codeluppi S, Yung YC, et al. A comparative strat- lel single-nucleus RNA-seq with DroNc-seq. Nat Methods egy for single-nucleus and single-cell transcriptomes con- 2017;14(10):955–8. firms accuracy in predicted cell-type expression from nu- 13. Macosko EZ, Basu A, Satija R, et al. Highly parallel genome- clear RNA. Sci Rep 2017;7(1):6031. wide expression profiling of individual cells using nanoliter 35. Satija R, Farrell JA, Gennert D, et al. Spatial reconstruc- droplets. Cell 2015;161(5):1202–14. tion of single-cell gene expression data. Nat Biotechnol 14. Svensson V, Natarajan KN, Ly LH, et al. Power analysis 2015;33(5):495–502. of single-cell RNA-sequencing experiments. Nat Methods 36. Butler A, Satija R. Integrated analysis of single cell tran- 2017;14(4):381–7. scriptomic data across conditions, technologies, and species. 15. Hashimshony T, Senderovich N, Avital G, et al. CEL-Seq2: bioRxiv 2017;p. 164889. sensitive highly-multiplexed single-cell RNA-seq. Genome 37. Tasic B, Menon V, Nguyen TN, et al. Adult mouse cortical cell Biol 2016;17(1):77. taxonomy revealed by single cell transcriptomics. Nat Neu- rosci 2016;19(2):335–46. 16. Petukhov V, Guo J, Baryawno N, et al. Accurate estimation of molecular counts in droplet-based single-cell RNA-seq ex- 38. The Tabula Muris Consortium, Quake SR, Wyss-Coray T, et al. periments. bioRxiv 2017;p. 171496. Single-cell transcriptomic characterization of 20 organs and 17. Soumillon M, Cacchiarelli D, Semrau S, et al. Characteriza- tissues from individual mice creates a Tabula Muris. bioRxiv tion of directed differentiation by high-throughput single- 2018;p. 237446. cell RNA-seq. bioRxiv 2014. 39. Han X, Wang R, Zhou Y, et al. Mapping the mouse cell atlas 18. Jaitin DA, Kenigsberg E, Keren-Shaul H, et al. Massively par- by microwell-seq. Cell 2018;172(5):1091–1107.e17. allel single-cell RNA-seq for marker-free decomposition of 40. Bagnoli JW, Ziegenhain C, Janjic A, et al. mcSCRB-seq: sen- tissues into cell types. Science 2014;343(6172):776–9. sitive and powerful single-cell RNA sequencing. bioRxiv 19. Klein AM, Mazutis L, Akartuna I, et al. Droplet barcoding for 2017;p. 188367. single-cell transcriptomics applied to embryonic stem cells. 41. Broad Institute Single Cell Portal. https://portals.broadinsti tute.org/single cell/study/dronc-seq-single-nucleus-rna-s Cell 2015;161(5):1187–1201. 20. Zilionis R, Nainys J, Veres A, et al. Single-cell barcod- eq-on-mouse-archived-brain. Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy059/5005022 by Ed 'DeepDyve' Gillespie user on 21 June 2018 Parekh et al. 9 42. Law CW, Chen Y, Shi W, et al. voom: precision weights un- .org/10.5524/100447. lock linear model analysis tools for RNA-seq read counts. 46. Grun ¨ D, Kester L, van Oudenaarden A. Validation of Genome Biol 2014;15(2):R29. noise models for single-cell transcriptomics. Nat Methods 43. Lun ATL, McCarthy DJ, Marioni JC. A step-by-step workflow 2014;11(6):637–40. for low-level analysis of single-cell RNA-seq data with Bio- 47. Tian L, Su S, Amann-Zalcenstein D, et al. scPipe: a flexible conductor. F1000Res 2016;5. data preprocessing pipeline for single-cell RNA-sequencing 44. Soneson C, Robinson MD. Bias, robustness and scalability data. bioRxiv 2017;p. 175927. in single-cell differential expression analysis. Nat Methods 48. Islam S, Zeisel A, Joost S, et al. Quantitative single-cell 2018;15(4):255–61. RNA-seq with unique molecular identifiers. Nat Methods 45. Parekh S, Ziegenhain C, Vieth B, et al. Supporting data for 2014;11(2):163–6. ‘zUMIs - A fast and flexible pipeline to process RNA sequenc- ing data with UMIs.’ GigaScience Database 2018;http://dx.doi Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy059/5005022 by Ed 'DeepDyve' Gillespie user on 21 June 2018 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png GigaScience Oxford University Press

zUMIs - A fast and flexible pipeline to process RNA sequencing data with UMIs

Free
9 pages

Loading next page...
 
/lp/ou_press/zumis-a-fast-and-flexible-pipeline-to-process-rna-sequencing-data-with-0X2PnZdO55
Publisher
BGI
Copyright
© The Author(s) 2018. Published by Oxford University Press.
eISSN
2047-217X
D.O.I.
10.1093/gigascience/giy059
Publisher site
See Article on Publisher Site

Abstract

Background: Single-cell RNA-sequencing (scRNA-seq) experiments typically analyze hundreds or thousands of cells after amplification of the cDNA. The high throughput is made possible by the early introduction of sample-specific bar codes (BCs), and the amplification bias is alleviated by unique molecular identifiers (UMIs). Thus, the ideal analysis pipeline for scRNA-seq data needs to efficiently tabulate reads according to both BC and UMI. Findings: zUMIs is a pipeline that can handle both known and random BCs and also efficiently collapse UMIs, either just for exon mapping reads or for both exon and intron mapping reads. If BC annotation is missing, zUMIs can accurately detect intact cells from the distribution of sequencing reads. Another unique feature of zUMIs is the adaptive downsampling function that facilitates dealing with hugely varying library sizes but also allows the user to evaluate whether the library has been sequenced to saturation. To illustrate the utility of zUMIs, we analyzed a single-nucleus RNA-seq dataset and show that more than 35% of all reads map to introns. Also, we show that these intronic reads are informative about expression levels, significantly increasing the number of detected genes and improving the cluster resolution. Conclusions: zUMIs flexibility makes if possible to accommodate data generated with any of the major scRNA-seq protocols that use BCs and UMIs and is the most feature-rich, fast, and user-friendly pipeline to process such scRNA-seq data. Keywords: single-cell RNA-sequencing; digital gene expression; unique molecular identifiers ; pipeline porate unique molecular identifiers (UMIs) to label individual Introduction cDNA molecules with a random nucleotide sequence before am- The recent development of increasingly sensitive protocols al- plification [ 7]. This enables the computational removal of am- lows for the generation of RNA-sequencing (RNA-seq) libraries plification noise and thus increases the power to detect expres- of single cells [1]. The throughput of such single-cell RNA-seq sion differences between cells [8, 9]. To increase the throughput, (scRNA-seq) protocols is rapidly increasing, enabling the pro- many protocols also incorporate sample-specific bar codes (BCs) filing of tens of thousands of cells [ 2, 3] and opening exciting to label all cDNA molecules of a single cell with a nucleotide se- possibilities to analyze cellular identities [4, 5]. As the required quence before library generation [10]. This allows for early pool- amplification from such small starting amounts introduces sub- ing, which further decreases amplification noise [ 6]. Addition- stantial amounts of noise [6], many scRNA-seq protocols incor- Received: 17 October 2017; Revised: 16 March 2018; Accepted: 15 May 2018 The Author(s) 2018. Published by Oxford University Press. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy059/5005022 by Ed 'DeepDyve' Gillespie user on 21 June 2018 2 zUMIs - RNA-seq with UMIs Table 1: Features of available UMI pipelines for the quantification of gene expression data. Compatible Open Quality UMI BC Down- UMI library Name Reference source filter collapsing Mapper detection Intron sampling protocols Cell Ranger [2] yes BC+UMI Hamming STAR A no yes [2] distance CEL-seq [15] yes BC+UMI Identity bowtie2 WL no no [15, 46] only dropEst [16] yes BC Frequency- TopHat2 or WL,top- yes no [2, 13, 19] based Kallisto n,EM Drop-seq- [13] no BC+UMI Hamming STAR WL,top-n no no [13, 15, 17] tools distance scPipe [47] yes BC+UMI Hamming subread WL,top-n no no [13, 17, 18, 46] distance umis [14] yes BC Frequency- Kallisto WL,top- no no [2, 13, 17–19, based n,EM 46, 48] UMI-tools [25] yes BC+UMI Network- BWA WL no no [17, 19] based zUMIs This work yes BC+UMI Hamming STAR A,WL,top-n yes yes [2, 3, 12, 13, distance 15, 17, 18, 21, 46, 48] We consider whether the pipeline is open source, has sequence quality filters for cell BCs and UMIs, mappers, UMI-collapsing options, options for BC de tection (A, automatically infer intact BCs; WL, extract only the given list of known BCs; top-n, order BCs according the number of reads and keep the top n BCs; EM, merge BCs with given edit distance), whether it can count intron mapping reads, whether it offers a utility to make varying library sizes more comparable via downsampling, and finally with which RNA-seq library preparation protocols is it compatible ally, for cell types such as primary neurons, it has been proven Implementation and Operation to be more feasible to isolate RNA from single nuclei rather than Filtering and mapping whole cells [11, 12]. This decreases mRNA amounts further so that it has been suggested to count intron mapping reads origi- The first step in our pipeline is to filter reads that have low- nating from nascent RNAs as part of single-cell expression pro- quality BCs according to a user-defined threshold (Fig. 1). This files [ 11]. However, the few bioinformatic tools that process RNA- step eliminates the majority of spurious BCs and thus greatly seq data with UMIs and BCs have limitations (Table 1). For ex- reduces the number of BCs that need to be considered for count- ample, the Drop-seq-tools is not an open source [13]. While Cell ing. Similarly, we also filter low-quality UMIs. Ranger is open, it is exceedingly difficult to adapt the code to The remaining reads are then mapped to the genome using new or unknown sample BCs and other library types. Other tools the splice-aware aligner STAR [22]. The user is free to customize are specifically designed to work with one mapping algorithm mapping by using the options of STAR. Furthermore, if the user and focus mainly on transcriptome references [14, 15]. Further- wishes to use a different mapper, it is also possible to provide more, the only other UMI-RNA-seq pipeline providing the utility zUMIs with an aligned bam file instead of the fastq file with the to also consider intron mapping reads, dropEst [16], is only appli- cDNA sequence, with the sole requirement that only one map- cable to droplet-based protocols. Here, we present zUMIs, a fast ping position per read is reported in the bam file. and flexible pipeline that overcomes these limitations. Transcript counting Findings Next, reads are assigned to genes. In order to distinguish exon zUMIs is a pipeline to process RNA-seq data that were mul- and intron counts, we generate two mutually exclusive an- tiplexed using cell BCs and also contain UMIs. Read-pairs are notation files from the provided gtf, one detailing exon posi- filtered to remove reads with low-quality BCs or UMIs based tions, the other introns. Based on those annotations, Rsubread on sequence and then mapped to a reference genome (Fig.1). featureCounts [23] is used to first assign reads to exons and af- Next, zUMIs generates UMI and read count tables for exon and terward to check whether the remaining reads fall into introns, exon+intron counting. We reason that very low input material in other words, if a read is overlapping with intronic and ex- such as from single nuclei sequencing might profit from in- onic sequences, it will be assigned to the exon only. The output cluding reads that potentially originate from nascent RNAs. An- is then read into R using data.table [24], generating count ta- other unique feature of zUMIs is that it allows for downsampling bles for UMIs and reads per gene per BC. We then collapse UMIs of reads before collapsing UMIs, thus enabling the user to as- that were mapped either to the exon or intron of the same gene. sess whether a library was sequenced to saturation or whether Note that only the processing of intron and exon reads together deeper sequencing is necessary to depict the full mRNA com- allows for properly collapse of UMIs that can be sampled from plexity. Furthermore, zUMIs is flexible with respect to the length the intronic as well as from the exonic part of the same nascent and sequences of the BCs and UMIs, supporting protocols that mRNA molecule. have both sequences in one read [2, 3, 12, 13, 15, 17, 18]aswell Per default, we only collapse UMIs by sequence identity. If as protocols that provide UMI and BC in separate reads [19–21]. there is a risk that a large proportion of UMIs remains under- This makes zUMIs the only tool that is easily compatible with all collapsed due to sequence errors, zUMIs provides the option to major UMI-based scRNA-seq protocols. collapse UMIs within a given Hamming distance. We compare Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy059/5005022 by Ed 'DeepDyve' Gillespie user on 21 June 2018 Parekh et al. 3 Figure 1: Schematic of the zUMIs pipeline. Each of the gray panels from left to right depicts a step of the zUMIs pipeline. First, fastq files are filtered according to user-defined bar code (BC) and unique molecular identifier (UMI) quality thresholds. Next, the remaining cDNA reads are mapped to the reference genome using STAR. Gene-wise read and UMI count tables are generated for exon, intron, and exon+intron overlapping reads. To obtain comparable library sizes, reads can be downsampled to a desired range during the counting step. In addition, zUMIs also generates data and plots for several quality measures, such as the number of detected genes/UMIs per BCe and distribution of reads into mapping feature categories. (A) (B) Drop−seq−tools vs zUMIs 1 < 17 UMI−tools vs zUMIs 1 < 17 UMI−tools vs zUMIs Hamming1 1.000 0.997 0.995 0.992 0.990 Drop−seq−tools UMI−tools zUMIs 1e+02 (C) (D) zUMIs − 1 < 17 zUMIs − Hamming1 UMI−tools − edit1 directional Drop−seq−tools − 1 < 17 1e+00 1e−02 Unfiltered zUMIs − 1 < 17 zUMIs − Hamming1 1e−04 Drop−seq−tools − 1 < 17 UMI−tools − edit1 directional 50000 100000 150000 200000 UMIs unfiltered Figure 2: Comparison of different UMI collapsing methods. We compared Drop-seq-tools and UMI-tools with zUMIs using our HEK dataset (227 mio reads). (A) Run time to count exonic UMIs using zUMIs (hamming distance = 0), UMI-tools (”unique” mode) and Drop-seq-tools (edit distance = 0). (B) Box plots of correlation coefficients of gene-wise UMI counts of the same cell generated with different methods. UMI counts generated using zUMIs (quality filter “1 base under phred 17” or ha mming distance = 1) were correlated to UMI counts generated using Drop-seq-tools (quality filter “1 base under phred 17” ) and UMI-tools (“directional adjacency” mode ). (C) Comparison of the total number of UMIs per cell derived from different counting methods to “unfiltered” counts. (D) Violin plots of gene-wise dispersion estimates with different quality filtering and UMI collapsing methods. Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy059/5005022 by Ed 'DeepDyve' Gillespie user on 21 June 2018 UMIs filtered Time taken (minutes) Dispersion Pearson's R 4 zUMIs - RNA-seq with UMIs the two zUMIs UMI-collapsing options to the recommended di- or UMIs enough to justify the extra cost. In our HEK-cell exam- rectional adjacency approach implemented in UMI-tools [25]us- ple dataset, the number of detected genes starts leveling off at ing our in-house example dataset (see Methods section). zUMIs 1 million reads. Sequencing double that amount would only in- identity collapsing yields nearly identical UMI counts per cell crease the number of detected genes from 9,000 to 10,600 when as UMI-tools, while Hamming distance yields increasingly fewer counting exon reads (Fig.3D). In line with previous findings [ 8, UMIs per cell with increasing sequencing depth (Fig.2C). Smith 14], the saturation curve of exon+intron counting runs parallel et al [25] suggest that edit distance collapsing without consid- to the one for exon counting, both indicating that a sequencing ering the relative frequencies of UMIs might indeed overreach depth of 1 million reads per cell is sufficient for these libraries. and overcollapse the UMIs. We suspect that this is indeed what happens in our example data, where we find that gene-wise dis- Output and statistics persion estimates appear suspiciously truncated as expected if several counts are unduly reduce to one, the minimal number zUMIs outputs three UMI and three read count tables: gene-wise after collapsing (Fig.2D). counts for traditional exon counting, one for intron and one for However, note that the above-described differences are mi- exon+intron counts. If a user chooses the downsampling option, nor. By and large, there is good agreement between UMI counts six additional count tables per target read count are provided. To obtained by UMI-tools [25], the Drop-seq pipeline [13], and zU- evaluate library quality, zUMIs summarizes the mapping statis- MIs. The correlation between gene-wise counts of the same cell tics of the reads. While exon and intron mapping reads likely is >0.99 for all comparisons (Fig. 2B). In light of this, we consider represent mRNA quantities, a high fraction of intergenic and un- the >3 times higher processing speed of zUMIs to be a decisive mapped reads indicates low-quality libraries. Another measure advantage (Fig.2A). of RNA-seq library quality is the complexity of the library, for which the number of detected genes and the number of iden- tified UMIs are good measures (Fig. 1). We processed 227 mil- Cell BC selection lion reads with zUMIs and quantified expression levels for exon In order to be compatible with well-based and droplet-based and intron counts on a Unix machine using up to 16 threads, scRNA-seq methods, zUMIs needs to be able to deal with known which took less than 3 hours. Increasing the number of reads as well as random BCs. As default behavior, zUMIs infers which increases the processing time approximately linearly, where fil- BCs mark good cells from the data (Fig.3A, 3B). To this end, we tering, mapping, and counting each take up roughly one third of fit a k-dimensional multivariate normal distribution using the the total time (Fig.3E). We also observed that the peak random R-package mclust [26, 27] for the number of reads/BC, where k access memory usage for processing datasets of 227, 500, and is empirically determined by mclust via the Bayesian informa- 1,000 million pairs was 42 Gb, 89 Gb, and 172 Gb, respectively. tion criterion. We reason that only the kth normal distribution Finally, zUMIs could process the largest scRNA-seq dataset re- with the largest mean contains BCs that identify reads originat- ported to date with around 1.3 million brain cells and 30 billion ing from intact cells. We exclude all BCs that fall in the lower read-pairs generated with 10xGenomics Chromium (see Meth- 1% tail of this kth normal distribution to exclude spurious BCs. ods section) on a 22-core processor in only 7 days. The HEK dataset used here contains 96 cells with known BCs and zUMIs identifies 99 BCs as intact, including all the 96 known Intron counting BCs. Also, for the single-nucleus RNA-seq from Habib et al. [12], zUMIs identified a reasonable number of cells; Habib et al. report Recently, it has been shown that intron mapping reads in RNA- 10,877 nuclei and zUMIs identified 11,013 intact nuclei. However, seq likely originate from nascent mRNAs and are useful for gene we recommend to always check the elbow plot generated by zU- expression estimates [31, 32]. Additionally, novel approaches MIs (Fig.3B) to confirm that the cutoff used by zUMIs is valid for a leverage the ratios of intron and exon mapping reads to infer given dataset. In cases where the number of BCs or BC sequences information on transcription dynamics and cell states [33]. To are known, it is preferable to use this information. If zUMIs is ei- address this new aspect of analysis, zUMIs also counts and col- ther given the number of expected BCs or is provided with a list lapses intron-only mapping reads as well as intron and exon of BC sequences, it will use this information and forgo automatic mapping reads from the same gene with the same UMI. To as- inference. sess the information gain from intronic reads to estimate gene expression levels, we analyzed a publicly available DroNc-seq dataset from mouse brain ([12]; see Methods section). For the Downsampling ∼11,000 single nuclei of this dataset, the fraction of intron map- scRNA-seq library sizes can vary by orders of magnitude, which ping reads of all reads goes up to 61%. Thus, if intronic reads complicates normalization [28, 29]. A straight-forward solution are considered, the mean number of detected genes per cell for this issue is to downsample overrepresented libraries [30]. increases from 1,041 for exon counts to 1,995 for exon+intron zUMIs has an built-in function for downsampling datasets to a counts. Next, we used the resulting UMI count tables to investi- user-specified number of reads or a range of reads. By default, gate whether exon+intron counting improves the identification zUMIs downsamples all selected BCs to be within three absolute of cell types, as suggested by Lake et al. [11]. The validity and ac- deviations from the median number of reads per BC (Fig.3C). Al- curacy of counting introns for single-nucleus sequencing meth- ternatively, the user can provide a target sequencing depth, and ods has recently been demonstrated [34]. Following the Seurat zUMIs will downsample to the specified read number or omit pipeline to cluster cells [35, 36], we find that using exon+intron the cell from the downsampled count table if fewer reads were counts discriminates 28 clusters, while we could only discrimi- present. Furthermore, zUMIs also allows the user to specify a nate 19 clusters using exon counts (Fig.4A, 4B). The larger num- multiple target read number at once for downsampling. This ber of clusters is not simply due to the increase in the counted feature is helpful if the user wishes to determine whether the UMIs and genes. When we permute the intron counts across RNA-seq library was sequenced to saturation or whether fur- cells and add them to the exon counts, the added noise actu- ther sequencing would increase the number of detected genes ally reduces the number of identifiable clusters (Fig. 4E). Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy059/5005022 by Ed 'DeepDyve' Gillespie user on 21 June 2018 Parekh et al. 5 1.5e+07 (A) (B) (C) 2.0e+08 0.9 1.0e+07 1.5e+08 0.6 1.0e+08 0.3 5.0e+06 5.0e+07 0.0 0.0e+00 0.0e+00 0 1000 2000 3000 4000 log10(Reads) Cell Index Cells (D) (E) Counting Filtering Mapping Summarizing total Exon Exon+Intron 0e+00 1e+06 2e+06 3e+06 200 400 600 800 1000 Sequencing Depth Million Read pairs Figure 3: Utilities of zUMIs. Each of the panels shows the utilities of zUMIs pipeline. The plots from A–D show the results from the example HEK dataset used here. (A) The plot shows a density distribution of reads per BC. Cell BCs with reads right of the blue line are selected. (B) The plot shows the cumulative read distribution in the example HEK dataset where the BCs in light blue are the selected cells. (C) The bar plot shows the number of reads per selected cell BC with the red lines showing upper and lower median absolute deviation (MAD) cutoffs for adaptive downsampling. Here, the cells below the lower MAD have very low coverage and are discarded in downsampled count tables. (D) Cells were downsampled to six depths from 100,000 to 3,000,000 reads. For each sequencing depth, the genes detected per cell are shown. (E) Runtime for three datasets with 227, 500, and 1,000 million read-pairs. The runtime is divided in the main steps of the zUMIs pipeline as follows: filtering, mapping, counting, and summarizing. Each dataset was processed using 16 threads (“-p 16”). We continue to further characterize the seven clusters that neurons [37], suggesting that the split of cluster 7 might be bio- were subdivided by the addition of intron counts (Fig.4D). First, logically meaningful (Fig.4E). we identify DE genes between the newly formed clusters. If we In order to evaluate the power gained by exon+intron count- count only exon reads, there appear to be, on average, only 10 ing in a more systematic way, we perform power simulations us- DE genes between the subgroups, while exon+intron counting ing empirical mean and dispersion distributions from the largest yields ∼10 times more DE genes, thus corroborating the signal and most uniform cluster (∼1,500 cells) [9]. For a fair comparison, found with clustering. The log2-fold changes of those additional we include all detected genes, which is equivalent to the num- DE genes estimated with either counting strategy are generally ber of genes detected with exon+intron counting. Also, since we in good agreement; especially large log2-fold changes are de- call a gene detected as soon as one count is associated, exon tected with both exon and exon+intron counting (Fig.4F). Genes counting is necessarily a subset of exon+intron. Thus, there are, that are detected as DE in only one of our counting strategies on average, 4 times more genes in the lowest expression quan- have small log2-fold changes, and there are more of these small tile for exon counting than for exon+intron counting (Fig.4H). For changes detected using exon+intron counting. those genes, expression is too spurious to be used for differen- Detecting more genes naturally increases the chance to also tial expression analysis; for exon+intron counting, we have, on detect more informative genes. Here, we cross-reference the average, 60% power to detect a DE gene in the first mean expres- gene list with marker genes for transcriptomic subtypes de- sion bin with a well-controlled false discovery rate (FDR) (Fig.4G). tected for major cell types of the mouse brain [37] and find that In summary, the increased power for exon+intron counting and ∼5% of the additional genes are also marker genes, which cor- probably also the larger number of clusters are due to better de- responds well to the general frequency of marker genes among tection of lowly expressed genes. Furthermore, we think that al- the detected genes (4%). In the same vein, we also detect propor- though potentially noisy, the large number of additionally de- tionally more DE genes with exon+intron counting compared to tected genes makes exon+intron counting worthwhile, espe- exon counting. Thus, including introns simply allows us to bet- cially for single-nuclei sequencing techniques that are enriched ter detect present transcripts, while leaving the proportions of for nuclear nascent RNA transcripts, such as DroNc-seq [12]. Ad- interest unaltered. Having a closer look at cluster 7, it was split ditionally, exon+intron counting may help in extracting as much into a bigger (7) and a smaller cluster (24) using exon+intron information as possible from low coverage data as generated counting (Fig.4A-C), we find one marker gene (Il1rapl2) to be in the context of high-throughput cell atlas efforts (e.g., 10,000– DE between the subclusters using exon+intron counting, while 20,000 reads/cell [38, 39]. Last, users should always exclude the Il1rapl2 had only spurious counts using exon counts. Il1rapl2 is possibility of intronic reads stemming from genomic DNA con- a marker for transcriptomic subtypes of GABAergic Pvalb-type tamination in the library preparation by confirming low inter- Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy059/5005022 by Ed 'DeepDyve' Gillespie user on 21 June 2018 Density Genes Detected Cumulative Reads Runtime (hours) Number of Reads 2 7 6 1 5 8 3 Exon+ Exon Intron 6 zUMIs - RNA-seq with UMIs (A) (B) (C) (D) Exon Exon+Intron 60 27 exPFC 7 GABA 7 CA 19 6 DG SMC 14 13 1 0 OPC 7777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777 7 ODC MG 9 ASC −25 −30 8 8 8 END 0.00 0.25 0.50 0.75 1.00 10 1000 −50 Cell Type Fraction −40 −20 0 20 40 −50 −25 0 25 50 Number of DE Genes (E) (F) Exon+fake−Intron reads Exon Exon+Intron (H) 7 2000 20 01234 Il1rapl2− normalized counts (I) 0.75 (G) 0.50 0.4 Both 0.25 Exon 0.0 0.00 Exon+Intron 0.15 −0.4 0.10 −0.8 −20 0.05 −1.2 0.00 −2 −1 0 Log2−Fold Change Exon only −20 0 20 strata of mean−expression Figure 4: Contribution of intron reads to biological insights. We analyzed published single-nucleus RNA-seq data from mouse prefrontal cortex (PFC) and hippocampus [12] to assess the utility of counting intron in addition to exon reads. We processed the raw data with zUMIs to obtain expression tables with exon reads as well as exon+intron reads and then used the R-package Seurat [35, 36] to cluster cells. With exon counts, we identified 19 clusters (A), and with exon+intron counts we identified 27 clusters (B). Clusters are represented as t-SNE plots and colored according to the most frequent cell-type assignment in the original article [12]: glutamatergic neurons from the prefrontal cortex (exPFC), GABAergic interneurons (GABA), pyramidal neurons from the hippocampal CA region (CA), granule neurons from the hippocampal dentate gyrus region (DG), astrocytes (ASC), microglia (MG), oligodendrocytes (ODC), oligodendrocyte precursor cells (OPC), neuronal stem cells (NSC), smooth muscle cells (SMC) and endothelial cells (END). Different shades of those clusters indicate that multiple clusters had the same major cell type assigned. If we randomly sample counts from the intron data and add them to the exon counting, the noise reduces the number of clusters and the Seurat pipeline can only identify 9–11 clusters (E). The composition of each cluster based on exon+intron is detailed in panel (C), and cells that were not assigned a cell type in [12] are displayed as empty. The boxes mark the clusters that were not split when using exon data only. For example, cluster 7 from exon counting, which mainly consists of GABAergic neurons, was split into clusters 7, 24 (506, 66 cells) when using exon+intron counting. In (D), we show the numbers of genes that were differentially expressed (DE) (limma p-adj <0.05) between the clusters found only with exon+intron counts. The panel numbers represent the exon counting cluster numbers and the y-axis the exon+intron counting cluster number. The log2-fold changes corresponding to these contrasts are also used in (G). Among the genes that were additionally detected to be DE by exon+intron −5 counting was the marker gene Il1rapl2 (limma p-adj = 10 ). In (F), we present a violin plot of the normalized counts for Il1rapl2 in cells of the GABAergic subclusters 7 and 24. Log2-fold changes calculated with exon+intron counts correlate well with exon counts (G). Note that for exon counting only, half as many genes could be evaluated as for exon+intron counting and thus only half of the exon+intron genes are depicted in (G). Large log fold changes (LFCs) are found to be significant with both counting strategies (purple points are close to the bisecting line). We conducted simulations based on mean and dispersion measured using exon cluster 0 (1,616 cells, ∼90% exPFC). In (I) we show the expected true positive rate and the false discovery rate for a scenario comparing 300 vs 300 cells. Results for exon and exon+intron counting were stratified into five quantiles according to the mean expression of genes, where stratum 1 contains lowly expressed genes and stratum 5 the most highly expressed genes. The numbers of genes falling into each of the bins using exon+intron and exon counting are depicted in (H). genic mapping fractions using the statistics output provided by Analyzed RNA-seq datasets zUMIs. HEK293T cells were cultured in DMEM high glucose with L- glutamine (Biowest) supplemented with 10% fetal bovine serum (Thermo Fisher) and 1% penicillin/streptomycin (Sigma-Aldrich) Conclusion in a 37 C incubator with 5% carbon dioxide. Cells were passaged zUMIs is a fast and flexible pipeline for processing raw reads to and split every 2 or 3 days. For single-cell RNA-seq, HEK293T cells were dissociated by incubation with 0.25% Trypsin (Sigma- obtain count tables for RNA-seq data using UMIs. To our knowl- edge, it is the only open source pipeline that has a BC and Aldrich) for 5 minutes at 37 C. The single-cell suspension was washed twice with phosphate-buffered saline, and dead cells UMI quality filter, allows intron counting, and has an integrated downsampling functionality. These features ensure that zUMIs were stained with Zombie Yellow (Biolegend) according to the is applicable to most experimental designs of RNA-seq data, manufacturer’s protocol. Single cells were sorted into DNA including single-nucleus sequencing techniques, droplet-based LoBind 96-well polymerase chain reaction (PCR) plates (Eppen- methods where the BC is unknown, as well as plate-based UMI- dorf) containing lysis buffer with a Sony SH-800 cell sorter in 3- methods with known BCs. Finally, zUMIs is computationally ef- drop purity mode using a 100-μmnozzle. Next, single-cell RNA- ficient, user-friendly, and easy to install. seq libraries were constructed from one 96-well plate using a slightly modified version of the mcSCRB-seq protocol. Reverse transcription was performed as described previously [40], with the only change being the use of KAPA HiFi HotStart enzyme for Methods Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy059/5005022 by Ed 'DeepDyve' Gillespie user on 21 June 2018 Log2−Fold Change Exon+Intron Exon+Intron cluster ID No. of FDR TPR Cluster ID (Exon+Intron) Genes Parekh et al. 7 PCR amplification of cDNA. Resulting libraries were sequenced ters. To illustrate that the additional clusters found by count- using an Illumina HiSeq1500 with 16 cycles in Read 1 to decode ing exon+intron reads are not spurious, we use intron-only UMI cell BCs (6 bases) and UMIs (10 bases) and 50 cycles in Read 2 to counts from the same data to add to the observed exon-only sequence into the cDNA fragment, obtaining ∼227 million reads. counts. More specifically, to each gene we add scran-size factor- Raw fastq files were processed using zUMIs, mapping to the hu- corrected intron counts from the same gene after permuting man genome (hg38) and Ensembl gene models (GRCh38.84). them across cells. We assessed the cluster numbers from 100 In addition, we analyzed data from 1.3 million mouse brain such permutations. cells generated on the 10xGenomics Chromium platform [2]. Se- quences were downloaded from the National Center for Biotech- Comparison of UMI collapsing strategies nology Information Sequence Read Archive under accession number SRP096558. The data consist of 30 billion read-pairs In order to validate zUMIs and compare different UMI collaps- from 133 individual samples. In these data, read 1 contains 16 bp ing methods, we used the HEK dataset described above. We ran for the cell BC and 10 bp for the UMI and read 2 contains 114 bp zUMIs (1) without quality filtering, (2) filtering for onebase un- of cDNA. zUMIs was run using default settings, and we allowed der Phred 17, and (3) collapsing similar UMI sequences within a 7 threads per job for a total of up to 42 threads on an Intel Xeon hamming distance of 1. To compare with other available tools, E5-2699 22-core processor. we ran the same dataset using the Drop-seq-tools version 1.13 Finally, we obtained mouse brain DroNc-seq read data [12] [13] and quality filter “1 base under Phred 17” without edit dis- from the Broad Institute Single Cell Portal [41]. This dataset con- tance collapsing. Last, the HEK dataset was used with UMI tools sists of ∼1,615 million read-pairs from ∼11,000 single nuclei. [25] in (1) “unique” and (2) “directional adjacency” mode with Read 1 contains a 12 bpcell BC and a 8 bpUMI and read 2 60 bpof edit distance set to 1. Also, we compared the output of zU- cDNA. MIs from the DroNc-seq dataset when using default parame- The two mouse datasets were mapped to genome version ters (“1 base under Phred 20”) to UMI-tools in (1) “unique,” (2) mm10 and applying Ensembl gene models (GRCm38.75). “directional adjacency,” and (3) “cluster” settings. For each set- ting and tool combination, we compared per-cell/per-nuclei UMI contents in a linear model fit. Power simulations and DE analysis We evaluated the power to detect differential expression with thehelpofthe powsimRpackage [9]. For the DroNc-seq dataset, Availability of source code and requirements we estimated the parameters of the negative binomial dis- Project name: zUMIs tribution from one of the identified clusters, namely, cluster Project home page: https://github.com/sdparekh/zUMIs 0, compromising 1,500 glutamatergic neuronal cells from the Operating system(s): UNIX prefrontal cortex (Fig.4D). Since we detect more genes with Programming language: shell, R, perl exon+intron counting (4,433 compared to 1,782), we included Other requirements: STAR >= 2.5.3a, R >= 3.4, Rsubread >= this phenomenon in our read count simulation by drawing mean 1.26.1, pigz >= 2.3 & samtools >= 1.1 expression values for a total of 4,433 genes. This means that License: GNU GPLv3.0 the table includes sparse counts for the exon counting. Log - Research Resource Identification Initiative ID: SCR 016139 fold changes were drawn from a gamma distribution with shape equal to 1 and scale equal to 2. In each of the 25 simulation it- erations, we draw an equal sample size of 300 cells per group Availability of supporting data and test for differential expression using limma-trend [42]on log counts per million (CPM) values with scran [43] library size All data that were generated for this project were submitted to correction. The true positive rate and FDR are stratified over the GEO under accession GSE99822. An archival copy of the source empirical mean expression quantile bins. code and test data are available via the GigaScience repository For the differential expression analysis between clusters, we GigaDB [45]. use the same DE estimation procedure as in the simulations: scran normalization followed by limma-trend DE-analysis (c.f. [44]). Abbreviations BC: barcode; DE: differentially expressed; FDR: false discovery Cluster identification rate; MAD: median absolute deviation; PCR: polymerase chain reaction; PFC: prefrontal cortex; scRNA-seq: single-cell RNA se- After processing the DroNc-seq data [12] with zUMIs as de- quencing; UMI: unique molecular identifier. scribed above, we cluster cells based on UMI counts derived from exons only and exons+introns reads using the Seurat pipeline [35, 36]. First, cells with fewer than 200 detected genes were fil- Competing interests tered out. The filtered data were normalized using the LogNor- malize function. We then scale the data by regressing out the The authors declare that they have no competing interests. effects of the number of transcripts and genes detected per cell using the ScaleData function. The normalized and scaled data are then used to identify the most variable genes by fitting a Funding relationship between mean expression (ExpMean) and disper- This work has been supported by the Deutsche Forschungsge- sion (LogVMR) using the FindVariableGenes function. The iden- tified variable genes are used for principle component analysis, meinschaft (DFG) through SFB1243 subprojects A14/A15. and the top 20 principle components are then used to find clus- ters using graph-based clustering as implemented in FindClus- Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy059/5005022 by Ed 'DeepDyve' Gillespie user on 21 June 2018 8 zUMIs - RNA-seq with UMIs ing and sequencing using droplet microfluidics. Nat Protoc Author contributions 2017;12(1):44–73. S.P. and C.Z. designed and implemented the pipeline. B.V. tested 21. Hochgerner H, Lonnerber ¨ g P, Hodge R, et al. STRT-seq-2i: the pipeline and helped in power simulations. All authors con- dual-index 5’ single cell and nucleus RNA-seq on an address- tributed to writing the manuscript. able microwell array. Sci Rep 2017;7(1):16327. 22. Dobin A, Davis CA, Schlesinger F, et al. STAR: ultrafast uni- versal RNA-seq aligner. Bioinformatics 2013;29(1):15–21. References 23. Liao Y, Smyth GK, Shi W. featureCounts: an efficient general 1. Sandberg R. Entering the era of single-cell transcriptomics in purpose program for assigning sequence reads to genomic biology and medicine. Nat Methods 2014;11(1):22–4. features. Bioinformatics 2014;30(7):923–30. 2. Zheng GXY, Terry JM, Belgrader P, et al. Massively parallel 24. Dowle M, Srinivasan A. data.table: Extension of ‘data.frame.” digital transcriptional profiling of single cells. Nat Commun 2017, https://CRAN.R-project.org/package=data.table,r package version 1.10.4. 2017;8:14049. 3. Rosenberg AB, Roco CM, Muscat RA, et al. Single-cell profiling 25. Smith TS, Heger A, Sudbery I. UMI-tools: modelling sequenc- of the developing mouse brain and spinal cord with split- ing errors in unique molecular identifiers to improve quan- pool barcoding. Science 2018;360(6385):176–82. tification accuracy. Genome Res. 2017. 4. Wagner A, Regev A, Yosef N. Revealing the vectors of cel- 26. Fraley C, Raftery AE. Model-based clustering, discrimi- lular identity with single-cell genomics. Nat Biotechnol nant analysis, and density estimation. J Am Stat Assoc 2016;34(11):1145–60. 2002;97(458):611–31. 5. Regev A, Teichmann SA, Lander ES, et al. The Human Cell 27. Fraley C, Raftery AE Enhanced Model-Based Clustering, Atlas. Elife, 2017; 6. Density Estimation and Discriminant Analysis Software: 6. Parekh S, Ziegenhain C, Vieth B, et al. The impact of ampli- MCLUST. , J. Classification , 2003, 20, 263–286. fication on differential expression analyses by RNA-seq. Sci 28. Vallejos CA, Risso D, Scialdone A, et al. Normalizing single- Rep 2016;6:25533. cell RNA sequencing data: challenges and opportunities. Nat 7. Kivioja T, Vah ¨ ar ¨ autio A, Karlsson K, et al. Counting absolute Methods 2017;14(6):565–71. 29. Evans C, Hardin J, Stoebel DM. Selecting between-sample numbers of molecules using unique molecular identifiers. Nat Methods 2012;9(1):72–4. RNA-seq normalization methods from the perspective of 8. Ziegenhain C, Vieth B, Parekh S, et al. Comparative anal- their assumptions. Brief Bioinform 2017. ysis of single-cell RNA sequencing methods. Mol Cell 30. Grun ¨ D, van Oudenaarden A. Design and analysis of single- 2017;65(4):631–43.e4. cell sequencing experiments. Cell 2015;163(4):799–810. 9. Vieth B, Ziegenhain C, Parekh S, et al. powsimR: power anal- 31. Hendriks GJ, Gaidatzis D, Aeschimann F, et al. Extensive os- ysis for bulk and single cell RNA-seq experiments. Bioinfor- cillatory gene expression during C. elegans larval develop- matics 2017;33(21):3486–3488. ment. Mol Cell 2014;53(3):380–92. 10. Ziegenhain C, Vieth B, Parekh, et al. Quantitative 32. Gaidatzis D, Burger L, Florescu M, et al. Analysis of intronic single-cell transcriptomics. Brief Funct Genomics 2018; and exonic reads in RNA-seq data characterizes transcrip- doi:10.1093/bfgp/ely009. tional and post-transcriptional regulation. Nat Biotechnol 11. Lake BB, Ai R, Kaeser GE, et al. Neuronal subtypes and diver- 2015;33(7):722–9. sity revealed by single-nucleus RNA sequencing of the hu- 33. La Manno G, Soldatov R, Hochgerner H, et al. RNA velocity in single cells. bioRxiv 2017;p. 206052. man brain. Science 2016;352(6293):1586–90. 12. Habib N, Avraham-Davidi I, Basu A, et al. Massively paral- 34. Lake BB, Codeluppi S, Yung YC, et al. A comparative strat- lel single-nucleus RNA-seq with DroNc-seq. Nat Methods egy for single-nucleus and single-cell transcriptomes con- 2017;14(10):955–8. firms accuracy in predicted cell-type expression from nu- 13. Macosko EZ, Basu A, Satija R, et al. Highly parallel genome- clear RNA. Sci Rep 2017;7(1):6031. wide expression profiling of individual cells using nanoliter 35. Satija R, Farrell JA, Gennert D, et al. Spatial reconstruc- droplets. Cell 2015;161(5):1202–14. tion of single-cell gene expression data. Nat Biotechnol 14. Svensson V, Natarajan KN, Ly LH, et al. Power analysis 2015;33(5):495–502. of single-cell RNA-sequencing experiments. Nat Methods 36. Butler A, Satija R. Integrated analysis of single cell tran- 2017;14(4):381–7. scriptomic data across conditions, technologies, and species. 15. Hashimshony T, Senderovich N, Avital G, et al. CEL-Seq2: bioRxiv 2017;p. 164889. sensitive highly-multiplexed single-cell RNA-seq. Genome 37. Tasic B, Menon V, Nguyen TN, et al. Adult mouse cortical cell Biol 2016;17(1):77. taxonomy revealed by single cell transcriptomics. Nat Neu- rosci 2016;19(2):335–46. 16. Petukhov V, Guo J, Baryawno N, et al. Accurate estimation of molecular counts in droplet-based single-cell RNA-seq ex- 38. The Tabula Muris Consortium, Quake SR, Wyss-Coray T, et al. periments. bioRxiv 2017;p. 171496. Single-cell transcriptomic characterization of 20 organs and 17. Soumillon M, Cacchiarelli D, Semrau S, et al. Characteriza- tissues from individual mice creates a Tabula Muris. bioRxiv tion of directed differentiation by high-throughput single- 2018;p. 237446. cell RNA-seq. bioRxiv 2014. 39. Han X, Wang R, Zhou Y, et al. Mapping the mouse cell atlas 18. Jaitin DA, Kenigsberg E, Keren-Shaul H, et al. Massively par- by microwell-seq. Cell 2018;172(5):1091–1107.e17. allel single-cell RNA-seq for marker-free decomposition of 40. Bagnoli JW, Ziegenhain C, Janjic A, et al. mcSCRB-seq: sen- tissues into cell types. Science 2014;343(6172):776–9. sitive and powerful single-cell RNA sequencing. bioRxiv 19. Klein AM, Mazutis L, Akartuna I, et al. Droplet barcoding for 2017;p. 188367. single-cell transcriptomics applied to embryonic stem cells. 41. Broad Institute Single Cell Portal. https://portals.broadinsti tute.org/single cell/study/dronc-seq-single-nucleus-rna-s Cell 2015;161(5):1187–1201. 20. Zilionis R, Nainys J, Veres A, et al. Single-cell barcod- eq-on-mouse-archived-brain. Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy059/5005022 by Ed 'DeepDyve' Gillespie user on 21 June 2018 Parekh et al. 9 42. Law CW, Chen Y, Shi W, et al. voom: precision weights un- .org/10.5524/100447. lock linear model analysis tools for RNA-seq read counts. 46. Grun ¨ D, Kester L, van Oudenaarden A. Validation of Genome Biol 2014;15(2):R29. noise models for single-cell transcriptomics. Nat Methods 43. Lun ATL, McCarthy DJ, Marioni JC. A step-by-step workflow 2014;11(6):637–40. for low-level analysis of single-cell RNA-seq data with Bio- 47. Tian L, Su S, Amann-Zalcenstein D, et al. scPipe: a flexible conductor. F1000Res 2016;5. data preprocessing pipeline for single-cell RNA-sequencing 44. Soneson C, Robinson MD. Bias, robustness and scalability data. bioRxiv 2017;p. 175927. in single-cell differential expression analysis. Nat Methods 48. Islam S, Zeisel A, Joost S, et al. Quantitative single-cell 2018;15(4):255–61. RNA-seq with unique molecular identifiers. Nat Methods 45. Parekh S, Ziegenhain C, Vieth B, et al. Supporting data for 2014;11(2):163–6. ‘zUMIs - A fast and flexible pipeline to process RNA sequenc- ing data with UMIs.’ GigaScience Database 2018;http://dx.doi Downloaded from https://academic.oup.com/gigascience/article-abstract/7/6/giy059/5005022 by Ed 'DeepDyve' Gillespie user on 21 June 2018

Journal

GigaScienceOxford University Press

Published: May 26, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off