Algorithm and implementation
D3E takes a read-count table as an input, with rows and columns corresponding to transcripts and cells, respectively. The user should split the columns into two or more groups by providing cell labels in the input file. If there are more than two groups of cells, they must be compared one pair at a time. D3E uses four steps to process the data. First, input data is normalised using the same algorithm as DESeq2 (see Implementation) and filtered by removing the genes that are not expressed in any of the cells. Second, the Cramér-von Mises test, the Kolmogorov-Smirnov (KS) test, or the likelihood ratio test [11] is used to identify the genes with a significant change in expression between the two samples of interest. Third, the transcriptional bursting model is fitted to the expression data for each gene in both samples using either the method of moments or a Bayesian method [17]. Fourth, the change in parameters between the two samples is calculated for each gene (Fig. 1
d).
A command-line version of D3E written in Python can be downloaded from GitHub (
https://hemberg-lab.github.io/D3E
), and the source code is available under the GPL licence. Furthermore, there is also a web-version available at
http://www.sanger.ac.uk/sanger/GeneRegulation_D3E
. Due to the time required to run D3E, the web version limits the number of genes and cells that may be analyzed, and it can only use the method of moments for estimating parameters.
DE analysis module
To compare distributions obtained from two different sets of cells, D3E uses either the Cramér-von Mises test, the KS test or the likelihood ratio test to quantify the difference in gene expression (see Implementation). The first two tests are non-parametric which is advantageous since it allows us to apply D3E to any single-cell dataset, not just the ones collected using RNA-seq. The null hypothesis for all three tests is that the two samples are drawn from the same distribution. The premise of D3E is that when two samples are drawn from the same population of cells, the test should result in a high p-value. On the other hand, if the cells are drawn from two populations with different transcript distributions, then the resulting p-value should be low.
We first evaluated D3E using synthetic data. Fortunately, there is a widely used, experimentally validated stochastic model available for single-cell gene expression [22]. We refer to this model as the transcriptional bursting model (Fig. 1
a), and it is characterized by three parameters: α, the rate of promoter activation; β, the rate of promoter inactivation; γ, the rate of transcription when the promoter is in the active state; and a transcript degradation rate λ. The stationary distribution of the transcriptional bursting model takes the form of a Poisson-Beta mixture distribution [17, 22]
$$\begin{array}{@{}rcl@{}} PB(n \; | \; \alpha, \beta, \gamma, \lambda) & = & \text{Poi}(n \; | \; \gamma x/\lambda) \bigwedge_{x} \text{Beta}(x \;|\; \alpha, \beta) \\ & = & \frac{\gamma^{n}e^{-\gamma/\lambda}\Gamma(\alpha/\lambda+n)\Gamma(\alpha/\lambda + \beta/\lambda)}{\lambda^{n}\Gamma(n+1)\Gamma(\alpha/\lambda + \beta/\lambda + n)\Gamma(\alpha/\lambda)}\\&&\Phi\left(\frac{\alpha}{\lambda}, \frac{\alpha}{\lambda} + \frac{\beta}{\lambda} + n; \frac{\gamma}{\lambda}\right), \end{array} $$
where n is the number of transcripts of a particular gene, x is an auxiliary variable, Γ is the Euler Gamma function, and Φ(a,b;c) is the confluent hypergeometric function.
An important feature of the Poisson-Beta distribution is that the three parameters α,β and γ are normalised by the rate of mRNA degradation λ. Consequently, when fitting the parameters for the Poisson-Beta distribution from a stationary sample, only three parameters can be estimated, and they are unique up to a common multiplicative factor, λ. Since a single-cell RNA-seq experiment corresponds to a snapshot of individual cells, it is often reasonable to assume that the samples are drawn from the stationary distribution. The inability to uniquely identify all four parameters from single-cell RNA-seq data means that it is only approprite to apply the transcriptional bursting model to DE analysis when λs are constant between the compared samples, or when the degradation rates are known for both samples. Without knowledge of λ it is impossible to unambiguously determine how the parameters have changed.
To evaluate the sensitivity of the Cramér-von Mises, the KS and the likelihood ratio test to changes in the parameters, we selected triplets of parameters (α,β,γ) from a range that is characteristic for single-cell RNA-seq data [15]. For each parameter triplet one of the parameters was varied, while fixing the remaining two, and a series of tests was carried out on the corresponding Poisson-Beta samples. For each combination of parameters, we assumed that there were 50 cells from each condition when generating the data. The results can be summarized by a set of matrices, where rows and columns correspond to values of the varied parameter, and the elements in the matrix are p-values from the tests (Fig. 2
c). Ideally we would like to find high p-values on the diagonal and low p-values away from the diagonal. We used a heuristic for characterizing the pattern of p-values, and for each matrix we obtained a single score, S (see Implementation). When S=0, high p-values are only found on the diagonal, suggesting that D3E has successfully identified genes where there was a change in one of the parameters.
The results suggest that all three tests are capable of accurately detecting changes in the parameters in certain regions of the parameter space (Fig. 2
a, Additional file 1). For all three tests, changes in β are the most difficult to detect while changes in γ are the easiest to identify. The methods perform poorly when γ is small and either β is large or α is small. In this regime, the Poisson-Beta distribution is similar to the Poisson distribution with a mean close to zero, and it is challenging to identify which parameter has changed, and by how much. From a biological perspective, when a transcription rate is small and a gene has a small duty cycle (small α or big β) there are almost no transcripts produced since the promoter spends most of its time in the inactive state. Therefore, changes in either of the three parameters will be difficult to distinguish. The performance of each method can be summarized by the average score across all combination of parameters and we find that the likelihood ratio test is the most accurate, followed by the Cramér-von Mises test and the KS test (Fig. 2
b, S1). The accuracy, however, is also mirrored by the computational costs; analyzing the data in Fig. 2
a takes about 18 s for the KS-test, 30 min for the Cramér-von Mises test and 8.3 h for the likelihood ratio test using a MacBook Pro laptop with a 3.2 GHz Intel i5 processor and 16 Gb of RAM.
We also considered the scenario when the two distributions are different, but the mean is identical. This is a situation where it is all but impossible for methods which only use the mean to reliably detect that there has been a change in the expression profile. In contrast, D3E is able to reliably identify a change in expression. Our results show that a change of α and β by a factor of 2, which is roughly equivalent to changing the variance by the same factor, is sufficient for the p-value to drop below.05 for a sample of 50 cells (Fig. 2
d).
A particular challenge for DE analysis is to determine the p-value threshold for when a change can be considered significant. The traditional approach is to use a fixed value, e.g..05, and then adjust for multiple hypothesis testing. For D3E two different methods are offered for selecting the critical p-value. On the one hand, the user can specify a false discovery rate, and the tool will use the Benjamini-Hochberg procedure [4] to identify the critical p-value. Alternatively, one may take an empirical approach whereby one of the two samples is first split into two parts. By definition, the two parts should have identical distributions for all genes, which means that it can be used as a negative control. D3E applies one of the three tests to the negative control, and records the lowest p-value identified, p
∗. When comparing the two original distributions, only genes with a p-value below a×p
∗ are considered significant, where the default value for the parameter a is.1. We advice the user that the reduced sample size of the negative control set is likely to result in a less stringent cut-off than what would be expected if the negative control had the same sample size as the original data. To evaluate the heuristic strategy, we generated 1,000 pair of samples with the same number of reads and cells, using identical parameter values for the samples in each pair. Using the Cramér-von Mises test we recorded the lowest observed p-value, p
∗, and we used.1×p
∗ as a threshold for when to call a test significant. For both the method of moments and the Bayesian method, we found that 97 % of the genes were detected as not DE. The control experiment demonstrates that D3E is capable of accurately distinguishing situations where the parameters are unchanged.
To further evaluate the performance of D3E relative to other DE methods, we generated additional synthetic data sets where one of the three parameters was varied while the other two were fixed as before. For each data set we designated genes as DE where the parameter had changed by at least a factor of 1.5, 2 or 4. The arbitrary decision of what constitutes a significant change allows us to define the calls of the DE algorithms as either true positive, false positive, true negative or false negative. The results can be summarized as a ROC curve, and again we find that changes in β are more difficult to detect compared to the other two parameters (Fig. 3). Importantly, we find that for larger parameter changes, D3E is always amongst the best performing methods (Fig. 3).
Parameter estimation module
The likelihood ratio test is a parametric test, and thus it requires estimates of the parameters α, β and γ. D3E has two approaches; the method of moments which is relatively fast and a more computationally costly Bayesian inference method (see Implementation) [17]. We evaluated the accuracy of the parameter estimates, both assuming perfect sampling as well as in the scenario when some transcripts are lost due to low starting levels of mRNA [16]. In our simulations, we first generated parameters for the Poisson-Beta model by drawing from a distribution similar to the one for the data reported by Islam et al. [15]. We then assumed that the dropout probability scaled with the mean expression level as \(p_{dropout} = 1 - e^{-\mu ^{2}/b}\) [23], where the parameter b controls the rate of dropouts. In our simulations we used b=10,100 and 200. Our simulations show that the relative error for the estimates of the parameters α, β and γ are relatively robust to dropout events (Fig. 4
a).
An important advantage of using the transcriptional bursting model (Fig. 1) is that it is possible to derive other quantities - the average burst size, the burst frequency, the mean expression level, and the proportion of time in the active state (duty cycle) - which are easier to measure and interpret biologically than the parameters α, β and γ. Importantly, the transcriptional bursting model allows us to learn more about how the expression level has changed between the two conditions. In the transcriptional bursting model, there are three different ways to increase the mean expression level; by decreasing the degradation rate, by increasing the burst frequency, or by increasing the burst size. By comparing the correlations between the change in mean expression levels and the change in burst size or burst frequency, it is possible to gain additional biological insights relating to what aspect has led to the change in gene expression. In contrast to the estimates of the rate parameters the correlation estimates are sensitive to dropout events (Fig. 4
b), and one must thus be careful when interpreting the results.
Application to experimental data
The tests on synthetic data suggest that D3E can reliably identify differentially expressed genes. A more useful test of the algorithm, however, involves experimental data which has been reliably validated. Unlike bulk data [27], unfortunately there are no gold-standard datasets available. Nonetheless, to further evaluate D3E, we considered the single-cell RNA-seq data from two and four-cell mouse embryos where qPCR data from the same cell-types was collected for 90 genes [6]. Unfortunately, the correlation of changes in gene expression between the qPCR and RNA-seq data (ρ
Δ
=.46) (Additional file 2: Figure S2) is even worse than the correlation of the individual samples (ρ
2=.6, ρ
4=.5). Thus, it does not come as a surprise that the overlap between the genes which are considered DE in the qPCR experiment has little overlap with genes which are considred DE from RNA-seq by any of the five algorithms that we compared (Additional file 3: Table S1). Even so, we find large differences in the number of genes identified as DE, ranging from 1 (edgeR) to 35 (DESeq2).
To further evaluate D3E, we applied it to the two datasets collected by Islam et al. [15] from 48 mouse embryonic stem cells and 44 mouse embryonic fibroblasts. To establish the p-value threshold, we first separated the stem cells into two groups, and compared the expression (see Implementation). We used this approach for determining the threshold for D3E, SCDE, edgeR and limma, while for DESeq2, we used the adjusted p-value reported by the software. When comparing the two cell-types, D3E identified 4716 genes as DE, DESeq2 identified 6360 genes, limma-voom identified 7245 genes, edgeR identified 1140 genes, and SCDE identified 1092 genes (Fig. 5
d). Surprisingly, the agreement between the five methods is quite low with only a core set of 380 genes identified by all three methods. If we require a gene to be identified of 4 out of 5 methods, then an additional 495 genes are idenitified as DE, suggesting that there is a set of around 900 genes which can confidently considered DE. To further evaluate the set of genes identified as DE by each method, we investigated the distiribution of fold-change values (Fig. 5
a). The distributions gives an indication of how large fold changes are required for detection, and we note most of the genes have a higher expression in fibroblasts compared to stem cells. Compared to DESeq2, SCDE and edgeR, we also notice that D3E is able to identify genes with a lower fold change. Indeed, there were several examples of genes where the change in mean expression was modest, but they were still identified by D3E as differentially expressed (Fig. 5
c).
Next, we took advantage of the transcriptional bursting model underlying D3E, and we fitted the parameters α, β, and γ for all genes. We found that for 85 % of the genes, at least one of the parameters changed by at least 2-fold, suggesting that there are substantial differences between the two cell-types. The results show that all three parameters follow log-normal distributions, spanning approximately one or two orders of magnitude in both cell-types (Fig. 5
b, Additional files 4 and 5). With the exception of the duty cycle which is constrained to be in the interval [0, 1], the derived quantities showed a similar distribution.
We calculated the three derived quantities for each condition for the 2105 genes where we were able to obtain degradation rates for both cell-types [31, 32]. Next, we compared the changes in degradation rate, burst frequency and burst size to the change in mean expression level (Fig. 5
e). The results clearly demonstrate that it is the change in burst size which underlies the change in mean expression levels (ρ=.91), suggesting that altering the burst size is the driving mechanism behind differences in mean expression between conditions. However, considering our simulation experiments (Fig. 4
b), it is likely that the true correlation is lower.
Another property of interest is the coefficient of variation (CV), defined as the standard deviation divided by the mean, which is used to quantify the gene expression noise. The CV is inversely correlated with the mean, and the transcriptional bursting model reveals that the change in CV is mainly correlated with the change in the duty cycle (ρ=.47), while the effect of a change in burst size is considerably smaller (ρ=.24, Additional file 6: Figure S3). To further demonstrate the use of the transcriptional bursting module, we also investigated changes in the auto-correlation times of each gene. The auto-correlation provides information about the time-scale of the noise, i.e. how quickly the gene expression level varies. The expected value of the autocorrelation, τ
c
, is given by (Implementation)
$$ \tau_{c} = \frac{\sigma^{2}}{\frac{\mu}{\lambda}+\frac{1}{\alpha + \beta}}. $$
((1))
Comparison of τ
c
and the change in the mean for the Islam et al. data reveals that the two quantities are strongly correlated (ρ=.87, Additional file 7: Figure S4). However, when investigating all the quantities on the right hand side of Eq. (1) the comparison shows that it is the change in variance which is most strongly correlated with the change in autocorrelation times (ρ=.90, Additional file 7: Figure S4).
We also applied the methodology to cells from the early and late blastocyst from mouse embryos [10]. Since we do not have access to degradation rates for these cell-types, not all of the correlations can be explored. Nevertheless, the correlations that can be calculated are largely similar to the ones observed for the Islam et al. data (Additional file 8: Figure S5–Additional file 9: Figure S6). Taken together, these results demonstrate that it is possible to generate testable hypotheses about how changes in the property of a gene has come about.
Discussion
DE analysis is one of the most common uses of bulk RNA-seq, and we expect that it will become an important application for single-cell RNA-seq as well. Here, we have presented D3E, a tool for analyszing DE for single-cell data. The main difference between D3E and other methods is that D3E compares the full distribution of each gene rather than just the first moment. Therefore, it becomes possible to identify genes where the higher moments have changed, with the mean remaining constant. To the best of our knowledge, D3E is the first method for DE analysis which considers properties other than the first and second moment of the distribution. Using synthetic data, we demonstrate that D3E can reliably detect when only the shape, but not the mean is changed.
One of the main challenges in developing a DE analysis method for single-cell RNA-seq data is that, unlike for bulk data, there are no gold-standards available [27]. Comparison of qPCR and RNA-seq data revealed only a modest correlation between the two methods, implying that the two methods are inconsistent. Thus, one must resort to synthetic data for evaluation. Fortunately, for single-cell gene expression, there is an analytically tractable transcriptional bursting model available which has been experimentally validated. Even with synthetic data, however, it is not obvious how one should define a change in expression. Consider the situation where one of the parameters changes by a small amount which is just sufficient to be detected given the limits of the technical noise, the read depth and the sample size. Then the question is whether or not the change is sufficient to be biologically meaningful. Another challenge stems from the difficulty of disentangling the technical and the biological noise. The transcriptional bursting model does not account for the technical noise in single-cell experiments which can be considerable [3, 7, 13, 28]. Since each cell can only be sequenced once, one cannot carry out technical replicates in the same way that can be done for bulk experiments.
Given the lack of ground truth and the difficulties involved in rigorously defining change in gene expression at the single-cell level, we recommend that more than one algorithm is used to identify DE genes. If the aim is to identify the genes which are most likely to have changed significantly, we believe that a consensus approach is the best one to use. Such a strategy will minimize the number of false positives with the drawback of increasing the number of false negatives.
D3E implements two different non-parametric methods and one parametric method for comparing probability distributions. The methods emphasize different aspects of the distributions, and they are also associated with different computational costs. An important future research question is to determine what method is the most appropriate for single-cell DE analysis. Since the results are sensitive to technical noise, such developments should ideally be carried out taking the specific details of the protocol into consideration.
We have shown that the transcriptional bursting model makes it possible to extract additional, biologically relevant results from the DE analysis. However, to be able to fully utilize the transcriptional bursting model, the mRNA degradation rates must be known, or assumed to be constant. Determining degradation rates directly remains experimentally challenging, and today they are only available for a handful of cell-types. However, alternative strategies have been proposed, whereby degradation rates are estimated from RNA-seq data using distribution of reads along the length of a gene [14, 37]. The RNA-seq based methods make it possible to estimate degradation rates without further experiments, and they could thus significantly expand the range of samples where the transcriptional bursting model can be applied. Another restriction of our method is that the groups of cells must be known in advance. In many cases, e.g. when comparing samples from a mutant and wild-type or when comparing different stages of development [6], it is straightforward to assign labels to cells. However, there are also situations, e.g. when analyzing a tissue-sample, when the cell-labels are unknown and in these scenarios D3E is no longer applicable.