dataesr/scientific-paragraphs-categorization huggingface dataset:
Harmonizing continuous trajectories While so far we considered datasets that have a clear stratification of cells into discrete subpopulations, 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  [49] , HEMATO-Paul  [50] ; 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.
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 ). 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 and Supplementary Figure  9 ). 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.
Rapid integration of multiple datasets To demonstrate the scalability of our framework in the context of harmonizing multiple (and possibly large) dataset, we ran scVI to integrate a cohort of 26 datasets spanning 105,476 cells from multiple tissues and technologies, which was made available by the authors of Scanorama (a method based on truncated singular value decomposition followed by nearest neighbor matching  [19] ). Using the hardware specified in the original paper  [19]  (Intel Xeon E5-2650v3 CPU limited to 10 cores with 384 GB of RAM), Seurat Alignment and MNN required over 24 hours, while Scanorama completed its run in 20 minutes. Using a simpler configuration (eight-core Intel i7-6820HQ CPU with 32 GB RAM) along with one NVIDIA Tesla K80 GPU (GK210GL; addressing 24 GB RAM), we found that scVI integrates all datasets and learns a common embedding in less than 50 minutes.
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 a k-nearest neighbors classifier (setting k 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 assignments q Φ (c | x) 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. Finally, we also applied the domain adaptation method Correlation Alignment for Unsupervised Domain Adaptation (CORAL,  [55] ). This method was not initially developed for single-cell analysis but is an insightful benchmark from the machine learning literature. We start by exploring the four dataset pairs in Figure  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).
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  4 ). 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 of n = 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 its k-nearest neighbors.
As a second measure we computed for each PBMC-CITE cell the overlap between its k-nearest PBMC-CITE neighbors in the harmonized mRNA-based space and in the protein space. We then report the average across all cells in Supplementary Figure  13. 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 ). 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 of CD19 expression in cells assigned with the "B cell" label; Figure  5d  and Supplementary Figure  13 ). 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 substantial 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  [56]  of the respective cell subsets (Supplementary Figure  13 ). Cell type annotation in a single dataset based on "seed" labels An important variant of the annotation problem lies within the context of an ab initio labeling of a single dataset where only a subset of the cells can be confidently annotated based on the raw data.
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; see Online Methods). Notably, several important biomarkers (CD4, CTLA4, and GITR) are detected in less than 5% of the cells, which renders their use for annotation not straightforward.
Furthehrmore, 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 ). 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.
Cell type taxonomy and hierarchical classification with scANVI Another subtle yet important variation of the annotation problem is when the labels are not mutually exclusive but rather form a taxonomy of cell types or states. To effectively annotate cells in this setting, we extended scANVI to perform hierarchical classification, which as before we carry out from first principles, relying on probabilistic graphical models (Online Methods). To demonstrate this extended version, we use a dataset of the mouse nervous system  [60]  that was annotated using a cell type taxonomy with several levels of granularity. At the lowest (most granular) level, the cells are stratified into 265 cell sub-types. At the second lowest level of granularity these 265 subtypes are grouped into 39 subsets, each corresponding to a more coarse definition of a cell type.
g. , increasing the number of classifier training epochs; see Supplementary Note 3). However, when we take the hierarchy of the labels into account, the performance of scANVI increases to 37%, thus outperforming the other methods by a significant margin. Notably, while we tested the extrapolation of seed labeling and the hierarchical mode only in the context of a single dataset, this variation of the scANVI model can also be directly applied in the context of multiple datasets (i. e. , transferring hierarchical annotations between datasets).
Hypotheses testing in harmonized datasets: the case of differential expression With their probabilistic representation of the data, scVI and scANVI each provide a natural way of performing various types of hypotheses testing (Online Methods). This is different from other approaches  [12, 13, 18, 19, 20]  where the dataset alignment procedures do not carry direct probabilistic interpretation, and the resulting harmonized data can thus not be directly used for these purposes. To demonstrate this, we focus on the problem of differential expression. As a first case study, we use two of the PBMC datasets (PBMC-8K and PBMC-68K) and looked for differentially expressed genes in two settings: comparing the B cells to dendritic cells, and similarly for CD4+ versus CD8+ T cells. For evaluation, we used reference sets of differentially expressed genes that were obtained from published bulk-level analysis of similar cell subsets (microarrays,  [61, 62] , as in  [29] ).
While this benchmark relies on real data, a clear caveat is the lack of a well defined ground truth. To address this, we used a second benchmark based on simulations with Symsim  [59]. The simulated data consists of five subpopulations of varying degrees of transcriptional distance, profiled in two different "batches" of different technical quality (Online Methods, Supplementary Figure  16 ). This framework allowed us to derive an exact log fold changes (LFC) between every pair of simulated subpopulations, which enable a more accurate evaluation of performance. In both benchmark studies, we assume that labels are only available for one of the two input batches or datasets (in the real data we assume that PBMC-8K is the annotated one). To apply scVI, we first harmonized the input pair of datasets and transferred labels using a k-nearest neighbors classifier on the joint latent space (k = 10).
We then consider these annotations (predicted and pre-labeled) as fixed and sample 100 cell pairs, each pair consisting of one cell from each population. For each cell pair we sample gene expression values from the variational posterior, while marginalizing over the different datasets, to compute the probability for differential expression in a dataset-agnostic manner. Aggregating across all selected pairs results in approximate Bayes factors that reflect the evaluated extent of differential expression (Online Methods). Since scANVI assigns posterior probability for associating any cell to any label, it enables a more refined scheme. Specifically, instead of sampling pairs of cells we are sampling pairs of points in the latent space, while conditioning on the respective label. This approach therefore does not assume a fixed label for each cell (or point in latent space) as in the scVI scheme, but rather a distribution of possible labels thus making it potentially more robust to mis-labeling.
As a reference, we look at the accuracy of differential expression analysis in each PBMC dataset separately (using their prior annotations to define the sets of cells we are comparing), which can computed with scVI (as in  [29] ) and edgeR. Reassuringly, we observe that the performance of scVI on the joint data is not lower than it is in either dataset in isolation. We also find that while scVI performs moderately better than scANVI, both methods compare favorably to edgeR in their accuracy. In our simulations, we considered differential expression between every possible pair out of the five simulated subpopulations. For evaluation, we computed the Spearman and Kendall rank correlation coefficients between the true LFC and the inferred Bayes factors (for scVI and scANVI) or estimated LFC (for edgeR).
Our results in Supplementary Figure  17ab  show that with this artificial, yet more clearly defined objective, scVI was substantially more accurate than edgeR and that in the harmonized data scANVI provided more exact and stable estimates than scVI (Supplementary Figure  17 ). Mislabeling of a certain proportion of cells in a dataset is a plausible scenario that may occur in any study. An important challenge is therefore to maintain the validity of downstream analysis despite such "upstream" annotation errors. To evaluate robustness in this setting, we repeated the simulation analysis, while introducing labeling errors at different rates. Specifically, prior to evaluating differential expression between two simulated sub-populations, we flip the labels of a certain proportion (up to 30%) of the respective cells in the annotated batch.
We then proceed as before and assign labels to cells in the unannotated batch by scVI or scANVI, followed by differential expression analysis. Our results (Supplementary Figure  18a ) suggest that scANVI is clearly more robust to this type of mislabeling than scVI (or edgeR, applied on the scVI-derived labels). Repeating the same analysis on the PBMC data (where the differential expression ground truth is obviously not available), we observe similar level of robustness in scANVI, albeit with not much difference compared to scVI and edgeR. Overall, our results demonstrate that both scVI and scANVI are capable of conducting differential expression effectively, while working directly on a harmonized dataset. Furthermore, we observe that both methods and especially scANVI are robust to mislabeling, providing further motivation for explicitly modeling label uncertainty.
Discussion In this study, we demonstrated that scVI provides a principled approach to harmonization of scRNA-seq data through joint probabilistic representation of multiple dataset, while accounting for technical hurdles such as 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; see Online 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 first shown that it performs well in the most prevalent application of transferring labels from a reference dataset to an unannotated one.
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 have shown 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. Notably, although scANVI achieves high accuracy when transferring labels from one dataset to another, it was not designed to automatically identify previously unobserved labels. Indeed, in Supplementary Figure  19  we demonstrate that increasing the number of labels in the model (C) to values beyond the number of observed labels does not alter the results much.
Nevertheless, we observed that unannotated cell populations that have an unobserved label are associated with low levels of mixing between the input datasets. We therefore advocate that clusters from an unannotated dataset that do not mix well should be inspected closely and, if appropriate, should be manually assigned with a new label. One concern in applying methods based on neural networks  [21, 65, 66, 67, 68]  in single-cell genomics and other domains is the robustness to hyperparameters choices  [69]. This concern has been addressed to some extent by recent progress in the field, proposing search algorithms based on held-out log-likelihood maximization  [67]. 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.
The development of scVI and scANVI required several modeling and implementation choices. In Supplementary Note 4, Supplementary Figures  20  and 21  we discuss the rationale behind the choice of a zero-inflated negative binomial (ZINB) distribution as well as robustness to choice of priors. Briefly, we find that exclusion of zero inflation from the model results in approximately similar performance, except for the case of harmonizing Smart-Seq2 and 10x data sets, in which ZINB performs significantly better. Such results might suggest that zero-inflation may be more suitable for certain technologies than others. Similarly, we investigate the prior on the library size which is defined per batch and show that computing the same prior for both the datasets (rather than each dataset individually, as we do by default) affects the performance only in the case of harmonizing the same pair of Smart-Seq2 and 10x datasets.
Since these datasets have very different sources of technical noise, this may suggest that it is indeed advisable to explicitly account for such differences during model fit. 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  [70]  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)  [71]. 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.
Equations to derive the evidence lower bound (ELBO) are derived in Supplementary Note 6. We optimize the sum ELBO = L + U over the neural networks parameters and the inverse-dispersion parameters (in a variational Bayesian inference fashion). Remarkably, the approximate posterior q Φ (c n | z n ) can be used as a classifier, assigning cells to cell types based on the location on the latent space. We sample from the variational posterior using the reparametrization trick  [30]  as well as "mini-batches" from the dataset to compute unbiased estimate of the objective gradients' with respect to the parameters. We use Adam  [78]  as a first-order stochastic optimizer to update the model parameters.
This set includes the C observed populations for which annotated cells are available, and C unobserved population that accounts for cell types for which an annotation is not available to the algorithm. Ad-hoc training procedures for scANVI inference are described in Supplementary Note 3. In the case of a one single dataset, we use a 2 fully-connected layers with 256 hidden units classifier with c = 1 epochs of classifier training in between each variational update. In the case of transfer of labels, we use Algorithm 2 with 1 fully-connected layers with 128 hidden units classifier, with c = 100 epochs of classifier training in between each variational update.
Bayesian differential expression scVI For each gene g and pair of cells (z a , z b ) with observed gene expression (x a ,x b ) and dataset identifier (s a ,s b ), we can formulate two mutually exclusive hypotheses: H g 1 := E s f g w (z a ,s) > E s f g w (z b ,s)v s. H g 2 := E s f g w (z a ,s)  E s f g w (z b ,s), where the expectation E s is taken with the empirical frequencies. Notably, we propose a hypothesis testing that do not calibrate the data to one batch but will find genes that are consistently differentially expressed. Evaluating which hypothesis is more probable amounts to evaluating a Bayes factor  [81]  (Bayesian generalization of the p-value) which is expressed as K = log e p(H g 1 | x a ,x b ) p(H g 2 | x a ,x b ). The sign of K indicates which of H g 1 and H g 2 is more likely.