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 '18 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 '18 at 8:45
  • $\begingroup$ What's your end goal after making a text file with the values? $\endgroup$
    – Devon Ryan
    Nov 26 '18 at 12:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.