DMRcaller.R
In DMRcaller: Differentially Methylated Regions caller

## ----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()