# bedtools single-nucleotide coverage in BED-specified regions for multiple BAMs

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.

• Do you want coverage for each sample separately, or summed up across all samples? – winni2k Nov 25 '18 at 8:37
• 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 Nov 25 '18 at 8:45
• What's your end goal after making a text file with the values? – Devon Ryan Nov 26 '18 at 12:13