## ---- echo=FALSE, results="hide", message=FALSE----------------------------
require(knitr)
opts_chunk$set(error=FALSE, message=FALSE, warning=FALSE)

## ----setup, echo=FALSE, message=FALSE--------------------------------------
library(batchelor)
set.seed(100)

## --------------------------------------------------------------------------
library(scRNAseq)
data(allen)

library(SingleCellExperiment)
sce1 <- as(allen, "SingleCellExperiment")
counts(sce1) <- assay(sce1)

library(scater)
sce1 <- sce1[1:2000,] # reducing the size for demo purposes.
sce1 <- normalize(sce1) # quick and dirty normalization.

## --------------------------------------------------------------------------
sce2 <- sce1
logcounts(sce2) <- logcounts(sce2) + rnorm(nrow(sce2), sd=2)

combined <- cbind(sce1, sce2)
combined$batch <- factor(rep(1:2, c(ncol(sce1), ncol(sce2))))

plotPCA(combined, colour_by="batch") # checking there is a batch effect.

## --------------------------------------------------------------------------
f.out <- fastMNN(A=sce1, B=sce2)
str(reducedDim(f.out, "corrected"))

## --------------------------------------------------------------------------
rle(f.out$batch)

## --------------------------------------------------------------------------
f.out2 <- fastMNN(combined, batch=combined$batch)
str(reducedDim(f.out2, "corrected"))

## --------------------------------------------------------------------------
plotReducedDim(f.out, colour_by="batch", use_dimred="corrected")

## --------------------------------------------------------------------------
cor.exp <- assay(f.out)[1,]
hist(cor.exp, xlab="Corrected expression for gene 1") 

## --------------------------------------------------------------------------
classic.out <- mnnCorrect(sce1, sce2)

## --------------------------------------------------------------------------
classic.out

## --------------------------------------------------------------------------
rescale.out <- rescaleBatches(sce1, sce2)
rescale.out

## --------------------------------------------------------------------------
# Just picking a bunch of random genes:
chosen.genes <- sample(nrow(combined), 1000)
f.out3 <- fastMNN(combined, batch=combined$batch,
    subset.row=chosen.genes)
str(reducedDim(f.out2, "corrected"))

## --------------------------------------------------------------------------
# Pretend the first X cells in each batch are controls.
restrict <- list(1:100, 1:200) 
rescale.out <- rescaleBatches(sce1, sce2, restrict=restrict)

## --------------------------------------------------------------------------
normed <- multiBatchNorm(A=sce1, B=sce2)
names(normed)

## --------------------------------------------------------------------------
pca.out <- multiBatchPCA(A=sce1, B=sce2)
names(pca.out)

## --------------------------------------------------------------------------
sessionInfo()