Assuming you have the coordinates for DMRs in each of the samples and can get them into R it's fairly simple using GenomicRanges.

Here is a simple example of how you could try to get it done starting with e.g. BED-like coordinates for your peaks. Might have some bugs but if you're lucky it works straight away by copy-pasteing into R but can hopefully get you started..

In general you get data into R by providing the path to your file on interest and using one of the built-in functions like 'read.table'.

data<-read.table(file="c:/path/to/my/data/mydata.txt",sep="\t")
head(data) #for peak at top 6 rows

Load/install required package (in R)

if('GenomicRanges' %in% rownames(installed.packages() ) ) { library(GenomicRanges) } else { source("https://bioconductor.org/biocLite.R") ; biocLite("GenomicRanges") ; library(GenomicRanges) }

Create mock DMRs (substitute with own coordinates by reading in your file so it looks something like this)

DMRs1<-data.frame(chr=c("chr1","chr1","chr2","chr2"),start=c(100,500,1000,2000),end=c(100,500,1000,2000)+100,stringsAsFactors=F)
DMRs2<-data.frame(chr=c("chr1","chr1","chr2","chr2"),start=c(300,525,1000,1800),end=c(300,525,1000,1800)+100,stringsAsFactors=F)
DMRs3<-data.frame(chr=c("chr1","chr1","chr2","chr2"),start=c(350,500,1000,2200),end=c(350,500,1000,2200)+100,stringsAsFactors=F)
head(DMRs1) ##check object, your files should look something like this once read in..

Make GRanges objects

r1<-GRanges(seqnames=DMRs1[,"chr"],ranges=IRanges(start=DMRs1[,"start"],end=DMRs1[,"end"]))
r2<-GRanges(seqnames=DMRs2[,"chr"],ranges=IRanges(start=DMRs2[,"start"],end=DMRs2[,"end"]))
r3<-GRanges(seqnames=DMRs3[,"chr"],ranges=IRanges(start=DMRs3[,"start"],end=DMRs3[,"end"]))

Get DMR overlaps using coverage and convert back to GRanges

myRanges<-coverage(c(r1,r2,r3))
myRanges<-as(myRanges,"GRanges")

Only keep regions where at least 2 DMRs have some degree of overlap

myRanges<-myRanges[score(myRanges)>1] # >2 would require all three experiments to have overlapping DMR

Merge contiguous regions with variable number of overlapping peaks

#myRanges<-myRanges[width(myRanges)>50] #only keep overlaps larger than e.g. 50 bp, remove first '#' to use
myRanges<-reduce(myRanges) #merge contiguous regions

Make container for map-back

elementMetadata(myRanges)$consensusRegions2dmr1Row<-NA
elementMetadata(myRanges)$consensusRegions2dmr2Row<-NA
elementMetadata(myRanges)$consensusRegions2dmr3Row<-NA

Map back to original peaks

myMatches<-as.matrix(findOverlaps(myRanges,r1))
elementMetadata(myRanges)$consensusRegions2dmr1Row[myMatches[,1]]<-myMatches[,2]

myMatches<-as.matrix(findOverlaps(myRanges,r2))
elementMetadata(myRanges)$consensusRegions2dmr2Row[myMatches[,1]]<-myMatches[,2]

myMatches<-as.matrix(findOverlaps(myRanges,r3))
elementMetadata(myRanges)$consensusRegions2dmr3Row[myMatches[,1]]<-myMatches[,2]

Output results matrix

out<-data.frame(chr=as.character(seqnames(myRanges)),start=start(myRanges),end=end(myRanges),as(elementMetadata(myRanges),"matrix"),stringsAsFactors=F)

write.table(out,file="peaksOverlapsOut.txt",sep="\t",quote=F,row.names=F) #open in e.g. excel and check