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.