Question: calculate average target coverage from bam file
gravatar for cmccabe
9 days ago by
cmccabe120 wrote:

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
ngs • 110 views
ADD COMMENTlink modified 9 days ago by Dan Gaston6.8k • written 9 days ago by cmccabe120
gravatar for Dan Gaston
9 days ago by
Dan Gaston6.8k
Dan Gaston6.8k wrote:

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.

ADD COMMENTlink written 9 days ago by Dan Gaston6.8k
gravatar for Kevin Blighe
9 days ago by
Kevin Blighe1.2k
Republic of Ireland
Kevin Blighe1.2k wrote:

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

ADD COMMENTlink modified 9 days ago • written 9 days ago by Kevin Blighe1.2k

Thank you both very much :).

ADD REPLYlink written 5 days ago by cmccabe120
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 567 users visited in the last hour