Commit bca95c80 authored by Safia Saci's avatar Safia Saci
Browse files

Update dmr_analysis.R

parent fe785cc7
# Packages ----------
library("dmrseq")
library("bsseq")
library("rtracklayer")
......@@ -53,7 +54,7 @@ cov_keep <- cov[keep, ]
mvalue_keep <- mvalue[keep, ]
# Construct the bsseq object --------
my_bsseq <- BSseq(M = mvalue_keep,
Cov = cov_keep,
......@@ -67,6 +68,8 @@ rm(cov, mvalue, windows_genome, cov_keep, mvalue_keep, window)
pData(my_bsseq) <- pData
#bsseq_test <- chrSelectBSseq(my_bsseq, seqnames = "12")
# dmrseq -------
dmr_line <- dmrseq(my_bsseq,
testCovariate = "line",
adjustCovariate = "condition",
......@@ -74,7 +77,7 @@ dmr_line <- dmrseq(my_bsseq,
smooth = FALSE)
dmr_cond <- dmrseq(bsseq_test,
dmr_cond <- dmrseq(my_bsseq,
testCovariate = "condition",
adjustCovariate = "line",
minNumRegion = 3,
......@@ -86,23 +89,25 @@ saveRDS(dmr_line, file = "dmrseq_result_line_.RDS")
saveRDS(dmr_cond, file = "dmrseq_result_cond.RDS")
# how many regions were significant at the FDR (q-value) threshold of 0.05
# Exploring results --------------
## how many regions were significant at the FDR (q-value) threshold of 0.05
sum(dmr_line$qval < 0.05 & !is.na(dmr_line$qval))
#[1] 67
sum(dmr_cond$qval < 0.05 & !is.na(dmr_cond$qval))
#[1] 0
# select just the regions below FDR 0.05 and place in a new data.frame
## select just the regions below FDR 0.05 and place in a new data.frame
sigRegions <- dmr_line[dmr_line$qval < 0.05 & !is.na(dmr_line$qval),]
# proportion of regions with hyper-methylation
## proportion of regions with hyper-methylation
sum(sigRegions$stat > 0) / length(sigRegions)
#[1] 0.4179104
#raw mean methylation differences
##raw mean methylation differences
rawDiff <- meanDiff(bs = my_bsseq,
dmrs = sigRegions,
testCovariate = "line")
......@@ -111,10 +116,10 @@ rawDiff <- meanDiff(bs = my_bsseq,
# [7] 0.53302195 0.56387873 -0.50272785 0.55070516
# sort the significant dmrs by the beta value
## sort the significant dmrs by the beta value
sigRegions$beta <- sort(sigRegions$beta, decreasing= TRUE)
# plot the significant dmrs
## plot the significant dmrs
for (i in 1:67) {
png(paste0("sig_dmr_", i, ".png"))
plotDMRs(my_bsseq,
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment