reference sequence coverage please

[sbsusers]

Hi all,

How to calculate how much percentage of reference sequence is covered/not covered (because of gaps) after alignment?

Thanks,

[Richard Finney]

1) figure out size of reference (should be about 3.1 GB for h. sapiens).

2) Figure out what's covered. Look at samtools depth command. Examine the output. Example is this :
chr2 89197150 13
chr2 89197151 13
chr2 89197152 13
chr2 89197153 13
chr2 89197154 7

note that there's coverage for many locations but sometimes it's zero so we can't just count the number of lines outputted so we have to do something like this ...

samtools depth yourfile.bam | awk '{if ($3>0) print $0}' | wc

3) Do a little algebra using answers from #1 and #2.