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

## ----style, echo=FALSE, results='asis'-----------------------------------
BiocStyle::markdown()

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

## ------------------------------------------------------------------------
ngenes <- 10000
ncells <- 200
mu <- 2^runif(ngenes, -1, 5)
gene.counts <- matrix(rnbinom(ngenes*ncells, mu=mu, size=10), nrow=ngenes)

## ------------------------------------------------------------------------
library(org.Mm.eg.db)
all.ensembl <- unique(toTable(org.Mm.egENSEMBL)$ensembl_id)
rownames(gene.counts) <- sample(all.ensembl, ngenes)

## ------------------------------------------------------------------------
nspikes <- 100
ncells <- 200
mu <- 2^runif(nspikes, -1, 5)
spike.counts <- matrix(rnbinom(nspikes*ncells, mu=mu, size=10), nrow=nspikes)
rownames(spike.counts) <- paste0("ERCC-", seq_len(nspikes))
all.counts <- rbind(gene.counts, spike.counts)

## ------------------------------------------------------------------------
library(scran)
sce <- SingleCellExperiment(list(counts=all.counts))
isSpike(sce, "MySpike") <- grep("^ERCC", rownames(sce))

## ------------------------------------------------------------------------
mm.pairs <- readRDS(system.file("exdata", "mouse_cycle_markers.rds", package="scran"))

## ------------------------------------------------------------------------
assigned <- cyclone(sce, pairs=mm.pairs)
head(assigned$scores)

## ------------------------------------------------------------------------
table(assigned$phases)

## ------------------------------------------------------------------------
sce <- computeSumFactors(sce)
summary(sizeFactors(sce))

## ------------------------------------------------------------------------
larger.sce <- SingleCellExperiment(list(counts=cbind(all.counts, all.counts, all.counts)))
clusters <- quickCluster(larger.sce)
larger.sce <- computeSumFactors(larger.sce, cluster=clusters)

## ------------------------------------------------------------------------
sce2 <- computeSpikeFactors(sce)
summary(sizeFactors(sce2))

## ------------------------------------------------------------------------
sce <- computeSpikeFactors(sce, general.use=FALSE)

## ------------------------------------------------------------------------
sce <- normalize(sce)

## ------------------------------------------------------------------------
fit <- trendVar(sce, parametric=TRUE)

## ------------------------------------------------------------------------
decomp <- decomposeVar(sce, fit)
top.hvgs <- order(decomp$bio, decreasing=TRUE)
head(decomp[top.hvgs,])

## ---- fig.cap=""---------------------------------------------------------
plot(decomp$mean, decomp$total, xlab="Mean log-expression", ylab="Variance")
o <- order(decomp$mean)
lines(decomp$mean[o], decomp$tech[o], col="red", lwd=2)
points(fit$mean, fit$var, col="red", pch=16)

## ------------------------------------------------------------------------
alt.fit <- trendVar(sce, use.spikes=FALSE) 
alt.decomp <- decomposeVar(sce, alt.fit)

## ------------------------------------------------------------------------
batch <- rep(c("1", "2"), each=100)
design <- model.matrix(~batch)
alt.fit2 <- trendVar(sce, design=design)
alt.decomp2 <- decomposeVar(sce, alt.fit)

## ------------------------------------------------------------------------
null.dist <- correlateNull(ncol(sce))
# Only using the first 200 genes as a demonstration.
cor.pairs <- correlatePairs(sce, subset.row=top.hvgs[1:200], null.dist=null.dist)
head(cor.pairs)

## ------------------------------------------------------------------------
null.dist2 <- correlateNull(design=design, iter=1e5) # fewer iterations, to speed it up.
cor.pairs2 <- correlatePairs(sce, subset.row=top.hvgs[1:200], 
    null.dist=null.dist2, design=design)

## ------------------------------------------------------------------------
cor.genes <- correlatePairs(sce, subset.row=top.hvgs[1:200], 
    null.dist=null.dist, per.gene=TRUE)

## ------------------------------------------------------------------------
y <- convertTo(sce, type="edgeR")

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