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

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


###################################################
### code chunk number 2: segmentSeq.Rnw:29-30
###################################################
  library(segmentSeq)


###################################################
### code chunk number 3: segmentSeq.Rnw:35-42
###################################################
if(require("parallel")) 
{
    numCores <- min(8, detectCores())
    cl <- makeCluster(numCores)
} else {
    cl <- NULL
}


###################################################
### code chunk number 4: segmentSeq.Rnw:45-46 (eval = FALSE)
###################################################
## cl <- NULL


###################################################
### code chunk number 5: segmentSeq.Rnw:49-50
###################################################
set.seed(1)


###################################################
### code chunk number 6: segmentSeq.Rnw:58-67
###################################################
datadir <- system.file("extdata", package = "segmentSeq")
libfiles <- c("SL9.txt", "SL10.txt", "SL26.txt", "SL32.txt")
libnames <- c("SL9", "SL10", "SL26", "SL32")
replicates <- c("AGO6", "AGO6", "AGO4", "AGO4")

aD <- readGeneric(files = libfiles, dir = datadir,
                  replicates = replicates, libnames = libnames,                   
                  polyLength = 10, header = TRUE, gap = 200)
aD


###################################################
### code chunk number 7: segmentSeq.Rnw:72-74
###################################################
sD <- processAD(aD, cl = cl)
sD


###################################################
### code chunk number 8: segmentSeq.Rnw:78-79
###################################################
if(nrow(sD) != 1452) stop("sD object is the wrong size (should have 1452 rows). Failure.")


###################################################
### code chunk number 9: segmentSeq.Rnw:88-89
###################################################
hS <- heuristicSeg(sD = sD, aD = aD, RKPM = 1000, largeness = 1e8, getLikes = TRUE, cl = cl)


###################################################
### code chunk number 10: segmentSeq.Rnw:92-94
###################################################
if(nrow(hS) != 507) stop("hS object is the wrong size (should have 507 rows). Failure.")
if(any(abs(colSums(exp(hS@locLikelihoods)) - c(88, 227)) > 2)) stop("hS object contains wrong number of loci. Likely failure.")


###################################################
### code chunk number 11: segmentSeq.Rnw:106-108
###################################################
cS <- classifySeg(sD = sD, aD = aD, cD = hS, cl = cl)      
cS


###################################################
### code chunk number 12: segmentSeq.Rnw:111-113
###################################################
if(abs(nrow(cS) - 64) > 2) stop("cS object is the wrong size (should have ~142 rows). Likely failure.")
if(any(abs(colSums(exp(cS@locLikelihoods)) - c(29,31)) > 2)) stop("cS object contains wrong number of loci. Likely failure.")


###################################################
### code chunk number 13: segPlot
###################################################
par(mfrow = c(2,1), mar = c(2,6,2,2))
plotGenome(aD, hS, chr = ">Chr1", limits = c(1, 1e5),
           showNumber = FALSE, cap = 50)
plotGenome(aD, cS, chr = ">Chr1", limits = c(1, 1e5),
           showNumber = FALSE, cap = 50)


###################################################
### code chunk number 14: figSeg
###################################################
par(mfrow = c(2,1), mar = c(2,6,2,2))
plotGenome(aD, hS, chr = ">Chr1", limits = c(1, 1e5),
           showNumber = FALSE, cap = 50)
plotGenome(aD, cS, chr = ">Chr1", limits = c(1, 1e5),
           showNumber = FALSE, cap = 50)


###################################################
### code chunk number 15: segmentSeq.Rnw:145-147
###################################################
loci <- selectLoci(cS, FDR = 0.05)
loci


###################################################
### code chunk number 16: segmentSeq.Rnw:153-154
###################################################
groups(cS) <- list(NDE = c(1,1,1,1), DE = c("AGO6", "AGO6", "AGO4", "AGO4"))


###################################################
### code chunk number 17: segmentSeq.Rnw:158-159
###################################################
cS <- getPriors(cS, cl = cl)


###################################################
### code chunk number 18: segmentSeq.Rnw:163-164
###################################################
cS <- getLikelihoods(cS, nullData = TRUE, cl = cl)


###################################################
### code chunk number 19: segmentSeq.Rnw:169-170
###################################################
topCounts(cS, NULL, number = 3)


###################################################
### code chunk number 20: segmentSeq.Rnw:175-176
###################################################
topCounts(cS, "NDE", number = 3)


###################################################
### code chunk number 21: segmentSeq.Rnw:181-182
###################################################
topCounts(cS, "DE", number = 3)


###################################################
### code chunk number 22: <stopCluster
###################################################
if(!is.null(cl))
    stopCluster(cl)


###################################################
### code chunk number 23: segmentSeq.Rnw:194-195
###################################################
sessionInfo()