inst/doc/baySeq.R

In baySeq: Empirical Bayesian analysis of patterns of differential expression in count data

### R code from vignette source 'baySeq.Rnw'

###################################################
### code chunk number 1: <style-Sweave
###################################################
BiocStyle::latex()


###################################################
### code chunk number 2: baySeq.Rnw:34-36
###################################################
set.seed(102)
options(width = 90)


###################################################
### code chunk number 3: baySeq.Rnw:40-41
###################################################
library(baySeq)


###################################################
### code chunk number 4: baySeq.Rnw:45-46
###################################################
if(require("parallel")) cl <- makeCluster(8) else cl <- NULL


###################################################
### code chunk number 5: baySeq.Rnw:51-53
###################################################
data(simData)
simData[1:10,]


###################################################
### code chunk number 6: baySeq.Rnw:57-59
###################################################
  replicates <- c("simA", "simA", "simA", "simA", "simA",
                  "simB", "simB", "simB", "simB", "simB")


###################################################
### code chunk number 7: baySeq.Rnw:68-70
###################################################
groups <- list(NDE = c(1,1,1,1,1,1,1,1,1,1),
               DE = c(1,1,1,1,1,2,2,2,2,2))


###################################################
### code chunk number 8: baySeq.Rnw:81-82
###################################################
CD <- new("countData", data = simData, replicates = replicates, groups = groups)


###################################################
### code chunk number 9: baySeq.Rnw:87-88
###################################################
  libsizes(CD) <- getLibsizes(CD)


###################################################
### code chunk number 10: plotMA
###################################################
plotMA.CD(CD, samplesA = "simA", samplesB = "simB",
          col = c(rep("red", 100), rep("black", 900)))


###################################################
### code chunk number 11: figPlotMA
###################################################
plotMA.CD(CD, samplesA = "simA", samplesB = "simB",
          col = c(rep("red", 100), rep("black", 900)))


###################################################
### code chunk number 12: baySeq.Rnw:111-112
###################################################
CD@annotation <- data.frame(name = paste("count", 1:1000, sep = "_"))


###################################################
### code chunk number 13: baySeq.Rnw:119-120
###################################################
CD <- getPriors.NB(CD, samplesize = 1000, estimation = "QL", cl = cl)


###################################################
### code chunk number 14: baySeq.Rnw:127-131
###################################################
CD <- getLikelihoods(CD, cl = cl, bootStraps = 3, verbose = FALSE)
CD@estProps
CD@posteriors[1:10,]
CD@posteriors[101:110,]


###################################################
### code chunk number 15: baySeq.Rnw:135-136
###################################################
CD@estProps[2]


###################################################
### code chunk number 16: baySeq.Rnw:143-144
###################################################
topCounts(CD, group = "DE")  


###################################################
### code chunk number 17: plotPosteriors
###################################################
plotPosteriors(CD, group = "DE", col = c(rep("red", 100), rep("black", 900)))


###################################################
### code chunk number 18: figPlotPosteriors
###################################################
plotPosteriors(CD, group = "DE", col = c(rep("red", 100), rep("black", 900)))


###################################################
### code chunk number 19: baySeq.Rnw:175-177
###################################################
data(mobData)
data(mobAnnotation)


###################################################
### code chunk number 20: baySeq.Rnw:190-192
###################################################
seglens <- mobAnnotation$end - mobAnnotation$start + 1
cD <- new("countData", data = mobData, seglens = seglens, annotation = mobAnnotation)


###################################################
### code chunk number 21: baySeq.Rnw:196-197
###################################################
libsizes(cD) <- getLibsizes(cD, estimationType = "quantile")


###################################################
### code chunk number 22: baySeq.Rnw:205-206
###################################################
cDPair <- cD[,1:4]


###################################################
### code chunk number 23: baySeq.Rnw:210-211
###################################################
replicates(cDPair) <- as.factor(c("D3/D3", "D3/D3", "WT/D3", "WT/D3"))


###################################################
### code chunk number 24: baySeq.Rnw:217-218
###################################################
NDE <- c(1,1,1,1)


###################################################
### code chunk number 25: baySeq.Rnw:223-224
###################################################
mobile <- c("non-mobile","non-mobile","mobile","mobile")


###################################################
### code chunk number 26: baySeq.Rnw:229-230
###################################################
groups(cDPair) <- list(NDE = NDE, mobile = mobile)


###################################################
### code chunk number 27: baySeq.Rnw:235-236
###################################################
cDPair <- getPriors.NB(cDPair, samplesize = 1e4, cl = cl)


###################################################
### code chunk number 28: plotPriors
###################################################
plotNullPrior(cDPair)


###################################################
### code chunk number 29: figPlotPriors
###################################################
plotNullPrior(cDPair)


###################################################
### code chunk number 30: baySeq.Rnw:265-266
###################################################
cDPair <- getLikelihoods(cDPair, nullData = TRUE, cl = cl)


###################################################
### code chunk number 31: baySeq.Rnw:272-273
###################################################
cDPair


###################################################
### code chunk number 32: baySeq.Rnw:278-279
###################################################
summarisePosteriors(cD)


###################################################
### code chunk number 33: baySeq.Rnw:286-287
###################################################
topCounts(cDPair, group = 2, normaliseData = TRUE)


###################################################
### code chunk number 34: baySeq.Rnw:292-293
###################################################
topCounts(cDPair, group = NULL, number = 500)


###################################################
### code chunk number 35: plotPairPosteriors
###################################################
plotPosteriors(cDPair, group = 2, samplesA = 1:2, samplesB = 3:4)


###################################################
### code chunk number 36: figPlotPairPosteriors
###################################################
plotPosteriors(cDPair, group = 2, samplesA = 1:2, samplesB = 3:4)


###################################################
### code chunk number 37: plotMAPost
###################################################
plotMA.CD(cDPair, samplesA = c(1,2), samplesB = c(3,4),
          col = rgb(red = exp(cDPair@posteriors[,2]), green = 0, blue = 0))


###################################################
### code chunk number 38: figPlotMAPost
###################################################
plotMA.CD(cDPair, samplesA = c(1,2), samplesB = c(3,4),
          col = rgb(red = exp(cDPair@posteriors[,2]), green = 0, blue = 0))


###################################################
### code chunk number 39: baySeq.Rnw:334-335
###################################################
cD@replicates <- as.factor(c("D3/D3", "D3/D3", "WT/D3", "WT/D3", "WT/WT", "WT/WT"))


###################################################
### code chunk number 40: baySeq.Rnw:341-342
###################################################
NDE <- factor(c(1,1,1,1,1,1))


###################################################
### code chunk number 41: baySeq.Rnw:347-348
###################################################
d3dep <- c("wtRoot","wtRoot","wtRoot","wtRoot","dicerRoot","dicerRoot")


###################################################
### code chunk number 42: baySeq.Rnw:353-354
###################################################
mobile <- c("dicerShoot","dicerShoot","wtShoot","wtShoot","wtShoot","wtShoot")


###################################################
### code chunk number 43: baySeq.Rnw:359-360
###################################################
groups(cD) <- list(NDE = NDE, d3dep = d3dep, mobile = mobile)  


###################################################
### code chunk number 44: baySeq.Rnw:367-369
###################################################
cD <- getPriors.NB(cD, cl = cl)
cD <- getLikelihoods(cD, nullData = TRUE, cl = cl)


###################################################
### code chunk number 45: baySeq.Rnw:373-374
###################################################
topCounts(cD, group = "mobile", normaliseData = TRUE)  


###################################################
### code chunk number 46: baySeq.Rnw:378-379
###################################################
topCounts(cD, group = "d3dep", normaliseData = TRUE)  


###################################################
### code chunk number 47: baySeq.Rnw:386-387
###################################################
if(!is.null(cl)) stopCluster(cl)


###################################################
### code chunk number 48: baySeq.Rnw:393-394
###################################################
sessionInfo()