################################################### ### 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)