## ----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)) ## ----correlation_plot, fig.width=6,fig.height=5,out.width='0.6\\textwidth', message=FALSE,cache=FALSE---- # compute the spatial correlation of methylation levels plotMethylationDataSpatialCorrelation(methylationDataList[["WT"]], distances = c(1,100,10000), regions = NULL, conditionsNames = c("WT"), context = c("CG"), labels = LETTERS, col = NULL, pch = c(1,0,16,2,15,17), lty = c(4,1,3,2,6,5), contextPerRow = FALSE, log = "x") ## ----correlation_compute, message=FALSE,cache=FALSE, eval=FALSE--------------- # # compute the coverage in the two contexts # correlation_CG_wt <- computeMethylationDataSpatialCorrelation(methylationDataList[["WT"]], # context="CG", # distances=c(1,10,100,1000,10000)) ## ----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") ## ----difference_methylation_proportions, fig.width=12,fig.height=5,out.width='0.95\\textwidth', message=FALSE,cache=FALSE---- # loading synthetic data data("syntheticDataReplicates") # create vector with colours for plotting cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7") # plotting the difference in proportions plot(start(methylationData), methylationData$readsM1/methylationData$readsN1, ylim=c(0,1), col=cbbPalette[2], xlab="Position in Chr3 (bp)", ylab="Methylation proportion") points(start(methylationData), methylationData$readsM2/methylationData$readsN2, col=cbbPalette[7], pch=4) points(start(methylationData), methylationData$readsM3/methylationData$readsN3, col=cbbPalette[3], pch=2) points(start(methylationData), methylationData$readsM4/methylationData$readsN4, col=cbbPalette[6], pch=3) legend(x = "topleft", legend=c("Treatment 1", "Treatment 2", "Control 1", "Control 2"), pch=c(1,4,2,3), col=cbbPalette[c(2,7,3,6)], bty="n", cex=1.0) ## ----computing biological replicates, fig.width=15,fig.height=5,out.width='0.95\\textwidth', message=FALSE,cache=FALSE---- # loading betareg library to allow using computeDMRsReplicates function library(betareg) # creating condition vector condition <- c("a", "a", "b", "b") # computing DMRs using the neighbourhood method DMRsReplicatesBins <- computeDMRsReplicates(methylationData = methylationData, condition = condition, regions = NULL, context = "CG", method = "bins", binSize = 100, test = "betareg", pseudocountM = 1, pseudocountN = 2, pValueThreshold = 0.01, minCytosinesCount = 4, minProportionDifference = 0.4, minGap = 0, minSize = 50, minReadsPerCytosine = 4, cores = 1) print(DMRsReplicatesBins) ## ----session_info,echo=TRUE--------------------------------------------------- sessionInfo()