---
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}
library(BiocStyle)
library(knitr)
opts_chunk$set(error=FALSE, message=FALSE, warning=FALSE)
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)
```

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.
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 through `quickCluster()` to avoid pooling together very different cells.

```{r}
library(scran)
set.seed(1000)    
clusters <- quickCluster(sce.gse81076, method="igraph", min.mean=0.1)
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)
```

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, min.mean=0.1, method="igraph")
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 <- 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.

# 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 matrix of corrected values for downstream analyses like clustering or visualization.

```{r}
set.seed(100) 
original <- list(
    GSE81076=logcounts(rescaled.gse81076)[chosen,],
    GSE85241=logcounts(rescaled.gse85241)[chosen,]
)

# Slightly convoluted call to avoid re-writing code later.
# Equivalent to fastMNN(GSE81076, GSE85241, k=20, d=50, approximate=TRUE)
mnn.out <- do.call(fastMNN, c(original, list(k=20, d=50, approximate=TRUE)))
dim(mnn.out$corrected)
```

Each row of the `corrected` matrix corresponds to a cell in one of the batches.
The `batch` field contains a run-length encoding object specifying the batch of origin of each row.

```{r}
mnn.out$batch
```

Advanced users may also be interested in the list of `DataFrame`s in the `pairs` 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}
mnn.out$pairs
```

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, 
at the cost of reducing accuracy by allowing MNN pairs to form between cells of different type.
- The order of the supplied batches does matter, as all batches are corrected towards the first.
We suggest setting the largest and/or most heterogeneous batch as the first.
This ensures that sufficient MNN pairs will be identified between the first and other batches for stable correction.
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.
Batches are chosen to maximize the number of MNN pairs in order to provide a stable correction.
- When `approximate=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.

# Examining the effect of correction

We create a new `SingleCellExperiment` object containing log-expression values for all cells, along with information regarding the batch of origin.
The MNN-corrected values are stored as dimensionality reduction results, befitting the principal components analysis performed within `fastMNN()`. 

```{r}
omat <- do.call(cbind, original)
sce <- SingleCellExperiment(list(logcounts=omat))
reducedDim(sce, "MNN") <- mnn.out$corrected
sce$Batch <- as.character(mnn.out$batch)
sce
```

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."}
set.seed(100)
# Using irlba to set up the t-SNE, for speed.
osce <- runPCA(sce, ntop=Inf, method="irlba")
osce <- runTSNE(osce, use_dimred="PCA")
ot <- plotTSNE(osce, colour_by="Batch") + ggtitle("Original")

set.seed(100)
csce <- runTSNE(sce, use_dimred="MNN")
ct <- plotTSNE(csce, 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."}
ct.gcg <- plotTSNE(csce, colour_by="ENSG00000115263") + 
    ggtitle("Alpha cells (GCG)")
ct.ins <- plotTSNE(csce, colour_by="ENSG00000254647") + 
    ggtitle("Beta cells (INS)")
ct.sst <- plotTSNE(csce, colour_by="ENSG00000157005") + 
    ggtitle("Delta cells (SST)")
ct.ppy <- plotTSNE(csce, colour_by="ENSG00000108849") + 
    ggtitle("PP cells (PPY)")
multiplot(ct.gcg, ct.ins, ct.sst, ct.ppy, cols=2)
```

# Using the corrected values in downstream analyses

For downstream analyses, the MNN-corrected values can be treated in the same manner as any other dimensionality reduction result.
For example, it is straightforward to use the MNN-corrected values directly used for clustering analyses, as shown below.

```{r}
snn.gr <- buildSNNGraph(sce, use.dimred="MNN")
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."}
csce$Cluster <- factor(clusters$membership)
plotTSNE(csce, colour_by="Cluster")
```

Differential expression analyses should be performed on the **original** log-expression values or counts.
We do not use the corrected values here (which no longer correspond to genes anyway) except to obtain the clusters or trajectories to be characterized.
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 with Stouffer's Z method to consolidate results across batches.
Intra-batch comparisons are robust to differences between batches but assume that each pair of clusters is present in at least one batch.

```{r}
m.out <- findMarkers(sce, clusters$membership, block=sce$Batch,
    direction="up")        
demo <- m.out[["4"]] # looking at cluster 4 (probably alpha cells).
demo <- demo[demo$Top <= 5,]

library(org.Hs.eg.db)
data.frame(row.names=rownames(demo),
    Symbol=mapIds(org.Hs.eg.db, keytype="ENSEMBL", 
        keys=rownames(demo), column="SYMBOL"),
    Top=demo$Top, FDR=demo$FDR)
```

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

Another approach is to define a design matrix containing the batch of origin as the sole factor.
`findMarkers()` will then fit a linear model to the log-expression values, similar to the use of `r Biocpkg("limma")` for bulk RNA sequencing data [@ritchie2015limma].
This handles situations where multiple batches contain unique clusters, whereby comparisons can be implicitly performed between shared cell types.
However, the use of a linear model makes some strong assumptions about the homogeneity of the batch effect across cell types and the equality of variances.

```{r}
# Setting up the design matrix (we remove intercept for full rank
# in the final design matrix with the cluster-specific terms).
design <- model.matrix(~sce$Batch)
design <- design[,-1,drop=FALSE]

m.alt <- findMarkers(sce, clusters$membership, design=design,
    direction="up")
demo <- m.alt[["4"]]
demo <- demo[demo$Top <= 5,]

data.frame(row.names=rownames(demo),
    Symbol=mapIds(org.Hs.eg.db, keytype="ENSEMBL", 
        keys=rownames(demo), column="SYMBOL"),
    Top=demo$Top, FDR=demo$FDR)
```

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.
Thus, it is tempting to use these corrected values directly for DE analyses.
This is inappropriate for various reasons:
    - 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.
    - 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.

# Concluding remarks

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