### 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:49-61
###################################################
chrlens <- c(1e6, 2e5)

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, 
                  chrs = c(">Chr1", ">Chr2"), chrlens = chrlens,
                  polyLength = 10, header = TRUE, gap = 200)
aD


###################################################
### code chunk number 5: segmentSeq.Rnw:66-68
###################################################
sD <- processAD(aD, gap = 100, cl = cl)
sD


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


###################################################
### code chunk number 7: segmentSeq.Rnw:90-92
###################################################
classSegs <- classifySeg(sD = sD, aD = aD, cD = hS, cl = cl)      
classSegs


###################################################
### code chunk number 8: 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, classSegs, chr = ">Chr1", limits = c(1, 1e5),
           showNumber = FALSE, cap = 50)


###################################################
### code chunk number 9: 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, classSegs, chr = ">Chr1", limits = c(1, 1e5),
           showNumber = FALSE, cap = 50)


###################################################
### code chunk number 10: segmentSeq.Rnw:123-125
###################################################
loci <- selectLoci(classSegs, FDR = 0.05)
loci


###################################################
### code chunk number 11: segmentSeq.Rnw:131-132
###################################################
groups(classSegs) <- list(NDE = c(1,1,1,1), DE = c("AGO6", "AGO6", "AGO4", "AGO4"))


###################################################
### code chunk number 12: segmentSeq.Rnw:136-137
###################################################
classSegs <- getPriors(classSegs, cl = cl)


###################################################
### code chunk number 13: segmentSeq.Rnw:141-142
###################################################
classSegs <- getLikelihoods(classSegs, nullData = TRUE, cl = cl)


###################################################
### code chunk number 14: segmentSeq.Rnw:147-148
###################################################
topCounts(classSegs, NULL, number = 3)


###################################################
### code chunk number 15: segmentSeq.Rnw:153-154
###################################################
topCounts(classSegs, "NDE", number = 3)


###################################################
### code chunk number 16: segmentSeq.Rnw:159-160
###################################################
topCounts(classSegs, "DE", number = 3)


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


###################################################
### code chunk number 18: segmentSeq.Rnw:172-173
###################################################
sessionInfo()