3
$\begingroup$

For example I have 3 bed files:

window.bed:

chr1    0   1000
chr1    1000    2000

lp1.bed:

chr1    20  50
chr1    30  40
chr1    1500    1900

lp2.bed:

chr1    68  786
chr1    899 987

How can I create this output:

chr1    0   1000    2   2
chr1    1000    2000    1   0

In this output columns 1, 2, 3 are the standard bed format cols, of the chromosome, start, and end of the features taken from window.bed above. Column 4 is the number of intersections between window.bed and lp1.bed per feature in window.bed. For example, the first element in column 4 is '2' because there are two intersections between window.bed and lp1.bed. Column 5 is the same but for the number of intersections between window.bed and lp2.bed for each feature in window.bed.

In the case that a feature is a subset of another feature in either lp1.bed, or lp2.bed, they should be counted as two intersections. For example, the second feature of lp1.bed is a subset of the first feature. This is reflected in the first element of column 4 of the output.

I suspect that it is doable with bedtools intersect, but I cannot work out the syntax.

$\endgroup$
0

2 Answers 2

2
$\begingroup$

The relevant bedtools option seems to be -c:

    -c  For each entry in A, report the number of overlaps with B.
        - Reports 0 for A entries that have no overlap with B.
        - Overlaps restricted by -f, -F, -r, and -s.

We can get the number of overlaps per region per file with:

$ bedtools intersect -c -a window.bed -b lp1.bed 
chr1    0   1000    2
chr1    1000    2000    1

$ bedtools intersect -c -a window.bed -b lp2.bed 
chr1    0   1000    2
chr1    1000    2000    0

So we're almost there, we just need to combine them, and that's a job for paste:


$ paste  <(bedtools intersect -c -a window.bed -b lp2.bed) \
         <(bedtools intersect -c -a window.bed -b lp1.bed)
chr1    0   1000    2   chr1    0   1000    2
chr1    1000    2000    0   chr1    1000    2000    1

But that gives us extra fields, so we can use cut to remove them:

$ paste  <(bedtools intersect -c -a window.bed -b lp1.bed) \
         <(bedtools intersect -c -a window.bed -b lp2.bed | cut -f 4)
chr1    0   1000    2   2
chr1    1000    2000    1   0

The paste command just joins lines from its two input files in the order it receives them, but the bedtools commands will always output the same order as the lines in window.bed so this should be fine.

The <(command) construct is called "process substitution" and is a way of passing the output of a command as a file to another command that expects files as input. Note that it isn't supported by all shells, but if you're using bash or zsh which are the most common these days, it should work.

$\endgroup$
2
$\begingroup$

Via bedmap --count, which gives the option to specify a custom value when regions do not map:

paste <(bedmap --echo --count --unmapped-val 0 --delim '\t' window.bed lp1.bed) <(bedmap --count --unmapped-val 0 window.bed lp2.bed) > answer.bed

Ref. https://bedops.readthedocs.io/en/latest/content/reference/statistics/bedmap.html#element-and-overlap-statistics

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.