How to split bedfile into non-overlapping regions and compute aggregation function on duplicate segments

Question

I have a bed-like file in which I have several overlapping regions with associated features. I want to split the overlapping regions into disjoint regions in order to process the file (e.g. with pandas), group by regions and compute an aggregate function.

The input data format is the following (overlapping are marked with *):

chr9    8855925    8986425    -1  *
chr9    8855965    8986325     0  *
chr9    8893857    9022105    -1

I want to split into disjoint intervals, then I want to groupby per-interval and compute an aggregation function (e.g. the mean of the last field):

chr9    8855925    8855965    -1
chr9    8855965    8893857  -0.5
chr9    8893857    8986325 -0.67
chr9    8986325    8986425    -1
chr9    8986425    9022105    -1

I was thinking of using bedtools and pandas to aggregate. However, I was unable to find the right tool to do this job in the bedtools/pybedtools documentation.

Answer 1

To make disjoint intervals, you could use BEDOPS bedops --partition, piping to bedmap --mean to get the mean signal over disjoint regions.

Starting with the input bedgraph file, convert it to five-column BED with GNU awk, putting the signal in the fifth column per UCSC convention:

$ awk -vOFS="\t" '{ print $1, $2, $3, ".", $4 }' /tmp/in.bedgraph | sort-bed - > /tmp/in.bed
$ cat in.bed
chr9    8855925 8986425 .   -1
chr9    8855965 8986325 .   0
chr9    8893857 9022105 .   -1

The disjoint set from these intervals would look like this:

$ bedops --partition /tmp/in.bed
chr9    8855925 8855965
chr9    8855965 8893857
chr9    8893857 8986325
chr9    8986325 8986425
chr9    8986425 9022105

We can pipe these intervals to bedmap --mean with the original five-column BED file as a "map" file, calculating the mean signal over the disjoint intervals:

$ bedops --partition /tmp/in.bed | bedmap --echo --mean --delim '\t' - /tmp/in.bed > /tmp/answer.bedgraph
$ cat /tmp/answer.bedgraph
chr9    8855925 8855965 -1.000000
chr9    8855965 8893857 -0.500000
chr9    8893857 8986325 -0.666667
chr9    8986325 8986425 -1.000000
chr9    8986425 9022105 -1.000000

The --mean option calculates the arithmetic mean. Other aggregation functions are available to process the signal or score (fifth) column in map BED files, like --stdev, --median, --mad, etc.

When calculating statistics, the default precision is six digits. The --prec <n> option can be used to specify desired precision in the score result.

See bedmap --help or the online documentation for a full description.

To avoid making and cleaning up intermediate files, bash process substitutions could be used to generate a one-liner:

$ bedops --partition <(sort-bed in.bedgraph | awk -vOFS="\t" '{ print $1, $2, $3, ".", $4 }') | bedmap --echo --mean --delim '\t' - <(sort-bed in.bedgraph | awk -vOFS="\t" '{ print $1, $2, $3, ".", $4 }') > answer.bedgraph

This can be useful for cleaner pipelines or with processing very large, whole-genome scale inputs on systems with limited disk space.