<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>

CC BY-SA 4.0
2018-11-25T08:37:35.207Z
Do you want coverage for each sample separately, or summed up across all samples?
winni2k
CC BY-SA 4.0
2018-11-25T08:45:04.187Z
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?
winni2k
CC BY-SA 4.0
2018-11-26T12:13:23.710Z
What's your end goal after making a text file with the values?
Devon Ryan
<sam><samtools><bed><bedtools><coverage>
bedtools single-nucleotide coverage in BED-specified regions for multiple BAMs
# 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

