Count 0, 5, 20 covered sites in bam

Question

[DVA]

Hello,

I wonder how to use a bam file to learn how many sites are not covered at all, how many are covered >=5 times and 20 times. We got such stats from our sequencing company, but want to see we can do it ourselves as well.

Thank you very much.

[GenoMax]

Are you looking to get the stats at individual base level or a interval window?

[DVA]

Thank you for the comment @genomax, I am looking at individual base level.

Answer 1

[Pierre Lindenbaum]

samtools depth -a in.bam | awk '{D=int($3);if(D<5) {D=0;} else if(D<20) {D=5;} else D=20; a[D]++;} END {for(x in a) printf("%s\t%d\n",x,a[x]);}'

5   14484
20  51
0   3771

[DVA]

Thanks for the help! Appreciate it.

Answer 2

[Devon Ryan]

You can get the relevant percentage of the genome with plotCoverage, though picard has some similar tools.

[DVA]

Thank you very much!