Expected allele frequency distribution of SNVs in real NGS data

Question

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.

Answer 1

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.