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

Update dmr_analysis.R

parent 600f2040
......@@ -67,14 +67,14 @@ rm(cov, mvalue, windows_genome, cov_keep, mvalue_keep, window)
pData(my_bsseq) <- pData
#bsseq_test <- chrSelectBSseq(my_bsseq, seqnames = "12")
result_line <- dmrseq(my_bsseq,
dmr_line <- dmrseq(my_bsseq,
testCovariate = "line",
adjustCovariate = "condition",
minNumRegion = 3,
smooth = FALSE)
result_cond <- dmrseq(bsseq_test,
dmr_cond <- dmrseq(bsseq_test,
testCovariate = "condition",
adjustCovariate = "line",
minNumRegion = 3,
......@@ -82,21 +82,45 @@ result_cond <- dmrseq(bsseq_test,
saveRDS(result_line, file = "dmrseq_result_line_.RDS")
saveRDS(result_cond, file = "dmrseq_result_cond.RDS")
saveRDS(dmr_line, file = "dmrseq_result_line_.RDS")
saveRDS(dmr_cond, file = "dmrseq_result_cond.RDS")
png("dmrseq_line.png")
plotDMRs(my_bsseq,
extend = 8000,
regions=result_line[1,],
testCovariate="line")
dev.off()
# 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
sigRegions <- dmr_line[dmr_line$qval < 0.05 & !is.na(dmr_line$qval),]
# proportion of regions with hyper-methylation
sum(sigRegions$stat > 0) / length(sigRegions)
#[1] 0.4179104
#raw mean methylation differences
rawDiff <- meanDiff(bs = my_bsseq,
dmrs = sigRegions,
testCovariate = "line")
# > head(rawDiff)
# [1] -0.06697107 0.65499397 -0.72611870 -0.64181925 -0.49965266 -0.45594677
# [7] 0.53302195 0.56387873 -0.50272785 0.55070516
# sort the significant dmrs by the beta value
sigRegions$beta <- sort(sigRegions$beta, decreasing= TRUE)
# plot the significant dmrs
for (i in 1:67) {
png(paste0("sig_dmr_", i, ".png"))
plotDMRs(my_bsseq,
extend = 2500,
regions=sigRegions[i,],
testCovariate="line")
dev.off()
}
png("dmr_test_cond.png")
plotDMRs(bsseq_test_cond,
extend = 8000,
regions=result_cond[1,],
testCovariate="condition")
dev.off()
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