i have to write a python program using pybedtools
I have three bedfiles which have results from a chip seq experiment and a refgenes text file, how do i determine the peaks which are bound in all three experiments. the peaks must overlap each other at least 20%..
1 answer
Here's a generic solution for finding three peaks that overlap from N input files, with more than 20% overlap between peak pairs.
$ bedops -u file1.bed file2.bed ... fileN.bed \
| bedmap --echo --echo-map --fraction-both 0.2 - \
| awk -F"|" '(split($2, a, ";") == 3)' \
> answer.bed
If peaks within one input file overlap by more than 20%, it is possible they will be in answer.bed if there is a third qualifying pair. But you can parse this file to filter those out, so that what is left are pairings between peaks from three different inputs. One easy way is to preprocess the inputs, appending an ID tag that ties an element to its input file, then use that ID after mapping to filter any within-ID pairs.
I'm honestly not too familiar with pybedtools. Maybe you can use subprocess to make system calls in Python: https://docs.python.org/2/library/subprocess.html
I now got the three peaks
Log in to answer this question.