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
-3/5 options to extract single-nucleotide read-end coverage.