4 2 1 a n a ly s i s Large-scale single-cell RNA sequencing (scRNA-seq) data sets that are produced in different laboratories and at different times contain batch effects that may compromise the integration and interpretation of the data. Existing scRNA-seq analysis methods incorrectly assume that the composition of cell populations is either known or identical across batches. We present a strategy for batch correction based on the detection of mutual nearest neighbors (MNNs) in the high-dimensional expression space. Our approach does not rely on predefined or equal population compositions across batches; instead, it requires only that a subset of the population be shared between batches. We demonstrate the superiority of our approach compared with existing methods by using both simulated and real scRNA-seq data sets. Using multiple droplet-based scRNA-seq data sets, we demonstrate that our MNN batch-effect-correction method can be scaled to large numbers of cells. The decreasing cost of scRNA-seq experiments [1] [2] [3] [4] has encouraged the establishment of large-scale projects such as the Human Cell Atlas, which profile the transcriptomes of thousands to millions of cells. For such large studies, logistical constraints inevitably dictate that data be generated separately i.e., at different times and with different operators. Data may also be generated in multiple laboratories using different cell dissociation and handling protocols, library-preparation technologies and/or sequencing platforms. All of these factors result in batch effects 5, 6 , in which the expression of genes in one batch differs systematically from that in another batch. Such differences can mask underlying biology or introduce spurious structure in the data; thus, to avoid misleading conclusions, they must be corrected before further analysis. Most existing methods for batch correction are based on linear regression. The limma package provides the removeBatchEffect function 7 , which fits a linear model containing a blocking term for the batch structure to the expression values for each gene. Subsequently, the coefficient for each blocking term is set to zero, and the expression values are computed from the remaining terms and residuals, thus yielding a new expression matrix without batch effects. The ComBat method 8 uses a similar strategy but performs an additional step involving empirical Bayes shrinkage of the blocking coefficient estimates. This procedure stabilizes the estimates in the presence of limited replicates by sharing information across genes. Other methods, such as RUVseq 9 and svaseq 10 , are also frequently used for batch correction, but their focus is primarily on identifying unknown factors of variation, for example, those due to unrecorded experimental differences in cell processing. After these factors are identified, their effects can be regressed out as described previously. Existing batch-correction methods were specifically designed for bulk RNA-seq. Thus, their application to scRNA-seq data is based on the assumption that the composition of the cell population within each batch is identical. Any systematic differences in mean gene expression between batches are attributed to technical differences that can be regressed out. However, in practice, the population composition is usually not identical across batches in scRNA-seq studies. Even if the same cell types are present in each batch, the abundance of each cell type in the data set can change depending upon subtle differences in procedures such as cell culture or tissue extraction, dissociation and sorting. Consequently, the estimated coefficients for the batch blocking factors are not purely technical but contain a nonzero biological component because of differences in composition. Batch correction based on these coefficients would thus yield inaccurate representations of the cellular expression profiles, and the results might potentially be worse than if no correction were performed. An alternative approach for data merging and comparison in the presence of batch effects uses a set of landmarks from a reference data set to project new data onto the reference 11, 12 . The rationale for this approach is that a given cell type in the reference batch will be most similar to cells of its own type in the new batch. Such projection strategies can be applied by using several dimensionality-reduction methods, such as principal component analysis (PCA) or diffusion maps, or by using force-based methods such as t-distributed stochastic neighbor embedding (t-SNE). This strategy depends on the selection of landmark points in high-dimensional space picked from the reference data set, which cover all cell types that might appear in the later batches. However, if the new batches include cell types that fall outside the transcriptional space explored in the reference batch, these cell types will not be projected to an appropriate position in the space defined by the landmarks (Supplementary Note 1) . 4 2 2 VOLUME 36 NUMBER 5 MAY 2018 nature biotechnology a n a ly s i s Here, we propose a new method for removal of discrepancies between biologically related batches according to the presence of MNNs between batches, which are considered to define the most similar cells of the same type across batches. The difference in expression values between cells in an MNN pair provides an estimate of the batch effect, which is made more precise by averaging across many such pairs. A correction vector is obtained from the estimated batch effect and applied to the expression values to perform batch correction. Our approach automatically identifies overlaps in population composition between batches and uses only the overlapping subsets for correction, thus avoiding the assumption of equal composition required by other methods. We demonstrate that our approach outperforms existing methods on a range of simulated and real scRNA-seq data sets involving different biological systems and technologies. RESULTS Our approach identifies cells that have mutually similar expression profiles between different experimental batches or replicates. We infer that any differences between these cells in the high-dimensional gene expression space are driven by batch effects (i.e., technical differences induced by the operator or other experimental artifacts) and do not represent the underlying biology of interest. We note that our definition of a batch effect may also incorporate some signal driven by biological features that are not of interest (for example, intersample differences due to genotype). After correction, multiple batches can be 'joined up' into a single data set (Fig. 1) . The first step of our method involves global scaling of the data through a cosine normalization. More precisely, if Y x is the expression vector for cell x, we define the cosine normalization as: Subsequently, we compute the Euclidean distance between the cosinenormalized expression profiles of pairs of cells. Calculating Euclidean distances on these normalized data is equivalent to using cosine distances on the original expression values (Supplementary Note 2). Cosine distances have been widely used for measuring cell similarities according to expression profiles 11, [13] [14] [15] and are appealing because they are scale independent 15 and thus robust to technical differences in sequencing depth and capture efficiency between batches. The next step involves identification of mutual nearest neighbors. Consider an scRNA-seq experiment consisting of two batches 1 and 2. For each cell i 1 in batch 1, we find the k cells in batch 2 with the smallest distances to i 1 , i.e., its k nearest neighbors in batch 2. We do the same for each cell in batch 2 to find its k nearest neighbors in batch 1. If a pair of cells from each batch is contained in each other's set of nearest neighbors, those cells are considered to be mutual nearest neighbors (Fig. 1) . We interpret these pairs as containing cells that belong to the same cell type or state despite being generated in different batches. Thus, any systematic differences in expression level between cells in MNN pairs should represent the batch effect. Our use of MNN pairs involves three assumptions: (i) there is at least one cell population that is present in both batches, (ii) the batch effect is almost orthogonal to the biological subspace, and (iii) the batch-effect variation is much smaller than the biological-effect variation between different cell types (more detailed discussion of these assumptions in Supplementary Note 3). The biological subspace refers to a set of basis vectors that represent biological processes; the length of each vector is equal to the number of genes. For example, some of these vectors may represent the cell cycle; some vectors may (1) (1) define expression profiles specific to each cell type; and other vectors may represent differentiation or activation states. The true expression profile of each cell can be expressed as the linear sum of these vectors. Meanwhile, the batch effect is represented by a vector of length equal to the number of genes, which is added to the expression profile for each cell in the same batch. Under our assumptions, it is straightforward to show that cells from the same population in different batches will form MNN pairs (Supplementary Note 4). This assumption can be more intuitively understood in that cells from the same population in different batches form parallel hyperplanes with respect to each other ( Fig. 1) . We also note that the orthogonality assumption is weak for a random one-dimensional batch-effect vector in high-dimensional data, especially given that local biological subspaces usually have much lower intrinsic dimensionality than the total number of genes in the data set. For each MNN pair, a pair-specific batch-correction vector is computed as the vector difference between the expression profiles of the paired cells. Although a set of biologically relevant genes (for example, highly variable genes) can facilitate identification of MNNs, the calculation of batch vectors does not need to be performed in the same space. Therefore, we can calculate the batch vectors for a different set of inquiry genes (Supplementary Note 5). A cell-specific © 2018 Nature America, Inc., part of Springer Nature. All rights reserved. nature biotechnology VOLUME 36 NUMBER 5 MAY 2018 4 2 3 a n a ly s i s batch-correction vector is then calculated as a weighted average of these pair-specific vectors, as computed with a Gaussian kernel. This approach stabilizes the correction for each cell and ensures that it changes smoothly between adjacent cells in the high-dimensional expression space. This Gaussian smoothing of batch vectors enables a locally linear batch correction; i.e., each MNN-pair batch vector contributes to the batch effect for cells in the neighborhood of the corresponding pair within each batch. Such locally linear correction of batch effects results in an overall correction that can tolerate nonconstant batch effects (Supplementary Fig. 1 ). We emphasize that this correction is performed for all cells, regardless of whether they participate in a MNN pair. Thus, correction can be performed on all cells in each batch, even if they do not have a corresponding cell type in the other batches. We generated simulated data for a simple scenario with two batches of cells, each consisting of varying proportions of three cell types (Online Methods). We applied each of three batch-correction methods-our MNN-based correction method, limma and ComBat-to the simulated data, then evaluated the results by inspecting t-SNE plots 16 (Online Methods). Proper removal of the batch effect should result in the formation of three clusters, one for each cell type, such that each cluster contains a mixture of cells from both batches. However, we observed this ideal result only after MNN correction (Fig. 2) . Expression data that were uncorrected or corrected with the other methods exhibited at least one cluster containing cells from only a single batch, thus indicating that the batch effect was not fully removed. This result is fully attributable to the differences in population composition, as discussed earlier. Repeating the simulation with identical proportions of all cell types in each batch yielded equivalent performance for all methods (Supplementary Fig. 2 ). To demonstrate the applicability of our method to real data, we considered two hematopoietic data sets generated in different laboratories through two different scRNA-seq protocols. In the first data set 12 , the authors used the SMART-seq2 protocol 17 to profile single cells from hematopoietic stem and progenitor cell populations in 12-week-old female mice. Using marker expression profiles from fluorescence-activated cell sorting (FACS), we retrospectively assigned known cell-type labels to cells (Online Methods). These labels included multipotent progenitors, lymphoid-primed multipotent progenitors, hematopoietic stem and progenitor cells, hematopoietic stem cells, common myeloid progenitors (CMPs), granulocyte-monocyte progenitors (GMPs) and megakaryocyte-erythrocyte progenitors (MEPs). In the second data set 18 , the authors used the massively parallel single-cell RNA-sequencing (MARS-seq) protocol to assess single-cell heterogeneity in myeloid progenitors from 6-to 8-week-old female mice. Again, indexed FACS was used to assign a cell-type label (MEP, GMP or CMP) to each cell. To assess performance, we performed t-SNE dimensionality reduction on the expression data for the highly variable genes, before and after correction with each of the three methods (MNN, limma and ComBat) (Fig. 3a-d and Online Methods). Only MNN correction correctly merged the cell types that were shared between batches, i.e., CMPs, MEPs and GMPs, while preserving the underlying differentiation hierarchy 12, 18 (Fig. 3e) . In contrast, the shared cell types still clustered by batch after correction with limma or ComBat, thus indicating that the batch effect had not been completely removed (coloring by batch in Supplementary Fig. 3 ). This result is attributable to the differences in cell-type composition between batches and is consistent with the simulation results. To ensure that these results were not due to an idiosyncrasy of the t-SNE method, we repeated our analysis with an alternative dimensionality-reduction approach (PCA) using only the cell types in common between the two batches ( Fig. 3f-i) . Among the methods, MNN correction was still the most effective at removing the batch effect. As a justification for the orthogonality of the batch effect to the biological hyperplane, we present a histogram of the angle between the batch vectors calculated by MNN and the first two singular value decomposition components of the reference batch used in MNN (i.e., the SMART-seq2 data set). Most angles are close to 90°, thus supporting the near-orthogonality assumption ( Supplementary  Fig. 3e) . A diffusion map 19 of the MNN-corrected data ( Supplementary  Fig. 3f-h) shows the same differentiation hierarchy of cell types as that in Figure 3e . Repeating the same analysis on a subset of randomly sampled genes (1,500 out of the total of 3,904 highly variable genes) yielded similar results, thus demonstrating the robustness of our analysis with respect to the input gene set (Supplementary Fig. 4 ). We further tested the ability of our method to combine more complex data sets generated through a variety of methods. Here, we focused on the pancreas because it is a highly heterogeneous tissue with several well-defined cell types. We combined scRNA-seq data on human pancreas cells from four different publicly available data sets 20-23 generated through two different scRNA-seq protocols (SMART-seq2 and scRNA-seq by multiplexed linear amplification (CEL-seq)/CELseq2). Cell-type labels were taken from the provided metadata or were derived according to the methodology described in the original publication (further details of data preprocessing in Online Methods). © 2018 Nature America, Inc., part of Springer Nature. All rights reserved. 4 2 4 VOLUME 36 NUMBER 5 MAY 2018 nature biotechnology a n a ly s i s We applied MNN, limma and ComBat to the combined data set and examined the corrected data. All three batch-correction methods improved the grouping of cells by their cell-type labels (Fig. 4a,b , Online Methods and Supplementary Fig. 5a-d) . This result is not surprising, because the discrepancy between cell-type composition in the four batches was modest (Supplementary Table 1) . However, even a small difference in composition was sufficient to cause ductal and acinar cells to be incorrectly separated after correction with limma or ComBat. By comparison, both cell types were coherently grouped across batches after MNN correction, in agreement with the simulation results. To determine the effect of correction on the quality of cell-type-based clustering, we assessed cluster separation by computing the average silhouette widths for each cell type ( Supplementary Fig. 5 and Online Methods). The average silhouette coefficient after MNN correction was significantly larger than those in the uncorrected and limma-and ComBat-corrected data (P < 0:05, two-sided Welch's t-test). Thus, MNN correction is able to decrease the between-batch variance within each cell type while preserving differences among cell types. We also computed the entropy of mixing (Online Methods) to quantify the extent of intermingling of cells from different batches. The data that were batch corrected with MNN showed higher entropy of mixing than did the uncorrected data and the data corrected with limma or ComBat ( Supplementary  Fig. 5 ). The improvement in the mixing of batches was observed in the reduced-dimension space obtained through either t-SNE or PCA (Supplementary Fig. 5e-l) . We again supported our assumption that batch effects are adequately removed when they lie orthogonally to the biological subspace ( Supplementary Fig. 5m-o) . The observed structure in the pancreas data was robust to the size of the input gene set, as demonstrated by random subsampling of the total highly variable gene set (Supplementary Fig. 6 ). After batch correction is performed, the corrected expression values can be used in routine downstream analyses such as clustering prior to differential gene expression identification. To provide a demonstration, we used the MNN-corrected expression matrix to simultaneously cluster cells from all four pancreas data sets. Our new cluster labels were in agreement with the previous cell-type assignments based on the individual batches, with an adjusted Rand index of 0.94 (a Rand index of 0 is equivalent to a random assignment, whereas a Rand index of 1 denotes a perfect match between previous and new assignments). Importantly, we obtained clusters for all batches in a single clustering step. This procedure ensured that the cluster labels were directly comparable between cells in different batches. In contrast, if clustering had been performed separately in each batch, there would have been no guarantee that a (weakly separated) cluster detected in one batch would have had a direct counterpart in another batch. © 2018 Nature America, Inc., part of Springer Nature. All rights reserved. nature biotechnology VOLUME 36 NUMBER 5 MAY 2018 4 2 5 a n a ly s i s We used our new clusters to perform a differential expression analysis between the delta islet cluster and the gamma islet cluster. Using cells from all batches, we detected 76 differentially expressed (DE) genes at a false discovery rate (FDR) of 5% (Fig. 4c) . This set included the marker genes for the cells included in the analysis (PPY and SST), genes involved in pancreatic islet cell development (PAX6) and genes recently implicated in delta islet function and the development of type 2 diabetes (CD9 and HADH) 22 . For comparison, we repeated the differential expression analysis by using only cells from each batch in which both cell types were present [20] [21] [22] . The results yielded only 12, 49 and 51 genes, respectively, at an FDR of 5%, which encompassed 14.5-57.9% of those detected when all cells were used (Fig. 4d) . Merging data sets is beneficial because it increases the number of cells without requiring additional experimental work, improves statistical power for downstream analyses such as differential gene expression and consequently provides additional biological insights. To this end, our MNN approach is critical because it ensures that merging is performed in a coherent manner. The advent of droplet-based cell capture, lysis, RNA reverse transcription and subsequent expression profiling by sequencing has allowed for single-cell expression experiments to be scaled up to tens and hundreds of thousands of cells 2, 3, 24 . These technologies are ideal for testing the scalability and applicability of our correction method to large scRNA-seq data sets. We specifically applied our MNN approach to two large data sets of droplet-based scRNA-seq derived from the commercial 10X Genomics Chromium platform 24 . We selected data sets comprising a mixture of cell identities and complexities, namely 68,000 peripheral blood mononuclear cells (PBMCs) and 4,000 T cells derived from different donors. PBMCs contain a milieu of peripheral adaptive and innate immune white blood cells as they circulate through the human vasculature, whereas peripheral T cells contain a mixture of naive and antigen-exposed lymphocytes involved in active immune surveillance. A naive merging of these two data sets without accounting for batch effects illustrated the separation of the T cells from their counterparts in the PBMC data (Fig. 5a,b) . The combination of these two data sets by using MNNs demonstrated that the separate peripheral T cells mapped to the T cell subsets within the PBMC mixture (Fig. 5c,d) . Importantly, other peripheral lymphocyte relationships were not distorted by the correction applied, despite the absence of MNNs in the T cell data set (Fig. 5c) . Specifically, 4, 446/4,459 (99.7%) of individual T cells mapped onto their appropriate counterparts in the PBMC data set (Fig. 5) . The remaining 13/4,459 (0.3%) mapped primarily to a small cluster of unknown ontogeny and to the edges of a large cluster of monocytes. In contrast, 14 non-T cells (0.3%; specifically monocytes) mapped to T cell clusters inappropriately. As the size of single-cell expression data sets increases, there will be a growing need for computational methods that can scale up to meet these requirements. To demonstrate the scalability of our method, we sampled different proportions of cells from the 68K PBMC data set, then corrected the batch effect between each subsample and the 4K T cell data. Within the range of 7,000 to 70,000 cells, we observed The y axis represents the −log 10 BenjaminiHochberg-adjusted P value (−log 10 P values >100 are censored at 100 for comparable scales), and the x axis is the log 2 fold change of expression in gamma cells over delta cells. Individual gene symbols are labeled when |log 2 fold change| > 3. More genes are consistently differentially expressed at an FDR of 5% in the combined data sets. (d) Venn diagrams representing the intersection of DE genes by using the cell-type labels after batch correction (blue circles) and using the original cell-type labels from each individual study (orange circles). Numbers in each segment are the total numbers of DE genes between delta and gamma islet cells in each batch. Each Venn diagram corresponds to a batch in which both cell types are present. © 2018 Nature America, Inc., part of Springer Nature. All rights reserved. 4 2 6 VOLUME 36 NUMBER 5 MAY 2018 nature biotechnology a n a ly s i s an approximately linear time increase (Fig. 5e) . Thus, our method is compatible with both the nature of droplet-technology-derived single-cell expression data and the scale of current and future data sets. Proper removal of batch effects is critical for valid data analysis and interpretation of results. This removal is especially pertinent as the scale and scope of scRNA-seq experiments increase, exceeding the capacity of data generation within a single batch. To answer the relevant biological questions, merging data from different batches-generated by different protocols, operators and/or platforms-is required. However, for biological systems that are highly heterogeneous, the composition of cell types and states is likely to differ across batches, owing to stochastic and uncontrollable biological variability. Existing batch-correction methods do not account for differences in cell composition between batches and fail to fully remove the batch effect in such cases. This failure can lead to misleading conclusions wherein batch-specific clusters are incorrectly interpreted as distinct cell types. By using both simulated data and real scRNA-seq data sets, we demonstrated that our MNN method is able to successfully remove the batch effect in the presence of differences in composition. Moreover, we demonstrated the MNN method's scalability on large droplet-based data sets. One prerequisite for our MNN method is that each batch must contain at least one shared cell population with another batch. This requirement is necessary for the correct identification of MNN pairs between batches. Batches without any shared structure are inherently difficult to correct, because the batch effects are completely confounded by biological differences. Such cases provide a motivation for using a 'cell control' , i.e., an easily reproducible cell population of known composition (from a cell line, for example) that is spiked into each sample for the purpose of removing batch effects across samples. A notable feature of our MNN correction method is that it adjusts for local variations in batch effects by using a Gaussian kernel. Our method is therefore able to accommodate differences in the size or direction of the batch effect between different cell subpopulations in high-dimensional space. Such differences are not easily handled by methods based on linear models (because they would require explicit modeling of predefined groupings of cells, which would defeat the purpose of using scRNA-seq to study population heterogeneity in the first place). Our results for the pancreas data set suggest that considering cell-type-specific batch effects (the default setting of MNN) rather than a globally constant batch effect for all cells improves batch-removal results (Supplementary Fig. 7 ). An important consequence is that a single control population might not suffice for accurate estimation of local batch effects. Instead, using an appropriately mixed population of cells to properly account for local variation may be necessary. We demonstrated in simulations and real data sets that MNN successfully combines cells with the same cell-type label, by bringing cells from different batches onto a common coordinate system that is defined by the first (reference) batch, such that all batches can be analyzed together. Therefore, MNN eliminates discrepancies between related batches without an analysis or interpretation of the origins and causes of batch effects (between each pair of batches). The study of the technical and biological origins of these discrepancies may also be interesting. For instance, one batch might contain cells from a geneknockout experiment, and the other batch might contain cells from a wild-type organism. In such cases, the correction vectors (provided as an output of the MNN algorithm) could potentially be examined to understand the differences between batches. Batch correction plays a critical role in the interpretation of scRNAseq data from both small studies, in which logistical constraints preclude the generation of data in a single batch, and large studies involving international consortia such as the Human Cell Atlas, in which scRNA-seq data are generated for a variety of related tissues at different times and by multiple laboratories. Our MNN method provides a superior alternative to existing methods for batch correction in the presence of compositional differences between batches. We anticipate that this method will improve the rigor of scRNA-seq data analysis and thus the quality of the biological conclusions.
Finally, in biomedical data analysis, mutual nearest neighbors (MNN) batch correction REF focuses on families of manifold deformations that are often encountered in biomedical data.
Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors
Nature Biotechnology
Nature Biotechnology
Biology Medicine