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

###################################################
### code chunk number 1: bsseq_analysis.Rnw:6-7
###################################################
options(width=70)


###################################################
### code chunk number 2: load
###################################################
library(bsseq)
library(bsseqData)


###################################################
### code chunk number 3: showData
###################################################
data(BS.cancer.ex)
BS.cancer.ex <- updateObject(BS.cancer.ex)
BS.cancer.ex
pData(BS.cancer.ex)


###################################################
### code chunk number 4: citation
###################################################
print(citation("bsseq"), style = "latex")


###################################################
### code chunk number 5: smooth (eval = FALSE)
###################################################
## BS.cancer.ex.fit <- BSmooth(BS.cancer.ex, mc.cores = 1, verbose = TRUE)


###################################################
### code chunk number 6: showDataFit
###################################################
data(BS.cancer.ex.fit)
BS.cancer.ex.fit <- updateObject(BS.cancer.ex.fit)
BS.cancer.ex.fit


###################################################
### code chunk number 7: cpgNumbers
###################################################
## The average coverage of CpGs on the two chromosomes
round(colMeans(getCoverage(BS.cancer.ex)), 1)
## Number of CpGs in two chromosomes
length(BS.cancer.ex)
## Number of CpGs which are covered by at least 1 read in all 6 samples
sum(rowSums(getCoverage(BS.cancer.ex) >= 1) == 6)
## Number of CpGs with 0 coverage in all samples
sum(rowSums(getCoverage(BS.cancer.ex)) == 0)


###################################################
### code chunk number 8: poisson
###################################################
logp <- ppois(0, lambda = 4, lower.tail = FALSE, log.p = TRUE)
round(1 - exp(6 * logp), 3)


###################################################
### code chunk number 9: smoothSplit (eval = FALSE)
###################################################
## ## Split datag
## BS1 <- BS.cancer.ex[, 1]
## save(BS1, file = "BS1.rda")
## BS2 <- BS.cancer.ex[, 2]
## save(BS1, file = "BS1.rda")
## ## done splitting
## 
## ## Do the following on each node
## 
## ## node 1
## load("BS1.rda")
## BS1.fit <- BSmooth(BS1)
## save(BS1.fit)
## save(BS1.fit, file = "BS1.fit.rda")
## ## done node 1
## 
## ## node 2
## load("BS2.rda")
## BS2.fit <- BSmooth(BS2)
## save(BS2.fit, file = "BS2.fit.rda")
## ## done node 2
## 
## ## join; in a new R session
## load("BS1.fit.rda")
## load("BS2.fit.rda")
## BS.fit <- combine(BS1.fit, BS2.fit)


###################################################
### code chunk number 10: keepLoci
###################################################
BS.cov <- getCoverage(BS.cancer.ex.fit)
keepLoci.ex <- which(rowSums(BS.cov[, BS.cancer.ex$Type == "cancer"] >= 2) >= 2 &
                     rowSums(BS.cov[, BS.cancer.ex$Type == "normal"] >= 2) >= 2)
length(keepLoci.ex)
BS.cancer.ex.fit <- BS.cancer.ex.fit[keepLoci.ex,]


###################################################
### code chunk number 11: BSmooth.tstat
###################################################
BS.cancer.ex.tstat <- BSmooth.tstat(BS.cancer.ex.fit, 
                                    group1 = c("C1", "C2", "C3"),
                                    group2 = c("N1", "N2", "N3"), 
                                    estimate.var = "group2",
                                    local.correct = TRUE,
                                    verbose = TRUE)
BS.cancer.ex.tstat


###################################################
### code chunk number 12: plotTstat
###################################################
plot(BS.cancer.ex.tstat)


###################################################
### code chunk number 13: dmrs
###################################################
dmrs0 <- dmrFinder(BS.cancer.ex.tstat, cutoff = c(-4.6, 4.6))
dmrs <- subset(dmrs0, n >= 3 & abs(meanDiff) >= 0.1)
nrow(dmrs)
head(dmrs, n = 3)


###################################################
### code chunk number 14: plotSetup
###################################################
pData <- pData(BS.cancer.ex.fit)
pData$col <- rep(c("red", "blue"), each = 3)
pData(BS.cancer.ex.fit) <- pData


###################################################
### code chunk number 15: plotDmr
###################################################
plotRegion(BS.cancer.ex.fit, dmrs[1,], extend = 5000, addRegions = dmrs)


###################################################
### code chunk number 16: plotManyRegions (eval = FALSE)
###################################################
## pdf(file = "dmrs_top200.pdf", width = 10, height = 5)
## plotManyRegions(BS.cancer.ex.fit, dmrs[1:200,], extend = 5000, 
##                 addRegions = dmrs)
## dev.off()


###################################################
### code chunk number 17: sessionInfo
###################################################
toLatex(sessionInfo())