inst/doc/scDD.R

In scDD: Mixture modeling of single-cell RNA-seq data to identify genes with differential distributions

## ----style-knitr, eval=TRUE, echo=FALSE, results="asis"--------------------
BiocStyle::latex()

## ----patternPlot, echo=FALSE, fig.show="hide", fig.width=7, fig.height=5----
### Note that the following code can be ignored for the purposes of
# analysis with scDD; it is simply used to generate the cartoon of 
# interesting DD patterns (illustration purposes only)
par(mfrow=c(2,2), tcl=-0.5, mai=c(0.4,0.4,0.5,0.3))
x <- seq(0, 6, by=0.05)

## traditional de 
# mu1 is 2 
# mu2 is 4
cord.x <- c(0,x,6) 
cord.y <- c(0,dnorm(x, 2, 0.75),0) 
curve(dnorm(x, 2 , 0.75),xlim=c(0,6),main="Traditional DE",
      xaxt="n", xlab="", ylab="", yaxt="n") 
polygon(cord.x,cord.y,col=rgb(0,0,1,1/4))
cord.x <- c(0,x,6) 
cord.y <- c(0,dnorm(x, 4, 0.75),0) 
lines(x, dnorm(x, 4 , 0.75))
polygon(cord.x,cord.y,col=rgb(1,0,0,1/4))
axis(side=1, at=c(2,4), labels=c(expression(mu[1]), expression(mu[2])),
     pos=0, cex.axis=1.5)
mtext("(A)", side = 3, line=0.5, adj=-0.1, cex=1.2, font=2)

x <- seq(0, 10, by=0.05)
## differential proportion
cord.x <- c(0,x,10) 
cord.y <- c(0,0.3*dnorm(x, 7, 1) + 0.7*dnorm(x, 3, 1),0) 
curve(0.3*dnorm(x, 7, 1) + 0.7*dnorm(x, 3, 1),xlim=c(0,10),main="DP",
      xaxt="n", xlab="", ylab="", yaxt="n") 
polygon(cord.x,cord.y,col=rgb(0,0,1,1/4))
cord.x <- c(0,x,10) 
cord.y <- c(0,0.3*dnorm(x, 3, 1) + 0.7*dnorm(x, 7, 1),0) 
lines(x, 0.3*dnorm(x, 3, 1) + 0.7*dnorm(x, 7, 1))
polygon(cord.x,cord.y,col=rgb(1,0,0,1/4))
axis(side=1, at=c(3,7), labels=c(expression(mu[1]), 
                                 expression(mu[2])), pos=0, cex.axis=1.5)
mtext("(B)", side = 3, line=0.5, adj=-0.1, cex=1.2, font=2)

## differential modes (DM)
cord.x <- c(0,x,6) 
cord.y <- c(0,dnorm(x, 2, 0.75),0) 
curve(dnorm(x, 2 , 0.75),xlim=c(0,6),main="DM", xaxt="n", 
      xlab="", ylab="", yaxt="n") 
polygon(cord.x,cord.y,col=rgb(0,0,1,1/4))
cord.x <- c(0,x,6) 
cord.y <- c(0,0.3*dnorm(x, 2, 0.6) + 0.7*dnorm(x, 4, 0.6),0) 
lines(x, 0.3*dnorm(x, 2, 0.6) + 0.7*dnorm(x, 4, 0.6))
polygon(cord.x,cord.y,col=rgb(1,0,0,1/4))
axis(side=1, at=c(2,4), labels=c(expression(mu[1]), 
                                 expression(mu[2])), pos=0, cex.axis=1.5)
mtext("(C)", side = 3, line=0.5, adj=-0.1, cex=1.2, font=2)

## Both DM and DP
cord.x <- c(0,x,10) 
cord.y <- c(0,0.5*dnorm(x, 2.5, 1) + 0.5*dnorm(x, 7.5, 1),0) 
curve(0.5*dnorm(x, 2.5, 1) + 0.5*dnorm(x, 7.5, 1),
      xlim=c(0,10),main="DB", xaxt="n", xlab="", ylab="", yaxt="n",
      ylim=c(0,max(0.5*dnorm(x, 2.5, 1) + 0.5*dnorm(x, 7.5, 1)))) 
polygon(cord.x,cord.y,col=rgb(0,0,1,1/4))
cord.x <- c(0,x,10) 
cord.y <- c(0,0.8*dnorm(x, 5, 2),0) 
lines(x, 0.8*dnorm(x, 5, 2))
polygon(cord.x,cord.y,col=rgb(1,0,0,1/4))
axis(side=1, at=c(2.5, 5, 7.5), labels=c(expression(mu[1]), 
                                         expression(mu[3]), 
                                         expression(mu[2])), 
     pos=0, cex.axis=1.5)
mtext("(D)", side = 3, line=0.5, adj=-0.1, cex=1.2, font=2)

## ----load lib, message=FALSE-----------------------------------------------
library(scDD)

## ----load exDat------------------------------------------------------------
data(scDatExSim)

## ----check class-----------------------------------------------------------
class(scDatExSim)
dim(scDatExSim)

## ----prior-----------------------------------------------------------------
prior_param=list(alpha=0.01, mu0=0, s0=0.01, a0=0.01, b0=0.01)

## ----main engine-----------------------------------------------------------
scDatExSim <- scDD(scDatExSim, prior_param=prior_param, testZeroes=FALSE)

## ----main results----------------------------------------------------------
RES <- results(scDatExSim)
head(RES)

## ----partition results-----------------------------------------------------
PARTITION.C1 <- results(scDatExSim, type="Zhat.c1")
PARTITION.C1[1:5,1:5]

## ----main engine perm------------------------------------------------------
scDatExSim <- scDD(scDatExSim, prior_param=prior_param, 
                 testZeroes=FALSE, permutations=100)

## ----load exdat2-----------------------------------------------------------
data(scDatEx)

## ----check class2----------------------------------------------------------
class(scDatEx)
dim(scDatEx)

## ----set num---------------------------------------------------------------
nDE <- 5
nDP <- 5
nDM <- 5
nDB <- 5
nEE <- 5
nEP <- 5
numSamples <- 100
seed <- 816

## ----simset----------------------------------------------------------------
SD <- simulateSet(scDatEx, numSamples=numSamples, 
                  nDE=nDE, nDP=nDP, nDM=nDM, nDB=nDB, 
                  nEE=nEE, nEP=nEP, sd.range=c(2,2), modeFC=4, plots=FALSE, 
                  random.seed=seed)

# load the SingleCellExperiment package to use rowData method
library(SingleCellExperiment)
head(rowData(SD))

## ----load summarizedexp, message=FALSE-------------------------------------
library(SingleCellExperiment)

## ----load exdatlist--------------------------------------------------------
data(scDatExList)

## ----create condition------------------------------------------------------
condition <- c(rep(1, ncol(scDatExList$C1)), rep(2,  ncol(scDatExList$C2)))

## ----rownames--------------------------------------------------------------
# Example of row and column names
head(rownames(scDatExList$C1))
head(colnames(scDatExList$C2))

names(condition) <- c(colnames(scDatExList$C1), colnames(scDatExList$C2))

## ----create sce------------------------------------------------------------
sce <- SingleCellExperiment(assays=list(normcounts=cbind(scDatExList$C1, 
                                                         scDatExList$C2)),
                                colData=data.frame(condition))
show(sce)

## ----load exdat3e----------------------------------------------------------
data(scDatEx)
show(scDatEx)

## ----preprocess------------------------------------------------------------
scDatEx <- preprocess(scDatEx, zero.thresh=0.9)
show(scDatEx)

## ----threshe---------------------------------------------------------------
scDatEx.scran <- preprocess(scDatEx, zero.thresh=0.5, scran_norm=TRUE)
show(scDatEx.scran)

## ----load exdat4-----------------------------------------------------------
data(scDatExSim)

## ----load sumExp, message=FALSE--------------------------------------------
library(SingleCellExperiment)

## ----plot DE, eval=TRUE, message=FALSE-------------------------------------
de <- sideViolin(normcounts(scDatExSim)[1,], scDatExSim$condition, 
           title.gene=rownames(scDatExSim)[1])

## ----plot DP, eval=TRUE, message=FALSE-------------------------------------
dp <- sideViolin(normcounts(scDatExSim)[6,], scDatExSim$condition, 
           title.gene=rownames(scDatExSim)[6])

## ----plot DM, eval=TRUE, message=FALSE-------------------------------------
dm <- sideViolin(normcounts(scDatExSim)[11,], scDatExSim$condition, 
           title.gene=rownames(scDatExSim)[11])

## ----plot DB, eval=TRUE, message=FALSE-------------------------------------
db <- sideViolin(normcounts(scDatExSim)[16,], scDatExSim$condition, 
           title.gene=rownames(scDatExSim)[16])

## ----plot EP, eval=TRUE, message=FALSE-------------------------------------
ep <- sideViolin(normcounts(scDatExSim)[21,], scDatExSim$condition, 
           title.gene=rownames(scDatExSim)[21])

## ----plot EE, eval=TRUE, message=FALSE-------------------------------------
ee <- sideViolin(normcounts(scDatExSim)[26,], scDatExSim$condition, 
           title.gene=rownames(scDatExSim)[26])

## ----plotGrid, fig.show='hide', fig.width=8.5, fig.height=11, message=FALSE----
library(gridExtra)
grid.arrange(de, dp, dm, db, ep, ee, ncol=2)

## ----sessionInfo, eval=TRUE------------------------------------------------
sessionInfo()