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.

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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.