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

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


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


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


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


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


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


###################################################
### code chunk number 7: baySeq.Rnw:67-69
###################################################
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:80-81
###################################################
CD <- new("countData", data = simData, replicates = replicates, groups = groups)


###################################################
### code chunk number 9: baySeq.Rnw:86-87
###################################################
  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:110-111
###################################################
CD@annotation <- data.frame(name = paste("count", 1:1000, sep = "_"))


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


###################################################
### code chunk number 14: baySeq.Rnw:126-130
###################################################
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:134-135
###################################################
CD@estProps[2]


###################################################
### code chunk number 16: baySeq.Rnw:142-143
###################################################
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:174-176
###################################################
data(mobData)
data(mobAnnotation)


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


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


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


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


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


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


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


###################################################
### code chunk number 27: baySeq.Rnw:234-235
###################################################
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:264-265
###################################################
cDPair <- getLikelihoods(cDPair, nullData = TRUE, cl = cl)


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


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


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


###################################################
### code chunk number 34: baySeq.Rnw:291-292
###################################################
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:333-334
###################################################
cD@replicates <- as.factor(c("D3/D3", "D3/D3", "WT/D3", "WT/D3", "WT/WT", "WT/WT"))


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


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


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


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


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


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


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


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


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