Nothing
## ----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()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.