This is a test version of Biostars. For the public version, visit https://www.biostars.org.
targets in a bed at or above a coverage

Is there a tool that will calculate percentage of targets in a bed file at or above 250x reads from a bam file? Currently I use bedtools to get all the reads that map to a target, followed by awk to calculate percentage above 250. However, since now our sequencing has increased there may be a more efficient way. Thank you :).

bed bam
• 2,926 views • link
updated by Ram Ram • written by bioguy24 bioguy24

2 answers

I'd recommend mosdepth. From the docs (especially relevant bullet in bold):

mosdepth can output:

  • per-base depth about 2x as fast samtools depth--about 25 minutes of CPU time for a 30X genome.
  • mean per-window depth given a window size--as would be used for CNV calling.
  • the mean per-region given a BED file of regions.
  • a distribution of proportion of bases covered at or above a given threshhold for each chromosome and genome-wide.
  • quantized output that merges adjacent bases as long as they fall in the same coverage bins e.g. (10-20)
  • threshold output to indicate how many bases in each region are covered at the given thresholds.
  • A summary of mean depths per chromosome and within specified regions per chromosome.

It should be fairly straightforward to go from this output to the percentage of regions above 250x.

• 0 views • link
written by bari.ballew bari.ballew

BEDtools will get you 90% there:

bedtools genomecov -ibam FILE.BAM -bg | awk '$4 > 249' > OUTPUT.BG
bedtools intersect -a TARGET.BED -b OUTPUT.BG -wao > OVERLAP

The first command reports the regions in your BAM with the desired coverage, the second reports the number of basepairs that overlap. Then it's a matter of dividing the overlap value by the size of the target.

• 0 views • link
written by harold.smith.tarheel harold.smith.tarheel

Log in to answer this question.

More posts like this