Delineating biological and technical variance in single cell expression data.

Single cell transcriptomics is becoming a common technique to unravel new biological phenomena whose functional significance can only be understood in the light of differences in gene expression between single cells. The technology is still in its early days and therefore suffers from many technical challenges. This review discusses the continuous effort to identify and systematically characterise various sources of technical variability in single cell expression data and the need to further develop experimental and computational tools and resources to help deal with it.

Introduction

Next generation sequencing (NGS) technologies have revolutionized the way of approaching molecular biology to advance our understanding of the working principles of biological systems, including identification of the building blocks. Genome sequencing is now widely used across diverse fields in biology, ranging from medicine and population studies to animal breeding. However, the information encoded in the genome is static, an ensemble of the cell’s potentialities manifest once the process of transcription is triggered. Therefore, studying the transcriptome is essential to understand how genome information is decoded in a particular cell specific context, as the cells ultimately constitute as adaptable and dynamic entities.

NGS evolved from a range of laboratory techniques developed for expression analysis over the years. Initial experimental approaches include the early Northern blotting (Alwine et al., 1977), which targets a single gene and measures its expression levels through hybridization of a labelled probe. Advances in increasing throughput of transcriptome studies came with microarrays (Schena et al., 1995), a technology that used a similar probing approach, but increased the number of quantified transcripts by using tens of thousands of probes on a chip, onto which the RNA sample is hybridized. Both approaches described above are limited by the fact that probe design requires previous knowledge of the transcript sequences. To this end, the use of sequencing technologies such as Sanger sequencing and its derivatives, including expressed sequence tag (EST), improved access to the diversity of the transcriptomic landscape by overcoming the probe design constraint. Currently, however, the most widely used application of NGS technologies to transcriptomics is RNA sequencing (RNA-Seq) (Mortazavi et al., 2008), by which -potentially- all mRNA molecules in a cell can be sequenced, and hence characterized and quantified.

The importance of RNA-Seq is not only founded in its ability to access unknown transcripts and spliced variants, but also to increase microarray’s dynamic range (i.e. the lowly expressed transcripts could be successfully detected) and sensitivity (i.e. the expression level measurements show higher accuracy). RNA-Seq has therefore become the technology of choice to provide a high-throughput and fully quantitative approach to studying the transcriptome of a broad range of species, including the ones lacking full genome sequence availability. The technology has therefore been widely applied, replacing microarrays for the analysis of gene expression profile differences among cell populations, comparative transcriptomics and disease biomarker identifications (Wang et al., 2009). However, it has become apparent that not all cells within a population behave similarly when it comes to gene expression or splicing and, in this context, bulk RNA-Seq fails to address some important questions (Sandberg, 2014).

Single cell expression technologies and applications

Over the years, single-cell approaches have been developed in combination with microscopy to visualize gene expression patterns in individual cells. For example, single-molecule RNA fluorescence in situ hybridization (RNA FISH) technology combines probe hybridization with fluorescent labelling to resolve the location of a target transcript (Lubeck and Cai, 2012). The main disadvantage of RNA FISH is that, although parallelizable, it only allows access to a limited subset of genes. The implementation of single-cell microarrays (Iscove et al., 2002) presented itself again as a high-throughput alternative to RNA FISH, and although it helps overcome this main limitation, it suffers the drawbacks of bulk microarrays. Furthermore, the limited amount of starting material and the relatively low sensitivity of microarrays enforced high levels of pre-amplification, which can introduce significant biases.

In the light of these limitations, RNA Sequencing was implemented at the single-cell level, theoretically enabling access to the transcriptome of every individual cell in a population (Ramsköld et al., 2012; Tang et al., 2010). Essentially, single-cell RNA-Seq requires the following steps: single cell isolation, mRNA capture and reverse transcription to cDNA, cDNA amplification to improve the low transcript yields rendered by single cells, and sequencing (Picelli et al., 2014).

Over the last few years, single-cell RNAseq has been proven useful to unravel biological phenomena that can only be understood in the light of differences in gene expression between single cells, including:

Studying early embryonic development: In early stages of embryonic development, only a few cells contribute to activating the molecular machinery for cell differentiation. The characterisation of transcription changes in individual inner cell mass (ICM) cells of blastocysts was proven crucial to understand the complex transition from ICMs to embryonic stem cells (ESCs) (Tang et al., 2010). This approach set a precedent for subsequent studies of later and more complex stages in the process of cell commitment and differentiation into specific lineages. In this context, a spatial-temporal profiling of gene expression in embryonic development in Caenorhabditis elegans was used to study the evolution of the germ layers. The authors noted that the gene expression program of the mesoderm is induced after those of the ectoderm and endoderm and strikingly, the endoderm gene expression program activates earlier than ectoderm expression program, a phenomenon that is conserved across many species (Hashimshony et al., 2014).

Measuring diversity in cell populations: Single cell analysis is the most powerful tool to study the diversity between individual cells treated as homogenous in a typical bulk RNA-seq experiment. It has proven potential of providing valuable insights in some of the key problems in biomedical field e.g. tumour heterogeneity, which poses substantial challenges in cancer treatment. For example, single cell analysis can unravel intra- and inter-tumour differences (Patel et al., 2014) as well as distinguishing between malignant and non-malignant cells (Tirosh et al., 2016).

Identification of new rare cell types: Complex tissues often contain previously unidentified cell types that cannot be studied using bulk RNA-Seq, as it provides only an estimate of expression influenced by the abundance of the different cell types present. Single cell transcriptomics provides a promise to address this underlying diversity in order to assess meaningful differences in phenotype. Using this strategy, authors identified and characterised a rare population of dormant neural cells which were activated upon brain injury (Llorens-Bobadilla et al., 2015). Another example is the development of a computational approach (scLVM) to identify subpopulations of cells using latent variable models to account for hidden factors such as cell cycle. Namely, different sub-populations of cells corresponding to the differentiation stages during naive T cells to T helper 2 cells were identified (Buettner et al., 2015). Identification of rare cells is of high relevance, particularly characterisation of progenitor cells to understand vertebrate development. To this end, single cell RNA-Seq has been used to unravel transcription heterogeneity and lineage commitment in myeloid progenitors, in order to further demonstrate how Cebpe deletion results into diminishing of certain myeloid lineages (Paul et al., 2015).

Mapping developmental hierarchies: transcription dynamics during development and disease can be studied in much greater details using single cell studies, as bulk RNA-seq, by averaging out signal from multiple cells, misses out on the signal from rare developmentally relevant cells. However, single cell transcriptome profiling over time is not feasible. Taking advantage of the fact that an experiment characterising hundreds of unsynchronised cells from a population typically provides a snapshot of cells at various stages during differentiation, various methods for pseudo-time inference form single cell RNA-seq data have recently been developed (Haghverdi et al., 2016; Reid and Wernisch, 2016; Trapnell et al., 2014) and reviewed (Bacher and Kendziorski, 2016). As an example of this, single cell expression data has successfully been used to reconstruct the developmental progression of cells and identify transient and terminal states together with the branching decisions (Treutlein et al., 2014).

Understanding diverse features of transcription control: Single cell transcriptomics has facilitated unravelling mechanistic details of transcription control such as kinetics and bimodality, as well as studying other features such as allelic biases and transcription networks. Even though single cell transcriptomics does not measure expression changes in one gene over time, an overall rate of transcription between individual cells can be acquired and approximately represent the stochasticity of expression of a vast number of genes, facilitating estimation of kinetics of gene expression (Kim and Marioni, 2013). Recent studies have unravelled the stochastic modes of gene expression, which were not apparent at the population level. The functional implications of this stochasticity (i.e. changes on the phenotype of seemingly identical cells) can be explained by variation in gene regulation processes across individual cells (Munsky et al., 2012). Allelic biases in gene expression have also been investigated, including stochastic allelic expression in early embryogenesis (Tang et al., 2011) as a particularly relevant example. Finally, single cell transcriptome data is successfully used to reconstruct gene regulatory networks (Moignard et al., 2015).

In summary, single cell analysis has a huge potential to bring new insights into diverse fields of biological research. In the next sections, we will put this in context by discussing the technical challenges currently faced by single cell analysis to extract the ‘biological’ or functionally relevant variability from the data, which hinder its theoretical potential.

Technical variability in single cells

Despite the promise held by the approach, single-cell RNA-Seq is not free from biases. Quite contrarily, the low availability of starting material (i.e. RNA extracted from an individual cell) introduces high technical variability, making single-cell RNA-Seq data analysis especially challenging (Stegle et al., 2015). This typically results into many missing values (technical) or true absence of expression (biological) in typically lowly expressed transcripts, and discriminating both, although important, is not currently feasible. Furthermore, the necessary amplification of starting material introduces additional biases, such as 3′ end enrichment of signal and preferential amplification of some transcripts and/or mRNA fragments. Reassuringly, bulk RNA-Seq experiments can be recapitulated in silico by pooling 30 or more single cell transcriptomes in silico (Marinov et al., 2014), and used to estimate technical variability.

The technical variability in single-cell RNA-Seq can be divided into two categories: Inter-cell variability and within cell variability (Fig. 1).

Fig. 1

Single-cell RNAseq data analysis workflow, including examples of computational methods available for each stage. Methods not developed specifically for single-cell RNAseq are marked with an asterisk. When reference is unspecified, see (Rostom et al., 2017) for full list.

Inter-cell variability

Inter-cell variability can appear as a result of the biological process under scrutiny, or can be due to unrelated phenomena, which can act as confounding factors. For example, the differences in cell cycle stage are responsible for a part of the gene expression variability between individual cells (Buettner et al., 2015), and differences in size affect the total amount of RNA extracted from individual cells (Stegle et al., 2015). In addition, cells can suffer stress or physical damage in cell capture systems used to isolate each single cell from the population, and capture sites can be occupied by more than one cell or by none. All of these events result in ‘low quality cells’, which have been estimated to be around 10–40% of the total number of cells in a dataset (Ilicic et al., 2016). The data from ‘low quality cells’ may be misleading and hinder meaningful biological interpretation. Hence, it is necessary to perform a cell filtering step in the data analysis pipeline.

Furthermore, RNA capture and amplification efficiency per cell is far from being uniform. This can be evaluated using external spike-ins as a control (Jiang et al., 2011, Stegle et al., 2015), a technique originally developed for bulk RNA-Seq. Simple measures include discarding cells with high proportion of reads mapping to mitochondrial genes or spike-in controls. As a complement to these, more complex statistical methods to eliminate confounding factors prior to data analysis have been developed. For example, variation due to differences in cell cycle stage was addressed by modelling cell-cycle variation as Gaussian processes, followed by linear regression, thus allowing the removal of noise caused by cell cycle (Buettner et al., 2015). Other quality control tools provide a model of the biological and technical features of low quality cells (Ilicic et al., 2016). The methods developed to identify sub-populations, including the popular principal component analysis (PCA), t-distributed stochastic neighbour embedding (t-SNE), or zero-inflated factor analysis (ZIFA) (Pierson and Yau, 2015) are also used to filter low quality cells. However, comprehensive quality controls with more than one metric are more effective for discarding true low quality samples, and are beginning to be implemented by R packages such as scone (McCarthy et al., 2016, Stegle et al., 2015). Even so, it should be noted that discarding multiple cells on the ground of technical quality can have a negative effect, since the number of singe cells studied in any experiment is generally limited due to high costs.

In addition to the biases listed above, batch effects are a source of technical variability that arises in the experimental stages required to generate sequencing data. Poor experimental design leads to completely confounded experiments in single-cell RNA-seq (that is, the observed variation cannot be separated into biological and technical), which can be avoided by introducing biological replicates of the different samples across all batches (Hicks et al., 2017).

Within cell variability

The main cause of technical noise across transcripts within a cell is the low amount of starting material available (mRNA) due to working with single cells (Brennecke et al., 2013). This noise component is dependent on the single-molecule capture efficiency, i.e. the fraction of mRNA molecules that are captured, amplified and subsequently sequenced from each cell (Marinov et al., 2014, Stegle et al., 2015). Notably, nano-litre volume sample preparation using microfluidics has been reported to improve capture efficiency (Wu et al., 2014) although more recent evaluations suggest that new library preparation methods, such as Smart-seq2, have outperformed these approaches when it comes to the number of molecules captured (Ziegenhain et al., 2017; Table 1).

Table 1

Qualitative comparison of library preparation methods for single-cell RNAseq. Note that full-length methods are not compatible with UMIs, and that all UMI methods capture only the 3′ end of the transcript. In-vitro transcription methods also include UMIs. Precision defined as reproducibility of the gene expression quantification. More the (+) signs, higher the costs or amplification bias or number of genes or cells (Ziegenhain et al., 2017).

	Technology	No. of cells	No. of genes	Amplification bias	Cost
Full-length cDNA methods	Smart-seq/C1	+	++	+++	++
	Smart-seq2	++	++++	+++	+++
UMI methods	SCRB-seq	+++	+++	++	+
	Drop-seq	++++	+	++	+
In-vitro transcription methods	MARS-seq	+++	+	+	+
	CEL-seq2/C1	++	++	+	++

Apart from capture efficiency, single cell analysis entails other inconveniences, namely the impossibility to use technical replicates and the PCR amplification biases, such as overrepresentation of the 3′ end of transcripts due to poly-A priming (Wu et al., 2014). Regarding uneven transcript coverage and length normalization, FPKM units (fragments per kilobase per million reads), although largely used in RNA-seq, have proven unreliable for cell-to-cell comparison when the total amount of RNA per cell differs largely (Lin et al., 2012). Such limitations, however, can be partially overcome by an alternative normalization using TPM units (transcripts per million) (Wagner et al., 2012), although whether TPM is the best unit of transcription readout is still disputed.

Specific experimental designs including external controls have been implemented to account for technical variation within each cell, including spike-in quantification of molecules of known abundance and sequence, whole-transcriptome spikes from a distantly related organism or a set of artificial spike-in mix (ERRC) (Jiang et al., 2011). These spike-ins, although used to detect low-quality cells, also serve as estimators of uneven reverse transcription and PCR amplification within each cell. To detect this, pool/split experiments are yet another way. In this case, RNA pooled from multiple cells is subsequently divided in separate reactions consisting of equal material to be used for the library construction. Variation in the pool/split experiment, as in the spike-in read counts, will be solely attributed to technical noise (Marinov et al., 2014). Using these additional controls demands further experimental and computational resources, takes over a relatively high number of the reads in the sequencing library, and not all experimental platforms can accommodate them (Bacher and Kendziorski, 2016). Furthermore, it has been recently shown that the amplification bias is stronger in endogenous genes than in ERCC spike-ins (Tung et al., 2017).

Another design for noise control is the use of unique molecular identifiers (UMIs). UMIs are short random DNA sequences (typically 6-8bp) linked to cDNA before the amplification step in order to estimate the absolute number of molecules in the cell lysate, while correcting for the amplification bias (Islam et al., 2014). During data analysis, reads with the same UMI (and mapping site) are interpreted as PCR duplicates and collapsed, preventing them to be interpreted as coming from genuine gene expression. However, UMIs will be present only in reads originating at the 3′ end of transcripts, and therefore are not compatible with isoform or allele specific expression studies, as they require full-length sequencing of the transcript. Given the limitations of both UMIs and spike-ins, a computational method that does not rely on any external controls has recently been developed; instead, CENSUS uses read counts and a generative model to estimate lysate (i.e. before reverse-transcription) transcript counts for each gene (Qiu et al., 2017).

Importantly, the low abundance of starting material has other implications, notably, not all transcripts are similarly affected by technical noise. To this extent, it is expected that transcripts with low levels of expression (resulting in low read counts) will be detected with less accuracy, as the sensitivity of methods available is limited, whereas expression of high-read count genes can be more reliably measured (Ramsköld et al., 2012). Brennecke et al. (2013) performed a single-cell RNA-Seq analysis of plant cells to establish the relationship between transcript level and sensitivity to technical noise using spike-ins as a control. They noted that, for spike-ins with an expected read count of up to 100 per cell, technical noise was maximal, making it impossible to assess whether inter-sample variability was biologically meaningful or solely due to noise. A more recent study by Kim et al. (2015) also used spike-ins to model technical noise and stated that only about 12% of variability in lowly expressed genes can be attributed to biological factors in contrast to over 55% of variability for highly expressed genes.

While technical improvements in single-cell RNA-Seq protocols are still in development, new computational tools are necessary to help discriminate genes where technical variability obfuscates the relevant biological one. We recently introduced a correlation-based method to identify the genes where biological variation was higher than technical, and create a high-quality subset that would be suitable for further analysis (Mantsoki et al., 2016). In addition, we noted that highly expressed genes generally have low coefficients of variation (CV), whereas lowly expressed genes present a wide range of CV values.

Future challenges and opportunities

There is an ongoing effort to identify and systematically characterise various sources of technical variability in single cell expression data and further develop computational tools and resources to help deal with it. Following the steps of development of UMIs and spike-ins to take into account the technical variability in single cell experiments, there is a constant need for development of experimental and computational protocols, as current ones are not exempt of limitations (Table 1), as it has been recently reviewed (Rostom et al., 2017, Ziegenhain et al., 2017).

Since technical variation can be introduced at various steps, ranging from pipetting variations during sample processing steps to amplification or sequencing bias, it is important to develop strategies that estimate and remove batch effects. Eliminating such biases in bulk RNA-seq can be done by performing technical replicates, which is constrained by the low amount of material in the case of single-cell RNA-seq. Hence, more research is necessary to allow correct discrimination of technical and biological variability. Furthermore, many computational challenges remain, including determining whether TPM is the best measure for transcript length normalization, or to what extent bulk RNA-seq analysis tools can be applied to single cell RNA-seq.

As more and more single cell datasets are generated, data integration will be inevitable. This poses a great challenge, as it was recently proven that the conversion of reads to molecules using UMIs is impacted by both biological and technical variation, indicating that UMI counts are not an unbiased estimator of gene expression levels (Tung et al., 2017). Thus, an additional challenge will be to account for biases in a way that is standardised enough to make datasets from different studies, which will be differently affected by noise, comparable.

In summary, single-cell RNA-Seq holds tremendous potential to unravel gene expression heterogeneity at an individual cell level. To fulfil this promise, the most urgent concern is that biological variation needs to be efficiently separated from technical variation. Faced with such a challenging endeavour, the joint efforts of both experimental and computational biologists are needed, as the increase in multi-disciplinary teams will help tackle different aspects of this issue and design comprehensive solutions. Apart from delineating technical and biological variation, the computational biology field can contribute in various ways to the ever growing single-cell analysis field, namely performing a systematic assessment of the many existing quality control and normalization methods, and developing statistical and computational tools for effectively discovering the biology behind high dimensional data generated by single cell RNA-seq.

Competing interests

The authors declare no completing interests.