I am trying to get the mean depth coverage for bam files with bedtools.

So far I've been using the command:

```bedtools coverage -hist -abam bamfile.bam -b targets_exome.bed > bamfile.bam.hist.all.txt```

But I find myself with a result similar to:

    1	9998	10030	HWI-H217:64:C48Y8ACXX:4:1110:12530:87823/1	1	-	9998	10030	0,0,0	1	32,	0,	0	32	32	1.0000000
    
    1	9998	10030	HWI-H217:64:C48Y8ACXX:4:1110:12530:87823/2	1	-	9998	10030	0,0,0	1	32,	0,	0	32	32	1.0000000
    
    1	9999	10096	HWI-H217:64:C48Y8ACXX:4:1303:11532:70370/1	0	-	9999	10096	0,0,0	1	97,	0,	0	97	97	1.0000000
    
    1	9999	10096	HWI-H217:64:C48Y8ACXX:4:1303:11532:70370/2	0	-	9999	10096	0,0,0	1	97,	0,	0	97	97	1.0000000

where the different features are overlapping (see #1 and #3).

I want to compute average coverage, defined as the number of reads covering the captured coding sequence of a haploid reference. How to do it in this case?
To get the mean depth of coverage for each interval specified in your BED file, use:

    bedtools coverage -a BED -b [BAM] -mean > MeanCoverageBED.bedgraph

If you wanted to get the overall mean, then just pipe this into awk and get the average of the average, something like:

    bedtools coverage -a BED -b [BAM] -mean | awk '{total+=$4} END {print total/NR}'


There are other commands with bedtools that I use as standard in my pipelines:

# output depth of coverage for all regions in the BAM file

*NB - sequential positions at the same read depth are merged into a single region*

    bedtools genomecov -ibam [BAM] -bga -split > CoverageTotal.bedgraph

# output the per base read depth for each region in the BED file

    bedtools coverage -a [BED] -b [BAM] -d > PerBaseDepthBED.bedgraph

# get percent genome covered

    zero=$(bedtools genomecov -ibam [BAM] -g [ref.fasta] -bga | awk '$4==0 {bpCountZero+=($3-$2)} {print bpCountZero}' | tail -1)
    nonzero=$(bedtools genomecov -ibam [BAM] -g [ref.fasta] -bga | awk '$4>0 {bpCountNonZero+=($3-$2)} {print bpCountNonZero}' | tail -1)
    percent=$(bc <<< "scale=6; ($nonzero / ($zero + $nonzero))*100")


