### R code from vignette source 'vignettes/bsseq/inst/doc/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())