Entering edit mode
7.0 years ago
reza
▴
300
"we calculated the mean and standard deviation read depth for all the variants"
The mentioned sentence is part of an paper and my question is: how can i do it in my own project? the following command estimate mean read depth for all variants for me?
samtools depth -a input.bam | awk '{c++;s+=$3}END{print s/c}'
Your command line is computing average read depth across all positions, not on variants only. You should have a variant VCF and do something like the following:
perl -ne 'print "$1\n" if /DP=(\d+)/' in.vcf | awk '$3<500{++x;y+=$3;z+=$3*$3}END{y/=x;z=sqrt(z/x-y*y);print y,z}'thanks for your answer. i checked your command and its output for my vcf file is: 0 0
It depends on where the depth information is kept in your vcf. Not all vcfs use the DP tag. For example, some use DP4, AD, ... You need to understand the basic of vcf first.
vcf file in my work created using samtools and has DP tag.
Sorry to open a very old chain. I am also trying to find a way to do this. However, I also got the same output (i.e. 0 0) as @reza? @reza did you ever find a method to do this? Thanks