This site is now read-only. You can find our new documentation site and support forum for posting questions here.
Be sure to read our welcome blog!
HET or HOM in RNA-seq calls with low allele frequency
I'm using the RNA-seq Best Practices pipeline to look at some of my data, and I'm seeing some calls that I'm having difficulty in understanding. An example site looks like this:
A G 5178.77 PASS AC=2;AF=1.00;AN=2;BaseQRankSum=-2.406;ClippingRankSum=-1.327;DP=220;ExcessHet=3.0103;FS=0.921;MLEAC=2;MLEAF=1.00;MQ=60.00;MQRankSum=1.513;QD=23.54;ReadPosRankSum=-0.752;SOR=0.590 GT:AD:DP:GQ:PL 1/1:18,202:220:82:5207,82,0
A G 4175.77 PASS AC=1;AF=0.500;AN=2;BaseQRankSum=-0.848;ClippingRankSum=0.482;DP=193;ExcessHet=3.0103;FS=1.648;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=-1.154;QD=21.64;ReadPosRankSum=3.068;SOR=0.950 GT:AD:DP:GQ:PL 0/1:30,163:193:99:4204,0,348
I have run the calling exactly as in your Best Practices (i.e. same hard filtering criteria). What I'm puzzled by is that an allele frequency of 18:202 and 30:163 yields HOM 1/1 and HET 0/1 in either case, respectively. It's a frequency of 8 % and 16 %, so it's relatively large difference (double), but my intuiting is that if a site that's only 16 % is OK to call HET, so would a 8 % site.
What's going on here? I assume that it's probably something with the mathematics underlying the pipeline, or my intuition is wrong in some manner of way. Is this one of those cases you mention in the Best Practices where low allele frequency can lead to wrong calls? If so, is there some way to change some part(s) of the pipeline to better accomodate such sites?
Thanks in advance!