SNP quality score in Samtools pileup

[wangzkai]

Hi,

I was examining the pileup by Samtools at a particular base of interest:

X 131016403 G G 103 0 60 53 T$T,,.t.....T.TT,,,..,.t,,...,t,t,tTT,.TT..T..Tt,,T,,t BFGGAGCEFEEE<B-GGGAGGFGGGGGG?GGFGGGGGBEGFGGGEGEFDGGEF

It looks like a clear heterozygous position with good coverage and decent base qualities, however it got a SNP quality score of 0 and a homozygous genotype call. Is there any possible explanation for this?

The data are from 75x2 PE reads and alignment was done using ELANDv2. Any help on this will be highly appreciated. Thanks!

[nilshomer]

Hi,

I was examining the pileup by Samtools at a particular base of interest:

X 131016403 G G 103 0 60 53 T$T,,.t.....T.TT,,,..,.t,,...,t,t,tTT,.TT..T..Tt,,T,,t BFGGAGCEFEEE<B-GGGAGGFGGGGGG?GGFGGGGGBEGFGGGEGEFDGGEF

It looks like a clear heterozygous position with good coverage and decent base qualities, however it got a SNP quality score of 0 and a homozygous genotype call. Is there any possible explanation for this?

The data are from 75x2 PE reads and alignment was done using ELANDv2. Any help on this will be highly appreciated. Thanks!

Maybe the mapping qualities for the variant reads are low?

[martian_bob]

Maybe the mapping qualities for the variant reads are low?

This is exactly what I found when I came across the same puzzling situation. Converting the data to BAM format and then visualizing it in IGV showed me that the apparently heterozygous SNP was getting all of its heterozygous bases from the low-quality ends of the reads - the SNP never once showed up in the beginning or middle of a read, only when it was within 7 nt of the end.

Odd? Yes. But if it were a true SNP, you'd expect to find it in half of the reads regardless of position.

[juan]

Yes that is the problem with SAMtools. The majority of variants are in the 2nd half of the read, hence you have lots of false positives.

[christophpale]

Does anyone have a code that can print out the positions within each of the reads where a given snp exist?

[pengchy]

Is there any more explanation?
I have found the following contrast examples:

scaffold2410 23912 G S 6 6 37 123 c,,,,,,,,,,,,,,,,,,,,ccc,,cc,,,,,,,,,,,,,ccc,,cc,cccc,,cc,,,,,,,cc,,c,cc,c,,c,,,,c,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, HHHHHHHHHHHHHHHHHFEHJCCHJEIHHHHHBHHHHHHFHHHHGHHHHHHHHHHCHHHHHHHJFFGJJJJHJHHHHHHHHHHHHHHHHHHHHH<HHHHHGHHHHGHHGHHHGHHGH<GJJIA

and :

scaffold12030   25942   A       R       37      44      25      44      gggggg,gg,g,,,,,,,ggggggggggggggggggggggggg,    HHHHHHGHHGHEHGHHHHHFHHHHHHHHHHGHBFHHHHHHHHHJ

Both of these two examples have similar reads quality and mapped on the reverse strand of reference, but with different "SNP quality", how these results produced?

Anyone who can give me any suggestions will be highly appreciated.

We have estimated the heterozygosis based on the results that filtered by VarFilter, obviously, we have under estimated the heterozygosis level.