###################################################
### chunk number 1: 
###################################################
set.seed(102)


###################################################
### chunk number 2: 
###################################################
library(baySeq)


###################################################
### chunk number 3: 
###################################################
if("snow" %in% installed.packages()[,1])
  {
    library(snow)
    cl <- makeCluster(4, "SOCK")
    clusterEvalQ(cl, library(baySeq))
  } else cl <- NULL


###################################################
### chunk number 4: 
###################################################
data(simCount)
data(libsizes)
simCount[1:10,]
libsizes


###################################################
### chunk number 5: 
###################################################
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))


###################################################
### chunk number 6: 
###################################################
CD <- new("countData", data = simCount, libsizes = libsizes, groups = groups)


###################################################
### chunk number 7: 
###################################################
CD@annotation <- data.frame(name = paste("count", 1:1000, sep = "_"))


###################################################
### chunk number 8: 
###################################################
CDP.Poi <- getPriors.Pois(CD, samplesize = 20, iterations = 1000, takemean = TRUE)


###################################################
### chunk number 9: 
###################################################
  CDP.Poi@priors


###################################################
### chunk number 10: 
###################################################
CDPost.Poi <- getLikelihoods.Pois(CDP.Poi, prs = c(0.5, 0.5), estimatePriors = TRUE, cl = cl)
CDPost.Poi@estProps
CDPost.Poi@posteriors[1:10,]
CDPost.Poi@posteriors[101:110,]


###################################################
### chunk number 11: 
###################################################
CDPost.Poi@estProps[2]


###################################################
### chunk number 12: 
###################################################
CDP.NBML <- getPriors.NB(CD, samplesize = 1000, estimation = "ML", cl = cl)


###################################################
### chunk number 13: 
###################################################
CDPost.NBML <- getLikelihoods.NBboot(CDP.NBML, prs = c(0.5,0.5), estimatePriors = TRUE, bootStraps = 2, cl = cl)
CDPost.NBML@estProps
CDPost.NBML@posteriors[1:10,]
CDPost.NBML@posteriors[101:110,]


###################################################
### chunk number 14: 
###################################################
CDPost.NBML@estProps[2]


###################################################
### chunk number 15: 
###################################################
topCounts(CDPost.NBML, group = 2)  


###################################################
### chunk number 16: 
###################################################
NBML.FPs <- 1:nrow(simCount) - getTPs(CDPost.NBML, 1, TPs = 1:100, decreasing = FALSE)
Poi.FPs <- 1:nrow(simCount) - getTPs(CDPost.Poi, 1, TPs = 1:100, decreasing = FALSE)
plot(x = NA, y = NA, xlim = c(0, 100), ylim = c(0, log(1000)), xlab = "Number of counts selected", ylab = "(Log) False Positives")
lines(x = 1:1000, y = log(NBML.FPs[1:1000]), type = "l", col = "red")
lines(x = 1:1000, y = log(Poi.FPs[1:1000]), type = "l", col = "blue")


###################################################
### chunk number 17: 
###################################################
if(!is.null(cl))
  stopCluster(cl)


###################################################
### chunk number 18: 
###################################################
set.seed(101)


###################################################
### chunk number 19: 
###################################################
library(baySeq)
if("snow" %in% installed.packages()[,1])
  {
    library(snow)
    cl <- makeCluster(4, "SOCK")
    clusterEvalQ(cl, library(baySeq))
  } else cl <- NULL


###################################################
### chunk number 20: 
###################################################
data(factData)
data(factlibsizes)


###################################################
### chunk number 21: 
###################################################
factgroups <- list(NDE = c(1,1,1,1,1,1,1,1),
                   DE.A.B = c(1,1,1,1,2,2,2,2),
                   DE.C.D = c(1,1,2,2,1,1,2,2))


###################################################
### chunk number 22: 
###################################################
CDfact <- new("countData", data = factCount, libsizes = factlibsizes, groups = factgroups)
CDfact@annotation <- data.frame(name = paste("count", 1:1000, sep = "_"))
CDfactP.NBML <- getPriors.NB(CDfact, samplesize = 1000, estimation = "ML", cl = cl)
CDfactPost.NBML <- getLikelihoods.NBboot(CDfactP.NBML, prs = c(0.2, 0.3,0.5), estimatePriors = TRUE, bootStraps = 2, cl = cl)
CDfactPost.NBML@estProps


###################################################
### chunk number 23: 
###################################################
topCounts(CDfactPost.NBML, group = 2)  


###################################################
### chunk number 24: 
###################################################
topCounts(CDfactPost.NBML, group = 3)