## ----pvalue, echo=TRUE, message=FALSE, cache=FALSE, eval=FALSE----------- # pValue <- 2*pnorm(-abs(zScore)) ## ----pvalue_adjust, echo=TRUE, message=FALSE, cache=FALSE, eval=FALSE---- # pValue <- p.adjust(pValue, method="fdr") ## ----load_presaved, message=FALSE,cache=FALSE---------------------------- library(DMRcaller) #load presaved data data(methylationDataList) ## ----load, echo=TRUE, message=FALSE, cache=FALSE, 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) ## ----pool_data, echo=TRUE, message=FALSE, cache=FALSE, eval=FALSE-------- # # load the data # methylationDataAll <- poolMethylationDatasets(methylationDataList) # # # In the case of 2 elements, this is equivalent to # methylationDataAll <- poolTwoMethylationDatasets(methylationDataList[[1]], # methylationDataList[[2]]) ## ----_read_pool_data, echo=TRUE, message=FALSE, cache=FALSE, eval=FALSE---- # # load the data # methylationDataAll <- readBismarkPool(c(file_wt, file_met13)) ## ----profile, fig.width=12,fig.height=4,out.width='0.95\\textwidth', message=FALSE,cache=FALSE---- 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")) ## ----profile_fine, echo=TRUE, message=FALSE, cache=FALSE, fig.width=5,fig.height=5,out.width='0.95\\textwidth', 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)) ## ----coverage, fig.width=6,fig.height=5,out.width='0.6\\textwidth', message=FALSE,cache=FALSE---- # 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) ## ----coverage_compute, message=FALSE,cache=FALSE, eval=FALSE------------- # # compute the coverage in the two contexts # coverageCGWT <- computeMethylationDataCoverage(methylationDataList[["WT"]], # context="CG", # breaks = c(1,5,10,15)) ## ----compute_DMRs_CG_noise_filter, message=FALSE,cache=FALSE------------- 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) ## ----compute_DMRs_CG_neighbourhood, message=FALSE,cache=FALSE------------ # 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) ## ----compute_DMRs_CG_bins, message=FALSE,cache=FALSE--------------------- # 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) ## ----compute_DMRs_CG_GE, message=FALSE,cache=FALSE----------------------- # 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) ## ----merge_DMRs_CG_noise_filter, message=FALSE,cache=FALSE--------------- 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) ## ----analyse_reads_inside_regions_for_condition, message=FALSE,cache=FALSE---- #retrive the number of reads in CHH context in WT in CG DMRs DMRsNoiseFilterCGreadsCHH <- analyseReadsInsideRegionsForCondition( DMRsNoiseFilterCGMerged, methylationDataList[["WT"]], context = "CHH", label = "WT") print(DMRsNoiseFilterCGreadsCHH) ## ----DMR_distribution, fig.width=15,fig.height=2.5,out.width='0.95\\textwidth', message=FALSE,cache=FALSE---- # compute the distribution of DMRs hotspots <- computeOverlapProfile(DMRsNoiseFilterCGMerged, chr_local, windowSize=5000, binary=TRUE) # plot the distribution of DMRs plotOverlapProfile(GRangesList("Chr3"=hotspots)) ## ----local_profile, fig.width=15,fig.height=5,out.width='0.95\\textwidth', message=FALSE,cache=FALSE---- # 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") ## ----session_info,echo=TRUE---------------------------------------------- sessionInfo()