samtools and awk depth calculation

[cmccabe]

Trying to calculate the average depth of a .bam file: Any ideas?

Essentialy I just want to get the ....X coverage.

 
samtools depth IonXpress_001.bam | awk '{sum+=$3} END { print "Average = ",sum/NR}' > output.txt
-bash: samtools: command not found
awk: cmd. line:1: fatal: division by zero attempted

Thank you .

EDIT: I also tried the below to calculate the genome of the .bam for 'NR', but it says samtools command not found.

 ./samtools view -H IonXpress_001.bam | grep -P '^@SQ' | cut -f 3 -d ':' | awk '{sum+=$1} END {print sum}'

= 3095693981

[cmccabe]

Figured it out, just forgot to add samtools to PATH.... oops, but the average is calculated to 0, should it not be $3 or is there somehing else? Thank you

 ./samtools depth IonXpress_001.bam | awk '{sum+=$1} END { print "Average = ",sum/3095693981}' > output.txt

[GenoMax]

Here is a thread with multiple options (including an awk solution): https://www.biostars.org/p/5165/

[cmccabe]

Thank you, I will try some of the other solutions in the post. Do you see why theawk is not working? I am trying to figure it out as our lab is doing more ngs and different data is often desired. Thank you .

[GenoMax]

Is that BAM directly from Ion software or did you make it yourself? Is the samtools depth part producing an output that looks reasonable?

[cmccabe]

The bam is from Ion software and I will post the output of samtools depth tomorrow. Thank you very much .

[cmccabe]

I found the perl script on that page incredibly useful.

./samtools mpileup in.bam | awk '{print $4}' | perl ~/coverage.pl

I am trying to modify the below to calculate the average coverage of a bed file, is that possible?

./samtools view -b in.bam <genomic region> | ./samtools mpileup - | awk '{print $4}' | perl ~/coverage.pl

Thank you .