### R code from vignette source 'vignettes/segmentSeq/inst/doc/methylationAnalysis.Rnw'

###################################################
### code chunk number 1: methylationAnalysis.Rnw:23-24
###################################################
  library(segmentSeq)


###################################################
### code chunk number 2: methylationAnalysis.Rnw:29-30 (eval = FALSE)
###################################################
## cl <- makeCluster(8)


###################################################
### code chunk number 3: methylationAnalysis.Rnw:34-35
###################################################
  cl <- NULL


###################################################
### code chunk number 4: methylationAnalysis.Rnw:40-49
###################################################
datadir <- system.file("extdata", package = "segmentSeq")
files <- c("short_18B_C24_C24_trim.fastq_CG_methCalls",
"short_Sample_17A_trimmed.fastq_CG_methCalls",
"short_13_C24_col_trim.fastq_CG_methCalls",
"short_Sample_28_trimmed.fastq_CG_methCalls")

mD <- readMeths(files = files, dir = datadir,
libnames = c("A1", "A2", "B1", "B2"), replicates = c("A","A","B","B"),
nonconversion = c(0.004777, 0.005903, 0.016514, 0.006134))


###################################################
### code chunk number 5: methDist
###################################################
par(mfrow = c(2,1))
distA <- plotMethDistribution(mD, samples = 1:2, main = "Distributions of methylation")
distB <- plotMethDistribution(mD, samples = 3:4, add = TRUE, col = "red")
plotMethDistribution(mD, samples = 1:2, subtract = distB[,2], main = "Differences between distributions")


###################################################
### code chunk number 6: figMethDist
###################################################
par(mfrow = c(2,1))
distA <- plotMethDistribution(mD, samples = 1:2, main = "Distributions of methylation")
distB <- plotMethDistribution(mD, samples = 3:4, add = TRUE, col = "red")
plotMethDistribution(mD, samples = 1:2, subtract = distB[,2], main = "Differences between distributions")


###################################################
### code chunk number 7: methylationAnalysis.Rnw:76-77
###################################################
sD <- processAD(mD, gap = 300, squeeze = 10, filterProp = 0.05, verbose = TRUE, strandSplit = TRUE, cl = cl)


###################################################
### code chunk number 8: methylationAnalysis.Rnw:86-88
###################################################
hS <- heuristicSeg(sD, mD, prop = 0.2, cl = cl, gap = 100, getLikes = FALSE)
hS


###################################################
### code chunk number 9: methylationAnalysis.Rnw:93-94
###################################################
hS <- mergeMethSegs(hS, mD, gap = 5000, cl = cl)


###################################################
### code chunk number 10: methylationAnalysis.Rnw:99-100
###################################################
hSL <- lociLikelihoods(hS, mD, cl = cl)


###################################################
### code chunk number 11: plotMeth (eval = FALSE)
###################################################
## plotMeth(mD, hSL, chr = "Chr1", limits = c(1, 50000), cap = 10)


###################################################
### code chunk number 12: figplotMeth
###################################################
plotMeth(mD, hSL, chr = "Chr1", limits = c(1, 50000), cap = 10)


###################################################
### code chunk number 13: methylationAnalysis.Rnw:139-140
###################################################
groups(hSL) <- list(NDE = c(1,1,1,1), DE = c("A", "A", "B", "B"))


###################################################
### code chunk number 14: methylationAnalysis.Rnw:145-146
###################################################
hSL <- getPriors.BB(hSL, cl = cl)


###################################################
### code chunk number 15: methylationAnalysis.Rnw:151-152
###################################################
hSL <- getLikelihoods.BB(hSL, cl = cl)


###################################################
### code chunk number 16: methylationAnalysis.Rnw:157-158
###################################################
topCounts(hSL, "DE")


###################################################
### code chunk number 17: methylationAnalysis.Rnw:161-162 (eval = FALSE)
###################################################
## stopCluster(cl)