inst/doc/methylationAnalysis.R

In segmentSeq: Methods for identifying small RNA loci from high-throughput sequencing data

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

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


###################################################
### code chunk number 2: methylationAnalysis.Rnw:29-30
###################################################
  library(segmentSeq)


###################################################
### code chunk number 3: methylationAnalysis.Rnw:35-42
###################################################
if(require("parallel")) 
{
    numCores <- min(8, detectCores())
    cl <- makeCluster(numCores)
} else {
    cl <- NULL
}


###################################################
### code chunk number 4: methylationAnalysis.Rnw:45-46 (eval = FALSE)
###################################################
## cl <- NULL


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

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 6: methDist
###################################################
par(mfrow = c(2,1))
dists <- plotMethDistribution(mD, main = "Distributions of methylation", chr = "Chr1")
plotMethDistribution(mD, 
                     subtract = rowMeans(sapply(dists, function(x) x[,2])), 
                     main = "Differences between distributions", chr = "Chr1")


###################################################
### code chunk number 7: figMethDist
###################################################
par(mfrow = c(2,1))
dists <- plotMethDistribution(mD, main = "Distributions of methylation", chr = "Chr1")
plotMethDistribution(mD, 
                     subtract = rowMeans(sapply(dists, function(x) x[,2])), 
                     main = "Differences between distributions", chr = "Chr1")


###################################################
### code chunk number 8: methylationAnalysis.Rnw:88-89
###################################################
sD <- processAD(mD, cl = cl)


###################################################
### code chunk number 9: methylationAnalysis.Rnw:92-93
###################################################
if(nrow(sD) != 249271) stop("sD object is the wrong size (should have 249271 rows). Failure.")


###################################################
### code chunk number 10: methylationAnalysis.Rnw:102-105
###################################################
thresh = 0.2
hS <- heuristicSeg(sD = sD, aD = mD, prop = thresh, cl = cl, gap = 100, getLikes = FALSE)
hS


###################################################
### code chunk number 11: methylationAnalysis.Rnw:108-109
###################################################
if(nrow(hS) != 2955) stop("hS object is the wrong size (should have 2955 rows). Failure.")


###################################################
### code chunk number 12: methylationAnalysis.Rnw:115-116
###################################################
hS <- mergeMethSegs(hS, mD, gap = 5000, cl = cl)


###################################################
### code chunk number 13: methylationAnalysis.Rnw:121-122 (eval = FALSE)
###################################################
## hSL <- lociLikelihoods(hS, mD, cl = cl)


###################################################
### code chunk number 14: methylationAnalysis.Rnw:125-126
###################################################
data(hSL)


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


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


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


###################################################
### code chunk number 18: methylationAnalysis.Rnw:170-171
###################################################
hSL <- methObservables(hSL)


###################################################
### code chunk number 19: methylationAnalysis.Rnw:175-176
###################################################
densityFunction(hSL) <- bbNCDist


###################################################
### code chunk number 20: methylationAnalysis.Rnw:180-181
###################################################
hSL <- getPriors(hSL, cl = cl)


###################################################
### code chunk number 21: methylationAnalysis.Rnw:186-187
###################################################
hSL <- getLikelihoods(hSL, cl = cl)


###################################################
### code chunk number 22: methylationAnalysis.Rnw:192-193
###################################################
topCounts(hSL, "DE")


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


###################################################
### code chunk number 24: methylationAnalysis.Rnw:206-207
###################################################
sessionInfo()