Dear all,
I got this answer on the google groups from Dr. Tim Triche and it worked for me. I add it here, in case any one of you had the same question.
...
yes, of course. And I say this having used both the DSS-backed N=1 DMR calling approach and also DMRcate's approach on the same pediatric samples. It's been an education -- the arrays are better for some things (e.g. CNV) and the WGBS for others (although maybe not the ones that would most immediately come to mind). Anyways...
Let's use the built-in example.
data(dmrcatedata)
myMs <- logit2(myBetas)
myMs.noSNPs <- rmSNPandCH(myMs, dist=2, mafcut=0.05)
patient <- factor(sub("-.*", "", colnames(myMs)))
type <- factor(sub(".*-", "", colnames(myMs)))
design <- model.matrix(~patient + type)
myannotation <- cpg.annotate("array", myMs.noSNPs, what="M", arraytype = "450K",
analysis.type="differential", design=design, coef=39)
dmrcoutput <- dmrcate(myannotation, lambda=1000)
Now we have an object of class "dmrcate.output" (try class(dmrcoutput) if you don't believe me) with an element named "result".
Let's use that to turn it into something that we can crank out BED files from.
dmrcoutput$gr <- as(dmrcoutput$results$coord, "GRanges")
dmrcoutput$gr$score <- dmrcoutput$results$meanbetafc
dmrcoutput$gr$p <- dmrcoutput$results$Stouffer
# in this example, it's useless, but just FYI:
pcutoff <- 0.05
dmrcoutput$gr <- subset(dmrcoutput$gr, p < pcutoff)
# now split by direction
hyperHypo <- split(dmrcoutput$gr, ifelse(dmrcoutput$gr$score > 0, "hyper", "hypo"))
sapply(hyperHypo, length)
# hyper hypo
# 336 407
# now export each
library(rtracklayer)
for (direction in names(hyperHypo)) {
export(hyperHypo[[direction]], paste0("/tmp/", direction, "DMRs.hg19.bed"))
}
And that's that. Similar logic for the output of DSS' callDMR() output.