This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How heterozygous GTs were defined in SNP calling vcf result

I use samtools-bcftools pipeline to call variations among a sample (30 individuals) with low coverage (~5X each). 

samtools mpileup -ugf genome.fa -b bam.list -t AD,INFO/AD,DP,SP  | bcftools call -vcO z -o out.vcf.gz

In the result, I found some strange heterozygous GTs:

(1) some heterozygous GTs has: DP=0, AD=0, SP=0

GT:PL:DP:SP:AD  0/1:0,0,0:0:0:0,0

(2) many heterozygous GTs have very low DPs:

GT:PL:DP:SP:AD 0/1:28,0,61:4:0:3,1

I don't think it is acceptable to deduce a heterozygous GT if the allele number of REF/ALT is only 3 and 1. The main problem is: How do I avoid such heterozygous GTs ? To add some parameters in samtools mpileup / bcftools call ? or use what software to filter them?

wgs snp samtools bcftools vcf
• 2,496 views • link
updated by Gabriel R. Gabriel R. • written by bill-zt • bill-zt

Well, I've used the traditional bcftools call -c model. When I switched to bcftools call -m, the GT:PL:DP:SP:AD 0/1:0,0,0:0:0:0,0 has been disappeared

• 0 views • link
written by bill-zt • bill-zt

1 answer

If you have a coverage of 4, there are two possibilities. Either you have a het site generating a (4 choose 3) or a homozygous site and the one base is an error (with probability that depends on the base quality). In that case, the het site is probably more likely if the base quality is now. The trick is the ratio of the likelihood. So 

log(homo ref)-log(het) = 28

based on the PL and the probability of error is 1/630 so it is quite high. You can use the PL values to make sure the desired state is many fold more likely than remaining models.

• 0 views • link
written by Gabriel R. Gabriel R.

Log in to answer this question.

More posts like this