###################################################
### 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")
  } else cl <- NULL


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


###################################################
### chunk number 5: 
###################################################
  replicates <- c(1,1,1,1,1,2,2,2,2,2)


###################################################
### chunk number 6: 
###################################################
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 7: 
###################################################
CD <- new("countData", data = simCount, replicates = replicates, libsizes = libsizes, groups = groups)


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


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


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


###################################################
### chunk number 11: 
###################################################
CDPost.Poi <- getLikelihoods.Pois(CDP.Poi, pET = "BIC", cl = cl)
CDPost.Poi@estProps
CDPost.Poi@posteriors[1:10,]
CDPost.Poi@posteriors[101:110,]


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


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


###################################################
### chunk number 14: 
###################################################
CDPost.NBML <- getLikelihoods.NB(CDP.NBML, pET = 'BIC', cl = cl)
CDPost.NBML@estProps
CDPost.NBML@posteriors[1:10,]
CDPost.NBML@posteriors[101:110,]


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


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


###################################################
### chunk number 17: 
###################################################
NBML.TPs <- getTPs(CDPost.NBML, group = 2, TPs = 1:100)
Poi.TPs <- getTPs(CDPost.Poi, group = 2, TPs = 1:100)


###################################################
### chunk number 18: FPsPlot
###################################################
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(1:1000 - NBML.TPs[1:1000]), type = "l", col = "red")
lines(x = 1:1000, y = log(1:1000 - Poi.TPs[1:1000]), type = "l", col = "blue")

legend(x = "topleft", lty = c(1,1), col = c("red", "blue"), legend = c("Negative Binomial", "Poisson-Gamma"))


###################################################
### chunk number 19: figFPs
###################################################
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(1:1000 - NBML.TPs[1:1000]), type = "l", col = "red")
lines(x = 1:1000, y = log(1:1000 - Poi.TPs[1:1000]), type = "l", col = "blue")

legend(x = "topleft", lty = c(1,1), col = c("red", "blue"), legend = c("Negative Binomial", "Poisson-Gamma"))


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


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


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


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


###################################################
### chunk number 24: 
###################################################
  replicates <- c(1,1,2,2,3,3,4,4)
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 25: 
###################################################
CDfact <- new("countData", data = factCount, replicates = replicates, libsizes = factlibsizes, groups = factgroups)
CDfact@annotation <- data.frame(name = paste("count", 1:1000, sep = "_"))
CDfactP.NBML <- getPriors.NB(CDfact, samplesize = 1000, estimation = "QL", cl = cl)
CDfactPost.NBML <- getLikelihoods.NB(CDfactP.NBML, pET = "BIC", cl = cl)
CDfactPost.NBML@estProps


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


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


###################################################
### chunk number 28: 
###################################################
data(simSeg)
data(libsizes)
simSeg[1:10,]
libsizes


###################################################
### chunk number 29: 
###################################################
replicates <- c(1,1,1,1,1,2,2,2,2,2)
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 30: 
###################################################
SD <- new("countData", data = simSeg[,-1], replicates = replicates, libsizes = libsizes, groups = groups, seglens = simSeg[,1])


###################################################
### chunk number 31: 
###################################################
SD@annotation <- data.frame(name = paste("gene", 1:1000, sep = "_"))


###################################################
### chunk number 32: 
###################################################
SDP.NBML <- getPriors.NB(SD, samplesize = 1000, estimation = "QL", cl = cl)
SDP.Pois <- getPriors.Pois(SD, samplesize = 20, cl = cl)


###################################################
### chunk number 33: 
###################################################
SDPost.Pois <- getLikelihoods.Pois(SDP.Pois, pET = "BIC", cl = cl)  
SDPost.NBML <- getLikelihoods.NB(SDP.NBML, pET = "BIC", cl = cl)


###################################################
### chunk number 34: 
###################################################
CSD <- new("countData", data = simSeg[,-1], replicates = replicates, libsizes = libsizes, groups = groups)
CSD@annotation <- data.frame(name = paste("gene", 1:1000, sep = "_"))

CSDP.NBML <- getPriors.NB(CSD, samplesize = 1000, estimation = "QL", cl = cl)
CSDPost.NBML <- getLikelihoods.NB(CSDP.NBML, pET = "BIC", cl = cl)
CSDP.Pois <- getPriors.Pois(CSD, samplesize = 20, cl = cl)
CSDPost.Pois <- getLikelihoods.Pois(CSDP.Pois, pET = "BIC", cl = cl)


###################################################
### chunk number 35: FPseglen
###################################################
SD.NBML.FPs <- 1:nrow(simSeg) - getTPs(SDPost.NBML, group = 2, TPs = 1:100)
CSD.NBML.FPs <- 1:nrow(simSeg) - getTPs(CSDPost.NBML, group = 2, TPs = 1:100)
SD.Pois.FPs <- 1:nrow(simSeg) - getTPs(SDPost.Pois, group = 2, TPs = 1:100)
CSD.Pois.FPs <- 1:nrow(simSeg) - getTPs(CSDPost.Pois, group = 2, TPs = 1:100)

plot(x = NA, y = NA, xlim = c(0, 100), ylim = c(0, log(100)), xlab = "Number of counts selected", ylab = "(Log) False Positives")
lines(x = 1:1000, y = log(CSD.NBML.FPs[1:1000]), type = "l", col = "red")
lines(x = 1:1000, y = log(SD.NBML.FPs[1:1000]), type = "l", col = "red", lty = 2)
lines(x = 1:1000, y = log(CSD.Pois.FPs[1:1000]), type = "l", col = "blue")
lines(x = 1:1000, y = log(SD.Pois.FPs[1:1000]), type = "l", col = "blue", lty = 2)
legend(x = "topleft", lty = c(1,2,1,2), col = c("red", "red", "blue", "blue"), legend = c("Negative Binomial (ignoring segment lengths)", "Negative Binomial (including segment lengths)", "Poisson-Gamma (ignoring segment lengths)", "Poisson-Gamma (including segment lengths)"))
                                                                               


###################################################
### chunk number 36: 
###################################################
NSDPost.NBML <- getLikelihoods.NB(SDP.NBML, pET = "BIC", nullData = TRUE, bootStraps = 1, cl = cl)
NCSDPost.NBML <- getLikelihoods.NB(CSDP.NBML, pET = "BIC", nullData = TRUE, bootStraps = 1, cl = cl)


###################################################
### chunk number 37: 
###################################################
topCounts(NSDPost.NBML, group = NULL)
topCounts(NCSDPost.NBML, group = NULL)


###################################################
### chunk number 38: FPNull
###################################################
NSD.FPs <- 1:nrow(simSeg) - getTPs(NSDPost.NBML, group = NULL, TPs = 101:200, decreasing = TRUE)
NCSD.FPs <- 1:nrow(simSeg) - getTPs(NCSDPost.NBML, group = NULL, TPs = 101:200, decreasing = TRUE)
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(NCSD.FPs[1:1000]), type = "l", col = "red")
lines(x = 1:1000, y = log(NSD.FPs[1:1000]), type = "l", lty = 2, col = "red")
legend(x = "topleft", lty = c(1,2), col = c("red", "red"), legend = c("Negative Binomial (ignoring segment lengths)", "Negative Binomial (including segment lengths)"))


###################################################
### chunk number 39: pltPriors
###################################################
plotPriors(SDP.NBML, group = 1)


###################################################
### chunk number 40: figPriors
###################################################
plotPriors(SDP.NBML, group = 1)