Single-cell RNA sequencing (scRNA-seq) has provided impressive temporal resolution to our understanding of continuous biological processes such as cell differentiation (Wagner et al., 2016; Tanay and Regev, 2017). It has uncovered hidden heterogeneity among cells and exposed the factors that determine each cell’s unique identity. One broad source of transcriptomic diversity is alternative splicing, and several studies have uncovered compelling evidence of changes in alternative splicing among single cells during differentiation (Welch et al., 2016; Qiu et al., 2017; Song et al., 2017; Huang and Sanguinetti, 2017).

A particularly surprising conclusion of several scRNA-seq studies was the observation that splicing was often bimodal among supposedly homogeneous cells (Shalek et al., 2013; Marinov et al., 2014; Song et al., 2017; Westoby et al., 2018). Individual cells had binary outcomes in splicing: some cells always spliced in a particular cassette exon, and some cells never spliced in the exon, but few cells showed truly intermediate inclusion within one cell. This unexpected result contrasted with previous single molecule imaging studies of several alternative exons that showed that cell-to-cell variability is minimized and tightly regulated by the splicing machinery in single cells (Waks et al., 2011; Maamar et al., 2013). Curiosity about this result led to investigations of mechanisms that might be responsible for stochastic splicing variability among apparently homogeneous cells, such as variation in DNA methylation (Linker et al., 2019).

We propose that the observed bimodality does not generally reflect a binary nature of splicing biology, but rather, that it exposes a technical limitation of the scRNA-seq data that have been collected so far. Because alternative isoforms of a gene share much of the same sequence, only the few RNA-seq reads mapping to the exact alternative splice junctions, or to the alternative exon itself, reveal its alternative splicing. When combined with the low mRNA capture efficiency of scRNA-seq and the PCR amplification of small amounts of starting material into a full-length sequencing library, these circumstances create the risk of bottlenecks that lose all but a few individual mRNAs of most genes in each cell.

The limitations of scRNA-seq are a known obstacle to studying splicing in single cells (Arzalluz-Luque and Conesa, 2018). Similar concerns have arisen with the use of scRNA-seq to infer allelic expression; a careful analysis showed that stochastic patterns resulted from technical noise (Kim et al., 2015). A recent study observed and modeled the high dropout rate of individual isoforms in scRNA-seq and advised that scRNA-seq is fundamentally unsuitable for measuring changes in alternative splicing (Westoby et al., 2020). Others have implemented workarounds, for example using sequence features to predict splicing outcomes in lieu of sufficient sequencing coverage (Huang and Sanguinetti, 2017), or attempting to identify excess variance beyond technical noise (Welch et al., 2016; Linker et al., 2019). These studies have identified true examples of differential splicing in single cells, but they fundamentally do not explain how scRNA-seq limitations have caused qualitative, not just quantitative, distortions in our understanding of alternative splicing.

Here, we show that scRNA-seq splicing data are consistent with a simple model (Figure 1). Consider a particular cassette exon whose true pattern of exclusion follows a unimodal distribution of isoform ratios across cells (i.e. most cells express both isoforms, with a ratio revolving around the same mean). This distribution can be distorted by extreme information loss during library preparation and sequencing, creating the illusion that individual cells only produce one isoform or the other. Our simulations make it clear that the reliability of splicing measurements is a function of the initial amount of mRNA, the efficiency of its recovery, the underlying splicing rate, and further distortions from PCR amplification of cDNA. These effects should be considered when interpreting previous studies that used qualitative changes in the observed distribution of the splicing rates (Song et al., 2017) or changes in their variance (Linker et al., 2019) as evidence for regulation of alternative splicing. Considering the true amount of information available for a cassette exon can allow for accurate observations of alternative splicing. Using a data normalization and filtering method to identify cassette exons with sufficient information, we are able to draw biologically relevant conclusions about alternative splicing in single cells.

Figure 1

Download asset Open asset

Our interest in splicing regulation led us to examine alternative splicing in several single cell differentiation datasets from mice and humans that were generated with methods that recover sequence from along the full length of mRNAs. To investigate the reported high variability of splicing between cells more closely, we began by examining the splicing of cassette exons in a high-coverage mouse scRNA-seq dataset (Chen et al., 2016), estimating their percent spliced-in as the fraction of splice junction reads that show exon inclusion (out of all reads that cover the junction). We use $\widehat{\mathrm{\Psi }}$ to denote these estimated rates, while $\mathrm{\Psi }$ denotes the actual rate as it is in the cell. For clarity, we define a single $\widehat{\mathrm{\Psi }}$ observation (which pertains to a specific cassette exon in an individual cell) as binary if it is close to 0 or 1 (i.e. the respective cell tends to express transcripts that either include the exon or exclude it, but not both). We then describe the distribution of an exon’s $\widehat{\mathrm{\Psi }}$ across cells as bimodal when its individual values are predominantly binary, where some cells have a $\widehat{\mathrm{\Psi }}$ close to 1 (most observed transcripts include the exon) and others have $\widehat{\mathrm{\Psi }}$ close to 0 (most observed transcripts do not include the exon). Strikingly, when we inspected several exons, we saw that they had more binary outcomes in cells with fewer reads covering their splice junctions, while cells with more reads were more likely to show non-binary $\widehat{\mathrm{\Psi }}$ values (Figure 2a). We realized that this effect of coverage may reflect a non-binary reality, since even if both isoforms are expressed in a certain cell, the likelihood of observing both isoforms is reduced as the number of captured mRNAs decreases. In contrast, if the underlying distribution were indeed bimodal with binary modes, as previously proposed (Shalek et al., 2013; Marinov et al., 2014; Song et al., 2017), then the read coverage would have little effect on the proportion of binary $\widehat{\mathrm{\Psi }}$ observations across cells.

To further explore this phenomenon, we extended our analysis to the full scRNA-seq datasets. In all cases, we found a strong effect of coverage on the observed binary $\widehat{\mathrm{\Psi }}$ in exons with intermediate splicing (average $\widehat{\mathrm{\Psi }}$ between 0.2 and 0.8). We consistently found that exons with low junction read coverage had more binary $\widehat{\mathrm{\Psi }}$ values and bimodal $\widehat{\mathrm{\Psi }}$ distributions (Figure 2b,c; Figure 2—figure supplement 1a,b). We found that the association between binary values and read coverage was not observed in exons that are binary but not bimodal (i.e. nearly constitutively excluded or included exons, with average $\widehat{\mathrm{\Psi }}$ close to 0 or 1; Figure 2c). Taken together, these observations suggest that the abundance of binary observations in exon inclusion patterns may reflect a distortion of an underlying unimodal splicing distribution (i.e. when cells in fact express both isoforms), rather than a truly bimodal splicing pattern in the analyzed cells.

Figure 2 with 1 supplement see all

Download asset Open asset

We then asked if the association of binary splicing outcomes with low read coverage could be due to some biological consequence of low transcription, or if it was a technical consequence of low sequencing coverage. We found that, among the cells in a single experiment, the cells with an overall higher number of splice junction reads also tended to have a smaller fraction of exons with binary values (Figure 2d, Figure 2—figure supplement 1b), suggesting an influence of technical coverage rather than transcription level. We further considered the effect of biological variations in transcription. Transcription at a single locus occurs in intermittent bursts of RNA synthesis (Raj and van Oudenaarden, 2008), suggesting a possible effect of size and frequency of transcriptional bursts on binary splicing. Using transcriptional bursting data from Larsson et al., 2019, we found that genes with either high burst frequency or large bursts did exhibit fewer binary splicing observations (Figure 2—figure supplement 1c), but a linear regression showed that burst frequency and size did not contribute to binary observations beyond the effect of read coverage (Figure 2—figure supplement 1d). This suggests that transcriptional bursting contributes to splicing variability only by determining the overall abundance of an mRNA, consistent with single-molecule fluorescence observations (Waks et al., 2011).

Simulations of RNA sequencing reveal technical sources of distortion of splicing estimates

A simple probabilistic exercise shows the potential loss of splicing information during sequencing. Single cell RNA-seq experiments that capture full-length transcripts have an estimated capture efficiency of only ∼10%, due to RNA degradation and inefficient reverse transcription (Marinov et al., 2014; Grün et al., 2014; Ziegenhain et al., 2017; Qiu et al., 2017). For instance, a gene that expresses 20 mRNA molecules in a cell might only have two mRNAs recovered, and if that gene is alternatively spliced with a true splicing rate $\mathrm{\Psi }$ of 0.5, there is approximately a 50% chance that those few recovered mRNAs will only represent one of the two isoforms that were originally present in the cell (Figure 3—figure supplement 1a,b). As many genes are expressed at just a few RNA molecules per cell, low recovery might affect many alternative splicing events (Shapiro et al., 2013; Zenklusen et al., 2008). Furthermore, while the empirically observed $\widehat{\mathrm{\Psi }}$ provides a maximum likelihood estimate for the true splicing rate, the uncertainty of this estimate (i.e. the range of alternative values with a nearly similar likelihood) decreases substantially with the number of observed molecules (Figure 3—figure supplement 1c and Materials and methods). As a result, the probability that the observed $\widehat{\mathrm{\Psi }}$ is close to the true underlying $\mathrm{\Psi }$ increases when more mRNA molecules are captured (Figure 3—figure supplement 1d).

Our theoretical reasoning above relied on a simple model where the number of observed mRNA molecules (rather than number of reads) is known and the only distorting factor is a limited capture efficiency. In practice, both of these assumptions are challenged due to additional factors, such as PCR amplification and variability in the capture efficiency across cells. To investigate the pertinent effects on $\widehat{\mathrm{\Psi }}$ distributions in this more complex setting, we designed a probabilistic simulator of alternative splicing in single cells (Figure 3—figure supplement 2). The model has two main components: we begin by simulating the underlying molecular content of each cell, by drawing gene expression levels and cassette exon splicing rates from a probabilistic model of cell state. We then simulate the technical process of extracting data from each cell using single-cell RNA sequencing with full transcript coverage. This part accounts for variability in capture rates, and the effects of PCR amplification, fragmentation and sequencing. It relies on SymSim, a simulation tool for single-cell RNA sequencing data (Zhang et al., 2019). The final product of our simulation is the number of splice junction reads that either span or skip each exon in each cell. These numbers are distorted in a way that reflects real nuisance factors. For instance, two reads could have originated from the same molecule due to amplification effects.

We used our simulator to investigate how the observed inclusion ($\widehat{\mathrm{\Psi }}$) of cassette exons differs from the underlying $\mathrm{\Psi }$, under different average capture rates, and setting the other technical parameters to values that are characteristic of Smart-seq2 datasets (see Materials and methods). We considered either a binary-bimodal regime of $\mathrm{\Psi }$ (i.e. both isoforms are expressed in the population, but rarely by the same cell; Figure 1a), or a non-binary regime (cells tend to express both isoforms; Figure 1b). We simulated splicing of cassette exons in 1500 genes, in a population of 300 single cells. In the bimodal simulation, 500 exons were modeled to have alternative splicing with a bimodal distribution across cells, and in the unimodal simulation, these 500 exons were modeled to have a unimodal distribution. To reflect the patterns seen in real data, we also simulated splicing of 500 exons that were nearly constitutively included and 500 that were nearly constitutively skipped.

As expected, in the binary-bimodal simulation, the observed $\widehat{\mathrm{\Psi }}$ reflected the underlying process well, independent of the average capture efficiency (Figure 3a). In contrast, when we modeled a non-binary splicing regime, the observed $\widehat{\mathrm{\Psi }}$ distributions were strikingly similar to the splicing distributions of cassette exons in real single-cell RNA-seq datasets (Figure 3b). Specifically, the loss of information due to mRNA recovery and library generation led many of the observed $\widehat{\mathrm{\Psi }}$ to become binary, and their observed distribution across cells to become bimodal. This tendency again correlated with coverage, whereby lowly covered exons showed the strongest effect, while exons with high coverage maintained a non-binary, unimodal distribution. Consistently, in this non-binary regime, the average of $\widehat{\mathrm{\Psi }}$ was similar to the true average of $\mathrm{\Psi }$, but the variance of $\widehat{\mathrm{\Psi }}$ increased (Figure 3—figure supplement 1e,f). Furthermore, as in the real data sets (Figure 2c), we also found that the dependency between read coverage and the chance of observing a binary $\widehat{\mathrm{\Psi }}$ is more pronounced in exons with an underlying $\mathrm{\Psi }$ that is far from binary (Figure 3c–f), highlighting again that such an association likely indicates an artifact.

Figure 3 with 2 supplements see all

Download asset Open asset

To address the extent of the distortion from mRNA recovery, we ran the simulator with varying capture efficiency and a fixed underlying $\mathrm{\Psi }=0.5$. We found that decreasing the average capture efficiency dramatically increased the number of binary $\widehat{\mathrm{\Psi }}$ observations, particularly for exons with low expression, although even highly expressed genes suffered great distortion in the observed $\widehat{\mathrm{\Psi }}$ when the average capture efficiency was very low (Figure 3g). These results reinforced the hypothesis that capture efficiency is a main technical factor that creates the appearance of bimodality in single-cell splicing.

Finally, we estimated the chance of observing only one type of isoform in any given cell (i.e. a binary $\widehat{\mathrm{\Psi }}$) as a function of the underlying $\mathrm{\Psi }$ and the number of transcripts that are present in the cell (Figure 3h). For this analysis, we set an average mRNA capture rate of 10%. Our results delineate the range of values in which an artifact is less expected. For instance, for an exon with 50% inclusion rate, a non-binary estimate is more likely if the respective gene has at least 50 transcripts in the cell. Notably, these estimates are more conservative than the theoretical analysis (Figure 3—figure supplement 1a,b), due to the effect of the technical nuisance factors we modeled in this simulation.

Accounting for mRNA recovery improves analysis of alternative splicing

We sought criteria that would identify reliable measurements of splicing in single cells and avoid distortions from low mRNA recovery. While ideally we should rely on the actual number of mRNA molecules recovered, full-length RNA-seq experiments generally do not report an absolute mRNA count. Previous studies assessed the quality of splicing observations based on the number of reads covering alternative splice junctions (Song et al., 2017; Linker et al., 2019), but the number of splice junction reads is influenced by the extent of PCR amplification and sequencing depth, and so it does not directly reflect the number of recovered mRNAs.

We realized that we could estimate the number of mRNA molecules that were captured into cDNA by using the Census normalization approach proposed by Qiu et al., 2017. This method infers a per-cell scaling factor between the relative abundance of each transcript, inferred from RNA-seq, and the actual number of mRNAs recovered. We found that some datasets with many reads per cell, such as the Song dataset (Song et al., 2017), nonetheless had few mRNAs recovered per cell (Figure 4a, Figure 4—figure supplement 1c), which may explain the extensive splicing bimodality in this dataset. The dataset with the highest recovery of mRNAs (Chen et al., 2016) indeed showed less binary splicing (Figure 2b).

Figure 4 with 4 supplements see all

Download asset Open asset

Next, we sought a threshold for mRNA recovery that would suppress spurious observations of binary splicing. Our simulations and our analytical calculations suggested that, for a capture efficiency of 10%, recovering 7–10 mRNA molecules would be more likely than not to capture both isoforms for exons with intermediate splicing (Figure 3h), and that the observed $\widehat{\mathrm{\Psi }}$ would most likely be within 0.1 of the real value (Figure 2—figure supplement 1d). In keeping with this, we saw in the real data that exons with an average of at least 10 mRNAs recovered per cell generally had substantially fewer binary observations (Figure 4b, Figure 4—figure supplement 1d; limited to exons with average $\widehat{\mathrm{\Psi }}$ between 0.05 and 0.95). Nonetheless, a subset of these exons still showed binary splicing in most cells. These exons had few splice junction reads relative to the estimated mRNA count. We expect that these represent anomalously low recovery of reads from the specific splice junctions of interest, perhaps due to annotation errors or poor recovery of fragments with particular sequence composition. We reasoned that a data-driven splice junction read criterion would prevent distortions arising from this low read recovery. To find an appropriate read minimum, we determined the number of splice junction reads expected to arise from a cassette exon in each cell in different datasets, based on that cell’s overall mRNA recovery (Figure 4c; see Materials and methods). This provided a second filtering criterion: we excluded exons with fewer splice junction reads than would be expected from 10 mRNAs, given that exon’s $\widehat{\mathrm{\Psi }}$ and the coverage rate in that particular cell. This metric is calculated for each exon separately, driven by the actual information in each cell, and it varies substantially between datasets.

We then selected the exons from each cell that met the combined requirement of 10 mRNAs and the splice junction reads that would arise from 10 mRNAs, and asked what effect this filter had on binary observations of splicing. As an example, our filter removed seemingly spurious binary observations for alternatively spliced exon 8 of Cadm1 (Figure 4d). The observations that remained after filtering showed a clear pattern of increased inclusion along pseudotime, with exon skipping in stem cells and exon inclusion in fully differentiated neurons (Spearman’s $r_s$ = 0.52, p=1.2 × 10⁷; Figure 4e). The correlation between splicing observations and differentiation pseudotime highlighted that cell type differences could contribute to observations of bimodality that are not the result of an artificial inflation of binary observations, and this was an important factor to consider in our analysis.

Bimodality in filtered exons is explained by differential splicing between cell types

Earlier papers had reported a surprising amount of bimodal splicing within a cell type, with extensive binary outcomes in individual cells, that is similar cells with different dominant isoforms. This suggested that some aspects of gene expression stochastically locked individual cells into producing primarily molecules of one or the other isoform. However, our results so far suggest that much of this observed variability may be an artifact of low mRNA recovery. In contrast, differences in splicing between cells of different types are consistent with conventional mechanisms such as regulation by cell-type-specific splicing factors. Therefore, when considering heterogeneous cell samples, large shifts in splicing between cell types may result in truly bimodal $\mathrm{\Psi }$. In contrast to the first scenario, we would expect the corresponding observations $\widehat{\mathrm{\Psi }}$ to remain bimodal even after low-quality data points were removed.

We therefore sought to distinguish between these two scenarios and estimate the extent of bimodal splicing within a cell type that is strongly supported by the data. Most of the datasets in our analysis came from differentiation time courses, and so to examine homogeneous cell subsets, we stratified the cells into groups indicative of their developmental stage. The groups were identified by clustering the cells using the normalized sum of reads from every gene (i.e. not directly representing their variation in splicing; Figure 4f) and labeled based on expression of known marker genes (Figure 4—figure supplement 1e).

We found that, prior to filtering, hundreds of exons with intermediate splicing (average $\widehat{\mathrm{\Psi }}$ between 0.2 and 0.8) seemed to have at least weakly bimodal splicing distributions (using a heuristic definition of at least 25% of cells with $\widehat{\mathrm{\Psi }}\le$ 0.25% and 25% with $\widehat{\mathrm{\Psi }}\ge$ 0.75), a rate that roughly matched previous descriptions (Song et al., 2017; Table 1). We then discarded the observations that did not meet our filtering criteria, and retained for every cluster only those exons with observations remaining in at least 50% of cells (Figure 4g). With that filtering we found that most remaining exons had fewer binary observations than the discarded exons (Figure 4—figure supplement 1f). Furthermore, cell clusters had no or very few remaining exons with bimodal splicing distributions (Figure 4h, Figure 4—figure supplement 1g,h, Table 1). For instance, in the Chen dataset, there were two remaining cases of bimodality, and both seem to be the result of sex-specific splicing (Figure 4—figure supplement 1i). This trend persisted when reanalyzing all datasets and clusters using a range of cutoffs for qualifying a splicing pattern as binary (Figure 4—figure supplement 1e).

Table 1

Dataset	Cell type cluster	Cells	Exons	Bimodal	% bimodal	Selected	Selected bimodal	% bimodal selected	p-val (adj)
Chen	ES	217	446	118	26%	94	0	0%	3.3e-14
Chen	Epi, early	98	402	107	27%	98	1	1%	9.3e-14
Chen	Epi, late	104	516	136	26%	76	1	1%	1.1e-09
Chen	Neuron, early	47	364	117	32%	43	0	0%	3.6e-08
Chen	Neuron, late	22	517	146	28%	61	0	0%	1.1e-09
Lescroart	Heart E6.75	172	286	77	27%	33	0	0%	2.0e-05
Lescroart	Heart E7.25	341	291	78	27%	36	0	0%	8.0e-06
Trapnell	Myoblast 00 hr	35	400	142	36%	41	0	0%	1.2e-08
Trapnell	Myoblast 24 hr	89	251	97	39%	27	0	0%	1.2e-06
Trapnell	Myoblast 48 hr	72	242	97	40%	32	1	3%	9.1e-07
Trapnell	Myoblast 72 hr	35	252	101	40%	24	1	4%	5.1e-05
Song	iPSC	62	616	269	44%	55	0	0%	3.3e-14
Song	NPC	73	212	92	43%	28	1	4%	1.2e-06
Song	MN	67	168	82	49%	19	4	21%	9.4e-03
Shalek	BMDC	13	149	51	34%	27	1	4%	6.6e-05

Filtering improves the identification of differentially spliced exons

Finally, we asked if our filtering criteria could help identify biologically relevant changes in splicing between cell types, by selecting exons whose biological variance is not overwhelmed by high technical variance. Focusing on three datasets with multiple cell type clusters (Chen, Trapnell, and Song), we examined the effect of retaining only exons that passed our filter in at least 50% of cells in each cluster. To avoid biases depending on the number of observations available for each remaining exon, we retained all data points for these exons, rather than removing individual cell data points that did not pass the filter. As a benchmark for comparison, we used a simpler criterion, akin to previous studies, of at least 10 splice junction reads in at least 50% of cells (Figure 4i). To focus on exons with substantial alternative splicing, we omitted exons from this analysis that were consistently excluded (average $\widehat{\mathrm{\Psi }}$ < 0.05) or included (average $\widehat{\mathrm{\Psi }}$ > 0.95) across all cells.

Taken together, these results challenge the idea that widespread bimodality among homogeneous cells arises from a binary nature of splicing in single cells, and instead suggest that splicing bimodality reflects biological differences between cell types or subtypes. Indeed, experimental observations of splicing bimodality support the importance of differences between cell types, rather than stochastic differences between homogeneous cells. The two cases of splicing bimodality confirmed with smFISH by Song et al., 2017 both showed bimodality between two characterized cell types. Similarly, Shalek et al., 2013 confirmed one case of splicing bimodality, in a gene whose expression they showed to be a marker of a hidden difference between mature and maturing cells, and so its splicing may reflect cell state difference at a wider context as well.

For each remaining exon in the two filtering schemes, we performed a Kruskal–Wallis one-way analysis of variance, which assigned it a p-value (Materials and methods). This test determines if the exon’s median $\mathrm{\Psi }$ changed significantly between any clusters. Both filtering schemes enriched substantially for exons that show significant changes across different clusters, compared with the unfiltered set of all exons with average $\widehat{\mathrm{\Psi }}$ between 0.05 and 0.95. Furthermore, the extent of enrichment increased with the strictness of the differential splicing test (Figure 4i–k; Figure 4—figure supplement 2). It is interesting to note that in the Chen dataset, the baseline (ten reads) filter was stricter than the combined filter, retaining 66 vs. 198 exons and thus achieving high precision at the expense of a lower recall. Indeed, due to the high mRNA recovery and low amplification in this dataset, our combined filter required only seven splice junction reads, rather than ten as in the baseline scheme. Conversely, in the Song dataset, which is predicted to have low mRNA recovery yet a large number of reads, the baseline filter retained almost all exons, while the combined filter retained only 104 of them, leading to significant over-representation of exons that are differentially spliced between clusters. The results show how differences in mRNA recovery and amplification between datasets may change the interpretation of an absolute number of reads. We corroborated these results using spatial autocorrelation as a way of identifying informative exons, rather than the clustering and differential splicing analysis. The autocorrelation test builds on the work of DeTomaso et al., 2019 to identify exons with a significantly high level of consistency in $\widehat{\mathrm{\Psi }}$ among transcriptionally similar cells (Figure 4—figure supplement 3).

Taken together, these results indicate that careful filtering of exons can help identify observations that are consistent with the gene expression space and the stratification of cells into types, and are thus likely more biologically meaningful (e.g. capturing known forms of regulation during differentiation) (Figure 4l, Figure 4—figure supplement 4; Liu et al., 2018). The concept of ‘meaningful’ is ascribed here to the observations, not the exons themselves, as there may be many additional exons with relevant variation (e.g. stable differences between cell types) that is not reliably measured due to low coverage. Indeed, the effects of the two filtering schemes and the differences between them exemplify how technical factors may distort our view of the data and consequently, our understanding of variation and regulation of splicing.

The surprising result that alternatively splicing is bimodal among single cells provoked curiosity and speculation. Bimodal outcomes might reflect hidden cell subtypes, but the bimodality was seen even among apparently homogeneous cells. Did splicing outcomes reflect some unknown, stochastic cell state?

We have shown here that the bimodal patterns could have an entirely different explanation: profound technical limitations of single cell RNA sequencing. A crucial limit on biologically meaningful splicing observations in a single cell is the number of mRNAs available to inform the measurement. This is determined both by the expression of the genes that contain the exons of interest, and by the capture efficiency of the experiment. It is important to note that the depth of a sequencing library does not necessarily reflect its quality. Along with low mRNA numbers, splicing observations are also distorted by uneven amplification efficiency and cDNA overamplification. Increasing PCR amplification cycles in an attempt to compensate for low capture efficiency has the risk of worsening the technical distortion. Indeed, in our analysis, the dataset with the highest read count per cell actually had quite low mRNA recovery and large technical distortion, creating an appearance of bimodal splicing (Song et al., 2017). Moreover, a qualitative change in the observed $\widehat{\mathrm{\Psi }}$ distribution of an exon between single cell subpopulations does not necessarily reflect a change in the underlying splicing rate, as changes in gene expression and mRNA recovery between samples can create the illusion of a splicing change.

Further developments in statistical analysis that carefully account for both missing and redundant information due to low capture efficiency could make splicing observations in single cells more reliable. We set the foundation for such analysis by proposing a probabilistic process that describes the biological and technical steps that generate single cell splicing data. We also introduced a simple approach that builds on the Census normalization (Qiu et al., 2017) to estimate the number of mRNAs recovered and the extent of artificial duplication of splicing information. This metric provides a practical filter for identifying exons with sufficient information to analyze. On the experimental side, improving the capture efficiency of scRNA-seq methods while moderating the extent of overamplification is crucial for increasing the subset of exons for which reliable observation can be made.

True biological insight into alternative splicing can indeed be found from high-quality scRNA-seq data, and we hope that new methods will allow better understanding of splicing regulation, cell-to-cell variation, and the importance of alternative splicing in defining cell fate (Hagemann-Jensen et al., 2020; Gupta et al., 2018). However, some limitations are inherent to the situation. Single cells express a limited number of mRNAs per gene; splicing observations in single cells will always be inherently noisy reflections of the underlying biology.

Analysis code is available at https://github.com/lareaulab/sc_binary_splicing. (Buen Abad Najar and Lareau, 2020; copy archived at https://github.com/elifesciences-publications/sc_binary_splicing).

All sequencing data reanalyzed in this study were acquired from GEO.

1. 1. Chen G
   2. Schell JP
   3. Benitez JA
   4. Petropoulos S
   5. Yilmaz M
   6. Reinius B
   7. Alekseenko Z
   8. Shi L
   9. Hedlund E
   10. Lanner F
   11. Sandberg R
   12. Deng Q

   (2016) NCBI Gene Expression Omnibus

   ID GSE74155. Single-cell analysis of allelic gene expression in pluripotency, differentiation and X-chromosome inactivation.

   https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE74155
2. 1. Lescroart F
   2. Wang X
   3. Lin X
   4. Swedlund B
   5. Gargouri S
   6. Sànchez-Dànes A
   7. Moignard V
   8. Dubois C
   9. Paulissen C
   10. Kinston S
   11. Göttgens B
   12. Blanpain C

   (2018) NCBI Gene Expression Omnibus

   ID GSE100471. Defining the early steps of cardiovascularlineage segregation by single cell RNA-seq.

   https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE100471
3. 1. Trapnell C
   2. Cacchiarelli D
   3. Grimbsby J
   4. Pokharel P
   5. Li S
   6. Morse M
   7. Mikkelsen T
   8. Rinn J

   (2014) NCBI Gene Expression Omnibus

   ID GSE52529. Pseudo-temporal ordering of individual cells reveals regulators of differentiation.

   https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE52529
4. 1. Song Y
   2. Botvinnik OB
   3. Lovci MT
   4. Kakaradov B
   5. Liu P
   6. Xu JL
   7. Yeo GW

   (2017) NCBI Gene Expression Omnibus

   ID GSE85908. Single-cell alternative splicing analysis with Expedition reveals splicing dynamics during neuron differentiation.

   https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE85908
5. 1. Fletcher RB
   2. Das D
   3. Gadye L
   4. Street KN
   5. Baudhuin A
   6. Wagner A
   7. Cole MB
   8. Flores Q
   9. Choi YG
   10. Yosef N
   11. Purdom E
   12. Dudoit S
   13. Risso D
   14. Ngai J

   (2017) NCBI Gene Expression Omnibus

   ID GSE95601. Olfactory stem cell differentiation: horizontal basal cell (HBC) lineage.

   https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE95601

1. 1. Arzalluz-Luque Á
   2. Conesa A

   (2018) Single-cell RNAseq for the study of isoforms—how is that possible?

   Genome Biology 19:110.

   https://doi.org/10.1186/s13059-018-1496-z

   • PubMed
   • Google Scholar
2. 1. Barash Y
   2. Blencowe BJ
   3. Frey BJ

   (2010) Model-based detection of alternative splicing signals

   Bioinformatics 26:i325–i333.

   https://doi.org/10.1093/bioinformatics/btq200

   • Google Scholar
3. Software

   1. Buen Abad Najar CF
   2. Lareau L

   (2020) sc_binary_splicing, version 5c2c825

   GitHub.

   https://github.com/lareaulab/sc_binary_splicing
4. 1. Busskamp V
   2. Lewis NE
   3. Guye P
   4. Ng AH
   5. Shipman SL
   6. Byrne SM
   7. Sanjana NE
   8. Murn J
   9. Li Y
   10. Li S
   11. Stadler M
   12. Weiss R
   13. Church GM

   (2014) Rapid neurogenesis through transcriptional activation in human stem cells

   Molecular Systems Biology 10:760.

   https://doi.org/10.15252/msb.20145508

   • PubMed
   • Google Scholar
5. 1. Chen G
   2. Schell JP
   3. Benitez JA
   4. Petropoulos S
   5. Yilmaz M
   6. Reinius B
   7. Alekseenko Z
   8. Shi L
   9. Hedlund E
   10. Lanner F
   11. Sandberg R
   12. Deng Q

   (2016) Single-cell analyses of X chromosome inactivation dynamics and pluripotency during differentiation

   Genome Research 26:1342–1354.

   https://doi.org/10.1101/gr.201954.115

   • PubMed
   • Google Scholar
6. 1. Cole MB
   2. Risso D
   3. Wagner A
   4. DeTomaso D
   5. Ngai J
   6. Purdom E
   7. Dudoit S
   8. Yosef N

   (2019) Performance assessment and selection of normalization procedures for Single-Cell RNA-Seq

   Cell Systems 8:315–328.

   https://doi.org/10.1016/j.cels.2019.03.010

   • PubMed
   • Google Scholar
7. 1. DeTomaso D
   2. Jones MG
   3. Subramaniam M
   4. Ashuach T
   5. Ye CJ
   6. Yosef N

   (2019) Functional interpretation of single cell similarity maps

   Nature Communications 10:4376.

   https://doi.org/10.1038/s41467-019-12235-0

   • PubMed
   • Google Scholar
8. 1. Dobin A
   2. Davis CA
   3. Schlesinger F
   4. Drenkow J
   5. Zaleski C
   6. Jha S
   7. Batut P
   8. Chaisson M
   9. Gingeras TR

   (2013) STAR: ultrafast universal RNA-seq aligner

   Bioinformatics 29:15–21.

   https://doi.org/10.1093/bioinformatics/bts635

   • PubMed
   • Google Scholar
9. 1. Faigenbloom L
   2. Rubinstein ND
   3. Kloog Y
   4. Mayrose I
   5. Pupko T
   6. Stein R

   (2015) Regulation of alternative splicing at the single-cell level

   Molecular Systems Biology 11:845.

   https://doi.org/10.15252/msb.20156278

   • PubMed
   • Google Scholar
10. 1. Fletcher RB
    2. Das D
    3. Gadye L
    4. Street KN
    5. Baudhuin A
    6. Wagner A
    7. Cole MB
    8. Flores Q
    9. Choi YG
    10. Yosef N
    11. Purdom E
    12. Dudoit S
    13. Risso D
    14. Ngai J

    (2017) Deconstructing olfactory stem cell trajectories at Single-Cell resolution

    Cell Stem Cell 20:817–830.

    https://doi.org/10.1016/j.stem.2017.04.003

    • PubMed
    • Google Scholar
11. 1. Grün D
    2. Kester L
    3. van Oudenaarden A

    (2014) Validation of noise models for single-cell transcriptomics

    Nature Methods 11:637–640.

    https://doi.org/10.1038/nmeth.2930

    • PubMed
    • Google Scholar
12. 1. Gupta I
    2. Collier PG
    3. Haase B
    4. Mahfouz A
    5. Joglekar A
    6. Floyd T
    7. Koopmans F
    8. Barres B
    9. Smit AB
    10. Sloan SA
    11. Luo W
    12. Fedrigo O
    13. Ross ME
    14. Tilgner HU

    (2018) Single-cell isoform RNA sequencing characterizes isoforms in thousands of cerebellar cells

    Nature Biotechnology 36:1197–1202.

    https://doi.org/10.1038/nbt.4259

    • Google Scholar
13. 1. Hagemann-Jensen M
    2. Ziegenhain C
    3. Chen P
    4. Ramsköld D
    5. Hendriks G-J
    6. Larsson AJM
    7. Faridani OR
    8. Sandberg R

    (2020) Single-cell RNA counting at allele and isoform resolution using Smart-seq3

    Nature Biotechnology 38:708–714.

    https://doi.org/10.1038/s41587-020-0497-0

    • Google Scholar
14. 1. Huang Y
    2. Sanguinetti G

    (2017) BRIE: transcriptome-wide splicing quantification in single cells

    Genome Biology 18:123.

    https://doi.org/10.1186/s13059-017-1248-5

    • Google Scholar
15. 1. Hubbard KS
    2. Gut IM
    3. Lyman ME
    4. McNutt PM

    (2013) Longitudinal RNA sequencing of the deep transcriptome during neurogenesis of cortical glutamatergic neurons from murine ESCs

    F1000Research 2:35.

    https://doi.org/10.12688/f1000research.2-35.v1

    • PubMed
    • Google Scholar
16. 1. Kim JK
    2. Kolodziejczyk AA
    3. Ilicic T
    4. Illicic T
    5. Teichmann SA
    6. Marioni JC

    (2015) Characterizing noise structure in single-cell RNA-seq distinguishes genuine from technical stochastic allelic expression

    Nature Communications 6:8687.

    https://doi.org/10.1038/ncomms9687

    • PubMed
    • Google Scholar
17. 1. Larsson AJM
    2. Johnsson P
    3. Hagemann-Jensen M
    4. Hartmanis L
    5. Faridani OR
    6. Reinius B
    7. Segerstolpe Åsa
    8. Rivera CM
    9. Ren B
    10. Sandberg R

    (2019) Genomic encoding of transcriptional burst kinetics

    Nature 565:251–254.

    https://doi.org/10.1038/s41586-018-0836-1

    • Google Scholar
18. 1. Lescroart F
    2. Wang X
    3. Lin X
    4. Swedlund B
    5. Gargouri S
    6. Sànchez-Dànes A
    7. Moignard V
    8. Dubois C
    9. Paulissen C
    10. Kinston S
    11. Göttgens B
    12. Blanpain C

    (2018) Defining the earliest step of cardiovascular lineage segregation by single-cell RNA-seq

    Science 359:1177–1181.

    https://doi.org/10.1126/science.aao4174

    • Google Scholar
19. 1. Li B
    2. Dewey CN

    (2011) RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome

    BMC Bioinformatics 12:323.

    https://doi.org/10.1186/1471-2105-12-323

    • PubMed
    • Google Scholar
20. 1. Linker SM
    2. Urban L
    3. Clark SJ
    4. Chhatriwala M
    5. Amatya S
    6. McCarthy DJ
    7. Ebersberger I
    8. Vallier L
    9. Reik W
    10. Stegle O
    11. Bonder MJ

    (2019) Combined single-cell profiling of expression and DNA methylation reveals splicing regulation and heterogeneity

    Genome Biology 20:30.

    https://doi.org/10.1186/s13059-019-1644-0

    • Google Scholar
21. 1. Liu J
    2. Geng A
    3. Wu X
    4. Lin RJ
    5. Lu Q

    (2018) Alternative RNA splicing associated with mammalian neuronal differentiation

    Cerebral Cortex 28:2810–2816.

    https://doi.org/10.1093/cercor/bhx160

    • PubMed
    • Google Scholar
22. 1. Maamar H
    2. Cabili MN
    3. Rinn J
    4. Raj A

    (2013) linc-HOXA1 is a noncoding RNA that represses Hoxa1 transcription in Cis

    Genes & Development 27:1260–1271.

    https://doi.org/10.1101/gad.217018.113

    • PubMed
    • Google Scholar
23. 1. Marinov GK
    2. Williams BA
    3. McCue K
    4. Schroth GP
    5. Gertz J
    6. Myers RM
    7. Wold BJ

    (2014) From single-cell to cell-pool transcriptomes: stochasticity in gene expression and RNA splicing

    Genome Research 24:496–510.

    https://doi.org/10.1101/gr.161034.113

    • PubMed
    • Google Scholar
24. 1. Qiu X
    2. Hill A
    3. Packer J
    4. Lin D
    5. Ma Y-A
    6. Trapnell C

    (2017) Single-cell mRNA quantification and differential analysis with Census

    Nature Methods 14:309–315.

    https://doi.org/10.1038/nmeth.4150

    • Google Scholar
25. 1. Raj A
    2. van Oudenaarden A

    (2008) Nature, nurture, or chance: stochastic gene expression and its consequences

    Cell 135:216–226.

    https://doi.org/10.1016/j.cell.2008.09.050

    • PubMed
    • Google Scholar
26. 1. Shalek AK
    2. Satija R
    3. Adiconis X
    4. Gertner RS
    5. Gaublomme JT
    6. Raychowdhury R
    7. Schwartz S
    8. Yosef N
    9. Malboeuf C
    10. Lu D
    11. Trombetta JJ
    12. Gennert D
    13. Gnirke A
    14. Goren A
    15. Hacohen N
    16. Levin JZ
    17. Park H
    18. Regev A

    (2013) Single-cell transcriptomics reveals bimodality in expression and splicing in immune cells

    Nature 498:236–240.

    https://doi.org/10.1038/nature12172

    • PubMed
    • Google Scholar
27. 1. Shapiro E
    2. Biezuner T
    3. Linnarsson S

    (2013) Single-cell sequencing-based technologies will revolutionize whole-organism science

    Nature Reviews Genetics 14:618–630.

    https://doi.org/10.1038/nrg3542

    • PubMed
    • Google Scholar
28. 1. Shen S
    2. Park JW
    3. Lu ZX
    4. Lin L
    5. Henry MD
    6. Wu YN
    7. Zhou Q
    8. Xing Y

    (2014) rMATS: robust and flexible detection of differential alternative splicing from replicate RNA-Seq data

    PNAS 111:E5593–E5601.

    https://doi.org/10.1073/pnas.1419161111

    • PubMed
    • Google Scholar
29. 1. Song Y
    2. Botvinnik OB
    3. Lovci MT
    4. Kakaradov B
    5. Liu P
    6. Xu JL
    7. Yeo GW

    (2017) Single-Cell Alternative Splicing Analysis with Expedition Reveals Splicing Dynamics during Neuron Differentiation

    Molecular Cell 67:148–161.

    https://doi.org/10.1016/j.molcel.2017.06.003

    • Google Scholar
30. 1. Street K
    2. Risso D
    3. Fletcher RB
    4. Das D
    5. Ngai J
    6. Yosef N
    7. Purdom E
    8. Dudoit S

    (2018) Slingshot: cell lineage and pseudotime inference for single-cell transcriptomics

    BMC Genomics 19:477.

    https://doi.org/10.1186/s12864-018-4772-0

    • PubMed
    • Google Scholar
31. 1. Tanay A
    2. Regev A

    (2017) Scaling single-cell genomics from phenomenology to mechanism

    Nature 541:331–338.

    https://doi.org/10.1038/nature21350

    • Google Scholar
32. 1. Trapnell C
    2. Cacchiarelli D
    3. Grimsby J
    4. Pokharel P
    5. Li S
    6. Morse M
    7. Lennon NJ
    8. Livak KJ
    9. Mikkelsen TS
    10. Rinn JL

    (2014) The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells

    Nature Biotechnology 32:381–386.

    https://doi.org/10.1038/nbt.2859

    • PubMed
    • Google Scholar
33. 1. Wagner A
    2. Regev A
    3. Yosef N

    (2016) Revealing the vectors of cellular identity with single-cell genomics

    Nature Biotechnology 34:1145–1160.

    https://doi.org/10.1038/nbt.3711

    • Google Scholar
34. 1. Waks Z
    2. Klein AM
    3. Silver PA

    (2011) Cell‐to‐cell variability of alternative RNA splicing

    Molecular Systems Biology 7:506.

    https://doi.org/10.1038/msb.2011.32

    • Google Scholar
35. 1. Welch JD
    2. Hu Y
    3. Prins JF

    (2016) Robust detection of alternative splicing in a population of single cells

    Nucleic Acids Research 44:e73.

    https://doi.org/10.1093/nar/gkv1525

    • PubMed
    • Google Scholar
36. 1. Westoby J
    2. Herrera MS
    3. Ferguson-Smith AC
    4. Hemberg M

    (2018) Simulation-based benchmarking of isoform quantification in single-cell RNA-seq

    Genome Biology 19:191.

    https://doi.org/10.1186/s13059-018-1571-5

    • PubMed
    • Google Scholar
37. 1. Westoby J
    2. Artemov P
    3. Hemberg M
    4. Ferguson-Smith A

    (2020) Obstacles to detecting isoforms using full-length scRNA-seq data

    Genome Biology 21:74.

    https://doi.org/10.1186/s13059-020-01981-w

    • PubMed
    • Google Scholar
38. 1. Xiong HY
    2. Barash Y
    3. Frey BJ

    (2011) Bayesian prediction of tissue-regulated splicing using RNA sequence and cellular context

    Bioinformatics 27:2554–2562.

    https://doi.org/10.1093/bioinformatics/btr444

    • PubMed
    • Google Scholar
39. 1. Zenklusen D
    2. Larson DR
    3. Singer RH

    (2008) Single-RNA counting reveals alternative modes of gene expression in yeast

    Nature Structural & Molecular Biology 15:1263–1271.

    https://doi.org/10.1038/nsmb.1514

    • PubMed
    • Google Scholar
40. 1. Zhang X
    2. Xu C
    3. Yosef N

    (2019) Simulating multiple faceted variability in single cell RNA sequencing

    Nature Communications 10:2611–2627.

    https://doi.org/10.1038/s41467-019-10500-w

    • PubMed
    • Google Scholar
41. 1. Ziegenhain C
    2. Vieth B
    3. Parekh S
    4. Reinius B
    5. Guillaumet-Adkins A
    6. Smets M
    7. Leonhardt H
    8. Heyn H
    9. Hellmann I
    10. Enard W

    (2017) Comparative analysis of Single-Cell RNA sequencing methods

    Molecular Cell 65:631–643.

    https://doi.org/10.1016/j.molcel.2017.01.023

    • PubMed
    • Google Scholar

Author details

1. Carlos F Buen Abad Najar

   Center for Computational Biology, University of California, Berkeley, Berkeley, United States

   Contribution

   Conceptualization, Software, Formal analysis, Investigation, Methodology, Writing - original draft

   Competing interests

   No competing interests declared

2. Nir Yosef

   1. Center for Computational Biology, University of California, Berkeley, Berkeley, United States
   2. Department of Electrical Engineering and Computer Science and the Center for Computational Biology, University of California, Berkeley, Berkeley, United States
   3. Ragon Institute of MGH, MIT, and Harvard, Cambridge, United States
   4. Chan Zuckerberg Biohub, San Francisco, San Francisco, United States

   Contribution

   Conceptualization, Formal analysis, Supervision, Methodology, Writing - original draft

   For correspondence

   niryosef@berkeley.edu

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0001-9004-1225

3. Liana F Lareau

   1. Center for Computational Biology, University of California, Berkeley, Berkeley, United States
   2. Department of Bioengineering, University of California, Berkeley, Berkeley, United States

   Contribution

   Conceptualization, Formal analysis, Supervision, Methodology, Writing - original draft

   For correspondence

   lareau@berkeley.edu

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0003-3223-3426

• 5,565

  views

• 576

  downloads

• 37

  citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.