## ---- 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) register(SerialParam()) set.seed(100) ## ------------------------------------------------------------------------ ngenes <- 10000 ncells <- 200 mu <- 2^runif(ngenes, 3, 10) gene.counts <- matrix(rnbinom(ngenes*ncells, mu=mu, size=2), 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, 3, 10) spike.counts <- matrix(rnbinom(nspikes*ncells, mu=mu, size=2), nrow=nspikes) rownames(spike.counts) <- paste0("ERCC-", seq_len(nspikes)) all.counts <- rbind(gene.counts, spike.counts) ## ------------------------------------------------------------------------ library(scran) sce <- newSCESet(countData=data.frame(all.counts)) isSpike(sce) <- rep(c(FALSE, TRUE), c(ngenes, nspikes)) ## ------------------------------------------------------------------------ sce <- computeSumFactors(sce) summary(sizeFactors(sce)) ## ------------------------------------------------------------------------ larger.sce <- newSCESet(countData=data.frame(cbind(all.counts, all.counts, all.counts))) clusters <- quickCluster(larger.sce) larger.sce <- computeSumFactors(larger.sce, cluster=clusters) ## ------------------------------------------------------------------------ sce2 <- computeSpikeFactors(sce) summary(sizeFactors(sce2)) ## ------------------------------------------------------------------------ mm.pairs <- readRDS(system.file("exdata", "mouse_cycle_markers.rds", package="scran")) ## ------------------------------------------------------------------------ assigned <- cyclone(sce, mm.pairs) head(assigned$scores) ## ------------------------------------------------------------------------ phase <- rep("S", ncol(sce)) phase[assigned$scores$G1 > 0.5] <- "G1" phase[assigned$scores$G2M > 0.5] <- "G2M" phase[assigned$scores$G1 > 0.5 & assigned$scores$G2M > 0.5] <- "unknown" table(phase) ## ------------------------------------------------------------------------ fit <- trendVar(sce) ## ------------------------------------------------------------------------ decomp <- decomposeVar(sce, fit) top.hvgs <- order(decomp$bio, decreasing=TRUE) head(decomp[top.hvgs,]) ## ------------------------------------------------------------------------ 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) ## ---- eval=FALSE--------------------------------------------------------- ## set.seed(102251) ## ncells <- 100 ## niters <- 20 ## averaged.GLM <- averaged.GLMM <- averaged.log <- 0 ## for (it in seq_len(niters)) { ## ## counter1 <- rnbinom(ncells, mu=5, size=1) ## counter1[1] <- 10000 ## counter2 <- rnbinom(ncells, mu=5, size=0.2) # Arguably, all more variable than above. ## counter3 <- rnbinom(ncells, mu=20, size=0.2) ## counter4 <- rnbinom(ncells, mu=100, size=0.2) ## ## # With GLMs. ## require(edgeR) ## y <- DGEList(rbind(counter1, counter2, counter3, counter4), lib.size=rep(1e6, ncells)) ## y <- estimateDisp(y, cbind(rep(1, ncells)), prior.df=0) ## averaged.GLM <- averaged.GLM + y$tagwise.dispersion ## ## # With GLMMs (random effect variance is logged). ## require(lme4) ## plate <- integer(ncells) ## plate[(ncells/2+1):ncells] <- 1 # Just to stop it complaining. ## averaged.GLMM <- averaged.GLMM + c( ## unlist(VarCorr(glmer(counter1 ~ (1|plate), family=negative.binomial(1)))), ## unlist(VarCorr(glmer(counter2 ~ (1|plate), family=negative.binomial(1)))), ## unlist(VarCorr(glmer(counter3 ~ (1|plate), family=negative.binomial(1)))), ## unlist(VarCorr(glmer(counter4 ~ (1|plate), family=negative.binomial(1))))) ## ## # With variances of log-values. ## averaged.log <- averaged.log + apply(cpm(y, log=TRUE), 1, var) ## ## } ## averaged.GLM/niters ## averaged.GLMM/niters ## averaged.log/niters ## ------------------------------------------------------------------------ null.dist <- correlateNull(ncol(sce)) cor.pairs <- correlatePairs(sce[top.hvgs[1:200],], null.dist=null.dist) head(cor.pairs) ## ------------------------------------------------------------------------ null.dist2 <- correlateNull(design=design) cor.pairs2 <- correlatePairs(sce[top.hvgs[1:200],], null.dist=null.dist2, design=design) ## ---- eval=FALSE--------------------------------------------------------- ## set.seed(1023423) ## ncells <- 100 ## null.dist <- correlateNull(ncells) ## all.p <- list() ## for (it in 1:10000) { ## x1 <- rpois(ncells, lambda=10) ## x2 <- rpois(ncells, lambda=20) ## rho2 <- cor(rank(x1, ties.method="random"), rank(x2, ties.method="random"), method="spearman") ## all.p[[it]] <- sum(null.dist >= rho2)/length(null.dist) ## } ## sum(unlist(all.p) <= 0.01)/10000 ## sum(unlist(all.p) <= 0.05)/10000 ## ------------------------------------------------------------------------ y <- convertTo(sce, type="edgeR") ## ------------------------------------------------------------------------ sessionInfo()