Question: calculate average target coverage from bam file
0
gravatar for cmccabe
9 weeks ago by
cmccabe120
Chicago
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 • 211 views
ADD COMMENTlink modified 9 weeks ago by Dan Gaston6.9k • written 9 weeks ago by cmccabe120
1
gravatar for Dan Gaston
9 weeks ago by
Dan Gaston6.9k
Canada
Dan Gaston6.9k 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 weeks ago by Dan Gaston6.9k
0
gravatar for Kevin Blighe
9 weeks ago by
Kevin Blighe6.5k
Republic of Ireland (√Čire)
Kevin Blighe6.5k 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 weeks ago • written 9 weeks ago by Kevin Blighe6.5k
1

Thank you both very much :).

ADD REPLYlink written 8 weeks ago by cmccabe120
Please log in to add an answer.

Help
Access

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