calculate average target coverage from bam file

Question

[bioguy24]

I am trying to calculate the average coverage for each target in a bed file. I tried the below samtools 1.2 command but it gives an error. Is there a better way? Thank you :)

targetbed (tab-delimited)
chr12   9264962 9265142 chr12:9264962-9265142   A2M
chr12   9265945 9266149 chr12:9265945-9266149   A2M
chr12   9268349 9268455 chr12:9268349-9268455   A2M
chr22   43088885    43089967    chr22:43088885-43089967 A4GALT
chr3    137843095   137843730   chr3:137843095-137843730    A4GNT
chr3    137849680   137850108   chr3:137849680-137850108    A4GNT
chr12   53701262    53701507    chr12:53701262-53701507 AAAS
chr12   53701618    53701723    chr12:53701618-53701723 AAAS

desired output (tab-delimited)
Targets                 Gene   Reads  --- header is only to show the field, it is not necessary ----
chr12:9264962-9265142   A2M     25
chr12:9265945-9266149   A2M     25
chr12:9268349-9268455   A2M     30
chr22:43088885-43089967 A4GALT     15
chr3:137843095-137843730    A4GNT     20
chr3:137849680-137850108    A4GNT     10
chr12:53701262-53701507 AAAS     22
chr12:53701618-53701723 AAA     35

samtools depth -a file.bam target.bed | awk '{c++;s+=$3}END{print $4,$5,s/c}' > targetcov.txt 
depth: invalid option -- 'a'
[E::sam_parse1] missing SAM header
[W::sam_read1] parse error at line 1
[bam_plp_destroy] memory leak: 2. Continue anyway.
awk: cmd. line:1: fatal: division by zero attempted

Answer 1

[Kevin Blighe]

BEDTools is your man (or woman) here:

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

Output is in bedgraph format

If you change -mean to -d, the per-base read depth in your target regions will be output

[bioguy24]

Thank you both very much :).

Answer 2

[DG]

While there are other tools to do the job, the issue is you appear to have 1) a mismatch between the installed version of samtools and where you are taking your appropriate command from and 2) your command line for specifying the input files is incorrect. The -a option is only found in the newer versions of samtools (>1.2). You are getting an error about that being an invalid option so you lilely have o.19 installed on the computer. Secondly you give regions with your BED file using the -b flag. Right now samtools is trying to read your target.bed file as a bam file.