Body: <p>I have a BED6 (BED + name, score, strand info) file that defines some regions of interest.  I also have a set of BAMs corresponding to different samples.  I would like to obtain output similar to <code>bedtools genomecov -d -5</code>, which generates counts of read start positions, but I would like it restricted to regions of interest.</p>
  
  <p>I can do this for individual BAMs and regions by filtering input with <code>samtools view</code>, e.g.,</p>
  
  <pre><code>bedtools genomecov -d -5 -strand '+' -ibam &lt;(samtools view -b sample1.bam "chr1:12345-23456")
  </code></pre>
  
  <p>but this requires extracting position and strand info from the BED file, and then eventually I have to aggregate all the outputs back into a single file. This ends up being very slow (e.g., hours for ~50 regions of interest X 12 samples).</p>
  
  <p><strong>Is there a better way to accomplish this?</strong> </p>
  
  <p>Ideally, something along the lines of how <code>bedtools multicov</code> works would be great:</p>
  
  <pre><code>bedtools multicov -bams /path/to/my/bams/*.bam -bed roi.bed
  </code></pre>
  
  <p>but it doesn't have the <code>-d</code> and <code>-3/5</code> options to extract single-nucleotide read-end coverage.</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-11-25T08:37:35.207Z
Comments.Text: Do you want coverage for each sample separately, or summed up across all samples?
Comments.UserDisplayName: winni2k
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-11-25T08:45:04.187Z
Comments.Text: 50 ~1kb regions (as hinted at in the example) should not take hours to compute. Can you give more meta-information on the regions and BAMs that you are trying to process?
Comments.UserDisplayName: winni2k
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-11-26T12:13:23.710Z
Comments.Text: What's your end goal after making a text file with the values?
Comments.UserDisplayName: Devon Ryan
Tags: <sam><samtools><bed><bedtools><coverage>
Title: bedtools single-nucleotide coverage in BED-specified regions for multiple BAMs
ThreadText: # bedtools single-nucleotide coverage in BED-specified regions for multiple BAMs
  
  (Asked by: username_0 on 2018-11-24)
  
  I have a BED6 (BED + name, score, strand info) file that defines some regions of interest. I also have a set of BAMs corresponding to different samples. I would like to obtain output similar to `bedtools genomecov -d -5`, which generates counts of read start positions, but I would like it restricted to regions of interest.
  
  I can do this for individual BAMs and regions by filtering input with `samtools view`, e.g.,
  
  ```
  bedtools genomecov -d -5 -strand '+' -ibam <(samtools view -b sample1.bam "chr1:12345-23456")
  
  ```
  
  but this requires extracting position and strand info from the BED file, and then eventually I have to aggregate all the outputs back into a single file. This ends up being very slow (e.g., hours for ~50 regions of interest X 12 samples).
  
  **Is there a better way to accomplish this?**
  
  Ideally, something along the lines of how `bedtools multicov` works would be great:
  
  ```
  bedtools multicov -bams /path/to/my/bams/*.bam -bed roi.bed
  
  ```
  
  but it doesn't have the `-d` and `-3/5` options to extract single-nucleotide read-end coverage.
  
  **Tags:** sam, samtools, bed, bedtools, coverage
  
  > **Comments:**
  > **username_1**: Do you want coverage for each sample separately, or summed up across all samples?  
  > **username_1**: 50 ~1kb regions (as hinted at in the example) should not take hours to compute. Can you give more meta-information on the regions and BAMs that you are trying to process?  
  > **username_2**: What's your end goal after making a text file with the values?  
  
  ---
  
  ## 0 answers
  
  