Question: calculate average target coverage from bam file
0
gravatar for cmccabe
8 months ago by
cmccabe160
Chicago
cmccabe160 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 • 655 views
ADD COMMENTlink modified 8 months ago by Dan Gaston7.0k • written 8 months ago by cmccabe160
1
gravatar for Dan Gaston
8 months ago by
Dan Gaston7.0k
Canada
Dan Gaston7.0k 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 8 months ago by Dan Gaston7.0k
0
gravatar for Kevin Blighe
8 months ago by
Kevin Blighe19k
University College London Cancer Institute
Kevin Blighe19k 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 8 months ago • written 8 months ago by Kevin Blighe19k
1

Thank you both very much :).

ADD REPLYlink written 8 months ago by cmccabe160
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: 1922 users visited in the last hour