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
2 answers
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
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.
Log in to answer this question.