### R code from vignette source 'MethylSeekR.Rnw' ################################################### ### code chunk number 1: install (eval = FALSE) ################################################### ## if (!requireNamespace("BiocManager", quietly=TRUE)) ## install.packages("BiocManager") ## BiocManager::install("BSgenome") ################################################### ### code chunk number 2: genomes (eval = FALSE) ################################################### ## library(BSgenome) ## available.genomes() ################################################### ### code chunk number 3: install (eval = FALSE) ################################################### ## if (!requireNamespace("BiocManager", quietly=TRUE)) ## install.packages("BiocManager") ## BiocManager::install("BSgenome.Hsapiens.UCSC.hg18") ################################################### ### code chunk number 4: MethylSeekR.Rnw:102-103 ################################################### library(MethylSeekR) ################################################### ### code chunk number 5: MethylSeekR.Rnw:112-113 ################################################### set.seed(123) ################################################### ### code chunk number 6: MethylSeekR.Rnw:128-130 (eval = FALSE) ################################################### ## system.file("extdata", "Lister2009_imr90_hg18_chr22.tab", ## package="MethylSeekR") ################################################### ### code chunk number 7: MethylSeekR.Rnw:140-143 ################################################### library("BSgenome.Hsapiens.UCSC.hg18") sLengths=seqlengths(Hsapiens) head(sLengths) ################################################### ### code chunk number 8: MethylSeekR.Rnw:159-163 ################################################### 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:189-191 (eval = FALSE) ################################################### ## system.file("extdata", "SNVs_hg18_chr22.tab", ## package="MethylSeekR") ################################################### ### code chunk number 10: MethylSeekR.Rnw:198-202 ################################################### 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:212-213 ################################################### meth.gr <- removeSNPs(meth.gr, snps.gr) ################################################### ### code chunk number 12: MethylSeekR.Rnw:236-238 ################################################### plotAlphaDistributionOneChr(m=meth.gr, chr.sel="chr22", num.cores=1) ################################################### ### code chunk number 13: MethylSeekR.Rnw:256-258 (eval = FALSE) ################################################### ## library(parallel) ## detectCores() ################################################### ### code chunk number 14: MethylSeekR.Rnw:270-273 ################################################### PMDsegments.gr <- segmentPMDs(m=meth.gr, chr.sel="chr22", seqLengths=sLengths, num.cores=1) head(PMDsegments.gr) ################################################### ### code chunk number 15: MethylSeekR.Rnw:292-295 ################################################### 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:333-335 (eval = FALSE) ################################################### ## savePMDSegments(PMDs=PMDsegments.gr, ## GRangesFilename="PMDs.gr.rds", TableFilename="PMDs.tab") ################################################### ### code chunk number 18: MethylSeekR.Rnw:376-382 ################################################### library(rtracklayer) session <- browserSession() genome(session) <- "hg18" query <- ucscTableQuery(session, table = "cpgIslandExt") CpGislands.gr <- track(query) genome(CpGislands.gr) <- NA ################################################### ### code chunk number 19: MethylSeekR.Rnw:392-394 ################################################### CpGislands.gr <- suppressWarnings(resize(CpGislands.gr, 5000, fix="center")) ################################################### ### code chunk number 20: MethylSeekR.Rnw:399-402 ################################################### stats <- calculateFDRs(m=meth.gr, CGIs=CpGislands.gr, PMDs=PMDsegments.gr, num.cores=1) stats ################################################### ### code chunk number 21: MethylSeekR.Rnw:425-430 ################################################### FDR.cutoff <- 5 m.sel <- 0.5 n.sel=as.integer(names(stats$FDRs[as.character(m.sel), ] [stats$FDRs[as.character(m.sel), ]