FAQ - How to get average depth of coverage for targeted seq regions using BedTools

Question
If I have a list of targeted resequencing regions and a bam file is
there any easy way to get the average depth of coverage for each region?
I could calculate it from all of the values in the histogram output from
coverageBed -hist but that seems messy. In addition is there any way to
adjust the minimum number of hits at a location for it to be included in
the fraction of bases in B that had non-zero coverage into bases in B
that had more-than-X coverage?

Reply
There is no option for directly computing the mean coverage. However, one can combine the -d option (reports the depth at each base in each interval) with the groupBy tool found in the filo package (https://github.com/arq5x/filo

) to compute the mean.  For example, here's a 5 column BED of the UCSC simple repeat track:

$ head simRep.chr22.bed
22      16052163        16052199        trf     65
22      16059465        16059496        trf     53
22      16060479        16060524        trf     90
22      16060479        16060592        trf     120
22      16060517        16060609        trf     129
22      16071925        16071962        trf     74
22      16078616        16078686        trf     86
22      16079642        16079731        trf     61
22      16079644        16079696        trf     59
22      16079648        16079724        trf     52

I can compute the coverage at each base in each of these intervals as follows (column 6 is the position in the interval, column 7 is the depth):

$ coverageBed -abam aln.bam -b simRep.22.bed -d | head
22      17825583        17825839        trf     372     1       2
22      17825583        17825839        trf     372     2       2
22      17825583        17825839        trf     372     3       2
22      17825583        17825839        trf     372     4       2
22      17825583        17825839        trf     372     5       2
22      17825583        17825839        trf     372     6       2
22      17825583        17825839        trf     372     7       3
22      17825583        17825839        trf     372     8       3
22      17825583        17825839        trf     372     9       3
22      17825583        17825839        trf     372     10      3

Now, I can compute the mean by grouping on columns 1-5 of the output (the orig. BED entry) and computing the mean on column 7 for each grouped interval (the last column in the output below is the mean):

$ coverageBed -abam aln.bam -b simRep.22.bed -d | groupBy -g 1,2,3,4,5 -c 7 -o mean | head
22      17825583        17825839        trf     372     2.96875
22      16506905        16558419        trf     20099   7.6543852156695271205
22      16506905        16558419        trf     23903   7.6543852156695271205
22      16506941        16537880        trf     14401   7.47002165551569241586
22      16515001        16516375        trf     832     4.88427947598253275885
22      17039134        17049375        trf     5493    6.38755980861244054836
22      17039134        17049375        trf     4544    6.38755980861244054836
22      17039134        17049375        trf     5071    6.38755980861244054836
22      18087713        18088319        trf     759     5.42904290429042912791
22      19529662        19529964        trf     469     5.91721854304635730415

Note that coverageBed does not maintain the order of the intervals in the original file, so you might want to resort the output first:

$ coverageBed -abam aln.bam -b simRep.22.bed -d | sort -k1,1 -k2,2n | groupBy -g 1,2,3,4,5 -c 7 -o mean

>  In addition is there any way to
> adjust the minimum number of hits at a location for it to be included in
> the fraction of bases in B that had non-zero coverage into bases in B
> that had more-than-X coverage?

There is no way to do this except by taking the raw output of the -d or -hist and computing this yourself.

Best,
Aaron