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

###################################################
### code chunk number 1: install (eval = FALSE)
###################################################
## source("http://www.bioconductor.org/biocLite.R")
## biocLite("BSgenome")


###################################################
### code chunk number 2: genomes (eval = FALSE)
###################################################
## library(BSgenome)
## available.genomes()


###################################################
### code chunk number 3: install (eval = FALSE)
###################################################
## source("http://www.bioconductor.org/biocLite.R")
## biocLite("BSgenome.Hsapiens.UCSC.hg18")


###################################################
### code chunk number 4: MethylSeekR.Rnw:95-96
###################################################
library(MethylSeekR)


###################################################
### code chunk number 5: MethylSeekR.Rnw:105-106
###################################################
set.seed(123)


###################################################
### code chunk number 6: MethylSeekR.Rnw:121-123 (eval = FALSE)
###################################################
## system.file("extdata", "Lister2009_imr90_hg18_chr22.tab",
## package="MethylSeekR")


###################################################
### code chunk number 7: MethylSeekR.Rnw:133-136
###################################################
library("BSgenome.Hsapiens.UCSC.hg18")
sLengths=seqlengths(Hsapiens)
head(sLengths)


###################################################
### code chunk number 8: MethylSeekR.Rnw:152-156
###################################################
methFname <- system.file("extdata", 
"Lister2009_imr90_hg18_chr22.tab", package="MethylSeekR")
meth.gr <- readMethylome(FileName=methFname, seqLengths=sLengths)
head(meth.gr)


###################################################
### code chunk number 9: MethylSeekR.Rnw:182-184 (eval = FALSE)
###################################################
## system.file("extdata", "SNVs_hg18_chr22.tab",
## package="MethylSeekR")


###################################################
### code chunk number 10: MethylSeekR.Rnw:191-195
###################################################
snpFname <- system.file("extdata", "SNVs_hg18_chr22.tab",
package="MethylSeekR")
snps.gr <- readSNPTable(FileName=snpFname, seqLengths=sLengths)
head(snps.gr)


###################################################
### code chunk number 11: MethylSeekR.Rnw:205-206
###################################################
meth.gr <- removeSNPs(meth.gr, snps.gr)


###################################################
### code chunk number 12: MethylSeekR.Rnw:229-231
###################################################
plotAlphaDistributionOneChr(m=meth.gr, chr.sel="chr22", 
num.cores=1)


###################################################
### code chunk number 13: MethylSeekR.Rnw:249-251 (eval = FALSE)
###################################################
## library(parallel)
## detectCores()


###################################################
### code chunk number 14: MethylSeekR.Rnw:263-266
###################################################
PMDsegments.gr <- segmentPMDs(m=meth.gr, chr.sel="chr22", 
seqLengths=sLengths, num.cores=1)
head(PMDsegments.gr)


###################################################
### code chunk number 15: MethylSeekR.Rnw:285-288
###################################################
plotAlphaDistributionOneChr(m=subsetByOverlaps(meth.gr, 
PMDsegments.gr[values(PMDsegments.gr)$type=="notPMD"]), chr.sel="chr22", 
num.cores=1)


###################################################
### code chunk number 16: myfig1
###################################################
plotPMDSegmentation(m=meth.gr, segs=PMDsegments.gr)


###################################################
### code chunk number 17: MethylSeekR.Rnw:326-328 (eval = FALSE)
###################################################
## savePMDSegments(PMDs=PMDsegments.gr, 
## GRangesFilename="PMDs.gr.rds", TableFilename="PMDs.tab")


###################################################
### code chunk number 18: MethylSeekR.Rnw:369-375
###################################################
library(rtracklayer)
session <- browserSession()
genome(session) <- "hg18"
query <- ucscTableQuery(session, "cpgIslandExt")
CpGislands.gr <- track(query, asRangedData = FALSE)
genome(CpGislands.gr) <- NA


###################################################
### code chunk number 19: MethylSeekR.Rnw:385-387
###################################################
CpGislands.gr <- 
suppressWarnings(resize(CpGislands.gr, 5000, fix="center"))


###################################################
### code chunk number 20: MethylSeekR.Rnw:392-395
###################################################
stats <- calculateFDRs(m=meth.gr, CGIs=CpGislands.gr,
PMDs=PMDsegments.gr, num.cores=1)
stats


###################################################
### code chunk number 21: MethylSeekR.Rnw:418-423
###################################################
FDR.cutoff <- 5 
m.sel <- 0.5 
n.sel=as.integer(names(stats$FDRs[as.character(m.sel), ]
[stats$FDRs[as.character(m.sel), ]<FDR.cutoff])[1])
n.sel


###################################################
### code chunk number 22: MethylSeekR.Rnw:429-434
###################################################
UMRLMRsegments.gr <- segmentUMRsLMRs(m=meth.gr, meth.cutoff=m.sel,
nCpG.cutoff=n.sel, PMDs=PMDsegments.gr, 
num.cores=1, myGenomeSeq=Hsapiens, 
seqLengths=sLengths)
head(UMRLMRsegments.gr)


###################################################
### code chunk number 23: myfig3
###################################################
plotFinalSegmentation(m=meth.gr, segs=UMRLMRsegments.gr,
PMDs=PMDsegments.gr,meth.cutoff=m.sel)


###################################################
### code chunk number 24: MethylSeekR.Rnw:499-501 (eval = FALSE)
###################################################
## saveUMRLMRSegments(segs=UMRLMRsegments.gr, 
## GRangesFilename="UMRsLMRs.gr.rds", TableFilename="UMRsLMRs.tab")