Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-05-15T19:06:55.030Z
Answers.LastActivityDate: 2023-05-15T19:06:55.030Z
Answers.UserDisplayName: Alex Reynolds
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-05-11T15:10:56.047Z
Answers.LastActivityDate: 2023-05-11T15:10:56.047Z
Answers.UserDisplayName: terdon
Body: <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>
  
Tags: <bedtools><bed>
Title: Create bed file where each column is the number of intersections with another bed file
ThreadText: # 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"](https://www.gnu.org/software/bash/manual/html_node/Process-Substitution.html) 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.
  
  
  ---
  
  