Dear all, I came across the following situation: I have three WGBS methylation data from the epigenomics roadmap:

E016_WGBS_FractionalMethylation_allCHROM.bedgraph E013_WGBS_FractionalMethylation_allCHROM.bedgraph E004_WGBS_FractionalMethylation_allCHROM.bedgraph

The values of methylation states on specific sites are:

grep "chrY 5926917" E016_WGBS_FractionalMethylation_allCHROM.bedgraph
    chrY    59269171    59269175    0.98 grep "chrY 5926917"
E004_WGBS_FractionalMethylation_allCHROM.bedgraph
    chrY    59269171    59269173    0.93 chrY   59269173    59269175    1 grep "chrY
5926917"  E013_WGBS_FractionalMethylation_allCHROM.bedgraph
    chrY    59269171    59269173    0.91 chrY   59269173    59269175    0.9

Notice that the site E016 state has a uniformed methylation state so the nucleotides at chrY 59269171 59269173 and chrY 59269173 59269175 are merged into one site which is not true for the rest of the sites. So if I perform multiIntersect:

grep 5926917 E016_E004_E013_WGBS_FractionalMethylation_allCHROM2_multiIntersect.bedgraph
    chrY    59269171    59269175    3   1,2,3   1   1   1
    chrY    59269175    59269209    0   none    0   0   0

All sites are merged into one If I do subsequent intersects

grep 59269173  E016_E004_E013_WGBS_FractionalMethylation_allCHROM2.bedgraph
    chrY    59269171    59269173    0.98    chrY    59269171    59269175    0.93    chrY    59269171    59269173    0.91    2   2

chrY 59269171 59269173 0.98 chrY 59269171 59269175 1 chrY 59269173 59269175 0.9 2 2 chrY 59269173 59269175 0.98 chrY 59269171 59269175 0.93 chrY 59269171 59269173 0.91 2 2 chrY 59269173 59269175 0.98 chrY 59269171 59269175 1 chrY 59269173 59269175 0.9 2 2

The multiintersect does not include the appropriate methylation states (not its function) The subsequent IntersectBed keep the full file but somehow assign the values wrongly, what I would have expect would be:

chrY    59269171    59269173    0.98    chrY    59269171    59269175    0.93    chrY    59269171    59269173    0.91    2   2
chrY    59269173    59269175    0.98    chrY    59269171    59269175    1   chrY    59269173    59269175    0.9 2   2

Is there a way to get the correct intersect output? Thank you.

There is a patch in BEDOPS 2.4.29 that fixes this issue, so that we get items in correct order.

Given three BED5 files, A.bed, B.bed, and C.bed (your epigenome samples, say, preformatted with awk with unique, lexicographically-ordered IDs in the ID field column):

$ more A.bed
chr1    10468   10470   setA    0.97
chr1    10470   10472   setA    1
chr1    10483   10485   setA    0.95

$ more B.bed
chr1    10468   10470   setB    0.79
chr1    10470   10472   setB    0.93
chr1    10483   10485   setB    0.88

$ more C.bed
chr1    10468   10470   setC    0.91
chr1    10470   10472   setC    0.97
chr1    10483   10485   setC    0.87

We can make a partition and a union of these three sets:

$ bedops --partition A.bed B.bed C.bed > partition.bed
$ bedops --everything A.bed B.bed C.bed > union.bed

We can then map the scores from union.bed to the elements in partition.bed:

$ bedmap --echo --echo-map-score --delim '\t' --multidelim '\t' --prec 2 partition.bed union.bed
chr1    10468   10470   0.97    0.79    0.91
chr1    10470   10472   1.00    0.93    0.97    
chr1    10483   10485   0.95    0.88    0.87

However, I don't know when we will put out this revision; it may be a week or two, perhaps longer, as there some other issues I need to resolve with this release that will take time to fix.

You can, if you like, check out the current development version of BEDOPS via:

$ git clone https://github.com/bedops/bedops.git
$ cd bedops
$ git checkout v2p4p29
$ make
$ make install

Then:

$ ./bin/bedops --partition ...
$ ./bin/bedmap --echo --echo-map-score ...

etc.

I'm sorry I don't have a production version ready for you to use, but if this is critical work for you, then the above should hopefully get you going.

If you have "holes" in your sets, you need to patch them.

Say we have modified our initial three sets:

$ more A.bed
chr1    10468   10470   setA    0.97
chr1    10470   10472   setA    1
chr1    10483   10485   setA    0.95

$ more B.bed
chr1    10468   10470   setB    0.79
chr1    10483   10485   setB    0.88

$ more C.bed
chr1    10468   10470   setC    0.91
chr1    10470   10472   setC    0.97
chr1    10483   10485   setC    0.87

Note here that B.bed is missing one of the three elements.

As you note, if we bedmap the partition against the union, we'll get gaps.

Make the partition set as before:

$ bedops --partition A.bed B.bed C.bed > partition.bed

For each input set, bedmap that set against the partition, using the --indicator function to mark where there are and are not overlaps. Use this value with awk to post-process output:

$ bedmap --indicator --echo --echo-map --delim '\t' partition.bed A.bed | awk -v label='setA' -v OFS='\t' '{ if ($1 == 1) { print $5, $6, $7, $8, $9 } else { print $2, $3, $4, label, "NAN" } }' > A.patched.bed

$ bedmap --indicator --echo --echo-map --delim '\t' partition.bed B.bed | awk -v label='setB' -v OFS='\t' '{ if ($1 == 1) { print $5, $6, $7, $8, $9 } else { print $2, $3, $4, label, "NAN" } }' > B.patched.bed

$ bedmap --indicator --echo --echo-map --delim '\t' partition.bed C.bed | awk -v label='setC' -v OFS='\t' '{ if ($1 == 1) { print $5, $6, $7, $8, $9 } else { print $2, $3, $4, label, "NAN" } }' > C.patched.bed

(This is easy to script for N arbitrary sets. I'm writing everything out here explicitly so you see how this works.)

You could use 0.0 instead of NAN, depending on your experiment and what these values mean in that context. The NAN value is a reserved IEEE 754 floating point value, so it can be used here when it is meaningful to use that value as a marker.

Incidentally, this is why bedmap returns nothing when there are no overlaps, as 0.0 and NAN are values that have very specific and different meanings, and it is not really possible for bedmap to guess correctly which value is appropriate to use when the mapping result is a proper empty set.

In any case, here is what our patched sets look like:

$ more A.patched.bed
chr1    10468   10470   setA    0.97
chr1    10470   10472   setA    1
chr1    10483   10485   setA    0.95

$ more B.patched.bed
chr1    10468   10470   setB    0.79
chr1    10470   10472   setB    NAN
chr1    10483   10485   setB    0.88

$ more C.patched.bed
chr1    10468   10470   setC    0.91
chr1    10470   10472   setC    0.97
chr1    10483   10485   setC    0.87

From here out, it's the same as before. Take the union of the patched sets:

$ bedops --everything *.patched.bed > union.bed

Then map as before:

$ bedmap --echo --echo-map-score --delim '\t' --multidelim '\t' --prec 2 partition.bed union.bed
chr1    10468   10470   0.97    0.79    0.91
chr1    10470   10472   1.00    nan     0.97
chr1    10483   10485   0.95    0.88    0.87

(Note: This requires bedmap v2.4.29 or greater.)

You could use BEDOPS bedops --partition to partition the genomic space and bedmap to map the unioned set against the partitioned space.

Abstractly:

$ bedops --partition A.bed B.bed C.bed > partitionedABC.bed
$ bedops --everything A.bed B.bed C.bed > unionABC.bed
$ bedmap --echo --echo-map --delim '\t' partitionedABC.bed unionABC.bed > answer.bed

You could use awk to re-columnize answer.bed to whatever special format you're aiming at.

Or, as a one-liner, without intermediate files to slow things down:

$ bedmap --echo --echo-map --delim '\t' <(bedops --partition A.bed B.bed C.bed) <(bedops --everything A.bed B.bed C.bed) > answer.bed

I will wait patiently. Thank you for your efforts.