I have an inbred mouse strain that I am sequencing and there should be little to NO heterozygosity. Yet with the default settings of UGT -heterozygosity (which is 0.001) many homs are being called as hets. When 230/250 reads are alternate and 20/250 are reference, it calls a het, even though it should be homozygous alternate.
What do you recommendations for this setting for inbred animals?
thanks, GATK is great!
Dear developers, I have looked into the forum for similar questions but I couldn't find any. I have several cases in which I get homozygous calls in positions with ~50% of reads calling the mutation (or less), please find here an example of a position validated by Sanger (as het) in which I have a high coverage (~400 reads in total) here is the results with UnifiedGenotyper:
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT L1 chr4 998101 . C T 7296.77 . AC=2;AF=1.00;AN=2;BaseQRankSum=16.344;DP=411;Dels=0.00;FS=0.000;HaplotypeScore=11.8258;MLEAC=2;MLEAF=1.00;MQ=59.94;MQ0=0;MQRankSum=-1.436;QD=17.75;ReadPosRankSum=-0.062 GT:AD:DP:GQ:PL 1/1:203,208:411:99:7325,581,0
can you help me in this case? I am really puzzled.
I have run it with v2.2, 2.4 and 2.5 and I always had the same genotype call (the excerpt here is from the 2.5, I have downloaded it just to check it was not caused by a bug already fixed). It's not a downsampling issue since I have high coverage samples (HaloPlex) and used higher dcov than the default.
I have two SNP calls as shown below:
chr1 86208198 . G C 903.25 PASS AC=2;AF=1.00;AN=2;BaseQRankSum=0.117;DP=33;Dels=0.00;FS=6.463;HRun=0;HaplotypeScore=8.5544;MQ=43.69;MQ0=0;MQRankSum=-0.506;QD=27.37;ReadPosRankSum=0.272;SB=-414.82 GT:AD:DP:GQ:PL 1/1:2,30:33:18.60:936,19,0
chr16 14895239 . C T 671.60 . AC=1;AF=0.50;AN=2;BaseQRankSum=-0.582;DP=33;Dels=0.00;FS=0.000;HRun=0;HaplotypeScore=1.5317;MQ=53.44;MQ0=0;MQRankSum=0.970;QD=20.35;ReadPosRankSum=-0.711;SB=-292.79 GT:AD:DP:GQ:PL 0/1:3,30:33:6.39:701,0,6
The first SNP is categorized as 1/1 and the second SNP as 0/1. For both the SNP's the variant allele ratio are 30/32=0.9375 and 30/33=0.909 which are approximately equal and above 0.9. On what criteria one SNP is determined as 0/1 and the other as 1/1?
As per my knowledge both the SNPs should be 1/1. Could anyone comment the reason for this discrepancy?
I have a VCF containing 7.4m SNPs over 70 individuals from an F2 intercross, called by the UnifiedGenotyper v2.3.6. I am trying to set appropriate thresholds for filtering these SNPs. The attached plots summarise the individual calls from this data set, with depth on the x-axis, genotype quality on the y-axis and frequency of particular DP+GQ combinations shown in greyscale. The first plot shows 0/1 (heterozygote) calls, the second shows 0/0 (homozyote) calls (the 1/1 plot looks similar to the 0/0 plot).
The homozygote plot shows a clear relationship between minimum depth and maximum GQ; it is impossible to get high GQs at low depth. However, this is not the case for heterozygotes. This makes intuitive sense to me - at low depth, one cannot be sure that a call really is homozygote; perhaps the other allele simply hasn't been sequenced. But we can have more confidence in a low depth heterozygote, as both alleles have been seen.
However, I am wondering what your recommendations for best practice are here; do you recommend using the same GQ thresholds for homozygote and heterozygote calls, or different thresholds? If the same thresholds, it seems like there will be a bias at low coverage; many (quite possibly real) homozygote calls will be excluded, which will make it appear that there is an excess of heterozygosity in low coverage individuals.
Also, there seems to be a periodicity in the homozygote (but not the heterozygote) GQ values; GQ values divisible by three have a different distribution to other GQ values. I assume this doesn't affect the results too much (after all, the scale is fairly arbitrary in the first place) but I'd be interested to know what causes this, if it is known.
Thanks for your help,