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

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


###################################################
### code chunk number 2: segmentSeq.Rnw:31-32
###################################################
  library(segmentSeq)


###################################################
### code chunk number 3: segmentSeq.Rnw:37-43
###################################################
if(FALSE) { # set to FALSE if you don't want parallelisation
    numCores <- min(8, detectCores())
    cl <- makeCluster(numCores)
} else {
    cl <- NULL
}


###################################################
### code chunk number 4: segmentSeq.Rnw:46-47
###################################################
set.seed(1)


###################################################
### code chunk number 5: segmentSeq.Rnw:55-64
###################################################
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 6: segmentSeq.Rnw:69-71
###################################################
sD <- processAD(aD, cl = cl)
sD


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


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


###################################################
### code chunk number 9: segmentSeq.Rnw:89-91
###################################################
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, 208)) > 2)) stop("hS object contains wrong number of loci. Likely failure.")


###################################################
### code chunk number 10: segmentSeq.Rnw:103-105
###################################################
cS <- classifySeg(sD = sD, aD = aD, cD = hS, cl = cl)      
cS


###################################################
### code chunk number 11: segmentSeq.Rnw:108-110
###################################################
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,36)) > 2)) stop("cS object contains wrong number of loci. Likely failure.")


###################################################
### code chunk number 12: 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 13: 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 14: segmentSeq.Rnw:142-144
###################################################
loci <- selectLoci(cS, FDR = 0.05)
loci


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


###################################################
### code chunk number 16: segmentSeq.Rnw:155-156
###################################################
cS <- getPriors(cS, cl = cl)


###################################################
### code chunk number 17: segmentSeq.Rnw:160-161
###################################################
cS <- getLikelihoods(cS, nullData = TRUE, cl = cl)


###################################################
### code chunk number 18: segmentSeq.Rnw:166-167
###################################################
topCounts(cS, NULL, number = 3)


###################################################
### code chunk number 19: segmentSeq.Rnw:172-173
###################################################
topCounts(cS, "NDE", number = 3)


###################################################
### code chunk number 20: segmentSeq.Rnw:178-179
###################################################
topCounts(cS, "DE", number = 3)


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


###################################################
### code chunk number 22: segmentSeq.Rnw:191-192
###################################################
sessionInfo()