### R code from vignette source 'vignettes/baySeq/inst/doc/baySeq.Rnw'

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


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


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


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


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


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


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


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


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


###################################################
### code chunk number 14: baySeq.Rnw:128-132
###################################################
CD <- getLikelihoods(CD, pET = 'BIC', cl = cl)
CD@estProps
CD@posteriors[1:10,]
CD@posteriors[101:110,]


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


###################################################
### code chunk number 16: baySeq.Rnw:144-145
###################################################
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:172-173
###################################################
  data(pairData)  


###################################################
### code chunk number 20: baySeq.Rnw:177-181
###################################################
pairCD <- new("countData", data = list(pairData[,1:4], pairData[,5:8]),
                 replicates = c(1,1,2,2),
                 groups = list(NDE = c(1,1,1,1), DE = c(1,1,2,2)),
              densityFunction = bbDensity)


###################################################
### code chunk number 21: baySeq.Rnw:185-186
###################################################
libsizes(pairCD) <- getLibsizes(pairCD)


###################################################
### code chunk number 22: baySeq.Rnw:191-192
###################################################
pairCD <- getPriors(pairCD, samplesize = 1000, cl = cl)


###################################################
### code chunk number 23: baySeq.Rnw:197-198
###################################################
  pairCD <- getLikelihoods(pairCD, pET = 'BIC', nullData = TRUE, cl = cl)


###################################################
### code chunk number 24: baySeq.Rnw:203-204
###################################################
  topCounts(pairCD, group = 2)


###################################################
### code chunk number 25: baySeq.Rnw:207-208
###################################################
  topCounts(pairCD, group = 1)


###################################################
### code chunk number 26: baySeq.Rnw:220-222
###################################################
data(mobData)
data(mobAnnotation)


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


###################################################
### code chunk number 28: baySeq.Rnw:241-242
###################################################
libsizes(cD) <- getLibsizes(cD, estimationType = "quantile")


###################################################
### code chunk number 29: baySeq.Rnw:250-251
###################################################
cDPair <- cD[,1:4]


###################################################
### code chunk number 30: baySeq.Rnw:255-256
###################################################
replicates(cDPair) <- as.factor(c("D3/D3", "D3/D3", "WT/D3", "WT/D3"))


###################################################
### code chunk number 31: baySeq.Rnw:262-263
###################################################
NDE <- c(1,1,1,1)


###################################################
### code chunk number 32: baySeq.Rnw:268-269
###################################################
mobile <- c("non-mobile","non-mobile","mobile","mobile")


###################################################
### code chunk number 33: baySeq.Rnw:274-275
###################################################
groups(cDPair) <- list(NDE = NDE, mobile = mobile)


###################################################
### code chunk number 34: baySeq.Rnw:280-281
###################################################
cDPair <- getPriors.NB(cDPair, samplesize = 1e4, cl = cl)


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


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


###################################################
### code chunk number 37: baySeq.Rnw:310-311
###################################################
cDPair <- getLikelihoods(cDPair, nullData = TRUE, cl = cl)


###################################################
### code chunk number 38: baySeq.Rnw:317-318
###################################################
cDPair


###################################################
### code chunk number 39: baySeq.Rnw:323-324
###################################################
summarisePosteriors(cD)


###################################################
### code chunk number 40: baySeq.Rnw:331-332
###################################################
topCounts(cDPair, group = 2, normaliseData = TRUE)


###################################################
### code chunk number 41: baySeq.Rnw:337-338
###################################################
topCounts(cDPair, group = NULL, number = 500)


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


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


###################################################
### code chunk number 44: 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 45: 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 46: baySeq.Rnw:377-378
###################################################
cD@replicates <- as.factor(c("D3/D3", "D3/D3", "WT/D3", "WT/D3", "WT/WT", "WT/WT"))


###################################################
### code chunk number 47: baySeq.Rnw:384-385
###################################################
NDE <- factor(c(1,1,1,1,1,1))


###################################################
### code chunk number 48: baySeq.Rnw:390-391
###################################################
d3dep <- c("wtRoot","wtRoot","wtRoot","wtRoot","dicerRoot","dicerRoot")


###################################################
### code chunk number 49: baySeq.Rnw:396-397
###################################################
mobile <- c("dicerShoot","dicerShoot","wtShoot","wtShoot","wtShoot","wtShoot")


###################################################
### code chunk number 50: baySeq.Rnw:402-403
###################################################
groups(cD) <- list(NDE = NDE, d3dep = d3dep, mobile = mobile)  


###################################################
### code chunk number 51: baySeq.Rnw:410-412
###################################################
cD <- getPriors.NB(cD, cl = cl)
cD <- getLikelihoods(cD, nullData = TRUE, cl = cl)


###################################################
### code chunk number 52: baySeq.Rnw:416-417
###################################################
topCounts(cD, group = "mobile", normaliseData = TRUE)  


###################################################
### code chunk number 53: baySeq.Rnw:421-422
###################################################
topCounts(cD, group = "d3dep", normaliseData = TRUE)  


###################################################
### code chunk number 54: baySeq.Rnw:429-430
###################################################
if(!is.null(cl)) stopCluster(cl)


###################################################
### code chunk number 55: baySeq.Rnw:436-437
###################################################
sessionInfo()