Pierre Boyeau,Jeffrey Regier,Adam Gayoso,Michael I. Jordan,Romain Lopez,Nir Yosef
An empirical Bayes method for differential expression analysis of single cells with deep generative models
differential expression,scRNA-seq,deep generative modeling
Significance Characterizing cell identities and the effects of perturbations through the analysis of expression differences between groups of cells is a crucial challenge for single-cell omics. Conducting such a differential expression (DE) analysis with single-cell RNA sequencing data has been challenging for several reasons, including low data sensitivity, noise, and bias which are, for example, induced by batch effects and differences between technologies. Our work offers lvm-DE—a principled framework enabling the use of deep generative models as a means for quantifying and assessing the significance of DE. lvm-DE provides interpretable discoveries and corrects for batch effects. We empirically validate lvm-DE performance on both synthetic and real data.
An empirical Bayes method for differential expression analysis of single cells with deep generative models

Characterizing cell identities and the effects of perturbations through the analysis of expression differences between groups of cells is a crucial challenge for single-cell omics. Conducting such a differential expression (DE) analysis with single-cell RNA sequencing data has been challenging for several reasons, including low data sensitivity, noise, and bias which are, for example, induced by batch effects and differences between technologies. Our work offers lvm-DE—a principled framework enabling the use of deep generative models as a means for quantifying and assessing the significance of DE. lvm-DE provides interpretable discoveries and corrects for batch effects. We empirically validate lvm-DE performance on both synthetic and real data.

Single-cell omics measurements promise to resolve the heterogeneity of cellular identities, characterizing subtle molecular changes between cells (1). In particular, the measurement of gene expression in single cells—with single-cell RNA-sequencing (scRNA-seq)—creates opportunities to characterize differences in complex scenarios such as rare cell-type detection or inter-cell-type/intercondition comparative analyses. However, technical factors such as the sequencing depth and batch effects obscure the underlying biological signal, requiring careful probabilistic modeling of the data to allow accurate detection of differentially expressed genes. Many differential expression (DE) methods have been developed to handle the specificity of scRNA-seq data (2–4), often building on earlier work on bulk RNA-seq (5, 6). However, there is a gap between the way these methods are applied and the structural characteristics of modern datasets. First, the performance of these approaches in modern datasets is highly dependent on their ability to correct for batch effects (7). Indeed, a large number of these approaches rely on generalized linear models (GLMs), while recent benchmarking studies in scRNA-seq data integration focus on nonlinear methods [e.g., MNNs (8)] for their evaluation (9). The linear method present in that benchmarking study, Combat (10), was found to be suitable only for simple scenarios. Second, most methods for DE analysis are unable to take advantage of the large sample sizes of current scRNA-seq experiments in learning their underlying gene expression model. Indeed, DE methods are often run on pairs of cell groups, subsampled from the original data. Gene expression level estimates could be improved by leveraging shared information across the whole dataset for model fitting, for example, to refine given gene expression levels from gene–gene correlation structures across batches. Recent advances in deep generative modeling (11) have led to the development of a range of methods for normalization and analysis of gene expression data at scale, able to deal with millions of cells and able to account for complex study designs (12). Most of these models [e.g., scSphere (13), scVI (14), scANVI (15), scGen (16) and CellBender (17)] represent the expected gene expression in cell n as a nonlinear function of a latent low-dimensional summary of the cell’s state zn and its batch identifier sn. Such approaches intrinsically lead to information sharing across cells and across genes and promise to address the aforementioned issues. However, little attention has been given to carrying out rigorous and accurate DE analyses with those models. Two exceptions are scVI (14) and scANVI (15) in which the uncertainty in the latent variable zn is used to approximate a Bayes factor and detect DE. However, the Bayesian approach studied in refs. 14 and 15 has three fundamental limitations. First, the underlying definition of DE genes is too sensitive and may tag genes with negligible differences in expression and often of limited biological relevance. Second, interpreting Bayes factors is challenging because of the difficulty in translating them to established multiple hypothesis testing metrics, for example, false discovery rate (FDR). Third, scVI proposes to average Bayes factors from randomly sampled pairs of cells from each group, which lacks theoretical justification. After presenting background information on DE and deep generative models (Section 1), we introduce lvm-DE, a general Bayesian framework for detecting DE that addresses the aforementioned issues (Section 2). lvm-DE takes as input a fitted deep generative model of scRNA-seq data, a pair of cell groups, and a target α and provides as output an estimate of the log fold change (LFC) for every gene as well as a list of DE genes, with FDR control at level α. Our Bayesian formulation of DE hypothesis testing uses a composite alternative, built from the LFC [as in DESeq2 (5)] to avoid detecting spurious or weakly expressed genes. We benchmark lvm-DE applied to scVI (14) and scSphere (13) against classical DE approaches on simulated and real-world datasets (Section 3). Our results demonstrate that lvm-DE is a valuable tool for understanding DE effects in scRNA. An open-source implementation of lvm-DE is available in beta as part of the scvi-tools Python package (18) (https://github.com/scverse/scvi-tools).

Before introducing our statistical framework, we describe two lines of work in scRNA-seq data analysis on which our work builds: DE and deep generative modeling.

DE analysis of transcriptomics data aims to determine which genes have significantly different expression levels between two groups of cells A and B (e.g., representing different cell types or cells obtained under different experimental conditions). It is known that the direct approach of comparing counts across groups is subject to noise and bias because of the various sources of technical factors that confound the data, including variation in library size and batch effects (19). To address this, all DE methods define (either implicitly or explicitly) a model with a latent variable hgc that represents the denoised or normalized expression level for each gene g and each condition c ∈ {A, B}. DE for gene g between condition A and B is then assessed by comparing some statistics measuring the differences between the distributions of hgA and hgB. Finally, the significance of these differences is corrected to account for the multiplicity of hypotheses (one per gene) by estimating the FDR. Several popular DE frameworks for scRNA-seq data were originally developed for bulk RNA-seq data. Those methods are still competitive although not explicitly designed for scRNA-seq (20). Two of the most prevalent methods, edgeR (6) and DESeq2 (5), model gene expressions as negative binomial distributions. Another widely used method, Limma-voom (21), uses log-transform counts instead. A common theme for all three methods is that they rely on GLM, making it convenient to represent and control for possible biases (e.g., due to batch effects or GC content). Although these GLMs are fit for every gene separately, it is important to note that their training procedures pool information across genes, relying on mean–variance relationships and empirical Bayes methods to improve the fit of some of the parameters. Another key contribution, featured in DESeq2, is to propose a composite null hypothesis in which the hypotheses are formulated in terms of the fold change being higher than some threshold. This results in exclusion of genes whose LFC is low, even if the differences in their expression are statistically significant. In practice, this procedure helps exclude genes with little biological interest. Another set of DE analysis frameworks have been designed specifically for scRNA-seq data, with popular methods including MAST (2), SCDE (3), and scDDboost (22). MAST uses a hurdle model that was designed to account for the abundance of so-called dropouts in scRNA-seq data; these are zero counts that are observed mostly with weakly expressed genes, due to limitations in RNA capture efficiencies. SCDE takes a conceptually similar approach and employs a Bayesian framework to model scRNA-seq data as a mixture of negative binomial and Poisson distributions. The negative binomial component of the mixture accounts for counts that are attributable to gene expression, while the Poisson component models dropout events (23). For within-cell-type DE analyses across conditions, scDDboost employs an empirical Bayes approach to detect DE from compositional changes within a cell type from subtype assignments obtained from clustering. In particular, it posits an empirical Bayes prior for subtype counts and captures differences in gene expression within subtypes from a negative binomial model. In the following, we compare our approach for DGM-based DE against this collection of DE methods, coming from both bulk (DESeq2, edgeR, and Limma-voom) and the scRNA-seq literature (MAST).

Given the steadily growing sizes of datasets (24), recent work has begun to adopt more flexible and nonlinear approaches for modeling gene expression in single cells. While these methods require more data for training, they are often complemented by scalable inference procedures. One popular line of work relies on advances in deep generative modeling (12), in which models make use of deep neural networks as link functions in their conditional distributions (e.g., refs. 13–17, 25–30). Most of these methods fit a hierarchical probabilistic model (i.e., a graphical model) that includes a low-dimensional latent variable embedding each cell in some cell state space. To ensure that these cell embeddings reflect biological variation as much as possible, and remain disentangled from technical variation, a common practice is to model nuisance factors as parameters or random variables within the generative model. One factor that is commonly accounted for is variation between cells due to sequencing depth. For instance, scVI (14) and scANVI (15) models include a random variable that represents the sequencing depth of each cell. This variable is used as a scaling factor, thus providing normalized expression levels. The hierarchical model also links the latent representation of the cell state directly to the normalized expression values, thus decoupling the cell embeddings from variations in sequencing depth. An alternative approach, scPhere (13), uses the sum of observed transcripts as a deterministic scaling factor but relies on an additional penalty term to ensure that the cell embeddings are decoupled from it. Deep generative models provide accommodating strategies to correct for one or more other nuisance factors present in single cells. Multiple approaches, including scVI, scANVI, or scPhere, model normalized gene expression as a function of both cell state embeddings and these nuisance factors. This function learned during the fitting process and parameterized as a multilayered perceptron can model intricate nonlinear nuisance effects on gene expression. Such nuisance factors can include batch or patient identity in multidonor datasets but have also been generalized to other covariates such as cell cycles (18). Distributional assumptions are other key components that can improve cell-state representation and gene expression modeling. For example, in scGen (16), the counts are normalized, log-transformed, and treated as samples from a normal distribution. Conversely, both scVI (14) and scPhere (13) model gene expression counts under a negative binomial likelihood. Other models (e.g., the original work of scVI) include the addition of an additional low-magnitude term, which aims to handle zero inflation in the data, due to limitations in sensitivity and the effects of transcription bursting (discussed in more depth in refs. 31 and 32). Another component of importance is the modeling of the cell embedding. DCA (27) relies on a deterministic embedding, treated as a parameter optimized using maximum likelihood estimation. Another line of work, variational autoencoders (VAEs), posits a prior distribution on the embedding that regularizes the embedding space and aims to learn more useful representations (11). While most unsupervised VAEs posit a noninformative isotropic Gaussian prior, alternative priors may improve cell representation in specific contexts. For instance, scANVI (15) leverages partial data annotations to learn more structured priors on the latent variable. scPhere embeds biological states in a hyperspherical space with a uniform prior, which could improve the modeling of rare cell types and hierarchical structures. While most applications of deep generative models in single cells revolve around the cell state space, few approaches consider normalized gene expression values and their uncertainty, that are central for DE. Aiming in this direction, DCA feeds denoised gene expression predicted by the model as inputs for DESeq2 instead of raw counts to better reduce the noise present in scRNA-seq datasets. This is unsuitable from a computational perspective because DESeq2 may be too long to run when dealing with groups of tens of thousands of cells, creating a bottleneck. Also, as DCA relies on an autoencoder, it does not provide uncertainty estimates on normalized gene expression and summarizes the data via its mean instead. This is known to potentially create a spurious signal in other downstream applications [e.g., gene–gene correlation estimation (33)]. scVI and scANVI (14, 15) both provide integrated tools to perform DE between two groups. First, expected expressions are sampled from the posteriors of representative cells. These values are then subject to Bayes factor computation to assess significance. This framework has several limitations. First, although Bayes factors are an appealing alternative to frequentist hypothesis testing for Bayesian models, they may be hard to interpret while analyzing real-world data and are difficult to combine with frequentist FDR control guarantees. Second, these approaches annotate DE genes as significant without considering the effect size (i.e., LFC), which can be extremely low for lowly expressed genes, which are relatively abundant in scRNA-seq data. Third, this approach is convenient since model fitting is performed only once using all cells, thus avoiding the need to refit a new model for each DE test (each considering different pairs of cell subsets). However, scVI does not explicitly account for cell type in its generative model. As a heuristic, scVI aggregates cells by averaging approximate Bayes factors across randomly sampled pairs of cells to compare cell populations. We provide more theoretical insight and an improved approach in this manuscript.

We introduce lvm-DE, a principled framework to enable the use of deep generative models as a means for quantifying and assessing the significance of differences in gene expression (Fig. 1A). lvm-DE takes as input a latent variable model, a scRNA-seq dataset with covariates (such as batch identifiers) and two groups of cells to be compared, as well as a significance level α. lvm-DE provides an estimate of the LFC of each gene between cells of the two groups. It also provides a list of significantly differentially expressed genes that control the FDR at level α. These estimates are derived while controlling for differences in library sizes and for variability due to observed technical confounders, such as the association of samples with different batch identifiers. lvm-DE comprises five steps (Fig. 1B). First, a model is fitted such that the observed expression values in each cell (X) are represented by a low-dimensional latent random variable (often denoted as z; in VAE-based models such as scVI, this amounts to learning an encoder network). This representation provides a summary of the state of each cell, while controlling for differences due to batch effects(s) and sequencing depth, but it does not use the annotation of cells into the two groups that we wish to compare (Fig. 1B, 1; Section A). Performing DE with an annotation-agnostic model is convenient as it does not require refitting the model when we want to compare different groups of cells. Second, annotations are used to build a posterior distribution of the latent representation of each of the two groups, providing a summary of their characteristic cell states (Fig. 1B, 2; Section B). Third, cell state posteriors are used to infer a distribution of normalized expression for every gene (Fig. 1B, 3) in each cell group. These high-dimensional distributions are generated by propagating the low-dimensional representations through the latent variable models (LVM). In VAE-based models, this amounts to applying the decoder network. Fourth, the normalized expression values are used to estimate the LFC and a posterior probability for the LFC to be larger than a given cutoff δ (Fig. 1B, 4; Section C) while correcting for batch effects (Section D). Finally, these posterior probabilities of DE are used to call genes whose significance is within a desired FDR α (Fig. 1B, 5; Section E).

Our formulation of lvm-DE is general, such that it applies broadly to many types of deep generative models, and can operate as a wrapper function for DE with these models. In the following, we provide the modeling prerequisites for using lvm-DE and then discuss additional considerations of modeling choices that can lead to better accuracy. Let xng be the observed number of transcripts in cell n that originate from each gene g, and let sn correspond to the batch annotation. A general form of generative process that qualifies for use by lvm-DE is as follows. First, a low-dimensional latent variable zn representing a cell’s biological state is sampled from a prior distribution. A neural network then takes a sample from zn and the batch information sn and maps it to another latent variable, hng, representing the underlying (i.e., normalized) expression level for gene g. Finally, the observations xng are drawn from model-specific conditional likelihood distributions, such as negative binomial or log-normal, with or without zero inflation (14, 16, 34). The likelihood of a single observation (cell) in these models decomposes as follows: Furthermore, lvm-DE assumes that there exists a tractable approximation to the posterior p(zn ∣ xn) (11, 35). While many published methods satisfy this basic prerequisite, here, we selected two such models—scVI and scPhere (13, 14)—both of which are based on a VAE but differ in their loss function and distributional assumptions. These two models employ variational approximations to the posterior. In particular, they fit distributions q(zn ∣ xn), whose parameters are functions of the input data xn and which are learned during model fitting. Beyond this basic prerequisite, we identified two considerations for improved performance (implementation details in SI Appendix, section B). The first consideration is the need to model the normalized expression levels hng, while controlling for nuisance variation in sequencing depth. In scVI and scPhere, this is addressed by modeling the gene expression count xng as a negative binomial distribution with mean parameter lnhng, where ln is a normalizing factor that can be interpreted as the inferred library size. The second consideration stems from the fact that lvm-DE exploits the variational distribution (i.e., an approximation of p(zn ∣ xn)) for Bayesian hypothesis testing. In this setting, it has been reported that VAE trained with the classical evidence lower bound (ELBO) may underestimate the variance of the variational posterior and thus underperform in tasks of decision-making, compared to models trained with alternate variational bounds (36). Consequently, we train scVI and scPhere using the importance-weighted evidence lower bound (IWELBO), which has been demonstrated to lead to improved performance in other application domains (37).

lvm-DE takes as input stratification of cells into groups (cn ∈ {1, …, K}, representing, for example, different cell types or different conditions under which cells were captured), along with a generative model that conforms with Eq. 1. To estimate the distribution of gene expression in every group of cells, we assume that there exists an oracle function Φ mapping a point in latent space z to the corresponding covariate c = Φ(z)∈{1, …, K}. This is a reasonable assumption because the low-dimensional representation given by scVI or other deep generative models often successfully stratifies cells according to their groups (14). Given a particular value C of the covariate c (e.g., a cell type or a condition), we denote by hgC the normalized expression level of gene g in condition C. We define hgC, conditioned on a given batch s, as where p(z ∣ Z ∈ Φ−1(C)) denotes the prior density p(z) truncated to the preimage of {C} under the oracle function Φ. For any latent code z and batch s, it is assumed possible to sample from p(hg ∣ s, z) using the generative model. Hence, our goal in the rest of this section is to define a strategy to sample from the left-hand side of the integral, p(z ∣ Z ∈ Φ−1(C)). Because Φ is inaccessible to us, we may estimate the prior density from the data, in the spirit of empirical Bayes methods. More precisely, we first approximate it using the posterior of the annotated data: Second, we estimate this quantity using the truncated dataset with cells belonging to condition C, yielding the mixture density where 𝒩C is the set of cell indices that belong to the condition C. To reduce the possible sensitivity to outlier cells, we developed a procedure to identify and remove these from the set 𝒩C. To do so, assuming that the latent mean representations within a cluster are normally distributed, we estimate the covariance matrix from posterior mean samples {𝔼p(z ∣ xi)[z]∣i∈𝒩C} and remove observations whose variational mean falls outside the 90%-confidence ellipse described by the covariance estimate. As demonstrated previously, the empirical prior appearing in Eq. 4 is intractable. A first solution is to replace the individual posterior distributions p(z ∣ xn) with their variational approximations qϕ(z ∣ xn), yielding the so-called plugin estimator A remaining caveat of this estimator is that it lends equal importance to each sample from the cell population and of the variational distribution. This is in a way similar to the original scVI, in which pairs of cells are sampled uniformly at random from the two populations and compared to each other. It is reasonable to expect, however, that in practice, some samples may be of significantly better quality, which explains the previously reported increased performance from using importance-sampling-based estimators (36). Based on this observation, we use the mixture prior from Eq. 4, approximated by self-normalized importance sampling (SNIS) with the mixture of Eq. 5 as the proposal distribution. In this formulation, we first sample , and the contribution of each sample zi is proportional to: where p(zi) denotes the prior density for zi, and the intractable evidence terms p(xn) are estimated with IWELBO estimators (K = 5, 000 particles). The final samples correspond to random samples with replacement from {zi, i ≤ N}, where each entry i of the set has weight (37), such that Consequently, samples from z that have high likelihood levels across a large pool of cells from 𝒩C will have higher weights in the posterior, thus further improving the robustness to outliers and to misassignment of cells into groups. In order to ensure that the importance sampling procedure is sufficiently stable, several diagnosis tools such as the Pareto smoothed importance sampling (PSIS) estimator (38) may be employed. Once we have samples from the distribution p(z ∣ Z ∈ Φ−1(C)), we may follow Eq. 2 to link these samples to the normalized expression levels hgC for DE analysis.

To access p(hgC ∣ s), it hence suffices to sample latent codes z from the variational posteriors qϕ(z ∣ x) to reweight them via Eq. 6 and to decode them using p(hg ∣ s, z) using the generative model of choice. In the case of scVI and scPhere, expression levels h and cell states z are deterministically linked by a neural network h = f(z, s).

A practical limitation of DE analysis is the identification of genes in which the differences in expression are significant and yet of low magnitude (low LFC) and thus of less or no interest. To address this, lvm-DE treats the LFC itself as a random variable and uses a Bayesian analysis to retain only genes in which the effect size is significantly larger than a cutoff. To quantify our confidence that the LFC is beyond a desirable level in a given batch s, we introduce the effect size random variable, where hgA ∼ p(hgA ∣ s),hgB ∼ p(hgB ∣ s) are approximated using the SNIS procedure described above. Here, ϵ is a small offset used for stability. When gene g has very low expression levels in both conditions, the pseudocount ϵ ensures that the associated LFC has a low magnitude (5, 39). Details regarding the choice of ϵ can be found in Appendix D. To test for DE in lvm-DE, we introduce three models for which gene g is either up-regulated (ℳ+g), down-regulated (ℳ−g), or not differentially expressed (ℳ0g) in population A compared to B. These three models are defined via the change variable rA, Bg as follows: where δ is a threshold ensuring that the selected genes are interesting in practice. We discuss heuristics for the choice of δ in Appendix E. The probability of each of those models may be estimated from the data, using the empirical Bayes method and the SNIS estimator. In particular, we note: the probability that gene g is DE (i.e., with absolute LFC greater than δ). The LFC is a widely used measure of effect size for expression measurements because of its interpretability (5, 40, 41). However, our framework could be extended to any difference measure between expression levels, including logit fold change for DNA methylation, as in scMET (42).

Our exposition in the previous section assumed that all the cells come from a single batch s. However, DE analyses are often required in multibatch or multidonor settings. Two scenarios, reflecting different perspectives on what batches refer to, appear when performing DE, depending on whether all batches contain cells from both conditions.

In the first scenario, the two compared subpopulations A and B have cells in overlapping batches. In the most natural setting, such as biological replicates, we expect a representation of all the compared conditions in all batches. This often is the case when batches refer to biological replicates or to different technologies; we then hope to detect DE genes reproducible across all the batches. In this case, let S be any set of batch indices for which there exists at least one observation of each condition. Following our empirical Bayes approach, we define batch-specific empirical priors relying on the subset of observations 𝒩Cs = {n∣sn=s,cn=C}. Then, the probability of each of the three models may be computed by averaging across batches (similar calculations were proposed in ref. 15). For example, the probability of upregulation is computed as A possible extension to continuous nuisance covariates is discussed in SI Appendix, section A.

In some setups, however, there may not be perfect overlap of each condition in all batches. (e.g., when comparing tumors or inflamed tissues to healthy controls) Particularly, when comparing tumors, or inflamed tissues to healthy controls, batches may denote cells from individual donors, and one aim is to detect differentially expressed genes for a cluster of cells between subsets of batches. In such a case, we redefine the LFC in Eq. 8 to make it independent of the batch assignments by sampling normalized gene expression in each cell using observed batches. In this more ambiguous scenario, we condition the posterior scales distributions with respect to observed batches; that is, p(hgC ∣ sn = snobs), snobs being the batch that cell n originates from.

Bayesian decision theory (43) shows that the optimal decision rule for our DE task consists of selecting the top genes ranked by pg (Eq. 10). However, selecting the number of “top” genes is not straightforward. The use of the Bayes factor to guide this choice carries the complication of not being interpretable when it comes to the expected error rates. A more intuitive frequentist solution is to set the significance threshold α over pg in a way that controls for the FDR. To this end, let dg be the (unknown) binary random variable denoting whether or not gene g is DE . Consider the multiple binary decision rule μk = (μgk, g ∈ G) that consists in calling DE the k genes with the highest DE probability pg. The false discovery proportion (FDP) of the decision rule μk can be expressed as We use the formulation of the posterior expected FDP, as introduced in ref. 44 and applied to genomics data in refs. 45 and 46, corresponding to Because this quantity is tractable after inference, a natural choice for k corresponds to the maximum value for which the posterior expected FDP is below the desired level α.

lvm-DE addresses two key tasks in DE analysis. The first corresponds to the robust estimation of the LFC in noisy and sparse data, including for rarely observed cell states. The second task is evaluating the statistical significance of the observed effect size and calling genes as differentially expressed for any target FDR. For both tasks, we apply the lvm-DE framework to scVI and scPhere. We denote scVI-lvm and scPhere-lvm the associated combinations, in contrast to scVI that natively provides a test for DE. We evaluate our models in comparison to the current state-of-the-art, including MAST (2), DESeq2 (5), edgeR (6), and Limma-voom (21) (referred to as voom). For each of these models, we follow the guidelines and implementation details specified in a recent comparative study of DE methods for scRNA-seq (20) (further details in Appendix F). We assess lvm-DE performance for LFC estimation and for the prediction of correct and calibrated DE gene sets using both synthetic and real datasets. In the real datasets, we focused the analysis on the top 3,000 highly variable genes using Seurat. Annotations were taken from the respective publications.

As a simulation scheme, we employ SymSim (47), which models biological and technical noise factors (e.g., transcriptional bursting and PCR amplification), as well as statistical dependencies between genes (relying on low-dimensional representation of the kinetics of promoters in each cell). We used SymSim to generate a hierarchy of cell states and focus our analysis on the comparison of two “sibling” types, denoted A and B (Fig. 2A). We generated a dataset of 10,000 cells and 1,000 genes across five cell types. The data correspond to two batches (5,000 cells each) with unique molecular identifiers (UMIs), which have similar properties to 10× Chromium data. To estimate the ground-truth LFC, we use the expected gene expression means, derived from the expected promoter-kinetic parameters in each gene and cell type (47). Ground-truth gene expression varies across clusters for all genes in SymSim; that is, there are no true negative genes. We restricted the analysis to genes that either had large (≥1 in absolute value) or very small (≤0.2 in absolute value) LFCs between the two compared clusters to define “true” DE genes without ambiguity. The resulting dataset contains a total number of 396 genes, including 127 differentially expressed genes. We selected these genes such that DE genes have nonambiguous high expression differences (absolute LFC higher than 1), while negatives have negligible expression differences (absolute LFC lower than 0.2). We first evaluated LFC estimates from the models. The lvm-DE-based models offered accurate estimates (Fig. 2B), on par with edgeR and DESeq2. This result suggests that lvm-DE properly aggregates individual cells to infer expression levels in the two subpopulations. An outlier in this analysis was MAST that exhibited a high mean-squared error range. This likely comes from the fact that this method relies on a hurdle model, whose LFC is difficult to compare to the ground-truth LFC. It was therefore discarded from further LFC estimation experiments. By definition, scVI-lvm and scPhere-lvm make use of all cells in the input data, including “external” cells that are not in the groups being compared. In our simulation, these are simulated cells that are not in group A or B (Fig. 2A). To test whether the inclusion of more external cells affects accuracy, we varied the number of external cells and conducted DE analysis between groups A and B (which remain unchanged). We find that our FDR control procedure is well calibrated and remains at the expected level in all cases. Furthermore, we find a steady increase in the true positive rate (Fig. 2C). These results indicate that an increase in the overall dataset sizes, but not necessarily in the compared populations, leads to an increase in power, thus supporting the notion that the DGMs can benefit from information on all cell types (e.g., gene–gene correlations) in a dataset. This is in contrast with the other benchmark methods listed here which do not absorb any information from cells outside the compared groups. The results in Fig. 2D support the suitability of lvm-DE for characterizing rare cell types, as capturing gene couplings can help alleviate the issue of small sample size. We tested this by leaving in the input data only a few cells of group B (representing a rare cell type), while keeping group A as before (representing a more prevalent background; see Fig. 2A). In particular, lvm-DE was retrained each time the rare subpopulation size was modified, such that all algorithms were trained on the same input data. For a rare subpopulation of size 50 and up, both lvm-DE methods provide gene lists that are accurate and control the FDR for several significance levels. In this experiment, we considered first a simple baseline, relying on a two-side t-test on log-counts-per-millions (CPM) counts, which fails at controlling the FDR. For the GLMs, we observed acceptable FDR control in small sample size setups, but interestingly, their FDR becomes higher than its expectation, with larger sizes of the rare subpopulation. An exception is the default DESeq2 procedure, which uses a composite null hypothesis (i.e., effect size threshold δ > 0). This method does not suffer from large FDR but loses power because of the difficulty in properly setting the effect size threshold in practice. Finally, we observe that the Bayes factor approach originally used in scVI yields high FDR levels. Consequently, we excluded this procedure from the next experiments. To further our confidence in these results, we expanded the analysis in two ways. In the first analysis, we modified our simulation process to generate additional and possibly more nuanced scenarios of DE, such as bimodality of the ground truth gene expression (48). Our procedure compares favorably to the other algorithms as it provides a high number of true positive genes and few false positives (SI Appendix, Fig. S2). In a second analysis, we used MUSCAT (49) as our simulation framework. Unlike SymSim, the data generated by MUSCAT are based on independent sampling from an underlying (here, negative binomial) distribution. Therefore, in this case, there is no benefit in using latent variable models, as no couplings between genes (which is to be expected in real datasets) can be leveraged. In such setups, the lvm-DE methods provided calibrated predictions, while having a similar or lower level of power compared to the GLMs, as expected (SI Appendix, Fig. S3).

To test the performance of lvm-DE with data, we used a 10X Chromium dataset of 9,432 peripheral blood mononuclear cells (PBMC) from a human donor (50) (Fig. 3A). These data provide an opportunity to establish the relative performance of lvm-DE in a configuration in which batch effects are not present. While there is no ground truth for DE in these data, several strategies can help evaluate the estimation of DE. First, we can ensure that lvm-DE does not detect differences in expression when comparing groups of cells of the same state. We therefore construct negative controls (Fig. 3B) by randomly splitting the B cell population (460 cells) into two groups (230 cells each) and predicting the LFCs with the different algorithms. In this configuration, we expect to measure negligible expression differences. lvm-DE compares favorably to its competitors for LFC estimation, since scVI-lvm along with scPhere-lvm predicts LFC of significantly lower amplitude than their GLMs’ counterparts. Another benchmark strategy is to compare the outcome of the DE analysis to results obtained from an independent dataset. To this end, we used estimations of LFCs and P-values from a bulk-RNA-seq dataset (51), in which DESeq2 was used for pairwise comparisons of monocytes, B cells, plasmacytoid, and myeloid dendritic cells. The bulk-level LFCs and significance scores were then compared against the results of applying the different algorithms on the single cell data (Fig. 3 C and D). We find that scVI-lvm and scPhere-lvm compare favorably to the GLMs in recapitulating both the LFC and the significance scores of the reference bulk data. Notably, the cell groups we analyzed had varying sizes (Fig. 3A), demonstrating that lvm-DE can capture state differences for abundant as well as less common cell types. A detailed comparison of LFC estimations for every method and for every pair of cell types is depicted in SI Appendix, Fig. S4.

lvm-DE can be employed to perform DE in large-scale multibatched data. As an illustration, we consider the PbmcBench data collection (52), consisting of two biological samples sequenced with seven different protocols, either droplet (drop-seq, 10Xv2, 10Xv3, inDrop, and CEL-seq) or well based (seq-well and Smart-Seq2). These data provide the opportunity to further validate the ability of our procedure to alleviate technical factors’ impact on DE. We compare the estimated LFC of scVI-lvm with bulk-RNA using the same reference as in the previous experiment. The inclusion of additional data sequenced with different technologies consistently increases the correlation between lvm-DE effect-size and significance predictions with bulk-RNA data (Fig. 4). In contrast, observing more replicates does not systematically improve the match of GLM’s predictions with bulk-RNA. This suggests that contrary to GLMs, lvm-DE succeeds in pooling information from heterogeneous experiments to better capture differences between cell states. In summary, lvm-DE benefits when either the number of observed cells increases or when more samples are available, even if the samples differ technically.

We now focus on a larger PBMC dataset, consisting of 44,721 cells from seven healthy donors and six patients with confirmed SARS-CoV2 infection (53) (Blish). The latent variable models were trained using donor ids as batches. This dataset provides an opportunity to evaluate procedures for DE, not only for comparing different cell types (to find their respective markers) but also comparing donor groups within each cell type. The latter comparison, revealing the influence of inflammation on different immune compartments, is usually more challenging with a more nuanced and less reproducible signal (54) (Fig. 5A). In this scenario, a popular approach consists in pseudobulking the data—that is, aggregating counts by sample prior to testing. We include this additional baseline for this experiment, relying on a two-sided t-test on log-CPM pseudobulk. Starting with between-cluster analysis, we divided the cells into types using annotations from ref. 53. We started by applying the algorithms on a negative control setting, obtained by evenly splitting dendritic cells as two clusters (203 cells each) and comparing the two clusters (Fig. 5B). As in the previous experiment, the LFC values estimated by scVI-lvm and scPhere-lvm in this negative control are significantly smaller than the rest of the algorithms. We next proceeded to a comparison between different cell types (irrespective of disease status). Consistent with our results above, we find that the DE analysis scVI-lvm and scPhere-lvm more closely match the bulk reference than the GLM-based methods, in terms of both effect sizes (LFC) and significance (ordering by FDR or P-values) estimates (Fig. 5C and SI Appendix, Fig. S5 for the individual plots). We also note that the lvm-based methods scale better in terms of computation time for this large dataset (Appendix K). We next moved to the more challenging task of DE between the patient and control population within each cell type. In this case, the lvm-DE approach uses the modified LFC to compare these populations, given that case and control subpopulations correspond to disjoint batches. SARS-CoV-2 infection has been demonstrated to lead to aberrant secretion of proinflammatory cytokines such as interferon gamma (55). As a reference to evaluate our DE analyses, we, therefore, used bulk transcriptome measurement of responses to different cytokines by different immune cell types. We first used a Microarray dataset (56), which provides reference LFC for cell-type-specific effects of interferon alpha. We find that scVI-lvm and scPhere-lvm better correlate with the reference microarray (significant under a t-test at significance level α = 0.05; Fig. 5D and SI Appendix, Fig. S6 for the individual plots). We repeated this analysis with another reference dataset of bulk RNA-sequencing measurement of responses to other cytokines (IL2, IL4, IL7, IL9, IL15, and IL21) and again found that the LFC and significance evaluations by lvm-DE had consistently high correlations with the reference (SI Appendix, Fig. S7). In summary, we find that also in this more challenging task of within-cell-type comparisons, lvm-DE provides more accurate estimations of gene expression compared to our benchmark methods.

While much effort has been made recently to design deep generative models for scRNA-seq and more generally multiomics data, there has been a paucity of research on leveraging these models for DE. In this work, we present lvm-DE, an empirical Bayes method for DE, taking full advantage of the flexibility and scalability of latent variable models. This flexibility permits nonlinear batch effects to be accommodated in the model, improving on global linear models for batch effect correction that assume that each batch has the same cell-type proportions, often an invalid assumption in scRNA-seq experiments (8). Furthermore, most DE pipelines, with few exceptions (57), do not take advantage of all available samples for between-type comparisons. Latent variable models that lvm-DE rely on can harness the learned shared gene-gene correlations to refine gene expression modeling and, hence, increase the number of detected DE genes. As such, lvm-DE helps characterize gene expression differences between cell populations in complex large-scale, batch-confounded datasets and provides calibrated DE gene sets relying on posterior expectations of the FDP. While we demonstrated broad applicability of lvm-DE, we discuss now two caveats that should be considered while implementing or using lvm-DE with a deep generative model. First, it is possible that lvm-DE may require slight modification of existing latent variable models, in order to ensure reasonable performance. In particular, we observed empirically that the library size modeling was key, for example, in DE comparisons involving two populations with different RNA content. In scVI, the library size is a latent variable, while in scPhere, it is deterministically assigned to the number of transcripts. We noticed that the latter scheme does not always successfully normalize counts and outlined what specific elements were required to unlock DE for scPhere (SI Appendix, section C). Second, lvm-DE may not be suitable to analyze very small datasets, containing less than a few thousand cells, as this setup may not be suitable to properly fit deep generative models. In such cases, it is more natural to use traditional methods, that i) scale well and ii) are not prone to overfitting. As ground-truth cell type labels are not available in most scRNA experiments, the success of DE heavily relies on a careful annotation strategy. Practical realities, including limited data availability or mislabeling, may deteriorate cell-type prediction accuracy and desired resolution, and hence the range of DE analyses. The emergence of high-quality, expert-annotated cell atlases (58, 59), along with scalable transfer-learning procedures (60), promises to provide straightforward and steady ways to annotate cells. As DE pipelines often follow the canonical approach of using the same gene counts to annotate (i.e., formulate hypotheses) and compare cell populations (testing these hypotheses), a core challenge is that the detected DE genes may reflect spurious differences, initially captured during the annotation process (61). Fortunately, emerging approaches may limit the circularity that is central to DE for scRNA (62). We evaluated how such a procedure, molecular cross-validation, could reduce spuriousness using synthetic data, showing marginal improvement over the canonical approach (SI Appendix, Fig. S8). As lvm-DE outlines a generic procedure to conduct DE for latent variable models, improving the LVM of choice can be a direction to improve the quality of the predictions. One idea could be to leverage available information on gene localization to improve the underlying representation of gene–gene correlations. For instance, hmmSeq (46) is a hierarchical Bayesian model for RNA-seq taking gene positioning into account to better capture coexpression of genes in the same chromosome. Designing models that handle the sparsity of scRNA-seq experiments is another way to improve DE for latent variable models. From a biological perspective, understanding the nature of zeros and the necessity to model their excess is an active topic of research (2, 31, 32). This question depends on the technology used for sequencing and has already been illustrated through the analysis of zero inflation for scRNA in latent variable models (15). Besides, our aim is to formulate a general framework for DE for such models, that is agnostic to such modeling choices. Still, we found that lvm-DE manages to provide meaningful and calibrated decisions for zero-inflated models (SI Appendix, section J and Fig. S9). Designing better optimization procedures for latent variable models in sparse, high-dimensional setups is another objective for scRNA-seq data. Sparsity can indeed hinder lvms’ ability to handle a larger number of genes. In scRNA-seq data, this limitation could be overcome in several ways for models resorting to amortized variational inference. The use of more refined training procedures, relying on tighter bounds to improve the fitness of the generative model (63) could unlock DE on larger datasets. The introduction of group-sparsity priors in the decoder could effectively leverage the fact that many genes share generalized biological functions to speed up inference (64). We hypothesize that scaling up latent variables models to larger datasets might not only help detect more DE genes but could also increase the resolution of cell population studies. While providing interpretability of decisions in terms of FDR was the main focus of our work, FDR control may not always be the most appropriate way to calibrate decisions for DE analyses. Predicting the direction of the LFC may be more important than detecting the presence of a DE in some applications, in which case False Sign Rate control (65) provides an alternative that should be considered. Another extension of our work is the formulation of differential analyses for latent variable models that integrate several omic modalities. scMET (42) is a Bayesian model for differential methylation. Leveraging existing correlations (66) between methylation and gene transcription in a unified latent variable model and formulating meaningful hypotheses could be a step forward to better understand complex differentiation scenarios, for example, for embryonic stem cell characterization (67).

Click here for additional data file.