# Chapter 1 Correcting batch effects

## 1.1 Motivation

Large single-cell RNA sequencing (scRNA-seq) projects usually need to generate data across multiple batches due to logistical constraints. However, the processing of different batches is often subject to uncontrollable differences, e.g., changes in operator, differences in reagent quality. This results in systematic differences in the observed expression in cells from different batches, which we refer to as “batch effects”. Batch effects are problematic as they can be major drivers of heterogeneity in the data, masking the relevant biological differences and complicating interpretation of the results.

Computational removal of batch-to-batch variation allows us to combine data across multiple batches for a consolidated downstream analysis. However, existing methods based on linear models (Ritchie et al. 2015; Leek et al. 2012) assume that the composition of cell populations are either known or the same across batches. To overcome these limitations, bespoke methods have been developed for batch correction of single-cell data (Haghverdi et al. 2018; Butler et al. 2018; Lin et al. 2019) that do not require a priori knowledge about the composition of the population. This allows them to be used in workflows for exploratory analyses of scRNA-seq data where such knowledge is usually unavailable.

## 1.2 Quick start

To demonstrate, we will use two separate 10X Genomics PBMC datasets generated in two different batches. Each dataset was obtained from the TENxPBMCData package and separately subjected to basic processing steps such as quality control and normalization. As a general rule, these upstream processing steps should be done within each batch where possible. For example, outlier-based QC on the cells is more effective when performed within a batch (Advanced Section 1.4), and we can more effectively model the mean-variance relationship on each batch separately (Basic Section 3.4).

#--- loading ---#
library(TENxPBMCData)
all.sce <- list(
pbmc3k=TENxPBMCData('pbmc3k'),
pbmc4k=TENxPBMCData('pbmc4k'),
pbmc8k=TENxPBMCData('pbmc8k')
)

#--- quality-control ---#
library(scater)
stats <- high.mito <- list()
for (n in names(all.sce)) {
current <- all.sce[[n]]
is.mito <- grep("MT", rowData(current)$Symbol_TENx) stats[[n]] <- perCellQCMetrics(current, subsets=list(Mito=is.mito)) high.mito[[n]] <- isOutlier(stats[[n]]$subsets_Mito_percent, type="higher")
all.sce[[n]] <- current[,!high.mito[[n]]]
}

#--- normalization ---#
all.sce <- lapply(all.sce, logNormCounts)

#--- variance-modelling ---#
library(scran)
all.dec <- lapply(all.sce, modelGeneVar)
all.hvgs <- lapply(all.dec, getTopHVGs, prop=0.1)

#--- dimensionality-reduction ---#
library(BiocSingular)
set.seed(10000)
all.sce <- mapply(FUN=runPCA, x=all.sce, subset_row=all.hvgs,
MoreArgs=list(ncomponents=25, BSPARAM=RandomParam()),
SIMPLIFY=FALSE)

set.seed(100000)
all.sce <- lapply(all.sce, runTSNE, dimred="PCA")

set.seed(1000000)
all.sce <- lapply(all.sce, runUMAP, dimred="PCA")

#--- clustering ---#
for (n in names(all.sce)) {
g <- buildSNNGraph(all.sce[[n]], k=10, use.dimred='PCA')
clust <- igraph::cluster_walktrap(g)$membership colLabels(all.sce[[n]]) <- factor(clust) } pbmc3k <- all.sce$pbmc3k
dec3k <- all.dec$pbmc3k pbmc3k ## class: SingleCellExperiment ## dim: 32738 2609 ## metadata(0): ## assays(2): counts logcounts ## rownames(32738): ENSG00000243485 ENSG00000237613 ... ENSG00000215616 ## ENSG00000215611 ## rowData names(3): ENSEMBL_ID Symbol_TENx Symbol ## colnames: NULL ## colData names(13): Sample Barcode ... sizeFactor label ## reducedDimNames(3): PCA TSNE UMAP ## mainExpName: NULL ## altExpNames(0): pbmc4k <- all.sce$pbmc4k
dec4k <- all.dec$pbmc4k pbmc4k ## class: SingleCellExperiment ## dim: 33694 4182 ## metadata(0): ## assays(2): counts logcounts ## rownames(33694): ENSG00000243485 ENSG00000237613 ... ENSG00000277475 ## ENSG00000268674 ## rowData names(3): ENSEMBL_ID Symbol_TENx Symbol ## colnames: NULL ## colData names(13): Sample Barcode ... sizeFactor label ## reducedDimNames(3): PCA TSNE UMAP ## mainExpName: NULL ## altExpNames(0): We then use the quickCorrect() function from the batchelor package to compute corrected values across the two objects. This performs all the steps to set up the data for correction (Section 1.3), followed by MNN correction to actually perform the correction itself (Section 1.6). Alternatively, we could use one of the other correction algorithms described in this chapter by modifying PARAM= appropriately. library(batchelor) quick.corrected <- quickCorrect(pbmc3k, pbmc4k, precomputed=list(dec3k, dec4k), PARAM=FastMnnParam(BSPARAM=BiocSingular::RandomParam())) quick.sce <- quick.corrected$corrected
quick.sce
## class: SingleCellExperiment
## dim: 31232 6791
## assays(1): reconstructed
## rownames(31232): ENSG00000243485 ENSG00000237613 ... ENSG00000198695
##   ENSG00000198727
## rowData names(1): rotation
## colnames: NULL
## colData names(1): batch
## reducedDimNames(1): corrected
## mainExpName: NULL
## altExpNames(0):

This yields low-dimensional corrected values for use in downstream analyses (Figure 1.1).

library(scater)
set.seed(00101010)
quick.sce <- runTSNE(quick.sce, dimred="corrected")
quick.sce$batch <- factor(quick.sce$batch)
plotTSNE(quick.sce, colour_by="batch")

## 1.3 Explaining the data preparation

The quickCorrect() function wraps a number of steps that are required to prepare the data for batch correction. The first and most obvious is to subset all batches to the common “universe” of features. In this case, it is straightforward as both batches use Ensembl gene annotation; more difficult integrations will require some mapping of identifiers using packages like org.Mm.eg.db.

universe <- intersect(rownames(pbmc3k), rownames(pbmc4k))
length(universe)
## [1] 31232
# Subsetting the SingleCellExperiment object.
pbmc3k <- pbmc3k[universe,]
pbmc4k <- pbmc4k[universe,]

# Also subsetting the variance modelling results, for convenience.
dec3k <- dec3k[universe,]
dec4k <- dec4k[universe,]

The second step is to rescale each batch to adjust for differences in sequencing depth between batches. The multiBatchNorm() function recomputes log-normalized expression values after adjusting the size factors for systematic differences in coverage between SingleCellExperiment objects. (Size factors only remove biases between cells within a single batch.) This improves the quality of the correction by removing one aspect of the technical differences between batches.

rescaled <- multiBatchNorm(pbmc3k, pbmc4k)
pbmc3k <- rescaled[[1]]
pbmc4k <- rescaled[[2]]

Finally, we perform feature selection by averaging the variance components across all batches with the combineVar() function. We compute the average as it is responsive to batch-specific HVGs while still preserving the within-batch ranking of genes. This allows us to use the same strategies described in Basic Section 3.5 to select genes of interest. In contrast, approaches based on taking the intersection or union of HVGs across batches become increasingly conservative or liberal, respectively, with an increasing number of batches.

library(scran)
combined.dec <- combineVar(dec3k, dec4k)
chosen.hvgs <- combined.dec$bio > 0 sum(chosen.hvgs) ## [1] 13431 When integrating datasets of variable composition, it is generally safer to err on the side of including more HVGs than are used in a single dataset analysis, to ensure that markers are retained for any dataset-specific subpopulations that might be present. For a top $$X$$ selection, this means using a larger $$X$$ (e.g., quickCorrect() defaults to 5000), or in this case, we simply take all genes above the trend. That said, many of the signal-to-noise considerations described in Basic Section 3.5 still apply here, so some experimentation may be necessary for best results. ## 1.4 No correction Before we actually perform any correction, it is worth examining whether there is any batch effect in this dataset. We combine the two SingleCellExperiments and perform a PCA on the log-expression values for our selected subset of HVGs. In this example, our datasets are file-backed and so we instruct runPCA() to use randomized PCA for greater efficiency - see Advanced Section 14.2.2 for more details - though the default IRLBA will suffice for more common in-memory representations. # Synchronizing the metadata for cbind()ing. # TODO: replace with combineCols when that comes out. rowData(pbmc3k) <- rowData(pbmc4k) pbmc3k$batch <- "3k"
pbmc4k$batch <- "4k" uncorrected <- cbind(pbmc3k, pbmc4k) # Using RandomParam() as it is more efficient for file-backed matrices. library(scater) set.seed(0010101010) uncorrected <- runPCA(uncorrected, subset_row=chosen.hvgs, BSPARAM=BiocSingular::RandomParam()) We use graph-based clustering on the components to obtain a summary of the population structure. As our two PBMC populations should be replicates, each cluster should ideally consist of cells from both batches. However, we instead see clusters that are comprised of cells from a single batch. This indicates that cells of the same type are artificially separated due to technical differences between batches. library(scran) snn.gr <- buildSNNGraph(uncorrected, use.dimred="PCA") clusters <- igraph::cluster_walktrap(snn.gr)$membership
tab.resid <- table(Cluster=clusters.resid, Batch=residuals$batch) tab.resid ## Batch ## Cluster 1 2 ## 1 478 2 ## 2 142 179 ## 3 22 41 ## 4 298 566 ## 5 340 606 ## 6 0 138 ## 7 404 376 ## 8 145 91 ## 9 2 636 ## 10 22 73 ## 11 6 51 ## 12 629 1110 ## 13 3 36 ## 14 91 211 ## 15 12 55 ## 16 4 8 ## 17 11 3 residuals <- runTSNE(residuals, dimred="corrected") residuals$batch <- factor(residuals$batch) plotTSNE(residuals, colour_by="batch") ## 1.6 MNN correction Consider a cell $$a$$ in batch $$A$$, and identify the cells in batch $$B$$ that are nearest neighbors to $$a$$ in the expression space defined by the selected features. Repeat this for a cell $$b$$ in batch $$B$$, identifying its nearest neighbors in $$A$$. Mutual nearest neighbors are pairs of cells from different batches that belong in each other’s set of nearest neighbors. The reasoning is that MNN pairs represent cells from the same biological state prior to the application of a batch effect - see Haghverdi et al. (2018) for full theoretical details. Thus, the difference between cells in MNN pairs can be used as an estimate of the batch effect, the subtraction of which yields batch-corrected values. Compared to linear regression, MNN correction does not assume that the population composition is the same or known beforehand. This is because it learns the shared population structure via identification of MNN pairs and uses this information to obtain an appropriate estimate of the batch effect. Instead, the key assumption of MNN-based approaches is that the batch effect is orthogonal to the biology in high-dimensional expression space. Violations reduce the effectiveness and accuracy of the correction, with the most common case arising from variations in the direction of the batch effect between clusters. Nonetheless, the assumption is usually reasonable as a random vector is very likely to be orthogonal in high-dimensional space. The batchelor package provides an implementation of the MNN approach via the fastMNN() function. (Unlike the MNN method originally described by Haghverdi et al. (2018), the fastMNN() function performs PCA to reduce the dimensions beforehand and speed up the downstream neighbor detection steps.) We apply it to our two PBMC batches to remove the batch effect across the highly variable genes in chosen.hvgs. To reduce computational work and technical noise, all cells in all batches are projected into the low-dimensional space defined by the top d principal components. Identification of MNNs and calculation of correction vectors are then performed in this low-dimensional space. # Again, using randomized SVD here, as this is faster than IRLBA for # file-backed matrices. We set deferred=TRUE for greater speed. set.seed(1000101001) mnn.out <- fastMNN(pbmc3k, pbmc4k, d=50, k=20, subset.row=chosen.hvgs, BSPARAM=BiocSingular::RandomParam(deferred=TRUE)) mnn.out ## class: SingleCellExperiment ## dim: 13431 6791 ## metadata(2): merge.info pca.info ## assays(1): reconstructed ## rownames(13431): ENSG00000239945 ENSG00000228463 ... ENSG00000198695 ## ENSG00000198727 ## rowData names(1): rotation ## colnames: NULL ## colData names(1): batch ## reducedDimNames(1): corrected ## mainExpName: NULL ## altExpNames(0): The function returns a SingleCellExperiment object containing corrected values for downstream analyses like clustering or visualization. Each column of mnn.out corresponds to a cell in one of the batches, while each row corresponds to an input gene in chosen.hvgs. The batch field in the column metadata contains a vector specifying the batch of origin of each cell. head(mnn.out$batch) 
## [1] 1 1 1 1 1 1

The corrected matrix in the reducedDims() contains the low-dimensional corrected coordinates for all cells, which we will use in place of the PCs in our downstream analyses.

dim(reducedDim(mnn.out, "corrected"))
## [1] 6791   50

A reconstructed matrix in the assays() contains the corrected expression values for each gene in each cell, obtained by projecting the low-dimensional coordinates in corrected back into gene expression space. We do not recommend using this for anything other than visualization (Chapter 3).

assay(mnn.out, "reconstructed")
## <13431 x 6791> matrix of class LowRankMatrix and type "double":
##                       [,1]       [,2]       [,3] ...    [,6790]    [,6791]
## ENSG00000239945 -2.522e-06 -1.851e-06 -1.199e-05   .  1.832e-06 -3.641e-06
## ENSG00000228463 -6.627e-04 -6.724e-04 -4.820e-04   . -8.531e-04 -3.999e-04
## ENSG00000237094 -8.077e-05 -8.038e-05 -9.631e-05   .  7.261e-06 -4.094e-05
## ENSG00000229905  3.838e-06  6.180e-06  5.432e-06   .  8.534e-06  3.485e-06
## ENSG00000237491 -4.527e-04 -3.178e-04 -1.510e-04   . -3.491e-04 -2.082e-04
##             ...          .          .          .   .          .          .
## ENSG00000198840 -0.0296508 -0.0340101 -0.0502385   . -0.0362884 -0.0183084
## ENSG00000212907 -0.0041681 -0.0056570 -0.0106420   . -0.0083837  0.0005996
## ENSG00000198886  0.0145358  0.0200517 -0.0307131   . -0.0109254 -0.0070064
## ENSG00000198695  0.0014427  0.0013490  0.0001493   . -0.0009826 -0.0022712
## ENSG00000198727  0.0152570  0.0106167 -0.0256450   . -0.0227962 -0.0022898

The most relevant parameter for tuning fastMNN() is k, which specifies the number of nearest neighbors to consider when defining MNN pairs. This can be interpreted as the minimum anticipated frequency of any shared cell type or state in each batch. Increasing k will generally result in more aggressive merging as the algorithm is more generous in matching subpopulations across batches. It can occasionally be desirable to increase k if one clearly sees that the same cell types are not being adequately merged across batches.

We cluster on the low-dimensional corrected coordinates to obtain a partitioning of the cells that serves as a proxy for the population structure. If the batch effect is successfully corrected, clusters corresponding to shared cell types or states should contain cells from multiple batches. We see that all clusters contain contributions from each batch after correction, consistent with our expectation that the two batches are replicates of each other.

library(scran)
snn.gr <- buildSNNGraph(mnn.out, use.dimred="corrected")
clusters.mnn <- igraph::cluster_walktrap(snn.gr)$membership tab.mnn <- table(Cluster=clusters.mnn, Batch=mnn.out$batch)
tab.mnn
##        Batch
## Cluster    1    2
##      1   337  606
##      2   289  542
##      3   152  181
##      4    12    4
##      5   517  467
##      6    17   19
##      7   313  661
##      8   162  118
##      9    11   56
##      10  547 1083
##      11   17   59
##      12   16   58
##      13  144   93
##      14   67  191
##      15    4   36
##      16    4    8

We can also visualize the corrected coordinates using a $$t$$-SNE plot (Figure 1.5). The presence of visual clusters containing cells from both batches provides a comforting illusion that the correction was successful.

library(scater)
set.seed(0010101010)
mnn.out <- runTSNE(mnn.out, dimred="corrected")

mnn.out$batch <- factor(mnn.out$batch)
plotTSNE(mnn.out, colour_by="batch")

See also Chapter 8 for a case study using MNN correction on a series of human pancreas datasets.

## 1.7 Further options

All of the batchelor functions can operate on a single SingleCellExperiment containing data from all batches. For example, if we were to recycle the uncorrected object from Section 1.4, we could apply MNN correction without splitting the object into multiple parts.

set.seed(10000)
single.correct <- fastMNN(uncorrected, batch=uncorrected$batch, subset.row=chosen.hvgs, BSPARAM=BiocSingular::RandomParam()) single.correct ## class: SingleCellExperiment ## dim: 13431 6791 ## metadata(2): merge.info pca.info ## assays(1): reconstructed ## rownames(13431): ENSG00000239945 ENSG00000228463 ... ENSG00000198695 ## ENSG00000198727 ## rowData names(1): rotation ## colnames: NULL ## colData names(1): batch ## reducedDimNames(1): corrected ## mainExpName: NULL ## altExpNames(0): It is similarly straightforward to simultaneously perform correction across >2 batches, either by having multiple levels in batch= or by providing more SingleCellExperiment objects (or even raw matrices of expression values). This is demonstrated below for MNN correction with an additional PBMC dataset (Figure 1.6). pbmc8k <- all.sce$pbmc8k
dec8k <- all.dec$pbmc8k quick.corrected2 <- quickCorrect(3k=pbmc3k, 4k=pbmc4k, 8k=pbmc8k, precomputed=list(dec3k, dec4k, dec8k), PARAM=FastMnnParam(BSPARAM=BiocSingular::RandomParam(), auto.merge=TRUE)) quick.sce2 <- quick.corrected2$corrected

set.seed(00101010)
quick.sce2 <- runTSNE(quick.sce2, dimred="corrected")
plotTSNE(quick.sce2, colour_by="batch")

In the specific case of MNN correction, we can also set auto.merge=TRUE to allow it to choose the “best” order in which to perform the merges. This is slower but can occasionally be useful when the batches involved have very different cell type compositions. For example, if one batch contained only B cells, another batch contained only T cells and a third batch contained B and T cells, it would be unwise to try to merge the first two batches together as the wrong MNN pairs would be identified. With auto.merge=TRUE, the function would automatically recognize that the third batch should be used as the reference to which the others should be merged.

## Session Info

R version 4.2.0 (2022-04-22)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.4 LTS

Matrix products: default
BLAS:   /home/biocbuild/bbs-3.15-bioc/R/lib/libRblas.so
LAPACK: /home/biocbuild/bbs-3.15-bioc/R/lib/libRlapack.so

locale:
[1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
[3] LC_TIME=en_GB              LC_COLLATE=C
[5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8       LC_NAME=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods
[8] base

other attached packages:
[1] scran_1.24.0                scater_1.24.0
[3] ggplot2_3.3.6               scuttle_1.6.2
[5] batchelor_1.12.0            SingleCellExperiment_1.18.0
[7] SummarizedExperiment_1.26.1 Biobase_2.56.0
[9] GenomicRanges_1.48.0        GenomeInfoDb_1.32.2
[11] HDF5Array_1.24.0            rhdf5_2.40.0
[13] DelayedArray_0.22.0         IRanges_2.30.0
[15] S4Vectors_0.34.0            MatrixGenerics_1.8.0
[17] matrixStats_0.62.0          BiocGenerics_0.42.0
[19] Matrix_1.4-1                BiocStyle_2.24.0
[21] rebook_1.6.0

loaded via a namespace (and not attached):
[1] bitops_1.0-7              filelock_1.0.2
[3] tools_4.2.0               bslib_0.3.1
[5] utf8_1.2.2                R6_2.5.1
[7] irlba_2.3.5               ResidualMatrix_1.6.0
[9] vipor_0.4.5               DBI_1.1.2
[11] colorspace_2.0-3          rhdf5filters_1.8.0
[13] withr_2.5.0               gridExtra_2.3
[15] tidyselect_1.1.2          compiler_4.2.0
[17] graph_1.74.0              cli_3.3.0
[19] BiocNeighbors_1.14.0      labeling_0.4.2
[21] bookdown_0.26             sass_0.4.1
[23] scales_1.2.0              stringr_1.4.0
[25] digest_0.6.29             rmarkdown_2.14
[27] XVector_0.36.0            pkgconfig_2.0.3
[29] htmltools_0.5.2           sparseMatrixStats_1.8.0
[31] highr_0.9                 fastmap_1.1.0
[33] limma_3.52.1              rlang_1.0.2
[35] DelayedMatrixStats_1.18.0 farver_2.1.0
[37] generics_0.1.2            jquerylib_0.1.4
[39] jsonlite_1.8.0            BiocParallel_1.30.2
[41] dplyr_1.0.9               RCurl_1.98-1.6
[43] magrittr_2.0.3            BiocSingular_1.12.0
[45] GenomeInfoDbData_1.2.8    ggbeeswarm_0.6.0
[47] Rcpp_1.0.8.3              munsell_0.5.0
[49] Rhdf5lib_1.18.2           fansi_1.0.3
[51] viridis_0.6.2             lifecycle_1.0.1
[53] stringi_1.7.6             yaml_2.3.5
[55] edgeR_3.38.1              zlibbioc_1.42.0
[57] Rtsne_0.16                grid_4.2.0
[59] ggrepel_0.9.1             parallel_4.2.0
[61] dqrng_0.3.0               crayon_1.5.1
[63] dir.expiry_1.4.0          lattice_0.20-45
[65] cowplot_1.1.1             beachmat_2.12.0
[67] locfit_1.5-9.5            CodeDepends_0.6.5
[69] metapod_1.4.0             knitr_1.39
[71] pillar_1.7.0              igraph_1.3.1
[73] codetools_0.2-18          ScaledMatrix_1.4.0
[75] XML_3.99-0.9              glue_1.6.2
[77] evaluate_0.15             BiocManager_1.30.18
[79] vctrs_0.4.1               purrr_0.3.4
[81] gtable_0.3.0              assertthat_0.2.1
[83] xfun_0.31                 rsvd_1.0.5
[85] viridisLite_0.4.0         tibble_3.1.7
[87] beeswarm_0.4.0            cluster_2.1.3
[89] bluster_1.6.0             statmod_1.4.36
[91] ellipsis_0.3.2           

### References

Butler, A., P. Hoffman, P. Smibert, E. Papalexi, and R. Satija. 2018. “Integrating single-cell transcriptomic data across different conditions, technologies, and species.” Nat. Biotechnol. 36 (5): 411–20.

Haghverdi, L., A. T. L. Lun, M. D. Morgan, and J. C. Marioni. 2018. “Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors.” Nat. Biotechnol. 36 (5): 421–27.

Leek, J. T., W. E. Johnson, H. S. Parker, A. E. Jaffe, and J. D. Storey. 2012. “The sva package for removing batch effects and other unwanted variation in high-throughput experiments.” Bioinformatics 28 (6): 882–83.

Lin, Y., S. Ghazanfar, K. Y. X. Wang, J. A. Gagnon-Bartsch, K. K. Lo, X. Su, Z. G. Han, et al. 2019. “scMerge leverages factor analysis, stable expression, and pseudoreplication to merge multiple single-cell RNA-seq datasets.” Proc. Natl. Acad. Sci. U.S.A. 116 (20): 9775–84.

Ritchie, M. E., B. Phipson, D. Wu, Y. Hu, C. W. Law, W. Shi, and G. K. Smyth. 2015. “limma powers differential expression analyses for RNA-sequencing and microarray studies.” Nucleic Acids Res. 43 (7): e47.