### R code from vignette source 'DMRcaller.Rnw' ################################################### ### code chunk number 1: pvalue (eval = FALSE) ################################################### ## pValue <- 2*pnorm(-abs(zScore)) ################################################### ### code chunk number 2: pvalue_adjust (eval = FALSE) ################################################### ## pValue <- p.adjust(pValue, method="fdr") ################################################### ### code chunk number 3: load_presaved ################################################### library(DMRcaller) #load presaved data data(methylationDataList) ################################################### ### code chunk number 4: load (eval = FALSE) ################################################### ## # specify the Bismark CX report files ## saveBismark(methylationDataList[["WT"]], ## "chr3test_a_thaliana_wt.CX_report") ## saveBismark(methylationDataList[["met1-3"]], ## "chr3test_a_thaliana_met13.CX_report") ## ## # load the data ## methylationDataWT <- readBismark("chr3test_a_thaliana_wt.CX_report") ## methylationDataMet13 <- readBismark("chr3test_a_thaliana_met13.CX_report") ## methylationDataList <- GRangesList("WT" = methylationDataWT, ## "met1-3" = methylationDataMet13) ################################################### ### code chunk number 5: pool_data (eval = FALSE) ################################################### ## # load the data ## methylationDataAll <- poolMethylationDatasets(methylationDataList) ## ## # In the case of 2 elements, this is equivalent to ## methylationDataAll <- poolTwoMethylationDatasets(methylationDataList[[1]], ## methylationDataList[[2]]) ################################################### ### code chunk number 6: _read_pool_data (eval = FALSE) ################################################### ## # load the data ## methylationDataAll <- readBismarkPool(c(file_wt, file_met13)) ################################################### ### code chunk number 7: profile ################################################### par(mar=c(4, 4, 3, 1)+0.1) plotMethylationProfileFromData(methylationDataList[["WT"]], methylationDataList[["met1-3"]], conditionsNames = c("WT","met1-3"), windowSize = 10000, autoscale = FALSE, context = c("CG")) ################################################### ### code chunk number 8: profile_fine (eval = FALSE) ################################################### ## regions <- GRanges(seqnames = Rle("Chr3"), ranges = IRanges(1,1E6)) ## ## # compute low resolution profile in 10 Kb windows ## profileCGWT <- computeMethylationProfile(methylationDataList[["WT"]], ## regions, ## windowSize = 10000, ## context = "CG") ## ## profileCGMet13 <- computeMethylationProfile(methylationDataList[["met1-3"]], ## regions, ## windowSize = 10000, ## context = "CG") ## ## profilesCG <- GRangesList("WT" = profileCGWT, "met1-3" = profileCGMet13) ## ## #plot the low resolution profile ## par(mar=c(4, 4, 3, 1)+0.1) ## par(mfrow=c(1,1)) ## plotMethylationProfile(profilesCG, ## autoscale = FALSE, ## labels = NULL, ## title="CG methylation on Chromosome 3", ## col=c("#D55E00","#E69F00"), ## pch = c(1,0), ## lty = c(4,1)) ################################################### ### code chunk number 9: coverage ################################################### # plot the coverage in the two contexts par(mar=c(4, 4, 3, 1)+0.1) plotMethylationDataCoverage(methylationDataList[["WT"]], methylationDataList[["met1-3"]], breaks = c(1,5,10,15), regions = NULL, conditionsNames=c("WT","met1-3"), context = c("CHH"), proportion = TRUE, labels=LETTERS, contextPerRow = FALSE) ################################################### ### code chunk number 10: coverage_compute (eval = FALSE) ################################################### ## # compute the coverage in the two contexts ## coverageCGWT <- computeMethylationDataCoverage(methylationDataList[["WT"]], ## context="CG", ## breaks = c(1,5,10,15)) ################################################### ### code chunk number 11: compute_DMRs_CG_noise_filter ################################################### chr_local <- GRanges(seqnames = Rle("Chr3"), ranges = IRanges(5E5,6E5)) # compute the DMRs in CG context with noise_filter method DMRsNoiseFilterCG <- computeDMRs(methylationDataList[["WT"]], methylationDataList[["met1-3"]], regions = chr_local, context = "CG", method = "noise_filter", windowSize = 100, kernelFunction = "triangular", test = "score", pValueThreshold = 0.01, minCytosinesCount = 4, minProportionDifference = 0.4, minGap = 0, minSize = 50, minReadsPerCytosine = 4, cores = 1) print(DMRsNoiseFilterCG) ################################################### ### code chunk number 12: compute_DMRs_CG_neighbourhood ################################################### # compute the DMRs in CG context with neighbourhood method DMRsNeighbourhoodCG <- computeDMRs(methylationDataList[["WT"]], methylationDataList[["met1-3"]], regions = chr_local, context = "CG", method = "neighbourhood", test = "score", pValueThreshold = 0.01, minCytosinesCount = 4, minProportionDifference = 0.4, minGap = 200, minSize = 1, minReadsPerCytosine = 4, cores = 1) print(DMRsNeighbourhoodCG) ################################################### ### code chunk number 13: compute_DMRs_CG_bins ################################################### # compute the DMRs in CG context with bins method DMRsBinsCG <- computeDMRs(methylationDataList[["WT"]], methylationDataList[["met1-3"]], regions = chr_local, context = "CG", method = "bins", binSize = 100, test = "score", pValueThreshold = 0.01, minCytosinesCount = 4, minProportionDifference = 0.4, minGap = 200, minSize = 50, minReadsPerCytosine = 4, cores = 1) print(DMRsBinsCG) ################################################### ### code chunk number 14: compute_DMRs_CG_GE ################################################### # load the gene annotation data data(GEs) #select the genes genes <- GEs[which(GEs$type == "gene")] # compute the DMRs in CG context over genes DMRsGenesCG <- filterDMRs(methylationDataList[["WT"]], methylationDataList[["met1-3"]], potentialDMRs = genes[overlapsAny(genes, chr_local)], context = "CG", test = "score", pValueThreshold = 0.01, minCytosinesCount = 4, minProportionDifference = 0.4, minReadsPerCytosine = 3, cores = 1) print(DMRsGenesCG) ################################################### ### code chunk number 15: merge_DMRs_CG_noise_filter ################################################### DMRsNoiseFilterCGMerged <- mergeDMRsIteratively(DMRsNoiseFilterCG, minGap = 200, respectSigns = TRUE, methylationDataList[["WT"]], methylationDataList[["met1-3"]], context = "CG", minProportionDifference = 0.4, minReadsPerCytosine = 4, pValueThreshold = 0.01, test="score") print(DMRsNoiseFilterCGMerged) ################################################### ### code chunk number 16: analyse_reads_inside_regions_for_condition ################################################### #retrive the number of reads in CHH context in WT in CG DMRs DMRsNoiseFilterCGreadsCHH <- analyseReadsInsideRegionsForCondition( DMRsNoiseFilterCGMerged, methylationDataList[["WT"]], context = "CHH", label = "WT") print(DMRsNoiseFilterCGreadsCHH) ################################################### ### code chunk number 17: DMR_distribution ################################################### # compute the distribution of DMRs hotspots <- computeOverlapProfile(DMRsNoiseFilterCGMerged, chr_local, windowSize=5000, binary=TRUE) # plot the distribution of DMRs plotOverlapProfile(GRangesList("Chr3"=hotspots)) ################################################### ### code chunk number 18: local_profile ################################################### # select a 20 Kb location on the Chr3 chr3Reg <- GRanges(seqnames = Rle("Chr3"), ranges = IRanges(510000,530000)) # create a list with all DMRs DMRsCGList <- list("noise filter" = DMRsNoiseFilterCGMerged, "neighbourhood" = DMRsNeighbourhoodCG, "bins" = DMRsBinsCG, "genes" = DMRsGenesCG) # plot the local profile par(cex=0.9) par(mar=c(4, 4, 3, 1)+0.1) plotLocalMethylationProfile(methylationDataList[["WT"]], methylationDataList[["met1-3"]], chr3Reg, DMRsCGList, conditionsNames = c("WT", "met1-3"), GEs, windowSize = 300, main="CG methylation") ################################################### ### code chunk number 19: session_info ################################################### sessionInfo()