Beyond comparisons of means: understanding changes in gene expression at the single-cell level

A statistical model to detect changes in expression patterns for scRNA-seq data sets

We propose a statistical approach to compare expression patterns between P pre-specified populations of cells. It builds upon BASiCS [8], a Bayesian model for the analysis of scRNA-seq data. As in traditional differential expression analyses, for any given gene i, changes in overall expression are identified by comparing population-specific expression rates ${\mu }_i^{\left(p\right)}$ (p=1,…,P), defined as the relative abundance of gene i within the cells in population p. However, the main focus of our approach is to assess differences in biological cell-to-cell heterogeneity between the populations. These are quantified through changes in population- and gene-specific biological over-dispersion parameters ${\delta }_i^{\left(p\right)}$ (p=1,…,P), designed to capture residual variance inflation (after normalization and technical noise removal) while attenuating the well-known confounding relationship between mean and variance in count-based data sets [9] (a similar concept was defined in the context of bulk RNA-seq by [10], using the term biological coefficient of variation). Importantly, such changes cannot be uncovered by standard differential expression methods, which are restricted to changes in overall expression. Hence, our approach provides novel biological insights by highlighting genes that undergo changes in cell-to-cell heterogeneity between the populations despite the overall expression level being preserved.

To disentangle technical from biological effects, we exploit spike-in genes that are added to the lysis buffer and thence theoretically present at the same amount in every cell (e.g., the 92 ERCC molecules developed by the External RNA Control Consortium [11]). These provide an internal control or gold standard to estimate the strength of technical variability and to aid normalization. In particular, these control genes allow inference on cell-to-cell differences in mRNA content, providing additional information about the analyzed populations of cells [12]. These are quantified through changes between cell-specific normalizing constants ${\varphi }_j^{\left(p\right)}$ (for the jth cell within the pth population). Critically, as described in Additional file 1: Note S1 and Fig. S1, global shifts in mRNA content between populations do not induce spurious differences when comparing gene-specific parameters (provided the offset correction described in ‘Methods’ is applied).

A graphical representation of our model is displayed in Fig. 1 (based on a two-group comparison). It illustrates how our method borrows information across all cells and genes (biological transcripts and spike-in genes) to perform inference. Posterior inference is implemented via a Markov chain Monte Carlo (MCMC) algorithm, generating draws from the posterior distribution of all model parameters. Post-processing of these draws allows quantification of supporting evidence regarding changes in expression patterns (mean and over-dispersion). These are measured using a probabilistic approach based on tail posterior probabilities associated with decision rules, where a probability cut-off is calibrated through the expected false discovery rate (EFDR) [13].

Our strategy is flexible and can be combined with a variety of decision rules, which can be altered to reflect the biological question of interest. For example, if the aim is to detect genes whose overall expression changes between populations p and p^′, a natural decision rule is $|log(\underset{i}{\overset{\left(p\right)}{\mu }}/{\mu }_i^{\left(p^{\prime }\right)}\left)\right|>{\tau }_0$, where τ₀≥0 is an a priori chosen biologically significant threshold for LFCs in overall expression, to avoid highlighting genes with small changes in expression that are likely to be less biologically relevant [6, 14]. Alternatively, changes in biological cell-to-cell heterogeneity can be assessed using $|log(\underset{i}{\overset{\left(p\right)}{\delta }}/{\delta }_i^{\left(p^{\prime }\right)}\left)\right|>{\omega }_0$, for a given minimum tolerance threshold ω₀≥0. This is the main focus of this article. As a default option, we suggest setting τ₀=ω₀=0.4, which roughly coincides with a 50 % increase in overall expression or over-dispersion in whichever group of cells has the largest value (this choice is also supported by the control experiments shown in this article). To improve the interpretation of the genes highlighted by our method, these decision rules can also be complemented by, e.g., requiring a minimum number of cells where the expression of a gene is detected.

More details regarding the model setup and the implementation of posterior inference can be found in ‘Methods’.

Alternative approaches for identifying changes in mean expression

To date, most differential expression analyses of scRNA-seq data sets have borrowed methodology from bulk RNA-seq literature (e.g., DESeq2 [6] and edgeR [5]). However, such methods are not designed to capture features that are specific to SC-level experiments (e.g., the increased levels of technical noise). Instead, BASiCS, SCDE [7] and MAST [15] have been specifically developed with scRNA-seq data sets in mind. SCDE is designed to detect changes in mean expression while accounting for dropout events, where the expression of a gene is undetected in some cells due to biological variability or technical artifacts. For this purpose, SCDE employs a two-component mixture model where negative binomial and low-magnitude Poisson components model amplified genes and the background signal related to dropout events, respectively. MAST is designed to capture more complex changes in expression, using a hurdle model to study both changes in the proportion of cells where a gene is expressed above background and in the positive expression mean, defined as a conditional value – given than the gene is expressed above background levels. Additionally, MAST uses the fraction of genes that are detectably expressed in each cell (the cellular detection rate or CDR) as a proxy to quantify technical and biological artifacts (e.g., cell volume). SCDE and MAST rely on pre-normalized expression counts. Moreover, unlike BASiCS, SCDE and MAST use a definition of changes in expression mean that is conceptually different to what would be obtained based on a bulk population (which would consider all cells within a group, regardless of whether a gene is expressed above background or not).

The performance of these methods is compared in Additional file 1: Note S2 using real and simulated data sets. While control of the false discovery rate (FDR) is not well calibrated for BASiCS when setting τ₀=0, this control is substantially improved when increasing the LFC threshold to τ₀=0.4 – which is the default option we recommend (Additional file 1: Table S1). Not surprisingly, the higher FDR rates of BASiCS lead to higher sensitivity. In fact, our simulations suggest that BASiCS can correctly identify more genes that are differentially expressed than other methods. While this conclusion is based on synthetic data, it is also supported by the analysis of the cell-cycle data set described in [16] (see Additional file 1: Fig. S2), where we observe that SCDE and MAST fail to highlight a large number of genes for which a visual inspection suggests clear changes in overall expression (Additional file 1: Figs. S3 and S4). We hypothesize that this is partly due to conceptual differences in the definition of overall expression and, for MAST, the use of CDR as a covariate.

Alternative approaches for identifying changes in heterogeneity of expression

To the best of our knowledge, BASiCS is the first probabilistic tool to quantify gene-specific changes in the variability of expression between populations of cells. Instead, previous literature has focused on comparisons based on the coefficient of variation (CV), calculated from pre-normalized expression counts (e.g., [17]), for which no quantitative measure of differential variability has been obtained. More recently, [9] proposed a mean-corrected measure of variability to avoid the confounding effect between mean expression and CV. Nonetheless, the latter was designed to compare expression patterns for sets of genes, rather than for individual genes.

Not surprisingly, our analysis suggests that a quantification of technical variability is critical when comparing variability estimates between cell populations (Additional file 1: Note S3 and Fig. S5). In particular, comparisons based on CV estimates can mask the biological signal if the strength of technical variability varies between populations.

A control experiment: comparing single cells vs pool-and-split samples

To demonstrate the efficacy of our method, we use the control experiment described in [17], where single mESCs are compared against pool-and-split (P&S) samples, consisting of pooled RNA from thousands of mESCs split into SC equivalent volumes. Such a controlled setting provides a situation where substantial changes in overall expression are not expected as, on average, the overall expression of SCs should match the levels measured in P&S samples. Additionally, the design of P&S samples should remove biological variation, leading to a homogeneous set of samples. Hence, P&S samples are expected to show a genuine reduction in biological cell-to-cell heterogeneity compared to SCs.

Here, we display the analysis of samples cultured in a 2i media. Hyper-parameter values for ${\mu }_i^{\left(p\right)}$’s and ${\delta }_i^{\left(p\right)}$’s were set to $a_{\mu }^2=a_{\delta }^2=0.5$, so that extreme LFC estimates are shrunk towards (−3,3) (see ‘Methods’). However, varying $a_{\mu }^2$ and $a_{\delta }^2$ leads to almost identical results (not shown), suggesting that posterior inference is in fact dominated by the data. In these data, expression counts correspond to the number of molecules mapping to each gene within each cell. This is achieved by using unique molecular identifiers (UMIs), which remove amplification biases and reduce sources of technical variation [18]. Our analysis includes 74 SCs and 76 P&S samples (same inclusion criteria as in [17]) and expression counts for 9378 genes (9343 biological and 35 ERCC spikes) defined as those with at least 50 detected molecules in total across all cells. The R code used to perform this analysis is provided in Additional file 2.

To account for potential batch effects, we allowed different levels of technical variability to be estimated in each batch (see Additional file 1: Note S4 and Fig. S6). Moreover, we also performed an independent analysis of each batch of cells. As seen in Additional file 1: Fig. S7, the results based on the full data are roughly replicated in each batch, suggesting that our strategy is able to remove potential artifacts related to this batch effect.

As expected, our method does not reveal major changes in overall expression between SCs and P&S samples as the distribution of LFC estimates is roughly symmetric with respect to the origin (see Fig. 2a) and the majority of genes are not classified as differentially expressed at 5 % EFDR (see Fig. 3b). However, this analysis suggests that setting the minimum LFC tolerance threshold τ₀ equal to 0 is too liberal as small LFCs are associated with high posterior probabilities of changes in expression (see Fig. 3a) and the number of differentially expressed genes is inflated (see Fig. 3b). In fact, counter-intuitively, 4710 genes (≈50 % of all analyzed genes) are highlighted to have a change in overall expression when using τ₀=0. This is partially explained by the high nominal FDR rates displayed in Additional file 1: Note S2.1 where, for τ₀=0, FDR is poorly calibrated when simulating under the null model. In addition, we hypothesize this heavy inflation is also due to small but statistically significant differences in expression that are not biologically meaningful. In fact, the number of genes whose overall expression changes is reduced to 559 (≈6 % of all analyzed genes) when setting τ₀=0.4. As discussed earlier, this minimum threshold roughly coincides with a 50 % increase in overall expression and with the 90th percentile of empirical LFC estimates when simulating under the null model (no changes in expression). Posterior inference regarding biological over-dispersion is consistent with the experimental design, where the P&S samples are expected to have more homogeneous expression patterns. In fact, as shown in Fig. 2b, the distribution of estimated LFCs in biological over-dispersion is skewed towards positive values (higher biological over-dispersion in SCs). This is also supported by the results shown in Fig. 3b, where slightly more than 2000 genes exhibit increased biological over-dispersion in SCs and almost no genes (≈60 genes) are highlighted to have higher biological over-dispersion in the P&S samples (EFDR = 5 %). In this case, the choice of ω₀ is less critical (within the range explored here). This is illustrated by the left panels in Fig. 3a, where tail posterior probabilities exceeding the cut-off defined by EFDR = 5 % correspond to similar ranges of LFC estimates.

mESCs across different cell-cycle stages

Our second example shows the analysis of the mESC data set presented in [16], which contains cells where the cell-cycle phase is known (G1, S and G2M). After applying the same quality control criteria as in [16], our analysis considers 182 cells (59, 58 and 65 cells in stages G1, S and G2M, respectively). To remove genes with consistently low expression across all cells, we excluded those genes with less than 20 reads per million (RPM), on average, across all cells. After this filter, 5,687 genes remain (including 5,634 intrinsic transcripts, and 53 ERCC spike-in genes). The R code used to perform this analysis is provided in Additional file 3.

As a proof of concept, to demonstrate the efficacy of our approach under a negative control, we performed permutation experiments, where cell labels were randomly permuted into three groups (containing 60, 60 and 62 samples, respectively). In this case, our method correctly infers that mRNA content as well as gene expression profiles do not vary across groups of randomly permuted cells (Fig. 4).

As cells progress through the cell cycle, cellular mRNA content increases. In particular, our model infers that mRNA content is roughly doubled when comparing cells in G1 vs G2M, which is consistent with the duplication of genetic material prior to cell division (Fig. 5a). Our analysis suggests there are no major shifts in expression levels between cell-cycle stages (Fig. 5b and upper triangular panels in Fig. 5d). Nonetheless, a small number of genes are identified as displaying changes in overall expression between cell-cycle phases at 5 % EFDR for τ₀=0.4 (Fig. 6). To validate our results, we performed gene ontology (GO) enrichment analysis within those genes classified as differentially expressed between cell-cycle phases (see Additional file 3). Not surprisingly, we found an enrichment of mitotic genes among the 545 genes classified as differentially expressed between G1 and G2M cells. In addition, the 209 differentially expressed genes between S and G2M are enriched for regulators of cytokinesis, which is the final stage of the cell cycle where a progenitor cell divides into two daughter cells [19].

Our method suggests a substantial decrease in biological over-dispersion when cells move from G1 to the S phase, followed by a slight increase after the transition from S to the G2M phase (see Fig. 5c and the lower triangular panels in Fig. 5d). This is consistent with the findings in [19], where the increased gene expression variability observed in G2M cells is attributed to an unequal distribution of genetic material during cytokinesis and the S phase is shown to have the most stable expression patterns within the cell cycle. Here, we discuss GO enrichment of those genes whose overall expression rate remains constant (EFDR = 5 %, τ₀=0.4) but that exhibit changes in biological over-dispersion between cell-cycle stages (EFDR = 5 %, ω₀=0.4). Critically, these genes will not be highlighted by traditional differential expression tools, which are restricted to differences in overall expression rates. For example, among the genes with higher biological over-dispersion in G1 with respect to the S phase, we found an enrichment of genes related to protein dephosphorylation. These are known regulators of the cell cycle [20]. Moreover, we found that genes with lower biological over-dispersion in G2M cells are enriched for genes related to DNA replication checkpoint regulation (which delays entry into mitosis until DNA synthesis is completed [21]) relative to G1 cells and mitotic cytokinesis when comparing to S cells. Both of these processes are likely to be more tightly regulated in the G2M phase. A full table with GO enrichment analysis of the results described here is provided in Additional file 3.

A statistical model to detect changes in expression patterns for scRNA-seq data sets

In this article, we introduce a statistical model for identifying genes whose expression patterns change between predefined populations of cells (given by experimental conditions or cell types). Such changes can be reflected via the overall expression level of each gene as well as through changes in cell-to-cell biological heterogeneity. Our method is motivated by features that are specific to scRNA-seq data sets. In this context, it is essential to normalize and remove technical artifacts appropriately from the data before extracting the biological signal. This is particularly critical when there are substantial differences in cellular mRNA content, amplification biases and other sources of technical variation. For this purpose, we exploit technical spike-in genes, which are added at the (theoretically) same quantity to each cell’s lysate. A typical example is the set of 92 ERCC molecules developed by the External RNA Control Consortium [11]. Our method builds upon BASiCS [8] and can perform comparisons between multiple populations of cells using a single model. Importantly, our strategy avoids stepwise procedures where data sets are normalized prior to any downstream analysis. This is an advantage over methods using pre-normalized counts, as the normalization step can be distorted by technical artifacts.

We assume that there are P groups of cells to be compared, each containing n_p cells (p=1,…,P). Let $X_{\mathit{\text{ij}}}^{\left(p\right)}$ be a random variable representing the expression count of a gene i (i=1,…,q) in the jth cell from group p. Without loss of generality, we assume the first q₀ genes are biological and the remaining q−q₀ are technical spikes. Extending the formulation in BASiCS, we assume that

with ${\mu }_i^{\left(p\right)}\equiv {\mu }_i$ for i=q₀+1,…,q and where CV stands for coefficient of variation (i.e., the ratio between standard deviation and mean). These expressions are the result of a Poisson hierarchical structure (see Additional file 1: Note S6.1). Here, ${\varphi }_j^{\left(p\right)}$’s act as cell-specific normalizing constants (fixed effects), capturing differences in input mRNA content across cells (reflected by the expression counts of intrinsic transcripts only). A second set of normalizing constants, $s_j^{\left(p\right)}$’s, capture cell-specific scale differences affecting the expression counts of all genes (intrinsic and technical). Among others, these differences can relate to sequencing depth, capture efficiency and amplification biases. However, a precise interpretation of the $s_j^{\left(p\right)}$’s varies across experimental protocols, e.g., amplification biases are removed when using UMIs [18]. In addition, θ_p’s are global technical noise parameters controlling the over-dispersion (with respect to Poisson sampling) of all genes within group p. The overall expression rate of a gene i in group p is denoted by ${\mu }_i^{\left(p\right)}$. These are used to quantify changes in the overall expression of a gene across groups. Similarly, the ${\delta }_i^{\left(p\right)}$’s capture residual over-dispersion (beyond what is due to technical artifacts) of every gene within each group. These so-called biological over-dispersion parameters relate to heterogeneous expression of a gene across cells. For each group, stable housekeeping-like genes lead to ${\delta }_i^{\left(p\right)}\approx 0$ (low residual variance in expression across cells) and highly variable genes are linked to large values of ${\delta }_i^{\left(p\right)}$. A novelty of our approach is the use of ${\delta }_i^{\left(p\right)}$ to quantify changes in biological over-dispersion. Importantly, this attenuates confounding effects due to changes in overall expression between the groups.

A graphical representation of this model is displayed in Fig. 1. To ensure identifiability of all model parameters, we assume that ${\mu }_i^{\left(p\right)}$’s are known for the spike-in genes (and given by the number of spike-in molecules that are added to each well). Additionally, we impose the identifiability restriction

Here, we discuss the priors assigned to parameters that are gene- and group-specific (see Additional file 1: Note S6.2 for the remaining elements of the prior). These are given by

Hereafter, without loss of generality, we simplify our notation to focus on two-group comparisons. This is equivalent to assigning Gaussian prior distributions for LFCs in overall expression (τ_i) or biological over-dispersion (ω_i). In such a case, it follows that

Hence, our prior is symmetric, meaning that we do not a priori expect changes in expression to be skewed towards either group of cells. Values for $a_{\mu }^2$ and $a_{\delta }^2$ can be elicited using an expected range of values for LFC in expression and biological over-dispersion, respectively. The latter is particularly useful in situations where a gene is not expressed (or very lowly expressed) in one of the groups, where, e.g., LFCs in overall expression are undefined (the maximum likelihood estimate of τ_i would be ±∞, the sign depending on which group expresses gene i). A popular solution to this issue is the addition of pseudo-counts, where an arbitrary number is added to all expression counts (in all genes and cells). This strategy is also adopted in models that are based on log-transformed expression counts (e.g., [15]). While the latter guarantees that τ_i is well defined, it leads to artificial estimates for τ_i (see Table 1). Instead, our approach exploits an informative prior (indexed by $a_{\mu }^2$) to shrink extreme estimates of τ_i towards an expected range. This strategy leads to a meaningful shrinkage strength, which is based on prior knowledge. Importantly – and unlike the addition of pseudo-counts – our approach is also helpful when comparing biological over-dispersion between the groups. In fact, if a gene i is not expressed in one of the groups, this will lead to a non-finite estimate of ω_i (if all expression counts in a group are equal to zero, the corresponding estimate of the biological over-dispersion parameters would be equal to zero). Adding pseudo-counts cannot resolve this issue, but imposing an informative prior for ω_i (indexed by $a_{\omega }^2$) will shrink estimates towards the appropriate range.

Generally, posterior estimates of τ_i and ω_i are robust to the choice of $a_{\mu }^2$ and $a_{\delta }^2$, as the data is informative and dominates posterior inference. In fact, these values are only influential when shrinkage is needed, e.g., when there are zero total counts in one of the groups. In such cases, posterior estimates of τ_i and ω_i are dominated by the prior, yet the method described below still provides a tool to quantify evidence of changes in expression. As a default option, we use $a_{\mu }^2=a_{\delta }^2=0.5$ leading to τ_i,ω_i∼ N(0,1). These default values imply that approximately 99 % of the LFCs in overall expression and over-dispersion are expected a priori to lie in the interval (−3,3). This range seems reasonable in light of the case studies we have explored. If a different range is expected, this can be easily modified by the user by setting different values for $a_{\mu }^2$ and $a_{\delta }^2$.

Posterior samples for all model parameters are generated via an adaptive Metropolis within a Gibbs sampling algorithm [25]. A detailed description of our implementation can be found in Additional file 1: Note S6.3.

Post hoc correction of global shifts in input mRNA content between the groups

The identifiability restriction in Eq. 3 applies only to cells within each group. As a consequence, if they exist, global shifts in cellular mRNA content between groups (e.g., if all mRNAs were present at twice the level in one population related to another) are absorbed by the ${\mu }_i^{\left(p\right)}$’s. To assess changes in the relative abundance of a gene, we adopt a two-step strategy where: (1) model parameters are estimated using the identifiability restriction in Eq. 3 and (2) global shifts in endogenous mRNA content are treated as a fixed offset and corrected post hoc. For this purpose, we use the sum of overall expression rates (intrinsic genes only) as a proxy for the total mRNA content within each group. Without loss of generality, we use the first group of cells as a reference population. For each population p (p=1,…,P), we define a population-specific offset effect:

and perform the following offset correction:

This is equivalent to replacing the identifiability restriction in Eq. 3 by

Technical details regarding the implementation of this post hoc offset correction are explained in Additional file 1: Note S6.4. The effect of this correction is illustrated in Fig. 7 using the cell-cycle data set described in the main text. As an alternative, we also explored the use of the ratio between the total intrinsic counts over total spike-in counts to define a similar offset correction based on

For the cell-cycle data set, both alternatives are equivalent. Nonetheless, the first option is more robust in cases where a large number of differentially expressed genes are present. Hereafter, we use ${\mu }_i^{\left(p\right)}$ and ${\varphi }_j^{\left(p\right)}$ to denote ${\tilde{\mu }}_i^{\left(p\right)}$ and ${\tilde{\varphi }}_j^{\left(p\right)}$, respectively.

A probabilistic approach to quantify evidence of changes in expression patterns

A probabilistic approach is adopted, assessing changes in expression patterns (mean and over-dispersion) through a simple and intuitive scale of evidence. Our strategy is flexible and can be combined with a variety of decision rules. In particular, here we focus on highlighting genes whose absolute LFC in overall expression and biological over-dispersion between the populations exceeds minimum tolerance thresholds τ₀ and ω₀, respectively (τ₀,ω₀≥0), set a priori. The usage of such minimum tolerance levels for LFCs in expression has also been discussed in [14] and [6] as a tool to improve the biological significance of detected changes in expression and to improve upon FDRs.

For a given probability threshold ${\alpha }_{{}_M}$ ($0.5<{\alpha }_{{}_M}<1$), a gene i is identified as exhibiting a change in overall expression between populations p and p^′ if

If τ₀→0, ${\pi }_i^M\left({\tau }_0\right)\to 1$ becoming uninformative to detect changes in expression. As in [26], in the limiting case where τ₀=0, we define

with

A similar approach is adopted to study changes in biological over-dispersion between populations p and p^′, using

for a fixed probability threshold ${\alpha }_{{}_D}$ ($0.5<{\alpha }_{{}_D}<1$). In line with Eqs. 11 and 12, we also define

with

Evidence thresholds ${\alpha }_{{}_M}$ and ${\alpha }_{{}_D}$ can be fixed a priori. Otherwise, these can be defined by controlling the EFDR [13]. In our context, these are given by

and

where I(A)=1 if event A is true, 0 otherwise. Critically, the usability of this calibration rule relies on the existence of genes under both the null and the alternative hypothesis (i.e., with and without changes in expression). While this is not a practical limitation in real case studies, this calibration might fail to return a value in benchmark data sets (e.g., simulation studies), where there are no changes in expression. As a default, if EFDR calibration is not possible, we set ${\alpha }_{{}_M}={\alpha }_{{}_D}=0.90$.

The posterior probabilities in Eqs. 10, 11, 13 and 14 can be easily estimated – as a post-processing step – once the model has been fitted (see Additional file 1: Note S6.5). In addition, our strategy is flexible and can be easily extended to investigate more complex hypotheses, which can be defined post hoc, e.g., to identify those genes that show significant changes in cell-to-cell biological over-dispersion but that maintain a constant level of overall expression between the groups, or conditional decision rules where we require a minimum number of cells where the expression of a gene is detected.

Software

Our implementation is freely available as an R package [27], using a combination of R and C++ functions through the Rcpp library [28]. This can be found in https://github.com/catavallejos/BASiCS, released under the GPL license.