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

