<p>Via <code>bedmap --count</code>, which gives the option to specify a custom value when regions do not map:</p>
<pre><code>paste &lt;(bedmap --echo --count --unmapped-val 0 --delim '\t' window.bed lp1.bed) &lt;(bedmap --count --unmapped-val 0 window.bed lp2.bed) &gt; answer.bed
</code></pre>
<p>Ref. <a href="https://bedops.readthedocs.io/en/latest/content/reference/statistics/bedmap.html#element-and-overlap-statistics" rel="nofollow noreferrer">https://bedops.readthedocs.io/en/latest/content/reference/statistics/bedmap.html#element-and-overlap-statistics</a></p>

CC BY-SA 4.0
2023-05-15T19:06:55.030Z
2023-05-15T19:06:55.030Z
Alex Reynolds
<p>The relevant <code>bedtools</code> option seems to be <code>-c</code>:</p>
<pre><code>    -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.
</code></pre>
<p>We can get the number of overlaps per region per file with:</p>
<pre><code>$ 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
</code></pre>
<p>So we're almost there, we just need to combine them, and that's a job for <code>paste</code>:</p>
<pre><code>
$ paste  &lt;(bedtools intersect -c -a window.bed -b lp2.bed) \
         &lt;(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
</code></pre>
<p>But that gives us extra fields, so we can use <code>cut</code> to remove them:</p>
<pre><code>$ paste  &lt;(bedtools intersect -c -a window.bed -b lp1.bed) \
         &lt;(bedtools intersect -c -a window.bed -b lp2.bed | cut -f 4)
chr1    0   1000    2   2
chr1    1000    2000    1   0
</code></pre>
<p>The <code>paste</code> 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 <code>window.bed</code> so this should be fine.</p>
<p>The <code>&lt;(command)</code> construct is called <a href="https://www.gnu.org/software/bash/manual/html_node/Process-Substitution.html" rel="nofollow noreferrer">&quot;process substitution&quot;</a> 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.</p>

CC BY-SA 4.0
2023-05-11T15:10:56.047Z
2023-05-11T15:10:56.047Z
terdon
<p>For example I have 3 bed files:</p>
<p>window.bed:</p>
<pre><code>chr1    0   1000
chr1    1000    2000
</code></pre>
<p>lp1.bed:</p>
<pre><code>chr1    20  50
chr1    30  40
chr1    1500    1900
</code></pre>
<p>lp2.bed:</p>
<pre><code>chr1    68  786
chr1    899 987
</code></pre>
<p>How can I create this output:</p>
<pre><code>chr1    0   1000    2   2
chr1    1000    2000    1   0
</code></pre>
<p>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.</p>
<p>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.</p>
<p>I suspect that it is doable with bedtools intersect, but I cannot work out the syntax.</p>

<bedtools><bed>
Create bed file where each column is the number of intersections with another bed file
# Create bed file where each column is the number of intersections with another bed file

(Asked by: username_0 on 2023-05-11)

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.

**Tags:** bedtools, bed

---

## 2 answers

### Answer by username_1 on 2023-05-15. Score: 2

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


---

### Answer by username_2 on 2023-05-11. Score: 2

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.


---

