---
title: Correcting batch effects in single-cell RNA-seq data
author: 
- name: Aaron T. L. Lun
  affiliation: Cancer Research UK Cambridge Institute, Li Ka Shing Centre, Robinson Way, Cambridge CB2 0RE, United Kingdom
- name: Michael D. Morgan
  affiliation: Wellcome Trust Sanger Institute, Wellcome Genome Campus, Hinxton, Cambridge CB10 1SA, United Kingdom
date: "`r Sys.Date()`"
vignette: >
  %\VignetteIndexEntry{05. Correcting batch effects}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}    
output: 
  BiocStyle::html_document:
    titlecaps: false
    toc_float: true
bibliography: ref.bib
---

```{r style, echo=FALSE, results='hide', message=FALSE, cache=FALSE}
library(BiocStyle)
library(knitr)
opts_chunk$set(error=FALSE, message=FALSE, warning=FALSE, cache=TRUE)
opts_chunk$set(fig.asp=1)
```

# Introduction

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 correction of these effects is critical for eliminating batch-to-batch variation, allowing data across multiple batches to be combined for valid downstream analysis.
However, existing methods such as `removeBatchEffect()` [@ritchie2015limma] assume that the composition of cell populations are either known or the same across batches.
This workflow describes the application of an alternative strategy for batch correction based on the detection of mutual nearest neighbours (MNNs) [@haghverdi2018batch].
The MNN approach does not rely on pre-defined or equal population compositions across batches, only requiring that a subset of the population be shared between batches.
We demonstrate its use on two human pancreas scRNA-seq datasets generated in separate studies.

# Processing the different datasets

## CEL-seq, GSE81076

### Loading in the data

This dataset was generated by @grun2016denovo using the CEL-seq protocol with unique molecular identifiers (UMIs) and ERCC spike-ins.
Count tables were obtained from the NCBI Gene Expression Omnibus using the accession number above.

```{r}
library(BiocFileCache)
bfc <- BiocFileCache("raw_data", ask = FALSE)
grun.fname <- bfcrpath(bfc, file.path("ftp://ftp.ncbi.nlm.nih.gov/geo/series",
    "GSE81nnn/GSE81076/suppl/GSE81076%5FD2%5F3%5F7%5F10%5F17%2Etxt%2Egz"))
```

We first read the table into memory.

```{r}
gse81076.df <- read.table(grun.fname, sep='\t', 
    header=TRUE, stringsAsFactors=FALSE, row.names=1)
dim(gse81076.df)
head(rownames(gse81076.df))
head(colnames(gse81076.df))
```

Unfortunately, the data and metadata are all mixed together in this file.
As a result, we need to manually extract the metadata from the column names.

```{r}
donor.names <- sub("^(D[0-9]+).*", "\\1", colnames(gse81076.df))
table(donor.names)
plate.id <- sub("^D[0-9]+(.*)_.*", "\\1", colnames(gse81076.df))
table(plate.id)
```

Another irritating feature of this dataset is that gene symbols were supplied, rather than stable identifiers such as Ensembl.
Stable identifiers are desirable as they facilitate reliable cross-referencing of gene identities between data sets^[The two data sets used in this workflow come from the same provider, who at least uses gene symbols consistently. So, technically, we could have skipped the conversion here. Nonetheless, we have chosen to perform it to demonstrate how one would deal with more heterogeneous data sources (e.g., mixtures of Ensembl identifiers and gene symbols, or multiple synonymous gene symbols).].
We convert all row names to Ensembl identifiers, removing `NA` or duplicated entries (with the exception of spike-in transcripts). 

```{r}
gene.symb <- gsub("__chr.*$", "", rownames(gse81076.df))
is.spike <- grepl("^ERCC-", gene.symb)
table(is.spike)

library(org.Hs.eg.db)
gene.ids <- mapIds(org.Hs.eg.db, keys=gene.symb, keytype="SYMBOL", column="ENSEMBL")
gene.ids[is.spike] <- gene.symb[is.spike]

keep <- !is.na(gene.ids) & !duplicated(gene.ids)
gse81076.df <- gse81076.df[keep,]
rownames(gse81076.df) <- gene.ids[keep]
summary(keep)
```

We create a `SingleCellExperiment` object to store the counts and metadata together.
This reduces the risk of book-keeping errors in later steps of the analysis.
Note that we re-identify the spike-in rows, as the previous indices would have changed after the subsetting.

```{r}
library(SingleCellExperiment)
sce.gse81076 <- SingleCellExperiment(list(counts=as.matrix(gse81076.df)),
	colData=DataFrame(Donor=donor.names, Plate=plate.id),
	rowData=DataFrame(Symbol=gene.symb[keep]))
isSpike(sce.gse81076, "ERCC") <- grepl("^ERCC-", rownames(gse81076.df)) 
sce.gse81076  
```

### Quality control and normalization

We compute quality control (QC) metrics for each cell [@mccarthy2017scater] and identify cells with low library sizes, low numbers of expressed genes, or high ERCC content.

```{r}
library(scater)
sce.gse81076 <- calculateQCMetrics(sce.gse81076, compact=TRUE)
QC <- sce.gse81076$scater_qc
low.lib <- isOutlier(QC$all$log10_total_counts, type="lower", nmad=3)
low.genes <- isOutlier(QC$all$log10_total_features_by_counts, type="lower", nmad=3)
high.spike <- isOutlier(QC$feature_control_ERCC$pct_counts, type="higher", nmad=3)
data.frame(LowLib=sum(low.lib), LowNgenes=sum(low.genes), 
	HighSpike=sum(high.spike, na.rm=TRUE))
```

Cells with extreme values for these QC metrics are presumed to be of low quality and are removed.
A more thorough analysis would examine the distributions of these QC metrics beforehand, but we will skip that step for brevity here.

```{r}
discard <- low.lib | low.genes | high.spike
sce.gse81076 <- sce.gse81076[,!discard]
summary(discard)
```

We compute size factors for the endogenous genes using the deconvolution method [@lun2016pooling].
This is done with pre-clustering by `quickCluster()` to avoid pooling together very different cells.

```{r}
library(scran)
library(BiocSingular)
set.seed(1000) # for irlba. 
clusters <- quickCluster(sce.gse81076, use.ranks=FALSE, BSPARAM=IrlbaParam())
table(clusters)
sce.gse81076 <- computeSumFactors(sce.gse81076, min.mean=0.1, clusters=clusters)
summary(sizeFactors(sce.gse81076))
```

We also compute size factors for the spike-in transcripts [@lun2017assessing].
Recall that we set `general.use=FALSE` to ensure that the spike-in size factors are only applied to the spike-in transcripts.

```{r}
sce.gse81076 <- computeSpikeFactors(sce.gse81076, general.use=FALSE)
summary(sizeFactors(sce.gse81076, "ERCC"))
```

We then compute normalized log-expression values for use in downstream analyses.

```{r}
sce.gse81076 <- normalize(sce.gse81076)
```

### Identifying highly variable genes

We identify highly variable genes (HVGs) using `trendVar()` and `decomposeVar()`, using the variances of spike-in transcripts to model technical noise.
We set `block=` to ensure that uninteresting differences between plates or donors do not inflate the variance.
The small discrepancy in the fitted trend in Figure \@ref(fig:var-gse81076) is caused by the fact that the trend is fitted robustly to the block-wise variances of the spike-ins, 
while the variances shown are averaged across blocks and not robust to outliers.

```{r var-gse81076, fig.cap="Variance of normalized log-expression values for each gene in the GSE81076 dataset, plotted against the mean log-expression. The blue line represents the mean-dependent trend fitted to the variances of the spike-in transcripts (red)."}
block <- paste0(sce.gse81076$Plate, "_", sce.gse81076$Donor)
fit <- trendVar(sce.gse81076, block=block, parametric=TRUE) 
dec <- decomposeVar(sce.gse81076, fit)

plot(dec$mean, dec$total, xlab="Mean log-expression", 
	ylab="Variance of log-expression", pch=16)
is.spike <- isSpike(sce.gse81076)
points(dec$mean[is.spike], dec$total[is.spike], col="red", pch=16)
curve(fit$trend(x), col="dodgerblue", add=TRUE)
```

We order genes by decreasing biological component, revealing some usual suspects such as insulin and glucagon.
We will be using this information later when performing feature selection prior to running `mnnCorrect()`.

```{r}
dec.gse81076 <- dec
dec.gse81076$Symbol <- rowData(sce.gse81076)$Symbol
dec.gse81076 <- dec.gse81076[order(dec.gse81076$bio, decreasing=TRUE),]
head(dec.gse81076)
```

```{r, echo=FALSE, results="hide"}
rm(gse81076.df)
gc()
```

## CEL-seq2, GSE85241

### Loading in the data

This dataset was generated by @muraro2016singlecell using the CEL-seq2 protocol with unique molecular identifiers (UMIs) and ERCC spike-ins.
Count tables were obtained from the NCBI Gene Expression Omnibus using the accession number above.

```{r}
muraro.fname <- bfcrpath(bfc, file.path("ftp://ftp.ncbi.nlm.nih.gov/geo/series",
    "GSE85nnn/GSE85241/suppl",
    "GSE85241%5Fcellsystems%5Fdataset%5F4donors%5Fupdated%2Ecsv%2Egz"))
```

We first read the table into memory.

```{r}
gse85241.df <- read.table(muraro.fname, sep='\t', 
    header=TRUE, row.names=1, stringsAsFactors=FALSE)
dim(gse85241.df)
head(rownames(gse85241.df))
head(colnames(gse85241.df))
```

We extract the metadata from the column names.

```{r}
donor.names <- sub("^(D[0-9]+).*", "\\1", colnames(gse85241.df))
table(donor.names)
plate.id <- sub("^D[0-9]+\\.([0-9]+)_.*", "\\1", colnames(gse85241.df))
table(plate.id)
```

Yet again, gene symbols were supplied instead of Ensembl or Entrez identifiers.
We convert all row names to Ensembl identifiers, removing `NA` or duplicated entries (with the exception of spike-in transcripts). 

```{r}
gene.symb <- gsub("__chr.*$", "", rownames(gse85241.df))
is.spike <- grepl("^ERCC-", gene.symb)
table(is.spike)

library(org.Hs.eg.db)
gene.ids <- mapIds(org.Hs.eg.db, keys=gene.symb, keytype="SYMBOL", column="ENSEMBL")
gene.ids[is.spike] <- gene.symb[is.spike]

keep <- !is.na(gene.ids) & !duplicated(gene.ids)
gse85241.df <- gse85241.df[keep,]
rownames(gse85241.df) <- gene.ids[keep]
summary(keep)
```

We create a `SingleCellExperiment` object to store the counts and metadata together.

```{r}
sce.gse85241 <- SingleCellExperiment(list(counts=as.matrix(gse85241.df)),
	colData=DataFrame(Donor=donor.names, Plate=plate.id),
	rowData=DataFrame(Symbol=gene.symb[keep]))
isSpike(sce.gse85241, "ERCC") <- grepl("^ERCC-", rownames(gse85241.df)) 
sce.gse85241  
```

### Quality control and normalization

We compute QC metrics for each cell and identify cells with low library sizes, low numbers of expressed genes, or high ERCC content.

```{r}
sce.gse85241 <- calculateQCMetrics(sce.gse85241, compact=TRUE)
QC <- sce.gse85241$scater_qc
low.lib <- isOutlier(QC$all$log10_total_counts, type="lower", nmad=3)
low.genes <- isOutlier(QC$all$log10_total_features_by_counts, type="lower", nmad=3)
high.spike <- isOutlier(QC$feature_control_ERCC$pct_counts, type="higher", nmad=3)
data.frame(LowLib=sum(low.lib), LowNgenes=sum(low.genes), 
	HighSpike=sum(high.spike, na.rm=TRUE))
```

Low-quality cells are defined as those with extreme values for these QC metrics and are removed.

```{r}
discard <- low.lib | low.genes | high.spike
sce.gse85241 <- sce.gse85241[,!discard]
summary(discard)
```

We compute size factors for the endogenous genes and spike-in transcripts, and use them to compute log-normalized expression values.

```{r}
set.seed(1000)
clusters <- quickCluster(sce.gse85241, use.ranks=FALSE, BSPARAM=IrlbaParam())
table(clusters)
sce.gse85241 <- computeSumFactors(sce.gse85241, min.mean=0.1, clusters=clusters)
summary(sizeFactors(sce.gse85241))
sce.gse85241 <- computeSpikeFactors(sce.gse85241, general.use=FALSE)
summary(sizeFactors(sce.gse85241, "ERCC"))
sce.gse85241 <- normalize(sce.gse85241)
```

### Identifying highly variable genes

We fit a trend to the spike-in variances as previously described, allowing us to model the technical noise for each gene (Figure \@ref(fig:var-gse85241)).
Again, we set `block=` to ensure that uninteresting differences between plates or donors do not inflate the variance.

```{r var-gse85241, fig.cap="Variance of normalized log-expression values for each gene in the GSE85241 dataset, plotted against the mean log-expression. The blue line represents the mean-dependent trend fitted to the variances of the spike-in transcripts (red)."}
block <- paste0(sce.gse85241$Plate, "_", sce.gse85241$Donor)
fit <- trendVar(sce.gse85241, block=block, parametric=TRUE) 
dec <- decomposeVar(sce.gse85241, fit)
plot(dec$mean, dec$total, xlab="Mean log-expression", 
	ylab="Variance of log-expression", pch=16)
is.spike <- isSpike(sce.gse85241)
points(dec$mean[is.spike], dec$total[is.spike], col="red", pch=16)
curve(fit$trend(x), col="dodgerblue", add=TRUE)
```

We order genes by decreasing biological component, as described above.

```{r}
dec.gse85241 <- dec
dec.gse85241$Symbol <- rowData(sce.gse85241)$Symbol
dec.gse85241 <- dec.gse85241[order(dec.gse85241$bio, decreasing=TRUE),]
head(dec.gse85241)
```

```{r, echo=FALSE, results="hide"}
rm(gse85241.df)
gc()
```

## Comments on additional batches

In @haghverdi2018batch, we originally performed batch correction across four separate pancreas scRNA-seq datasets.
For simplicity, we will only consider the two batches generated using CEL-seq(2) and ignore those generated using Smart-seq2 [@segerstolpe2016singlecell;@lawlor2017singlecell].
As one might expect, batch correction is easiest when dealing with data generated from the same technology,
as fewer systematic differences are present that can interfere with the biological structure.
Nonetheless, it is often possible to obtain good results when applying MNN correction to batches of data generated with different technologies.

It is also worth pointing out that both of the CEL-seq(2) batches above contain cells from multiple donors.
Each donor could be treated as a separate batch in their own right, 
reflecting (presumably uninteresting) biological differences between donors due to genotype, age, sex or other factors that are not easily controlled when dealing with humans.
For simplicity, we will ignore the donor effects within each study and only consider the removal of the batch effect between the two studies.
However, we note that it is possible to apply the MNN correction between donors in each batch and then between the batches - see `?fastMNN` for details.

# Feature selection across batches

To obtain a single set of features for batch correction, we compute the average biological component across all batches.
We then take all genes with positive biological components to ensure that all interesting biology is retained, equivalent to the behaviour of `denoisePCA()`.
However, the quality of the correction can often be sensitive to technical noise, which means that some discretion may be required during feature selection.
Users may prefer to take the top 1000-5000 genes with the largest average components, or to use `combineVar()` to obtain combined $p$-values for gene selection.

```{r}
universe <- intersect(rownames(dec.gse85241), rownames(dec.gse81076))
mean.bio <- (dec.gse85241[universe,"bio"] + dec.gse81076[universe,"bio"])/2
chosen <- universe[mean.bio > 0]
length(chosen)
```

We also 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.
(Keep in mind that the previously computed 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.

```{r}
rescaled <- batchelor::multiBatchNorm(
    sce.gse85241[universe,], 
    sce.gse81076[universe,]
)
rescaled.gse85241 <- rescaled[[1]]
rescaled.gse81076 <- rescaled[[2]]
```

**Comments from Aaron:**

- Technically, we should have performed variance modelling and feature selection _after_ calling `multiBatchNorm()`.
This ensures that the variance components are estimated from the same values to be used in the batch correction.
In practice, this makes little difference, and it tends to be easier to process each batch separately and consolidate all results in one step as shown above.
- The `batchelor::` prefix avoids ambiguity during the migration of `multiBatchNorm()` from `r Biocpkg("scran")` to `r Biocpkg("batchelor")`.
This will not be necessary in the next release.

# Performing MNN-based correction

Consider a cell $a$ in batch $A$, and identify the cells in batch $B$ that are nearest neighbours to $a$ in the expression space defined by the selected features.
Repeat this for a cell $b$ in batch $B$, identifying its nearest neighbours in $A$.
Mutual nearest neighbours are pairs of cells from different batches that belong in each other's set of nearest neighbours.
The reasoning is that MNN pairs represent cells from the same biological state prior to the application of a batch effect - see @haghverdi2018batch 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 can yield batch-corrected values.

We apply the `fastMNN()` function to the three batches to remove the batch effect, using the genes in `chosen`.
To reduce computational work and technical noise, all cells in all cells 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.
The function returns a `SingleCellExperiment` object containing low-dimensional corrected values for downstream analyses like clustering or visualization.

```{r}
set.seed(100) 
unc.gse81076 <- logcounts(rescaled.gse81076)[chosen,]
unc.gse85241 <- logcounts(rescaled.gse85241)[chosen,]

mnn.out <- batchelor::fastMNN(
    GSE81076=unc.gse81076, GSE85241=unc.gse85241,
    k=20, d=50, BSPARAM=IrlbaParam(deferred=TRUE)
)
mnn.out
```

Each column of `mnn.out` corresponds to a cell in one of the batches, while each row corresponds to an input gene in `chosen`.
The `corrected` matrix in the `reducedDims` slot contains the low-dimensional corrected coordinates for all cells.

```{r}
dim(reducedDim(mnn.out, "corrected"))
```

The `batch` field in the column metadata contains an object specifying the batch of origin of each cell. 

```{r}
# Using an Rle for pretty-printing of batch IDs
# (as all cells from the same batch are consecutive).
Rle(mnn.out$batch) 
```

Advanced users may also be interested in the list of `DataFrame`s in the `pairs` metadata field.
Each `DataFrame` describes the MNN pairs identified upon merging of each successive batch.
This may be useful for checking the identified MNN pairs against known cell type identity, e.g., to determine if the cell types are being paired correctly.

```{r}
metadata(mnn.out)$merge.info$pairs[[1]]
```

As previously mentioned, we have only used two batches here to simplify the workflow.
However, the MNN approach is not limited to two batches, and inclusion of more batches is as simple as adding more `SingleCellExperiment` objects to the `fastMNN()` call.

**Comments from Aaron:**

- The `k=` parameter specifies the number of nearest neighbours to consider when defining MNN pairs.
This should be interpreted as the minimum frequency of each cell type or state in each batch.
  - Larger values will improve the precision of the correction by increasing the number of MNN pairs.
    It also provides some robustness to violations of the assumption that the batch vector is orthogonal to the biological subspace [@haghverdi2018batch],
    by allowing the neighbour search to ignore biological variation in each batch to identify the correct MNN pairs.
  - However, larger values of `k` can also reduce accuracy by allowing incorrect MNN pairs to form between cells of different types.
    Thus, we suggest starting with the default `k` and increasing it if one is confident that the same cell types are not adequately merged across batches.
    This is better than starting with a large `k` as incorrect merging is much harder to diagnose than insufficient merging.
- When `BSPARAM=IrlbaParam(deferred=TRUE)`, `fastMNN()` uses methods from the `r CRANpkg("irlba")` package to perform the principal components analysis quickly.
While the run-to-run differences should be minimal, it does mean that `set.seed()` is required to obtain fully reproducible results.
The `deferred=` argument instructs `fastMNN()` to sacrifice some numerical precision for greater speed.

```{r, echo=FALSE, results="hide"}
gc()
```

# Examining the effect of correction

## By visualization

We examine the batch correction with some _t_-SNE plots.
Figure~\@ref(fig:tsne-batch) demonstrates how the cells separate by batch of origin in the uncorrected data.
After correction, more intermingling between batches is observed, consistent with the removal of batch effects.
Note that the E-MTAB-5601 dataset still displays some separation, which is probably due to the fact that the other batches are UMI datasets.

```{r tsne-batch, fig.width=10, fig.asp=0.6, fig.cap="t-SNE plots of the pancreas datasets, before and after MNN correction. Each point represents a cell and is coloured by the batch of origin."}
# Adding uncorrected values.
sce <- mnn.out
assay(sce, "original") <- cbind(unc.gse81076, unc.gse85241)

# Using irlba to set up the t-SNE, for speed.
set.seed(100)
osce <- runPCA(sce, exprs_values="original", ntop=Inf, BSPARAM=IrlbaParam())
osce <- runTSNE(osce, use_dimred="PCA")
ot <- plotTSNE(osce, colour_by="batch") + ggtitle("Original")

# Corrected.
set.seed(100)
sce <- runTSNE(sce, use_dimred="corrected")
ct <- plotTSNE(sce, colour_by="batch") + ggtitle("Corrected")

multiplot(ot, ct, cols=2)
```

We colour by the expression of marker genes for known pancreas cell types to determine whether the correction is biologically sensible.
Cells in the same visual cluster express the same marker genes (Figure \@ref(fig:tsne-markers)), indicating that the correction maintains separation of cell types.

```{r tsne-markers, fig.width=10, fig.height=10, fig.cap="t-SNE plots after MNN correction, where each point represents a cell and is coloured by its corrected expression of key marker genes for known cell types in the pancreas."}
# Replacing the row names for easier reference.
rowData(sce)$ENSEMBL <- rownames(sce)    
rowData(sce)$SYMBOL <- mapIds(org.Hs.eg.db, keytype="ENSEMBL", 
    keys=rownames(sce), column="SYMBOL")
rownames(sce) <- uniquifyFeatureNames(rownames(sce), rowData(sce)$SYMBOL)

ct.gcg <- plotTSNE(sce, by_exprs_values="reconstructed", colour_by="GCG") 
ct.ins <- plotTSNE(sce, by_exprs_values="reconstructed", colour_by="INS") 
ct.sst <- plotTSNE(sce, by_exprs_values="reconstructed", colour_by="SST") 
ct.ppy <- plotTSNE(sce, by_exprs_values="reconstructed", colour_by="PPY") 

multiplot(ct.gcg + ggtitle("Alpha cells"),
    ct.ins + ggtitle("Beta cells"),
    ct.sst + ggtitle("Delta cells"),
    ct.ppy + ggtitle("PP cells"),
    cols=2)
```

## With diagnostics

One useful diagnostic is the proportion of variance within each batch that is lost during MNN correction.
Specifically, this refers to the within-batch variance that is removed during orthogonalization with respect to the average correction vector at each merge step. 
This is returned via the `lost.var` field in the metadata of `mnn.out`, which contains a matrix of the variance lost in each batch (column) at each merge step (row).

```{r}
metadata(mnn.out)$merge.info$lost.var
```

Large proportions of lost variance suggest that correction is removing genuine biological heterogeneity.
This would occur due to violations of the assumption of orthogonality between the batch effect and the biological subspace [@haghverdi2018batch].
In this case, the proportion of lost variance is small, indicating that non-orthogonality is not a major concern.

```{r, echo=FALSE, results="hide"}
gc()
```

# Controlling the merge order

## Manual specification 

The order of the supplied batches will affect the result as the first batch is used to define the reference space to which all other batches are corrected.
Specifically, the first batch is defined as the reference batch;
the second batch is corrected to and merged with the current reference batch, yielding a new reference batch;
and so on for all batches in the supplied order, with an increasingly large reference batch at each step.
The use of a merged reference ensures that information from batches other than the first are used to identify MNN pairs in later batches.

The order of batches to merge can be manually specified in the `auto.order=` argument to `fastMNN()`.
In the example shown below, batch 2 is treated as the reference as it is the first specified batch in `auto.order=`.
The first batch is then corrected to the second batch to obtain a new reference batch - and so on, if more than two batches were present.

```{r}
mnn.out2 <- batchelor::fastMNN(
    GSE81076=unc.gse81076, GSE85241=unc.gse85241,
    k=20, d=50, auto.order=c(2,1), BSPARAM=IrlbaParam(deferred=TRUE)
)
metadata(mnn.out2)$merge.order # batch 2 (GSE85241) is first in the order.
metadata(mnn.out2)$merge.info$pairs[[1]] # 'first' now refers to GSE85241.
```

Using `auto.order=` will change the merge order without requiring a change to the supplied order of batches in `fastMNN()`.
Similarly, the order of batches (and cells) in the output will not be altered.
This makes it easy to explore different merge orders without altering the surrounding code.

```{r}
Rle(mnn.out2$batch) # same as mnn.out$batch
```

If very different batches (in terms of cell composition) are present, we suggest setting the largest, most heterogeneous batch as the first.
This ensures that sufficient MNN pairs will be identified between the first and other batches for stable correction.
Conversely, if two small batches without shared populations are supplied first, the wrong MNN pairs will be detected and the result of the merge will be incorrect.

```{r, echo=FALSE, results="hide"}
gc()
```

## Hierarchical merging 

### Motivation  

In more complex experiments, we may know beforehand that certain sets of batches are more similar to each other.
We might then want to merge those similar batches before attempting the more difficult merges involving batches with different cell type composition and/or gene expression.
Examples include:

- Merging batches that represent replicate experiments from the same condition, prior to merging across conditions.
- Merging batches generated with the same scRNA-seq technology prior to merging across technologies.

This strategy limits encourages detection of correct MNN pairs as similar batches should have more shared populations, 
By comparison, performing the more difficult merges first is more likely to introduce errors whereby distinct subpopulations are incorrectly placed together.
This unnecessarily propagates the error to later steps as the initial merge is used as a reference for subsequent merges.

### Setting up the inputs

Hierarchical merging can be achieved through multiple calls to `fastMNN()` with progressively merged batches.
To illustrate, assume that we want to remove `Donor` effects within each batch prior to merging across batches.
We split up the cells in GSE85241 according to the donor of origin, creating one `SingleCellExperiment` object for each donor.

```{r}
all.donors <- unique(rescaled.gse85241$Donor)
table(rescaled.gse85241$Donor)
by.donor.85241 <- vector("list", length(all.donors))
names(by.donor.85241) <- sort(all.donors)
for (x in all.donors) {
    by.donor.85241[[x]] <- rescaled.gse85241[,rescaled.gse85241$Donor==x]
}
```

We repeat this process for GSE81076.
For demonstration purposes, we will aggregate some of the donors together to ensure that there are enough cells in each level for MNN detection.

```{r}
adj.donors <- c(D101="A", D102="A", D10631="A",
    D17="B", D1713="B", D172444="B",
    D2="C", D3="C",
    D71="D", D72="D", D73="D", D74="D")[rescaled.gse81076$Donor]
table(adj.donors)

all.donors <- unique(adj.donors)
by.donor.81076 <- vector("list", length(all.donors))
names(by.donor.81076) <- sort(all.donors)
for (x in all.donors) {
    by.donor.81076[[x]] <- rescaled.gse81076[,adj.donors==x]
}
```

We use the `multiBatchPCA()` function to perform a PCA across _all_ batches to be merged.
This ensures that all cells are placed onto the same coordinate space, which would obviously not be possible if a PCA was performed for each batch separately.
Specifically, `multiBatchPCA()` performs a modified PCA to ensure that each supplied matrix contributes equally to the definition of the PC space.
This avoids problems with imbalances in the number of cells across batches, meanining that smaller batches (possibly with unique cell types) are not ignored.

```{r}
all.batches <- c(by.donor.85241, by.donor.81076)

# Cosine normalizing for consistency with fastMNN() defaults.
all.logcounts <- lapply(all.batches, logcounts)
scaled <- lapply(all.logcounts, batchelor::cosineNorm) 

set.seed(1000) # for irlba.
pc.all <- do.call(batchelor::multiBatchPCA, c(scaled, 
    list(d=50, BSPARAM=IrlbaParam(deferred=TRUE))
))
names(pc.all)
```

**Comments from Aaron:**

- Here, we have applied `multiBatchPCA()` to the batch-level inputs for convenience.
It is also possible to supply donor-level matrices to equalize contributions across donors.
- Recall that we ran `multiBatchNorm()` earlier to generate the `rescaled.gse*` objects.
As a result, cosine normalization is not technically necessary, as all batches should be on the same scale already (see `?cosineNorm` for a discussion of this).
Nonetheless, we run it here for consistency with our previous `fastMNN()` call where cosine normalization is turned on by default.

### Performing progressive merges 

We pass the matrices corresponding to the donors in GSE85241 to `fastMNN()`, setting `pc.input=TRUE` to indicate that dimensionality reduction has already been performed.
This uses the first donor to define the reference space^[Which can be changed with `auto.order=`, if so desired.] and merges cells from all other donors to the first.
In this manner, we remove donor effects within the GSE85241 batch.

```{r}
pcs.85241 <- pc.all[seq_along(by.donor.85241)]
mnn.out.85241 <- do.call(batchelor::fastMNN, c(pcs.85241, list(pc.input=TRUE)))
Rle(mnn.out.85241$batch)
```

Note that when `pc.input=TRUE`, a `DataFrame` is returned instead of a `SingleCellExperiment`.
This reflects the fact that per-gene identities are lost when PCs are supplied instead of per-gene expression vectors^[Though it is entirely possible to reconstruct these using rotation vectors, see below.].
We repeat this for the donors in the GSE81076 batch.

```{r}
pcs.81076 <- tail(pc.all, length(by.donor.81076))
mnn.out.81076 <- do.call(batchelor::fastMNN, c(pcs.81076, list(pc.input=TRUE)))
Rle(mnn.out.81076$batch)
```

The next step is to merge the two batches together.
To do this, we simply repeat the `fastMNN()` call with the donor-corrected values for each batch.
Again, we need to set `pc.input=TRUE` to prevent the function from unnecessarily (and incorrectly) repeating the cosine normalization and PCA steps on the corrected values.
We use a larger `k` in the final `fastMNN()` call to improve the robustness of MNN detection to outlier cells on the edges of each cluster.

```{r}
mnn.out3 <- batchelor::fastMNN(
    GSE81076=mnn.out.81076,
    GSE85241=mnn.out.85241, 
    pc.input=TRUE, k=100 # see comments below.
) 

Rle(mnn.out3$batch) # by dataset
Rle(c(mnn.out.81076$batch, mnn.out.85241$batch)) # by donor
```

This yields a final corrected expression matrix where both within-batch donor effects and batch effects have been corrected.
We examine the quality of each of the merge steps with t-SNE plots (Figure \@ref(fig:tsne-hmerge)).
Within each batch, the donors are generally well-mixed, and the final merge is consistent with Figure \@ref(fig:tsne-batch). 

```{r tsne-hmerge, fig.width=10, fig.asp=0.5, fig.cap="t-SNE plots after correcting for donor effects within each data set, and after correcting for batch effects between data sets (final). Each point represents a cell that is coloured according to its donor of origin (left, middle) or the data set of origin (right)."}
set.seed(1000)
par(mfrow=c(1,3))

library(Rtsne)
tout.85241 <- Rtsne(mnn.out.85241$corrected, pca=FALSE)
plot(tout.85241$Y[,1], tout.85241$Y[,2], main="GSE85241 donors",
    col=as.factor(mnn.out.85241$batch), xlab="tSNE1", ylab="tSNE2")

tout.81076 <- Rtsne(mnn.out.81076$corrected, pca=FALSE)
plot(tout.81076$Y[,1], tout.81076$Y[,2], main="GSE81076 donors",
    col=as.factor(mnn.out.81076$batch), xlab="tSNE1", ylab="tSNE2")

tout.all <- Rtsne(mnn.out3$corrected, pca=FALSE)
plot(tout.all$Y[,1], tout.all$Y[,2], main="Final",
    col=as.factor(mnn.out3$batch), xlab="tSNE1", ylab="tSNE2")
```

**Comments from Aaron:**

- Here, we have applied `multiBatchPCA()` to the batch-level inputs for convenience.
It is also possible to supply donor-level matrices to equalize contributions across donors.
- One might ask why a larger `k` was not needed in Figure \@ref(fig:tsne-batch) as well.
This was probably because the outliers were masked by inter-donor heterogeneity when the t-SNE plots were generated without removing the donor effects.
- In this specific example, cells from the same donor will occupy contiguous rows in the `mnn.out3$corrected` matrix.
However, this may not have been the case for the original ordering of cells in each `SingleCellExperiment`.
This requires some extra account-keeping to match up the final corrected matrix to the original ordering, e.g., when cross-referencing to metadata.

    ```{r}
    original.plate <- unlist(lapply(rescaled, "[[", i="Plate"))
    original.names <- unlist(lapply(rescaled, colnames))
    
    # Needs unique names: trigger error otherwise.
    stopifnot(anyDuplicated(original.names)==0L)

    m <- match(rownames(mnn.out3$corrected), original.names)
    new.plate <- original.plate[m]
    ```

```{r, echo=FALSE, results="hide"}
gc()
```

## Automatic specification

In situations where the nature of each batch is unknown, users can set `auto.order=TRUE` to allow `fastMNN()` to empirically choose which batches to merge at each step.
The first merge is performed between the pair of batches with the most MNN pairs.
Progressive merges are performed with the remaining batch that has the most MNN pairs with the current reference batch.
The aim is to maximize the number of MNN pairs at each step to provide a stable correction.
We demonstrate below using the within-donor subbatches in the GSE81076 data set.

```{r}
mnn.out.auto <- do.call(batchelor::fastMNN, c(pcs.81076, 
    list(pc.input=TRUE, auto.order=TRUE)))
names(by.donor.81076) # supplied order 
metadata(mnn.out.auto)$merge.order # automatically defined order
```

As with manual specification, the merge order does not affect the output order of cells.
This makes it easy to try different merge orderings without having to deal with a re-ordering of output cells.

```{r}
Rle(mnn.out.auto$batch) 
```

The obvious cost of this approach is that of computation time. 
Nearest-neighbour searches need to be performed between all pairs of batches, and then between each remaining batch and the reference at each merge step.
As such, we prefer manual definition of a merge order that makes better use of prior knowledge about the experiment design.

```{r, echo=FALSE, results="hide"}
gc()
```

# Using the corrected values in downstream analyses

## For cell-based analyses

The low-dimensional corrected values can be used in any procedure that involves computing and comparing cell-cell (Euclidean) distances.
Recall that the aim of batch correction is to bring together related cells from different batches while preserving biological differences between cells within each batch.
The exact values of the corrected coordinates may not be interpretable, but this is not a problem if we are interested in the relative magnitudes of the distances.
For example, the code below directly uses the MNN-corrected values for clustering. 

```{r}
snn.gr <- buildSNNGraph(sce, use.dimred="corrected")
clusters <- igraph::cluster_walktrap(snn.gr)
table(clusters$membership, sce$batch)
```

Figure \@ref(fig:tsne-cluster) shows strong correspondence between the cluster labels and separation in _t_-SNE space.

```{r tsne-cluster, fig.cap="t-SNE plot after MMN correction, where each point represents a cell and is coloured by its cluster identity."}
sce$Cluster <- factor(clusters$membership)
plotTSNE(sce, colour_by="Cluster")
```

The corrected values can be used in any procedure that operates on cell-cell distances.
This includes nearest-neighbor searches, non-linear visualization like _t_-SNE and trajectory inference.

## For gene-based analyses

For gene-based procedures like differential expression (DE) analyses or gene network construction, it is desirable to use the **original** log-expression values or counts.
The corrected values are only used to obtain cell-level results such as clusters or trajectories.
Batch effects are handled explicitly using blocking terms or via a meta-analysis across batches. 
We do not use the corrected values directly in gene-based analyses, for various reasons:

- It is usually inappropriate to perform DE analyses on batch-corrected values, due to the failure to model the uncertainty of the correction.
This usually results in loss of type I error control, i.e., more false positives than expected.
- The correction does not preserve the mean-variance relationship.
Applications of common DE methods like `r Biocpkg("edgeR")` or `r Biocpkg("limma")` are unlikely to be valid.
- Batch correction may (correctly) remove biological differences between batches in the course of mapping all cells onto a common coordinate system.
Returning to the uncorrected expression values provides an opportunity for detecting such differences if they are of interest.
Conversely, if the batch correction made a mistake, the use of the uncorrected expression values provides an important sanity check.

Indeed, in the specific case of `fastMNN()`, the batch-corrected values no longer correspond to per-gene expression values anyway.
This means that they cannot be directly used in gene-based analyses^[Though one can circumvent this, see below.].

Users should generally aim to avoid using batch-corrected values for per-gene analyses when within-batch alternatives are available. 
To illustrate, we perform a DE analysis on the uncorrected expression data using the clusters identified previously.
To model the batch effect, we set the batch of origin as the `block=` argument in `findMarkers()`.
This will perform all comparisons between clusters _within_ each batch, and then combine the $p$-values to consolidate results across batches.

```{r}
m.out <- findMarkers(sce, clusters$membership, block=sce$batch,
    direction="up", assay.type="original")
demo <- m.out[["10"]] # probably alpha cells.
demo <- demo[demo$Top <= 5,]
as.data.frame(demo[,1:3]) # only first three columns for brevity.
```

```{r, echo=FALSE, results="hide"}
# Checking that cluster 10 is what we think it is.
stopifnot(all(sapply(demo["GCG",-(1:3)], sign)==1))
```

Other approaches for handling batch effects during marker gene detection are discussed `r simpleSingleCell:::.link("de", "Blocking on uninteresting factors of variation", "elsewhere")`.
It is similarly possible to perform these analyses with standard Bioconductor packages for DE analysis such as `r Biocpkg("edgeR")` or `r Biocpkg("limma")`.
Note that the use of `block=` is roughly similar to the use of a batch-cluster interaction model and testing whether the average log-fold change across batches is equal to zero.

**Comments from Aaron:**

- Users of the older `mnnCorrect()` function will note that the function returned corrected expression values.
It is tempting to use these corrected values directly for DE analyses, but this would likely be inappropriate.
In addition to the reasons discussed above, the default parameters of `mnnCorrect()` do not return corrected values on the log-scale, but rather a cosine-normalized log-scale.
This makes it difficult to interpret the effect size of DE analyses based on the corrected values.

## Obtaining per-gene corrected values

When applied on gene expression values, `fastMNN()` will also return a matrix of corrected per-gene expression values in the `reconstructed` assay.
This is obtained by taking the cross-product of the corrected low-dimensional values with the rotation vectors from the initial PCA,
which effectively reverses the initial projection into a low-dimensional space during `multiBatchPCA()`.
The example below extracts the corrected expression values for insulin.

```{r}
assay(sce, "reconstructed")
summary(assay(sce)["INS",])
```

If `fastMNN()` was run on low-dimensional inputs, only the low-dimensional output will be reported.
Nonetheless, users can obtain per-gene corrected values by manually computing the cross-product using the PCA rotation vectors.
For example, the code below obtains corrected expression values for _GCG_ from our hierarchical merge.

```{r}
rotations <- metadata(pc.all)$rotation
cor.exp <- tcrossprod(mnn.out3$corrected,
    rotations["ENSG00000115263",,drop=FALSE])
summary(cor.exp)
```

Explicit calculation of all per-gene corrected values is probably ill-advised as this would involve the construction of a dense matrix.
This may be prohibitively memory-consuming for large data sets that are otherwise representable as sparse matrices.
Rather, corrected values can be computed for specific genes as they are needed, e.g., using the `LowRankMatrix` class.

Per-gene corrected values can be readily used for visualization, e.g., in Figure \@ref(fig:tsne-markers).
This can be more aesthetically pleasing than uncorrected expression values that may contain large shifts on the colour scale between cells in different batches.
However, use of the corrected values in any quantitative procedure should be treated with extreme caution.
If they must be used, they should be backed up by similar results from an analysis on the uncorrected values.

# Concluding remarks

We save the `SingleCellExperiment` object for use elsewhere.
This avoids the need to repeat all of the processing steps described above.

```{r}
saveRDS(file="pancreas_data.rds", sce)
```

All software packages used in this workflow are publicly available from the Comprehensive R Archive Network (https://cran.r-project.org) or the Bioconductor project (http://bioconductor.org).
The specific version numbers of the packages used are shown below, along with the version of the R installation.

```{r}
sessionInfo()
```

# References