This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Count 0, 5, 20 covered sites in bam

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.

wgs bam
• 2,039 views • link
updated by Pierre Lindenbaum Pierre Lindenbaum • written by DVA DVA

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

• 1 views • link
written by GenoMax GenoMax

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

• 1 views • link
written by DVA DVA

2 answers 

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
• 1 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Thanks for the help! Appreciate it.

• 1 views • link
written by DVA DVA

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

• 1 views • link
written by Devon Ryan Devon Ryan

Thank you very much!

• 1 views • link
written by DVA DVA

Log in to answer this question.

More posts like this