## ---- 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) sce1 <- ZeiselBrainData() sce1 sce2 <- TasicBrainData()[,1:500] sce2 ## ----------------------------------------------------------------------------- set.seed(19999) sce1 <- sce1[,sample(ncol(sce1), 500, replace=TRUE)] sce2 <- sce2[,sample(ncol(sce2), 500, replace=TRUE)] ## ----------------------------------------------------------------------------- universe <- intersect(rownames(sce1), rownames(sce2)) sce1 <- sce1[universe,] sce2 <- sce2[universe,] ## ----------------------------------------------------------------------------- # Ignoring spike-ins with use_altexps=FALSE. out <- multiBatchNorm(sce1, sce2, norm.args=list(use_altexps=FALSE)) sce1 <- out[[1]] sce2 <- out[[2]] ## ----------------------------------------------------------------------------- library(scran) dec1 <- modelGeneVar(sce1) dec2 <- modelGeneVar(sce2) combined.dec <- combineVar(dec1, dec2) chosen.hvgs <- combined.dec$bio > 0 summary(chosen.hvgs) ## ----------------------------------------------------------------------------- library(scater) combined <- correctExperiments(A=sce1, B=sce2, PARAM=NoCorrectParam()) combined <- runPCA(combined, subset_row=chosen.hvgs) plotPCA(combined, colour_by="batch") ## ----------------------------------------------------------------------------- library(batchelor) f.out <- fastMNN(A=sce1, B=sce2, subset.row=chosen.hvgs) str(reducedDim(f.out, "corrected")) ## ----------------------------------------------------------------------------- rle(f.out$batch) ## ----------------------------------------------------------------------------- f.out2 <- fastMNN(combined, batch=combined$batch, subset.row=chosen.hvgs) str(reducedDim(f.out2, "corrected")) ## ----------------------------------------------------------------------------- plotReducedDim(f.out, colour_by="batch", dimred="corrected") ## ----------------------------------------------------------------------------- cor.exp <- assay(f.out)[1,] hist(cor.exp, xlab="Corrected expression for gene 1", col="grey80") ## ----------------------------------------------------------------------------- # Using fewer genes as it is much slower. fewer.hvgs <- head(order(combined.dec$bio, decreasing=TRUE), 1000) classic.out <- mnnCorrect(sce1, sce2, subset.row=fewer.hvgs) ## ----------------------------------------------------------------------------- classic.out ## ----------------------------------------------------------------------------- rescale.out <- rescaleBatches(sce1, sce2) rescale.out ## ----------------------------------------------------------------------------- rescale.out <- runPCA(rescale.out, subset_row=chosen.hvgs, exprs_values="corrected") plotPCA(rescale.out, colour_by="batch") ## ----------------------------------------------------------------------------- regress.out <- regressBatches(sce1, sce2) assay(regress.out) ## ----------------------------------------------------------------------------- # 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, norm.args=list(use_altexps=FALSE)) names(normed) ## ----------------------------------------------------------------------------- # Using the same BSPARAM argument as fastMNN(), for speed. pca.out <- multiBatchPCA(A=sce1, B=sce2, subset.row=chosen.hvgs, BSPARAM=BiocSingular::IrlbaParam(deferred=TRUE)) names(pca.out) ## ----------------------------------------------------------------------------- sessionInfo()