I have a huge amount of ~20x human WGS samples, aligned, and all SNVs that were called with GATK under standard germline parameters set.

What I need to do is to model SNVs Allele Frequency (AF) for different underlying Copy Numbers. I'd better provide a toy example. For particular genomic region X:

If X is presented by 2 copies for the particular samples, we expect AF to be super-close to 1 or to 0.5.

If X is presented by 4 copies, I expect any particular AF to be close to 0.25, 0.5, 0.75 or 1.

Of course, I can use Binomial Distribution for these purposes. However, as we know, the distribution is not exactly Binomial due to alignment/sequencing biases and the median AF for all heterozygous SNVs is more close to 0.48 but not to 0.5 as we would expect. Another thing: for high copy numbers we expect higher coverages. And GATK use several filters so I suppose that we will not see SNVs with AF like 0.125 (in case if the segment has ploidy 8) - despite the super high coverage there GATK may reject this "weird" AF.

I have read several papers that model SNVs AFs (and I agree that Beta Binomial Distribution may be quite accurate), however, I was not convinced enough that I should use the particular modelling. From your experience (in case if you do SNVs calling), which probabilistic distribution should I use? How should I estimate parameters for each of them (should I expect for CN4 AF=0.5 more frequent than AF=0.75 or vice versa, how to estimate this from data)?

UPD: For simplicity we can say that we have a lot of previously identified regions with ploidy different from CN2, and I can take these coordinates from here. So I can use more or less "supervised" learning for parameters' estimation.


I don't have enough experience to answer which probabilistic distribution should be used.

However, this questions also also asks how to estimate parameters of the distributions. If a binomial distribution is chosen, then Heng Li's paper titled "A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data" 1 is probably the definitive one. Section 2.3.1 of that paper describes an EM algorithm for estimating allele frequencies from multiple samples under the assumption of Hardy-Weinberg equilibrium for arbitrary but constant ploidy.

  • $\begingroup$ Ideally even most popular tools like GATK also uses EM algorithm for the record to estimate the likelihood while samtools uses both EM and Brent's method. The real problems with the setting up of ploidy apriori but am not aware as of now about tools that estimates ploidy and then uses that information for SNV calls and finds out the AF. I know tools like ABSOLUTE have the power to estimate the ploidy but then you have to use them as an input for your later SNP calls. But I reckon most tools usually works with the assumption with ploidy set apriori. This is what is coming to my mind as of now. $\endgroup$ – ivivek_ngs Jun 6 '17 at 10:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.