Abstract

As single-cell transcriptomics becomes a mainstream technology, the natural next step is to integrate the accumulating data in order to achieve a common ontology of cell types and states. However, owing to various nuisance factors of variation, it is not straightforward how to compare gene expression levels across data sets and how to automatically assign cell type labels in a new data set based on existing annotations. In this manuscript, we demonstrate that our previously developed method, scVI, provides an effective and fully probabilistic approach for joint representation and analysis of cohorts of single-cell RNA-seq data sets, while accounting for uncertainty caused by biological and measurement noise. We also introduce single-cell ANnotation using Variational Inference (scANVI), a semi-supervised variant of scVI designed to leverage any available cell state annotations — for instance when only one data set in a cohort is annotated, or when only a few cells in a single data set can be labeled using marker genes. We demonstrate that scVI and scANVI compare favorably to the existing methods for data integration and cell state annotation in terms of accuracy, scalability, and adaptability to challenging settings such as a hierarchical structure of cell state labels. We further show that different from existing methods, scVI and scANVI represent the integrated datasets with a single generative model that can be directly used for any probabilistic decision making task, using differential expression as our case study. scVI and scANVI are available as open source software and can be readily used to facilitate cell state annotation and help ensure consistency and reproducibility across studies.

Introduction

Recent technological improvements in microfluidics and low volume sample handling [1] have enabled the emergence of single-cell transcriptomics [2,3] as a popular tool for analyzing biological systems [4,5,6]. This growing popularity along with a continued increase in the scale of the re-spective assays [7] has resulted in massive amounts of publicly available data and motivated large scale community efforts such as the Human Cell Atlas [8], Tabula Muris [9] and the BRAIN Initiative Cell Census Network [10]. The next natural step in the evolution of this field is therefore to integrate many available datasets from related tissues or disease models in order to increase statistical robustness [11], achieve consistency and reproducibility among studies [12,13], and ultimately converge to a common ontology of cell states and types [8,14].

A fundamental step toward the ideal of a common ontology is dataharmonization, namely integration of two or more transcriptomics datasets into a single dataset on which any downstream analysis can be applied. We use the term harmonization rather thanbatch effect correction in order to emphasize that the input datasets may come from very different sources (e.g., technology, laboratory), and from samples with a different composition of cell types. A wide range of methods have already been developed for this fundamental problem, initially for Microarrays and later on for bulk RNA sequencing, such as ComBat [15] and limma [16] which rely on generalized linear models with empirical Bayes shrinkage to avoid over-correction. More recently, similar methods have been proposed specifically for single-cell RNA sequencing (scRNA-seq), such as ZINB-WaVE [17], which explicitly accounts for the overabundance of zero entries in the data. However, because of their linear assumptions, these approaches may not be appropriate when provided with a heterogeneous sample that includes different cell states, each of which may be associated with a different sample-to-sample bias [12]. With these limitations in mind, the next generation of methods turned to non-linear strategies. Broadly speaking, each of these methods includes a combination of two components: (i) joint factorization of the input matrices (each corresponding to a different dataset) to learn a joint low-dimensional latent representation. This is usually done with well established numerical methods, such as integrative non-negative matrix factorization (LIGER [18]), singular value decomposition (Scanorama [19]), or canonical correlation analysis (Seurat Alignment [13]); (ii) additional non-linear transformation of the resulting latent representations so as to optimally “align” them onto each other. This is usually done using heuristics, such as alignment of mutual nearest neighbors (MNN [12], Scanorama [19] and Seurat Anchors [20]), dynamic time warping (Seurat Alignment [13]) or quantile normalization (LIGER [18]). While this family of methods has been shown to effectively overlay different datasets, it suffers from two important limitations. First, an explicit alignment procedure may be difficult to tune in a principled manner and consequently result in over-normalization, especially in challenging cases where the cell type composition is different between datasets and when technical differences between samples are confounded with biological differences of interest. Second, the alignment is done in an ad hoc manner and lacks probabilistic interpretability. Consequently, the resulting harmonized dataset is of limited use and cannot be directly applied for differential expression and other probabilistic decision-making tasks.

Another recent line of work makes use of neural networks to learn a joint representation of multiple datasets (SAUCIE [21]) or project one dataset into another (maximum mean discrepancy [MMD] ResNet [22]). These methods rely on an explicit non-parametric measure of discrepancy between probability distributions (MMD) to match either the latent spaces or directly the gene expression values from pairs of datasets. However, using the MMD with a universal kernel explicitly assumes that the cell type proportion is similar in all the datasets, which may be less suitable in the general case of data harmonization.

Besides harmonization, another important and highly related problem is that of automatedannotation of cell state. In principle, there are two ways to approach this problem. The first isab initio labeling of cells based on marker genes or gene signatures [13,23,24]. While this approach is intuitive and straightforward, its performance may be affected in the plausible case where marker genes are absent due to limitations in sensitivity. The second approach is to “transfer” annotations between datasets. In the simplest scenario, we have access to one dataset where states have been annotated eitherab initio, or using additional experimental measurements (e.g., protein expression [3,25] or lineage tracing [26]) and another, unannotated dataset from a similar condition or tissue. The goal is to use the labeled data to derive similar annotations for the second dataset, whenever applicable. This task is often complicated by factors such as differences in technology (e.g., using Smart-Seq2 data to annotate 10x Chromium data), partial overlap in cell type composition (i.e., not all labels should be transferred and not all unannotated cells should be assigned a label), complex organization of the labels (e.g., hierarchy of cell types and sub-types [27], continuum along phenotypic or temporal gradients), partial labeling (i.e., only a subset of cells from the “annotated” dataset can be assigned a label confidently), and the need to handle multiple (more than 2) datasets in a principled and scalable manner. One way to address the annotation problem with this approach is learning a classifier [27,28] in order to predict a fixed stratification of cells. Another, more flexible approach is to transfer annotations by first harmonizing the annotated and unannotated datasets, thus also gaining from the benefits of having a single dataset that can be subject to additional, joint, downstream analysis.

In this paper, we propose a strategy to address several of the outstanding hurdles in both of the harmonization and annotation problems. We first demonstrate that single-cell Variational Inference (scVI) [29] a deep generative model we previously developed for probabilistic representation of scRNA-seq data — performs well in both harmonization and harmonization-based annotation, going beyond its previously demonstrated capacity to correct batch effects. We then introduce single-cell ANnotation using Variational Inference (scANVI), a new method that extends scVI and provides a principled way to address the annotation problem probabilistically while leveraging any available label information. scANVI uses a semi-supervised generative model, which can be utilized for both approaches to the annotation problem. In the first scenario, we are concerned with a single dataset in which only a subset of cells can be confidently labeled (e.g., based on expression of marker genes) and annotations should then be transferred to other cells, when applicable. In the second scenario, annotated datasets are harmonized with unannotated datasets and then used to assign labels to the unannotated cells.

The inference procedure for both of the scVI and scANVI models relies on neural networks, stochastic optimization and variational inference [30,31] and scales to large numbers of cells and datasets. Furthermore, both methods provide a complete probabilistic representation of the data, which non-linearly controls not only for sample-to-sample bias but also for other technical factors of variation such as over-dispersion, library size discrepancies and zero-inflation. As such, each method provides a single probabilistic model that underlies the harmonized gene expression values (and the cell annotations, for scANVI), and can be used for any type of downstream hypotheses testing. We demonstrate the latter point through a differential expression analysis on harmonized data. Furthermore, through a comprehensive analysis of performance in various aspects of the harmonization and annotation problems and in various scenarios, we demonstrate that scVI and scANVI compare favorably to current state-of-the-art methods.

Results

In the following we demonstrate that our framework compares favorably to the existing state of the art in the harmonization and annotation problems in terms of accuracy, scalability, and adaptability to various settings. The first part of the paper focuses on the harmonization problem and covers a range of scenarios, including harmonization of datasets with varying levels of biological overlap, handling cases where the data is governed by a continuous (e.g., pseudotime) rather than discrete (cell types) form of variation, and processing multiple (> 20) datasets. While we demonstrate that scVI performs well in these scenarios, we also demonstrate that the latent space leaned by scANVI provides a proper harmonized representation of the input datasets — a property necessary for guaranteeing its performance in the annotation problem. In the second part of this manuscript we turn to the annotation problem and study its two main settings, namely transferring labels between datasets andab-inito labeling. In the first setting we consider the cases of datasets with a complete or partial biological overlap and use both experimentally- and computationally-derived labels to evaluate our performance. In the second settings, we demonstrate how scANVI can be used effectively to annotate a single dataset by propagating high confidence seed labels (based on marker genes) and by leveraging a hierarchical structure of cell state annotations. Finally, we demonstrate that the generative models inferred by scANVI and scVI can be directly applied for hypotheses testing, using differential expression as a case study.

Joint modeling of scRNA-seq datasets

We considerK different scRNA-seq datasets (Figure 1a). A dataset indexed byk ∈ {1, …,K} consists of aN_k × G_k matrix where each entry records the number of transcripts observed for each ofG_k genes in each ofN_k cells. In this work, we use a standard heuristic to filter the genes and generate a common (possibly large) gene set of sizeG (Online Methods), and thus format our data as a unique (Σ_kN_k)× G matrix whose individual entries are notedx_ng — namely the expression of geneg in celln. We uses_n ∈ {1, …,K} to represent the dataset from which each celln was generated. Furthermore, a subset of the cells may be associated with a cell state annotationc_n, which can describe either discrete cell types or hierarchical cell types. More complex structures over labels such as gradients are left as a future research direction.

• Download figure
• Open in new tab

Figure 1:

Harmonization of scRNA-seq datasets with generative models. (a) Schematic diagram of the variational inference procedure in both of the scVI and scANVI models. We show the order in which random variables in the generative model are sampled and how these variables can be used to derive biological insights. (b) The graphical models of scVI and scANVI. Vertices with black edges represent variables in both scVI and scANVI, and vertices with red edges are unique to scANVI. Shaded vertices represent observed random variables. Semi-shaded vertices represent variables that can be either observed or random. Empty vertices represent latent random variables. Edges signify conditional dependency. Rectangles (“plates”) represent independent replication. The complete model specification and definition of internal variables is provided in theOnline Methods

Since the problem of data harmonization of single-cell transcriptomics is difficult and can potentially lead to over-correction (Figure 1b) [32], we propose a fully-generative method as a robust and principled approach to address it. In our previous work [29], we built single-cell Variational Inference (scVI), a deep generative model where the expression levelx_ng is zero-inflated negative binomial (ZINB) when conditioned on the dataset identifier (s_n), and two additional latent random variables. The first, which we denote byl_n, is a one-dimensional Gaussian accounting for the variation in capture efficiency and sequencing depth. The second isz_n is a low dimensional Gaussian vector that represents the remaining variability (Figure 1c). This vector is expected to reflect biological differences between cells, and can be effectively used for visualization, clustering, pseudotime inference and other tasks. Since the scVI model explicitly conditions on the dataset identifier, it provides an effective way of controlling for technical sample-to-sample variability. However, scVI is unsupervised and does not make use of the available annotationsc_n, which can further guide the inference of an informative latent representationz_n. To this end, we present a more refined hierarchical structure forz_n. We drawz_n as a mixture conditioned on the cell annotationc_n and another latent variableu_n, accounting for further biological variability within a cell type (Figure 1d,Online Methods). We name the resulting approach single-cell ANnotation using Variational Inference (scANVI).

The variablesz_n, inferred either with scVI or scANVI, provide an embedding of all cells in a single, joint latent space. Since this latent space is inferred while controlling for the dataset of origin (s_n), it inherently provides a way to address the harmonization problem. The annotation of unlabeled cells can therefore be conducted with scVI using their proximity to annotated cells in the joint latent space (e.g., using majority vote over thek-nearest neighbors). The scANVI model provides a more principled way to annotate cells, namely through a Bayesian semi-supervised approach. Once fitted, the model is able to provide posterior estimates for the unobserved cell statec_n, which can be particularly useful when labels cannot be entirely trusted. Because the marginal distributionp(x_ng, c_n | s_n) ifc_n observed (resp.p(x_ng | s_n) otherwise) is not amenable to exact Bayesian computation, we use variational inference parameterized by neural networks to approximate it [30] (Online Methods).

Notably, scANVI and scVI both have a certain number of hyperparameters. In the following evaluations, conducted on different datasets and different scenarios, we use the exact same set of hyperparameters in order to demonstrate that our methods can be applied with a minimal requirement of hyperparameter tuning (Online Methods). We provide a robustness study for hyperparameters in the context of harmonization inSupplementary Figure 1. We further discuss the underlying assumptions of our framework in the context of competing harmonization methods inSupplementary Note 1.

Harmonizing pairs of datasets with a discrete population structure

We conducted a comparative study of harmonization algorithms on four different instances, each consisting of a pair of datasets. The first pair (PBMC-CITE [25], PBMC-8K [33]) represents the simplest case, in which the two datasets come from very similar biological settings (i.e., PBMCs) and are generated by the same technology (i.e., 10x) but in different labs (i.e., akin to batch correction). A second scenario is that of similar tissue but different technologies, which we expect to be more challenging as each technology comes with its own characteristics and biases [39]. For instance, some methods (10x, CEL-Seq2) profile the end of the transcript and use Unique Molecular Identifier (UMI) to mitigate inflation in counting, whereas others (e.g., most applications of Smart-Seq2) consider the full length of the transcript without controlling for this potential bias. Additionally, some protocols (e.g., Smart-Seq2) tend to have higher sensitivity and capture more genes per cell compared to others. Finally, studies using droplet based protocols tend to produce much larger numbers of cells compared to plate-based methods. We explore three such cases, including a bone marrow 10x and Smart-Seq2 pair from the Tabula Muris project (MarrowTM-10x, MarrowTM-ss2 [9]), a pancreas inDrop and CEL-Seq2 pair (Pancreas-InDrop, Pancreas-CEL-Seq2 [40]), and a dentate gyrus 10x and Fluidigm C1 pair (DentateGyrus-10x, DentateGyrus-C1 [41]).

Successful harmonization should satisfy two somewhat opposing criteria (Supplementary Figure 2). On the one hand, cells from the different datasets should be well mixed; namely, the set ofk-nearest neighbors (kNN) around any given cell (computed e.g., using euclidean distance in the harmonized latent space) should be balanced across the different datasets. For a fixed value ofk, this property can be evaluated using the entropy of batch mixing [12], which is akin to evaluating a simplek-nearest neighbors classifier for the batch identifier (Online Methods). While this property is important, it is not sufficient, since it can be achieved by simply randomizing the data. Therefore, in our evaluations we also consider the extent to which the harmonized data retains the original structure observed at each dataset in isolation. Here, we expect that the set ofk-nearest neighbors of any given cell in its original dataset should remain sufficiently close to that cell after harmonization. This property can be evaluated using a measure we callk-nearest neighbors purity (Online Methods). Clearly, this criterion can be simply optimized by not changing the input datasets, which will result in poor performance with respect to our first measure. Our evaluation therefore relies on both of these measures.

• Download figure
• Open in new tab

Figure 2:

Benchmarking of scRNA-seq harmonization algorithms (a)Visualization of the latent space of scVI using Uniform Manifold Approximation and Projection (UMAP [42]) for the bone marrow dataset pair [9], colored by cell type and batches. (b) Visualization of the Seurat Alignment latent space of the same dataset. (c) Quantitative comparison of the performance of scVI and Seurat Alignment by two measures. One the left retainment of structure for the benchmark algorithms on the four pairs of dataset evaluated via thek-nearest neighbors purity, defined as a jaccard index of neighbor sets (Online Methods). On the right dataset mixing for the benchmark algorithms on the four pairs of datasets evaluated by entropy of batch mixing.

Since our results depend on the neighborhood sizek, we consider a range of values -from a high resolution (k = 10) to a coarse (k = 500) view of the data. We compare scVI to several methods, including MNN [12], Seurat Alignment [13], ComBat [15] and principal component analysis (PCA). For each algorithm and pair of datasets, we report embeddings computed via a Uniform Manifold Approximation and Projection (UMAP) [42] (Figure 2a, Supplementary Figure 3 - 6) as well as the two evaluation metrics (Figure 2b-c). Overall, we observed that scVI performs well in terms of mixing, while comparing favorably to the other methods in terms of retainment of the original structure, for a wide range of neighborhood sizes and across all dataset pairs. As an example, consider the results of applying scVI and Seurat Alignment to the Tabula Muris bone marrow datasets (Figure 2a). As can be observed from both the entropy of batch mixing values and the visualization of the harmonized data, both methods perform well when it comes to mixing the input datasets. However, we observe that scVI performs better in retaining the original, continuous, structure of the data. Indeed, using cell type labels provided by the original publication [9], we observe that scVI captures the continuous process of myeloid cell development (from hematopoietic precursors to erythrocytes and monocytes) well, while Seurat Alignment groups the two trajectories (erythrocytes and monocytes) into separate clusters.

• Download figure
• Open in new tab

Figure 3:

Harmonizing datasets with different cellular composition. (a − d) show the case when no cell type is shared, where PBMC-8K contains all cells other than cell typec₀ and PBMC-CITE contains only cell typec₀. (a − b) UMAP visualization for the case wherec₀ corresponds to natural killer cells. (c − d) Harmonization statistics, aggregating the six experiments (settingc₀ to a different cell type in each experiment). (e − i) show the results of cell typec₀ is removed from the PBMC-8K dataset, while the PBMC-CITE dataset remains the same. (e − f) UMAP visualization for the case wherec₀ corresponds to CD4+ T cells. (g) batch entropy mixing for the cell type that was removed. (h) batch entropy mixing for each cell type in the 5 experiments where they were not removed. (i) Retainment of the original structure, aggregating all 6 experiments. Red arrows indicate the desired direction for each performance measure (e.g., low batch entropy is desirable in (d)). The boxplots are standard Tukey boxplots where the box is delineated by the first and third quartile and the whisker lines are the first and third quartile plus minus 1.5 times the box height. The dots are outliers that fall above or below the whisker lines.

While scANVI was designed for the problem of cell state annotation, we also wanted to evaluate its ability to harmonize datasets, which can be seen as a prerequisite. To evaluate this, we consider each dataset pair twice, each time using labels from one of the datasets (exploiting the semi-supervision framework of scANVI). Reassuringly, we found that scANVI is capable of effectively harmonizing the datasets, with a similar performance to that of scVI in terms of entropy of batch mixing andk-nearest neighbors purity (Figure 2b-c). We further explore the performance of scANVI in the annotation problem in the subsequent sections.

Harmonizing continuous trajectories

While so far we considered datasets that have a clear stratification of cells into discrete sub-populations, a conceptually more challenging case is harmonizing datasets in which the major source of variation forms a continuum, which inherently calls for accuracy at a higher level of resolution.

To explore this, we use a pair of datasets that provide a snapshot of hematopoiesis in mice (HEMATO-Tusi [43], HEMATO-Paul [44];Figure 4). These datasets consist of cells along the transition from common myeloid progenitor cells (Figure 4a-b; middle) through two primary differentiation trajectories myeloblast (top) and erythroblast-megakaryocyte (bottom). Notably, the HEMATO-Tusi dataset contains cells that appear to be more terminally differentiated, which are located at the extremes of the two primary branches. This can be discerned by the expression of marker genes (Figure 4e). For instance the HEMATO-Tusi unique erythroid cell population expressesHba-a2 (hemoglobin subunit) andAlas2 (erythroid-specific mitochondrial 5-aminolevulinate synthase) that are known to be present in reticulocytes [45,46]. At the other end, the granulocyte subset that is captured only by HEMATO-Tusi expressesItgam andS100a8.S100a8 is a neutrophil specific gene predicted by Nano-dissection [47] and is associated with GO processes such as leukocyte migration associated with inflammation and neutrophil aggregation.Itgam is not expressed in granulocyte-monocyte progenitor cells but is highly expressed in mature monocytes, mature eosinophils and macrophages [48]. We therefore do not expect mixing to take place along the entire trajectory. To account for this, we evaluated the extent of batch entropy mixing at different points along the harmonized developmental trajectory. As expected, we find that in most areas of the trajectory the two datasets are well mixed, while at the extremes, the entropy reduces significantly, using either scVI or Seurat Alignment (Figure 4c). Overall, we observe that scVI compares well in terms of both mixing the differentiation trajectories in each dataset and preserving their original, continuous, structure (Figure 4a-d).

• Download figure
• Open in new tab

Figure 4:

Harmonizing developmental trajectories. (a − b) UMAP visualization of the scVI latent space, with cells colored by the original labels from either the HEMATO-Paul (a) or HEMATO-Tusi (b) studies. The cells from the other dataset are colored in gray. (c) Batch entropy mixing of 20 of bins HEMATO-Tusi cells ordered by the potential of each cell derived from the PBA algorithm used in the original publication. Potential is a pseudotime measure that describes how differentiated a cell is. (center: common myeloid progenitors; moving left: erythrocyte branch; moving right: granulocyte branch). (d)k-nearest neighbors purity for scVI, Seurat, and scANVI. (e) Expression of marker genes that help determine the identity of batch-unique cells.

To validate scANVI in this context as well, we provided it with the categorical labels of cells along the two developmental trajectories, indicating their cell state (Figure 4c-d andSupplementary Figure 7). Even though this labeling scheme does not explicitly account for the ordering between states, we observe that scANVI is capable of mixing the two datasets, while retaining their original structure, achieving a level of accuracy comparable to that of scVI and better than that of Seurat Alignment.

Transferring cell type annotations between datasets

We next turned to evaluate scVI and scANVI in the context of harmonization-based annotation. Here, we test the extent to which annotations from a previously annotated dataset can be used to automatically derive annotations in a new unannotated dataset. For scVI and Seurat Alignment, we derive the annotations by first harmonizing the input datasets and then running ak-nearest neighbors classifier (settingk to 10) on the joint latent space, using the annotated cells to assign labels to the unannotated ones. Conversely, scANVI harmonizes the input datasets while using any amount of available labels. The prediction of unobserved labels is then conducted using the approximate posterior assignmentsq_Φ(c | z) of cell types (Online Methods), which can be derived directly from the model. An alternative approach, which we also include in our benchmark was taken by scmap-cluster [28], which instead of harmonizing, directly builds a classifier based on the labeled cells and then applies this classifier to the unlabeled cells.

We start by exploring the four dataset pairs inFigure 2, which have been annotated in their respective studies. In each experiment, we “hide” the cell type annotations from one dataset and transfer the second dataset labels to the first one. As a measure of performance, we report the weighted accuracy, which is the percent of cells that were correctly assigned to their correct (hidden) label, averaging over all labels (Online Methods). Importantly, the annotations in this first set of case studies were derived computationally. For example, by first clustering the cells, looking for marker genes expressed by each cluster and then assigning labels to the clusters accordingly. This level of annotation therefore makes the prediction problem relatively easy, and indeed, while we find that overall scANVI predicts unobserved labels more accurately, the differences between the methods are mild (Supplementary Figures 9 and 10).

To evaluate the accuracy of annotations without the need for computationally-derived labels, we turned to the PBMC-CITE dataset which includes measurements of ten key marker proteins in addition to mRNA [25], and the PBMC-sorted dataset [3], where cells were collected from bead purifications for eleven cell types (Supplementary Table 3). We applied scVI and scANVI to harmonize and annotate these two datasets along with a third dataset of PBMC (PBMC-68K [3]). Our analysis contains a combined set ofn = 169, 850 cells from the three datasets altogether. To generate a realistic scenario of cell type annotation, we only provide access to the experimentally-based labels from the PBMC-sorted dataset (Figure 5a-c). As an additional benchmark, we also evaluate Seurat Alignment, which was tested after removal of a randomly selected subset (40%) of the two large datasets (PBMC-68K and PBMC-sorted) due to scalability issues. Considering our harmonization performance measures (i.e., retainment of the original structure and batch mixing), we observe as before that scVI and scANVI perform similarly and compare favorably to Seurat Alignment. We then evaluated the accuracy of assigning unobserved labels, focusing on the PBMC-CITE dataset. Instead of using the labels from the original PBMC-CITE study as ground truth (which were computationally derived), we used the protein data, which provides an experimentally-derived proxy for cell state. To this end, we quantified the extent to which the similarity between cells in the harmonized mRNA-based latent space is consistent with their similarity at the protein level (Online Methods). We first computed the average discrepancy (sum of squared differences) between the protein measurements in each cell and the average over itsk-nearest neighbors. As a second measure we computed for each PBMC-CITE cell the overlap between itsk-nearest PBMC-CITE neighbors in the harmonized mRNA-based space and in the protein space. We then report the average across all cells inSupplementary Figure 11. Evidently, scANVI outperformed both scVI and Seurat Alignment for a wide range of neighborhood sizes, providing a representation for the mRNA data that is more consistent with the protein data (Figure 5c).

• Download figure
• Open in new tab

Figure 5:

Validation of cell type annotations using additional metadata. (a-b) UMAP plot of the scANVI latent space inferred for three harmonized datasets: PBMC-CITE [25], PBMC-sorted [3], and PBMC-68K [3]. Cells are colored by the batch they come from (a) and the PBMC-sorted labels(b). Cells from the PBMC-Cite and PBMC-68K are colored gray in (b). (c) The consistency of the harmonized PBMC-CITE mRNA data with the respective protein measurements, evaluated by mean squared error and for different neighborhood size. Lower values indicate higher consistency. (d) UMAP plot of the scANVI latent space, where cells are colored by normalized protein measurement. Only PBMC-CITE cells are displayed. (e) UMAP plot of the scANVI latent space, with cells from the PBMC-68k dataset colored according to their original label. For clarity of presentation, only cells originally labeled as dendritic cells or natural killer cells are colored. Evidently, a large number of these cells are mapped to a cluster of T-cells (right side of the plot).

As further support for the validity of these results, we confirmed that the labels inferred by scANVI in the PBMC-CITE dataset are consistent with the protein expression values (e.g., observing a uniquely high level ofCD19 expression in cells assigned with the “B cell” label;Figure 5d andSupplementary Figure 11). Interestingly, when considering the latent space learned by scANVI, we observed a certain amount of possible mislabeling in the original study of the PBMC-68K dataset (Figure 5e). In that study, the PBMC-68K cells were assigned with labels by taking the maximum correlation with cell subsets from the experimentally-annotated PBMC-sorted data. However this approach might be error prone, most likely due to low sensitivity and the influence of genes that are less relevant to cell type classification. Specifically, we observe a substantual amount of cells that are labeled by scANVI as T cells, but originally labeled as dendritic cells or natural killer cells. This re-annotation as T cells is supported by the expression of marker genes [49] of the respective cell subsets (Supplementary Figure 11).

Cell type annotation in a single dataset based on “seed” labels

An important variant of the annotation problem lies within the context of anab initio labeling of a single dataset where only a subset of the cells can be confidently annotated based on the raw data. This increasingly prevalent scenario may result from limited sensitivity of the scRNA-seq assay, where marker genes may only be confidently observed in a small subset of cells. One common way to address this problem is to compute some form of a distance metric between cells (e.g., after embedding with scVI or using Seurat PCA) and then assign labels based on proximity to annotated cells [3]. To benchmark our methods, we consider two such predictors: the first is clustering the cells and taking a majority vote inside each cluster, and the second is taking the majority vote of thek-nearest neighbors around each unannotated cell (k = 10). While these approaches are quite straightforward, their accuracy might suffer when the data does not form clear clusters [43], or when differences between labels are too subtle to be captured clearly by a transcriptome-wide similarity measure. To address these issues, scANVI takes an alternative approach, namely learning a latent embedding that is guided by the available labels, and then producing posterior probabilities for assigning labels to each cell.

As a case study, we compiled a dataset consisting of several experimentally sorted and labeled subsets of T cells from the PBMC-sorted dataset, including CD4 memory, CD4 naive, CD4 regulatory and CD8 naive. To make our analysis more realistic, we assume that the labels are completely unknown to us and therefore begin by trying to assign each T cell to its respective subset using marker genes (12 altogether; seeOnline Methods). Notably, several important biomarkers (CD4,CTLA4, andGITR) are detected in less than 5% of the cells, which renders their use for annotation not straightforward. Notably, many of these biomarkers are sparsely expressed to the extent that they are likely to be filtered out in the gene selection step of most harmonization procedures (Figure 6a).

• Download figure
• Open in new tab

Figure 6:

Cell type annotation in a single dataset using “seed” labeling. (a) discrepancies between marker genes that can be used to confidently label cells and highly variable genes in scRNA-seq analysis. (b − d) UMAP plot of the scVI latent space. (b) Seed cells are colored by their annotation (using known marker genes). (c) PBMC-sorted cell type labels from the original study based on marker-based sorting (d) The posterior probability of each cell being one of the four T cell subtypes obtained with scANVI.

To analyze this dataset, we first computed a signature score for each cell and for each label (i.e., T cell subset) using the scaled raw expression values of the respective marker genes (Online Methods). We then designated the top 50 scoring cells in each subset as the seed set of cells that are confidently annotated for that subset (Figure 6b). Reassuringly, this partial annotation is in agreement with the experimentally derived cell type labels available for this dataset (Figure 6c). However, this dataset does not form clear clusters, and in particular the seed sets of cells are not well separated. Such an observation makes clustering-based approaches potentially less precise. Indeed, usingk-means clustering on the scVI and Seurat Alignment latent space, we find that 74% and72% of the cells were assigned with their correct label. Similar analysis with two additional popular clustering algorithms (DBSCAN [50] and PhenoGraph [51]) further emphasizes the challenge of a cluster-based approach on this data. Specifically, DBSCAN does not partition the data into more than one cluster (scanning through a large number of parameter values;Online Methods), and PhenoGraph predicts 9 clusters and achieves an accuracy of 41% (Supplementary Figure 12).

Consistent with these results, the application of ak-nearest neighbors classifier resulted in a similar level of accuracy in the Seurat PCA latent space (71%), which is slightly improved when replacing it with the scVI latent space (73%;Supplementary Figure 12). Conversely, after fitting the scANVI model based on this partial labeling, the annotation posteriorq_Φ(c | z) (Figure 6d) provides a substantially more accurate cell type assignment, with 84% of cells annotated correctly. The scANVI latent space also has more distinct cell type clusters than both scVI and Seurat Alignment (Supplementary Figure 12).

Discussion

In this study, we demonstrated that scVI provides a principled approach to scRNA-seq data harmonization through joint probabilistic representation of multiple dataset, while accounting for technical hurdles such as over-dispersion, variable library size and limited sensitivity. We have demonstrated that scVI compares favorably to other methods in its accuracy and that it scales well, not only in terms of the number of cells (as in [29]) but also the number of input datasets (as opposed to other methods that work in a pairwise fashion and therefore scale quadratically with dataset size [19]). We have also shown that the harmonization step of scVI provides an effective baseline for automated transfer of cell type labels, from annotated datasets to new ones.

While the performance of scVI in the annotation problem compares favorably to other algorithms, it does not make use of any existing cell state annotations during model training, but rather after the latent space has been learned. To make better use of these annotations (which may be available for only some of the input datasets or only some cells within a dataset), we developed scANVI, a semi-supervised variant of scVI. While the latent space of scVI is defined by a Gaussian vector with diagonal unit variance, scANVI uses a mixture model, which enables it to directly represent the different cell states (each corresponding to a mixture component; seeOnline Methods) and provide a posterior probability of associating each cell with each label. We have demonstrated that similar to scVI, scANVI is capable of harmonizing datasets effectively. In addition, scANVI provides a way to address a number of variants of the annotation problem. Here, we have shown that it performs well in the most prevalent application of transferring labels from a reference dataset to an unannotated one, and that spurious annotations (when datasets do not overlap) are associated with low probability. We then demonstrated that scANVI can be used in the context of a single unannotated dataset, where high confidence (“seed”) labels are first inferred for a few cells (using marker genes) and then propagated to the remaining cells. Finally, we showed that scANVI is especially useful in the challenging case where the differences between cell states are too subtle to be captured clearly by a transcriptome-wide similarity measure, as well as in the case where the labels are organized in a hierarchy.

One concern in applying methods based on neural networks [21,57,58,59,60] in single-cell genomics and other domains is the robustness to hyperparameters choices [61]. This concern has been addressed to some extent by recent progress in the field, proposing search algorithms based on held-out log-likelihood maximization [59]. In this manuscript, we used an alternative approach that is more conducive for direct and easy application of our methods - namely we fix the hyperparameters and achieve state-of-the-art results on a substantial number of datasets and case studies.

Another important consideration while designing statistical models is the trade-off between goodness of fit and interpretability, which is still an open topic in machine learning research. Simple models such as the latent Dirichlet allocation [62] might not be particularly suited for scRNA-seq noise but are yet of interest, e.g., due to the immediate interpretation of word-level (i.e., gene-level) contributions to topics (i.e., cell types) [63]. Further effort in the use of deep generative models for applications in computational biology should come with attempts to perform model interpretation. For instance, SAUCIE [21] experimentally proposes to add an entropy regularization to a hidden layer of its denoising auto-encoder in order to infer clustering. Future principled efforts may focus on putting a suitable prior such as sparsity on neural networks weights (e.g., as in [64]). That way, individual neurons of the last hidden layer of the generative model would correspond to individual gene modules, directly readable from the weight sparsity motifs. Finally, as recent preprints propose proof of concepts for integrating single cell data across different data modalities such as Single molecule fluorescent in situ hybridization (smFISH), RNA-seq, ATAC-seq and DNA methylation [18,20], further work can utilize probabilistic graphical models that quantify measurement uncertainties in each assay, as well as the uncertainties of transferring information between modalities (e.g., predicting unmeasured gene expression in smFISH data).

An important distinguishing feature of both scVI and scANVI is that they rely on a fully-probabilistic model, thus providing a way to directly propagate uncertainties to any downstream analysis. While we have demonstrated this for differential expression analysis and cell type annotation, this can be incorporated to other tasks, such as differential abundance of sub-populations in case-control studies, correlation between genes and more. We therefore expect scVI, scANVI and similar tools to be of much interest as the field moves toward the goal of increasing reproducibility and consistency between studies and converging on to a common ontology of cell types. Both scVI and scANVI are available as an open-source software athttps://github.com/YosefLab/scVI. The code for reproducing the results in this manuscript has been deposited athttps://doi.org/10.5281/zenodo.2529945.

Author contributions

RL, EM, JR and NY conceived the statistical model. EM developed the software. CX, RL and EM applied the software to real data analysis. CX, RL, JR, NY, and MIJ wrote the manuscript. NY and MIJ supervised the work.

Competing interests

The authors declare no competing interests.

Online Methods

Supplementary Figures

• Download figure
• Open in new tab

Figure S1:

Robustness analysis for harmonization of the pair of datasets MarrowMT-10x / MarrowMT-ss2 with scVI. (a − c) We augment the number of hidden layers in the neural networkf_w and track acrossn = 5 random initializations for the batch entropy mixing (a), the held-out log likelihood (b) and the weighted accuracy of a nearest neighbor classifier on the latent space (c). (d − f) We increase the number of dimensions for the latent variablez and track acrossn = 5 random initialization the batch entropy mixing (d), the held-out log likelihood (e) and the weighted accuracy of a nearest neighbor classifier on the latent space (f)

• Download figure
• Open in new tab

Figure S2:

Schematic of the data harmonization problem. We are provided with two datasets (orange and blue), each consisting of two cell types (red and green). Our evaluation for the harmonization problem consists of two objectives: (1) mixing the two datasets well and (2) retaining the original structure in each dataset. Scenario 1 (top) is the case of under correction where objective (2) is achieved while objective (1) is not. Scenario 2 (middle) is the case of over correction where objective (1) is improved while objective (2) becomes worse. The bottom panel shows the desired scenario of mixing the datasets well while retaining the biological signal.

• Download figure
• Open in new tab

Figure S3:

Visualization of the benchmark PBMC-8K / PBMC-CITE. all positions for the scatter plots are derived using UMAP on the latent space of interest.

• Download figure
• Open in new tab

Figure S4:

Visualization of the benchmark MarrowMT 10x / ss2. all positions for the scatter plots are derived using UMAP on the latent space of interest.

• Download figure
• Open in new tab

Figure S5:

Visualization of the benchmark Pancreas InDrop / CEL-Seq2. all positions for the scatter plots are derived using UMAP on the latent space of interest.

• Download figure
• Open in new tab

Figure S6:

Visualization of the benchmark Dentate Gyrus 10X - Fluidigm C1. all positions for the scatter plots are derived using UMAP on the latent space of interest.

• Download figure
• Open in new tab

Figure S7:

(a − cb) Continuous trajectory obtained by the Seurat Alignment procedure for the HEMATO-Tusi and the HEMATO-Paul datasets. (c − d) Continuous trajectory obtained by the scANVI using the Tusi cell type labels for semi-supervision. (e − f) Continuous trajectory obtained by the scANVI using the Paul cluster labels for semi-supervision. All locations for scatter plots are computed via UMAP in their respective latent space.

• Download figure
• Open in new tab

Figure S8:

Large Scale data integration with scVI. (a − b) UMAP visualization of the scVI latent space colored by datasets (a) and by cell types (b). (c) accuracy of a nearest neighbor classifier based on scVI latent space

• Download figure
• Open in new tab

Figure S9:

Annotation results for the PBMC-8K / PBMC-CITE (a − b), the MarrowMT-10X / MarrowMT-SS2 datasets datasets (c d), Pancreas InDrop-CITESeq (e − f) and Dentate Gyrus 10X / Fluidigm C1 (g − h). Accuracies for transferring annotations from one dataset to another from ak-nearest neighbor classifier on Seurat Alignment,k-nearest neighbor classifier on scVI latent space, scANVI classifier and SCMAP classifier are shown. The aggregated results across for cell types that are shared between the two datasets is shown in box plots.

• Download figure
• Open in new tab

Figure S10:

Annotation results for the PBMC-8K / PBMC-CITE (a − b), the MarrowMT-10X / MarrowMT-SS2 datasets datasets (c − d), Pancreas InDrop-CITESeq (e − f) and Dentate Gyrus 10X / Fluidigm C1 (g − h). Accuracies for transferring annotations from one dataset to another from ak-nearest neighbors classifier on Seurat Alignment,k-nearest neighbors classifier on scVI latent space, scANVI classifier and SCMAP classifier are shown. The prediction accuracy for each cell type that is shared between the two datasets is shown on the y-axis and the size of the dots are proportional to the proportion of a cell type in the total population.

• Download figure
• Open in new tab

Figure S11:

(a)k-nearest neighbors purity of the merged latent space on the protein expression space as a function of the size of the neighborhood. (b) Protein expression heatmap showing consistency of PBMC-sorted labels and protein expression in PBMC-CITE. The protein expression per cell type is based onk-nearest neighbors imputation from the harmonized latent space obtained from scANVI trained with pure population labels. (c) We select individual cells that were labeled as Dendritic cells or Natural Killer cells in the original publication of the respective datasets, and compare the raw transcript count from cells inside the scANVI T cells cluster (DC*, NK*) against cells outside the T cells cluster (DC, NK). The expression of marker genes suggest that DC* and NK* is more likely to be T cells and thus the scANVI latent space is more accurate. (d) The batch entropy mixing of the three datasets in scVI, scANVI and Seurat Alignment merged space.

• Download figure
• Open in new tab

Figure S12:

Other methods of classifying T-cell subsets of the PBMC-Pure dataset. Coordinates for the scatter plots are derived from UMAP embedding based on the latent space of scANVI. (a) Ground truth labels from the purified PBMC populations (b)k-nearest neighbors classification labels when applied on scVI latent space from the seed set of cells (c)k-nearest neighbors classification labels when applied on Seurat Alignment latent space (d)k-means clustering based labels when applied to scVI latent space (e) DBSCAN clustering based labels when applied to scVI latent space. DBSCAN returns only one cluster but return some cells as unclassified. (f) PhenoGraph clusters on scVI latent space

• Download figure
• Open in new tab

Figure S13:

Differential Expression on mulitple datasets with scVI. Evaluation of consistency with the AUROC metric is shown for comparisons of CD4 vs CD8 T cells (a) and B cells vs Dendritic cells (b) on the PBMC-8K only (A), the PBMC-68k only (B) and the merged PBMC-8K / PBMC-68K (A+B) for scVI and edgeR. Error bars are obtained by multiple subsampling of the data to show robustness. (c − d) Same plots with the Kendall Tau metric.

Supplementary Note 1Statistical trade-offs in scRNA-seq harmonization

Harmonization is a hard and ill-defined problem. Especially, it can be difficult to formulate exactly its objective and at which level of granularity the “harmonization” is expected. Let us take the example of two PBMCs datasets that are exact biological replicates with same experimental conditions. In that case, the problem is well defined since we have a non-confounder assumption. Removing in a principled way all the variation in the data that corresponds to batches is reasonable. However, let’s take a slightly trickier case, one patient and one control dataset. If we look at the exact same metrics for harmonization, the resulting algorithm could discard any interesting biology that happens to the cells that get frozen. There is a confounder and the problem is much more difficult (ill-posed).

Therefore, it is extremely important to state how different models perform harmonization and what are their underlying hypothesis and modeling capabilities. For example, scmap [28] proposes a new gene filtering method to select gene that are claimed to be invariant to batch-effects. However, it is clear that over filtering can lead to ignoring biological information. MNNs [12] on its end assumes that the topology of cell types can be easily resolved by ak-nearest neighbor graph using a 𝕃² normalization and distance scheme. This allows MNNs [12] and all the neighbor-matching-based approaches [20,19] to remarkably merge batches with the risk of also merging cell types in the case where they are not detectable with this normalization scheme.

scVI and scANVI perform harmonization by learning a common generative model for a collection of gene expression probability distributions [p(x | z, s)]_{s∈{1, …,K}} indexed by the dataset-identifiers. The statistical richness of the collection of conditional variables [p(x | z, s)]_{s∈{1,…, K}} dictates the flexibility of our model towards integrating datasets.

One notable factor that sensibly contributes to harmonization capabilities is the prior for cell-specific scaling factorl_n that is dataset-specific. This helps probabilistically removing library-size caused discrepancies in the measurements and is more principled than normalization of the raw data [74]. Another important parameter is the parametrization of the neural network that maps variables (z, γ) to the expected frequencies 𝔼[w | z, s] =f_w(z, s). Since our functionf_w is now potentially non-linear, our model can benefit from more flexibility and fit batch-specific effects locally for each cell types or phenotypical condition. Especially, depending on how one designs the neural architecture off_w, it is possible to more flexibly correct dataset-specific effects. More specifically, we treatf_w as a feed-forward neural network for which at each layer we concatenate the batch-identifier with the hidden activations. A consequence of this design is that with more hidden layers inf_w, less parameters are shared and the family of distributions [p(x | z, s)]_{s∈{1, …,K}} has more flexibility to fit batch-specific effects (Supplementary Figure 1a-c). Conversely, when the latent dimension grows, the generative network has less incentive to use the dataset covariate and might mildly duplicate the information in its latent space (Supplementary Figure 1d-f). Throughout the paper, we fix those parameters (Online Methods and show competitive performance for all our datasets.

Another insight comes from how the variational distributionq(z | x, s) is parametrized. Our neural networks plays the role of an explicit stochastic mapping from the gene expressionx_n of a single-celln to a location in a latent spacez_n (a standard theme in scRNA-seq analysis, e.g LIGER [18], Seurat [13], scVI [29]). Via this map, we match an empirical distribution in gene expression space (i.e a dataset)p_data(x, s) = Ʃ_nδ_{(xn, sn)}(z), with a certain distribution latent spacep_data(x, s) = Ʃ_nδ_{Φ(xn, sn)}(z). In In certain cases, even though we designed this latent space to represent biological signal, the transformed dataset might still be confounded by technical effects [75]. In particular, this effect is susceptible to be severe when the generative model is not flexible enough to fit the dataset-specific effects or has a model misfit (as in, to some extent, the case with the CCA of [13] that assumes the data is log-normal). In this case, the go-to method is to empirically constrain the mapping (i.e the variational network in our case or the latent space itself in the case of SEURAT) to match the collections of variablesp_data(z, s)_{s∈{1,2026;;,K}}. However, each of these correction approaches has underlying modeling assumptions. For exam-ple, SAUCIE [21] and MMD ResNet [22] both propose to correct the variational distribution by adding on their objective function a non-parametric measure of distance between distributions (maximum mean discrepancy [76]). This approach specifically assumes that each dataset has the same cell type composition and the same biological signal and which is not reasonable in other cases than biological replicates. Seurat Alignment [13] on its end relies on a milder assumption that there is a common signal exactly reproducible between the two datasets and that CCA capture most of the biological variation, which is not obvious true considering limited suitability of linear Gaussian model for scRNA-seq. Seurat anchors [20] relies on CCA and suffer from the same problems as MNNs with its specific normalization scheme. Finally, the recent LIGER [18] method at first sight seems like a non-probabilistic version of scVI since it also learns a degenerate conditional distribution via its integrative non-negative matrix factorization [18] (NMF is a noisy-less version of a Gamma-Poisson generative model). However, it has a further quantile normalization of the latent spaces within clusters. First, the output of the clustering algorithm is not necessarily correct and could perturb downstream analyses. Second, if a cell type would be slightly different from one condition to another, this information would be lost in the final latent space. Overall, all these correction methods can therefore potentially lead to over-correction and statistical artifacts [32].

Supplementary Note 2 Evidence Lower Bound decomposition

For notational convenience, we do not write here that all the probabilities are conditioned on the batch identifiers. We also report the evidence lower bound (ELBO) only for one sample (i.e one cell) and drop the index notations by substituting {x_n, z_n, u_n, c_n, l_n}by {x, z, u, c, l}. This is without loss of generality since all the cells are independent and identically distributed under our model. Bayesian inference aims at maximizing the marginal likelihood of our datapoints. We also drop the parameters Θ (resp. Φ) of the generative model (resp. the variational distribution) for notational convenience. We will derive the ELBO in the case wherec is not observed (almost same calculations resolve the case wherec is observed). Notably, we derive the ELBO using Jensen inequality weighted by the variational distributionq(z, u, l, c | x). Similar derivations can be found in the variational autoencoder literature (eg. [65]). We assume our variational distribution factorizes as:

Then we simplify each individual term of the ELBO in 1 by recognizing KL divergence terms. In particular, we use subscript notation 𝕂𝕃_G (.||.), and 𝕂𝕃_M (.||.) to denote Gaussian and multinomial KL divergences.

Supplementary Note 3 Training Information

In order to train scanVI properly, several options are possible to train all the parameters (θ from the generative model,φ from the variational distribution andθ^C from the labels’ posterior exclusively). In all cases, parametersθ andφ should minimize the evidence lower bound (decomposed over the labeled samples 𝓛 and unlabeled samples 𝓤). Furthermore, [65] suggests to optimize jointly a modified objective that penalize the ELBO by a classification lossC on the labeled data so that the parametersθ^C also benefits from the learned data. In [77], we interpreted this extension of the ELBO as a correction term for noisy labels. In this manuscript, we propose a specific training procedure for scANVI, based on alternating minimization on the parameters (θ, φ) andθ^C. A joint training, conform to the one suggested in [65] is described inAlgorithm 1. Instead, an alternate training procedure would perform alternating minimization as described inAlgorithm 2. The advantage of joint training is presumably that it shapes the latent space directly, through the modification of the encoder’s weights. The advantage of the alternate training, is that the learning of the cell type edge converges properly at the end of each epoch, and improve indirectly the latent space quality, through the next steps of optimization. We have noticed in practice that using the joint training in the case of transferring labels from one dataset to another might break down the mixing in the latent space, as the lossC has only contributions from a unique dataset. Therefore, throughout the paper, we consistently used the alternating procedure ofAlgorithm 2.

• Download figure
• Open in new tab

Supplementary Note 4 Hierarchical classification

In this note, we explain the modeling details for an extension of scANVI where the cell types annotation can follow a two-level hierarchical structure. Let us note that this can in principle be adapted to any arbitrary tree representation of cell types taxonomy. This taxonomy needs to be hard-coded and knowna priori.

We do not need to formally modify the graphical model but only how we structure the variablec_n. Notably, we formally pose and consider the original Bayesian network as an hypergraph. In the generative model, the conditional distributionp(z | u, c) =p(z | u, y). Since the cell type group indicator variabley^g is not as refined than the cell type itselfy.

However, what sensibly changes is the parametrization of the variational distributionq(c | z) =q(y, y^g | z). The prior taxonomy knowledge encapsulates whether the assignment (y^g, y) = (i, j) is biologically possible (i.e cell typei is a sub-population of group cell typej). We encode this biological compatibility into a parent functionπ: {0, …,K} → {0, …,K^g} that maps a cell type to its parent in the hierarchy. We note for simplicity. We then use two neural networksf andf_g (with softmax non-linearities) that map the latent space to the joint approximate posteriorq(y, y^g | z) with the following rules:

Then, we can derive the marginal probability over finer cell types classes using the chain rule and Bayes rule:wherec(π_i) denotes the set of children of nodei children.

Acknowledgments

CX, RL and NY were supported by grant U19 AI090023 from NIH-NIAID. We thank Maxime Langevin, Yining Liu and Jules Samaran for helpful discussions and work on the scVI codebase as well as Allon Wagner and Chao Wang for their help on the choice for high specificity genes in the T cell study.

References

1. [1].↵
   A.Tanay andA.Regev, “Scaling single-cell genomics from phenomenology to mechanism,”Nature, vol.541, no.7637, pp.331–338,2017.

   OpenUrlCrossRefPubMed
2. [2].↵
   E.Macosko,A.Basu,R.Satija,J.Nemesh,K.Shekharet al., “Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets,”Cell, vol.161, no.5, pp.1202–1214,2017.

   OpenUrl
3. [3].↵
   G. X.Zheng,J. M.Terry,P.Belgrader,P.Ryvkin,Z. W.Bent,R.Wilson,S. B.Ziraldo,T. D.Wheeler,G. P.McDermott,J.Zhuet al., “Massively parallel digital transcriptional profiling of single cells,”Nature communications, vol.8, p.14049,2017.

   OpenUrl
4. [4].↵
   S.Semrau,J. E.Goldmann,M.Soumillon,T. S.Mikkelsen,R.Jaenisch, andA.Van Oudenaarden, “Dynamics of lineage commitment revealed by single-cell transcriptomics of differentiating embryonic stem cells,”Nature Communications, vol.8, no.1, p.1096,2017.

   OpenUrl
5. [5].↵
   J. T.Gaublomme,N.Yosef,Y.Lee,R. S.Gertneret al., “Single-cell genomics unveils critical regulators of Th17 cell pathogenicity,”Cell, vol.163, no.6, pp.1400–1412,2015.

   OpenUrlCrossRefPubMed
6. [6].↵
   A. P.Patel,I.Tirosh,J. J.Trombetta,A. K.Shalek,S. M.Gillespie,H.Wakimoto,D. P.Cahill,B. V.Nahed,W. T.Curry,R. L.Martuza,D. N.Louis,O.Rozenblatt-Rosen,M. L.Suv`a,A.Regev, andB. E.Bernstein, “Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma,”Science, vol.344, no.6190, pp.1396–1401, jun2014.

   OpenUrlAbstract/FREEFull Text
7. [7].↵
   P.Angerer,L.Simon,S.Tritschler,F. A.Wolf,D.Fischer, andF. J.Theis, “Single cells make big data: New challenges and opportunities in transcriptomics,”Current Opinion in Systems Biology, vol.4, pp.85 –91,2017.

   OpenUrl
8. [8].↵
   A.Regev,S. A.Teichmann,E. S.Lander,I.Amit,C.Benoistet al., “The human cell atlas,”eLife, vol.6, p.e27041, dec2017.

   OpenUrlCrossRefPubMed
9. [9].↵
   S. R.Quake,T.Wyss-Coray,S.Darmanis,T. M.Consortiumet al., “Single-cell transcriptomic characterization of 20 organs and tissues from individual mice creates a tabula muris,”bioRxiv, p.237446,2018.
10. [10].↵
    “Brain initiative cell census network (biccn),”2018. [Online]. Available:https://biccn.org/data/
11. [11].↵
    L.Wen andF.Tang, “Boosting the power of single-cell analysis,”Nature Biotechnology, vol.36, p.408, May2018.

    OpenUrl
12. [12].↵
    L.Haghverdi,A. T. L.Lun,M. D.Morgan, andJ. C.Marioni, “Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors,”Nature Biotechnology, vol.36, no.5, pp.421–427,2018.

    OpenUrlCrossRefPubMed
13. [13].↵
    A.Butler,P.Hoffman,P.Smibert,E.Papalexi, andR.Satija, “Integrating single-cell transcriptomic data across different conditions, technologies, and species,”Nature Biotechnology, vol.36, p.411, Apr2018.

    OpenUrlCrossRefPubMed
14. [14].↵
    A.Wagner,A.Regev, andN.Yosef, “Revealing the vectors of cellular identity with single-cell genomics,”Nature Biotechnology, vol.34, no.11, pp.1145–1160,2016.

    OpenUrlCrossRefPubMed
15. [15].↵
    W. E.Johnson,C.Li, andA.Rabinovic, “Adjusting batch effects in microarray expression data using empirical Bayes methods,”Biostatistics, vol.8, no.1, pp.118–127,2007.

    OpenUrlCrossRefPubMedWeb of Science
16. [16].↵
    M. E.Ritchie,B.Phipson,D.Wu,Y.Hu,C. W.Law,W.Shi, andG. K.Smyth, “limma powers differential expression analyses for rna-sequencing and microarray studies,”Nucleic Acids Res, vol.43, no.7, pp.e47–e47, Apr2015.

    OpenUrlCrossRefPubMed
17. [17].↵
    D.Risso,F.Perraudeau,S.Gribkova,S.Dudoit, andJ.-P.Vert, “A general and flexible method for signal extraction from single-cell RNA-seq data,”Nature Communications, vol.9, p.284,2018.

    OpenUrl
18. [18].↵
    J.Welch,V.Kozareva,A.Ferreira,C.Vanderburg,C.Martin, andE.Macosko, “Integrative inference of brain cell similarities and differences from single-cell genomics,”bioRxiv,2018. [Online]. Available:https://www.biorxiv.org/content/early/2018/11/02/459891
19. [19].↵
    B. L.Hie,B.Bryson, andB.Berger, “Panoramic stitching of heterogeneous single-cell transcriptomic data,”bioRxiv,2018. [Online]. Available:https://www.biorxiv.org/content/early/2018/07/17/371179
20. [20].↵
    T.Stuart,A.Butler,P.Hoffman,C.Hafemeister,E.Papalexi,W. M.Mauck,M.Stoeckius,P.Smibert, andR.Satija, “Comprehensive integration of single cell data,”bioRxiv,2018. [Online]. Available:https://www.biorxiv.org/content/early/2018/11/02/460147
21. [21].↵
    M.Amodio,K.Srinivasan,D.van Dijk,H.Mohsen,K.Yim,R.Muhle,K. R.Moon,S.Kaech,R.Sowell,R.Montgomery,J.Noonan,G.Wolf, andS.Krishnaswamy, “Exploring Single-Cell Data with Multitasking Deep Neural Networks,”bioRxiv,2017. [Online]. Available:https://www.biorxiv.org/content/early/2017/12/19/237065
22. [22].↵
    U.Shaham,K. P.Stanton,J.Zhao,H.Li,K.Raddassi,R.Montgomery, andY.Kluger, “Removal of batch effects using distribution-matching residual networks,”Bioinformatics, vol.33, no.16, pp.2539–2546,2017.

    OpenUrlCrossRef
23. [23].↵
    D.DeTomaso,M.Jones,M.Subramaniam,T.Ashuach,C. J.Ye, andN.Yosef, “Functional interpretation of single-cell similarity maps,”bioRxiv,2018. [Online]. Available:https://www.biorxiv.org/content/early/2018/08/29/403055
24. [24].↵
    D.DeTomaso andN.Yosef, “FastProject: A tool for low-dimensional analysis of single-cell RNA-seq data,”BMC Bioinformatics, no.1, p.315, dec2016.
25. [25].↵
    M.Stoeckius,C.Hafemeister,W.Stephenson,B.Houck-Loomis,P. K.Chattopadhyay,H.Swerdlow,R.Satija, andP.Smibert, “Simultaneous epitope and transcriptome measurement in single cells,”Nature methods, vol.14, no.9, p.865,2017.

    OpenUrl
26. [26].↵
    C.Weinreb,A. E.Rodriguez-Fraticelli,F. D.Camargo, andA. M.Klein, “Lineage tracing on transcriptional landscapes links state to fate during differentiation,”bioRxiv,2018. [Online]. Available:https://www.biorxiv.org/content/early/2018/11/11/467886
27. [27].↵
    F.Wagner andI.Yanai, “Moana: A robust and scalable cell type classification framework for single-cell RNA-Seq data,”bioRxiv,2018. [Online]. Available:https://www.biorxiv.org/content/early/2018/10/30/456129
28. [28].↵
    V. Y.Kiselev,A.Yiu, andM.Hemberg, “scmap: projection of single-cell RNA-seq data across data sets,”Nature Methods, vol.15, p.359, Apr2018.

    OpenUrl
29. [29].↵
    R.Lopez,J.Regier,M. B.Cole,M.Jordan, andN.Yosef, “Deep Generative Modeling for Single-cell Transcriptomics,”Nature Methods,2018.
30. [30].↵
    D. P.Kingma andM.Welling, “Auto-encoding variational Bayes,”in The International Conference on Learning Representations,2014.
31. [31].↵
    C.Louizos,K.Swersky,Y.Li,M.Welling, andR.Zemel, “The Variational Fair Autoencoder,”in International Conference on Learning Representations,2016.
32. [32].↵
    V.Nygaard,E. A.Rødland, andE.Hovig, “Methods that remove batch effects while retaining group differences may lead to exaggerated confidence in downstream analyses,”Biostatistics, vol.17, no.1, pp.29–39, Jan2016.

    OpenUrlCrossRefPubMed
33. [33].↵
    “10x genomics,”2017. [Online]. Available:https://support.10xgenomics.com/single-cell-gene-expression/datasets
34. [34].↵
    S.Picelli,O. R.Faridani,Å. K.Björklund,G.Winberg,S.Sagasser, andR.Sandberg, “Full-length rna-seq from single cells using smart-seq2,”Nature protocols, vol.9, no.1, p.171,2014.

    OpenUrl
35. [35].↵
    Y.Xin,J.Kim,M.Ni,Y.Wei,H.Okamoto,J.Lee,C.Adler,K.Cavino,A. J.Murphy,G. D.Yancopouloset al., “Use of the fluidigm c1 platform for rna sequencing of single mouse pancreatic islet cells,”Proceedings of the National Academy of Sciences, vol.113, no.12, pp.3293–3298,2016.

    OpenUrlAbstract/FREEFull Text
36. [36].↵
    D. A.Jaitin,E.Kenigsberg,H.Keren-Shaul,N.Elefant,F.Paul,I.Zaretsky,A.Mildner,N.Cohen,S.Jung,A.Tanayet al., “Massively parallel single-cell rna-seq for marker-free decomposition of tissues into cell types,”Science, vol.343, no.6172, pp.776–779,2014.

    OpenUrlAbstract/FREEFull Text
37. [37].↵
    A. M.Klein,L.Mazutis,I.Akartuna,N.Tallapragada,A.Veres,V.Li,L.Peshkin,D. A.Weitz, andM. W.Kirschner, “Droplet barcoding for single-cell transcriptomics applied to embryonic stem cells,”Cell, vol.161, no.5, pp.1187–1201,2015.

    OpenUrlCrossRefPubMed
38. [38].↵
    T.Hashimshony,N.Senderovich,G.Avital,A.Klochendler,Y.de Leeuw,L.Anavy,D.Gennert,S.Li,K. J.Livak,O.Rozenblatt-Rosenet al., “Cel-seq2: sensitive highly-multiplexed single-cell rna-seq,”Genome biology, vol.17, no.1, p.77,2016.

    OpenUrlCrossRefPubMed
39. [39].↵
    C.Ziegenhain,B.Vieth,S.Parekh,B.Reinius,A.Guillaumet-Adkins,M.Smets,H.Leonhardt,H.Heyn,I.Hellmann, andW.Enard, “Comparative analysis of single-cell rna sequencing methods,”Molecular Cell, vol.65, no.4, pp.631–643.e4, Feb2017.

    OpenUrlCrossRefPubMed
40. [40].↵
    M.Baron,A.Veres,S. L.Wolock,A. L.Faust,R.Gaujoux,A.Vetere,J. H.Ryu,B. K.Wagner,S. S.Shen-Orr,A. M.Kleinet al., “A single-cell transcriptomic map of the human and mouse pancreas reveals inter-and intra-cell population structure,”Cell systems, vol.3, no.4, pp.346–360,2016.

    OpenUrl
41. [41].↵
    H.Hochgerner,A.Zeisel,P.Lönnerberg, andS.Linnarsson, “Conserved properties of dentate gyrus neurogenesis across postnatal development revealed by single-cell rna sequencing,”Nature neuroscience, vol.21, no.2, p.290,2018.

    OpenUrlCrossRefPubMed
42. [42].↵
    L.McInnes andJ.Healy, “UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction,”arXiv, Feb.2018.
43. [43].↵
    B. K.Tusi,S. L.Wolock,C.Weinreb,Y.Hwang,D.Hidalgo,R.Zilionis,A.Waisman,J. R.Huh,A. M.Klein, andM.Socolovsky, “Population snapshots predict early haematopoietic and erythroid hierarchies,”Nature, vol.555, no.7694, pp.54–60, 22018.

    OpenUrl
44. [44].↵
    F.Paul,Y.Arkin,A.Giladi,D.Jaitin,E.Kenigsberg,H.Keren-Shaul,D.Winter,D.Lara-Astiaso,M.Gury,A.Weiner,E.David,N.Cohen,F.Lauridsen,S.Haas,A.Schlitzer,A.Mildner,F.Ginhoux,S.Jung,A.Trumpp,B.Porse,A.Tanay, andI.Amit, “Transcriptional heterogeneity and lineage commitment in myeloid progenitors,”Cell, vol.163, no.7, pp.1663 –1677,2015.

    OpenUrlCrossRefPubMed
45. [45].↵
    S.-H.Goh,M.Josleyn,Y. T.Lee,R. L.Danner,R. B.Gherman,M. C.Cam, andJ. L.Miller, “The human reticulocyte transcriptome,”Physiological genomics, vol.30, no.2, pp.172–178,2007.

    OpenUrlCrossRefPubMedWeb of Science
46. [46].↵
    “Transcription profiling by high throughput sequencing of murine hematopoietic progenitors and lineage cells.” [Online]. Available:https://www.ebi.ac.uk/gxa/experiments/E-MTAB-3079/Results
47. [47].↵
    W.Ju,C. S.Greene,F.Eichinger,V.Nair,J. B.Hodgin,M.Bitzer,Y.-s.Lee,Q.Zhu,M.Kehata,M.Liet al., “Defining cell-type specificity at the transcriptional level in human disease,”Genome research, pp. gr–155697,2013.
48. [48].↵
    I.Papatheodorou,N. A.Fonseca,M.Keays,Y. A.Tang,E.Barrera,W.Bazant,M.Burke,A.Füllgrabe,A. M.-P.Fuentes,N.Georgeet al., “Expression atlas: gene and protein expression across multiple studies and organisms,”Nucleic acids research, vol.46, no.D1, pp.D246–D251,2017.

    OpenUrl
49. [49].↵
    T.S. P.Heng,M. W.Painter,T. I. G. P.Consortium,K.Elpek,V.Lukacs-Kornek,N.Mauermann,S. J.Turley,D.Koller,F. S.Kim,A. J.Wagers,N.Asinovski,S.Davis,M.Fassett,M.Feuerer,D. H. D.Gray,S.Haxhinasto,J. A.Hill,G.Hyatt,C.Laplace,K.Leatherbee,D.Mathis,C.Benoist,R.Jianu,D. H.Laidlaw,J. A.Best,J.Knell,A. W.Goldrath,J.Jarjoura,J. C.Sun,Y.Zhu,L. L.Lanier,A.Ergun,Z.Li,J. J.Collins,S. A.Shinton,R. R.Hardy,R.Friedline,K.Sylvia, andJ.Kang, “The immunological genome project: networks of gene expression in immune cells,”Nature Immunology, vol.9, p.1091, Oct2008.

    OpenUrlCrossRefPubMedWeb of Science
50. [50].↵
    M.Ester,H.-P.Kriegel,J.Sander, andX.Xu, “A density-based algorithm for discovering clusters a density-based algorithm for discovering clusters in large spatial databases with noise,”in Proceedings of the Second International Conference on Knowledge Discovery and Data Mining,1996, pp.226–231.
51. [51].↵
    J. H.Levine,E. F.Simonds,S. C.Bendall,K. L.Davis,E. ad D.Amir,M. D.Tadmor,O.Litvin,H. G.Fienberg,A.Jager,E. R.Zunder,R.Finck,A. L.Gedman,I.Radtke,J. R.Downing,D.Pe’er, andG. P.Nolan, “Data-driven phenotypic dissection of AML reveals progenitor-like cells that correlate with prognosis,”Cell, vol.162, no.1, pp.184–197, jul2015.

    OpenUrlCrossRefPubMed
52. [52].↵
    A.Zeisel,H.Hochgerner,P.Lönnerberg,A.Johnsson,F.Memic,J.van der Zwan,M.Häring,E.Braun,L. E.Borm,G.La Manno,S.Codeluppi,A.Furlan,K.Lee,N.Skene,K. D.Harris,J.Hjerling-Leffler,E.Arenas,P.Ernfors,U.Marklund, andS.Linnarsson, “Molecular architecture of the mouse nervous system,”Cell, Aug2018.
53. [53].↵
    H. I.Nakaya,J.Wrammert,E. K.Lee,L.Racioppiet al., “Systems biology of vaccination for seasonal influenza in humans,”Nature Immunology, vol.12, no.8, pp.786–795, 72011.

    OpenUrlCrossRefPubMed
54. [54].↵
    G.Görgün,T. A. W.Holderried,D.Zahrieh,D.Neuberg, andJ. G.Gribben, “Chronic lymphocytic leukemia cells induce changes in gene expression of CD4 and CD8 T cells,”The Journal of Clinical Investigation, vol.115, no.7, pp.1797–805,2005.

    OpenUrlCrossRefPubMed
55. [55].↵
    M. D.Robinson,D. J.McCarthy, andG. K.Smyth, “edger: a bioconductor package for differential expression analysis of digital gene expression data,”Bioinformatics, vol.26, no.1, pp.139–140,2010.

    OpenUrlCrossRefPubMedWeb of Science
56. [56].↵
    C.Soneson andM. D.Robinson, “Bias, robustness and scalability in single-cell differential expression analysis,”Nature methods, vol.15, no.4, p.255,2018.

    OpenUrl
57. [57].↵
    J.Ding,A.Condon, andS. P.Shah, “Interpretable dimensionality reduction of single cell transcriptome data with deep generative models,”Nature Communications, vol.9, no.1, p.2002, dec2018.

    OpenUrl
58. [58].↵
    D.Wang andJ.Gu, “VASC: Dimension reduction and visualization of single cell RNA sequencing data by deep variational autoencoder,”bioRxiv,2017. [Online]. Available:https://www.biorxiv.org/content/early/2017/10/06/199315
59. [59].↵
    G.Eraslan,L. M.Simon,M.Mircea,N. S.Mueller, andF. J.Theis, “Single cell RNA-seq denoising using a deep count autoencoder,”bioRxiv,2018. [Online]. Available:https://www.biorxiv.org/content/early/2018/04/13/300681
60. [60].↵
    C. H.Grønbech,M. F.Vording,P. N.Timshel,C. K.Sønderby,T. H.Pers, andO.Winther, “scVAE: Variational auto-encoders for single-cell gene expression data,”bioRxiv,2018. [Online]. Available:https://www.biorxiv.org/content/early/2018/05/16/318295
61. [61].↵
    Q.Hu andC. S.Greene, “Parameter tuning is a key part of dimensionality reduction via deep variational autoencoders for single cell rna transcriptomics,”bioRxiv,2018. [Online]. Available:https://www.biorxiv.org/content/early/2018/09/20/385534
62. [62].↵
    D. M.Blei,A. Y.Ng, andM. I.Jordan, “Latent dirichlet alocation,”Journal of Machine Learning Research, vol.3, pp.993–1022, Mar.2003.

    OpenUrlCrossRefWeb of Science
63. [63].↵
    D. A.duVerle,S.Yotsukura,S.Nomura,H.Aburatani, andK.Tsuda, “CellTree: an R/bioconductor package to infer the hierarchical structure of cell populations from single-cell RNA-seq data,”BMC Bioinformatics, Sep2016.
64. [64].↵
    A. J.Riesselman,J. B.Ingraham, andD. S.Marks, “Deep generative models of genetic variation capture the effects of mutations,”Nature Methods, vol.15, no.10, pp.816–822,2018. [Online]. Available:https://doi.org/10.1038/s41592-018-0138-4

    OpenUrl
65. [65].↵
    D. P.Kingma,S.Mohamed,D.Jimenez Rezende, andM.Welling, “Semi-supervised Learning with Deep Generative Models,”in Advances in Neural Information Processing Systems27,2014, pp.3581–3589.

    OpenUrl
66. [66].↵
    D. P.Kingma andJ.Ba, “Adam: A method for stochastic optimization,”in International Conference on Learning Representations,2015.
67. [67].↵
    C. K.Sønderby,T.Raiko,L.Maaløe,S. K.Sønderby, andO.Winther, “Ladder variational autoencoders,”in Advances in Neural Information Processing Systems,2016.
68. [68].↵
    S.Ioffe andC.Szegedy, “Batch normalization: accelerating deep network training by reducing internal covariate shift,”in International Conference on Machine Learning,2015.
69. [69].↵
    L.Held andM.Ott, “On p-values and Bayes factors,”Annual Review of Statistics and Its Application, vol.5, no.1, pp.393–419,2018.

    OpenUrl
70. [70].↵
    R. E.Kass,A. E.Raftery,S.Association, andN.Jun, “Bayes factors,”Journal of the American Statistical Association, vol.90, no.430, pp.773–795,1995.

    OpenUrlCrossRefPubMedWeb of Science
71. [71].↵
    T. S.Andrews andM.Hemberg, “Dropout-based feature selection for scRNASeq,”bioRxiv,2018. [Online]. Available:https://www.biorxiv.org/content/early/2018/05/17/065094
72. [72].↵
    M. B.Cole,D.Risso,A.Wagner,D.DeTomaso,J.Ngai,E.Purdom,S.Dudoit, andN.Yosef, “Performance assessment and selection of normalization procedures for single-cell RNA-seq,”bioRxiv,2017. [Online]. Available:https://www.biorxiv.org/content/early/2018/05/18/235382
73. [73].↵
    M. J.Muraro,G.Dharmadhikari,D.Grün,N.Groen,T.Dielen,E.Jansen,L.van Gurp,M. A.Engelse,F.Carlotti,E. J.de Koninget al., “A single-cell transcriptome atlas of the human pancreas,”Cell systems, vol.3, no.4, pp.385–394,2016.

    OpenUrl
74. [74].↵
    C. A.Vallejos,D.Risso,A.Scialdone,S.Dudoit, andJ. C.Marioni, “Normalizing single-cell RNA sequencing data: challenges and opportunities,”Nature Methods, pp.565–571,2017.
75. [75].↵
    R.Lopez,J.Regier,M. I.Jordan, andN.Yosef, “Information Constraints on Auto-Encoding Variational Bayes,”Advances in Neural Information Processing Systems,2018.
76. [76].↵
    A.Gretton,K. M.Borgwardt,M. J.Rasch,B.Schölkopf, andA.Smola, “A Kernel Two-Sample Test,”Journal of Machine Learning Research, vol.13, pp.723–773,2012.

    OpenUrl
77. [77].↵
    M.Langevin,E.Mehlman,J.Regier,R.Lopez,M. I.Jordan, andN.Yosef, “A deep generative model for semi-supervised classification with noisy labels,”arXiv, vol.abs/1809.05957,2018. [Online]. Available:http://arxiv.org/abs/1809.05957