Allele-Frequency Estimation With Uncertain Genotype Calls

Why Keep Genotype Uncertainty?

A genotype call is not always certain. Sequencing reads or other assay data may be compatible with more than one genotype, especially when coverage is limited or the evidence is noisy. Reducing that evidence to a single best-guess genotype discards information: a call with probability 0.51 is treated just like one with probability 0.99. For estimating an allele frequency, we can instead retain the probabilities assigned to each possible genotype and use their expected allele counts.

This tutorial considers one biallelic locus in a sample of diploid individuals. The alternate allele is counted as 1 and the reference allele as 0. The goal is to estimate the population frequency $p$ of the alternate allele from uncertain genotype evidence. The estimation method depends on what the reported probabilities mean, so it is important to distinguish genotype likelihoods from posterior genotype probabilities.

Definition: For individual $i$, let $G_i\in\{0,1,2\}$ be the number of alternate alleles in the diploid genotype: 0 for reference homozygous, 1 for heterozygous, and 2 for alternate homozygous. A genotype likelihood $L_i(g)$ is the probability of the observed data $D_i$ given genotype $g$: $L_i(g)=P(D_i\mid G_i=g)$. A posterior genotype probability is $P(G_i=g\mid D_i)$; it combines the likelihood with a prior probability for the genotype.

Estimating Frequency From Posterior Genotype Probabilities

If posterior probabilities are available, the expected alternate-allele count for individual $i$ is the probability-weighted average of the three possible counts. Summing those expected counts gives the estimated total number of alternate alleles; dividing by the $2n$ allele copies in $n$ diploid individuals gives an allele-frequency estimate.

E(G_i\mid D_i)=\sum_{g=0}^{2} g\,P(G_i=g\mid D_i), \qquad \hat p=\frac{1}{2n}\sum_{i=1}^{n}E(G_i\mid D_i).

This is the posterior-expected allele count estimator for the supplied genotype probabilities. Its value depends on how those probabilities were obtained. For example, they may incorporate a frequency prior from a reference panel, or a model-based prior for the population being analyzed. The formula averages the uncertainty represented by the probabilities; it does not remove the assumptions or prior information used to produce them.

Worked Example: Averaging Posterior Genotype Probabilities

Suppose three hypothetical diploid individuals have posterior probabilities for $G=0,1,2$, respectively: $(0.70,0.25,0.05)$, $(0.10,0.60,0.30)$, and $(0.20,0.50,0.30)$. Each row sums to 1.

For individual 1, the expected alternate-allele count is $0(0.70)+1(0.25)+2(0.05)=0.35$. For individual 2, it is $0(0.10)+1(0.60)+2(0.30)=1.20$. For individual 3, it is $0(0.20)+1(0.50)+2(0.30)=1.10$.

The total expected count is $0.35+1.20+1.10=2.65$ alternate alleles among $2(3)=6$ allele copies. Thus

\hat p=\frac{2.65}{6}=0.4417.

For comparison, the most probable genotype for each person is 0, 1, and 1: the third individual’s largest probability is 0.50 for $G=1$, not 0.30 for $G=2$. Those hard calls contain $0+1+1=2$ alternate alleles, giving a hard-call estimate $2/6=0.3333$. The difference illustrates why the uncertain probabilities should not be silently replaced with most-probable calls.

From Genotype Likelihoods to an EM Update

Genotype likelihoods alone are not posterior genotype probabilities: they describe how compatible the observed data are with each genotype, without weighting those genotypes by their population frequencies. To estimate $p$ from likelihoods, we need a model for genotype frequencies. Here we use the standard Hardy–Weinberg model for a biallelic locus: allele copies are independently sampled with alternate-allele probability $p$, so genotype probabilities are $(1-p)^2$, $2p(1-p)$, and $p^2$ for $G=0,1,2$. We also assume individuals are independent draws from the same population, and that the genotype likelihoods correctly represent the data conditional on each genotype.

Under these assumptions, the probability of the observed data for individual $i$, conditional on $p$, is the sum over possible genotypes of likelihood times genotype probability. Assuming the individuals’ data are independent conditional on $p$, the sample likelihood is the product of these sums.

P(D_i\mid p)=\sum_{g=0}^{2}L_i(g)\binom{2}{g}p^g(1-p)^{2-g}, \qquad \mathcal{L}(p)=\prod_{i=1}^{n}P(D_i\mid p).

The $G_i$ are unobserved, which is why this is a latent-variable problem. The expectation-maximization (EM) algorithm alternates two steps. In the expectation step, calculate each genotype’s conditional probability using the current frequency estimate. In the maximization step, update $p$ to the expected total alternate-allele count divided by $2n$. Repeating these steps produces a sequence of frequency estimates.

EM update: Given a current estimate $p^{(t)}$, calculate $$ \tau_{ig}^{(t)} =\frac{L_i(g)\binom{2}{g}(p^{(t)})^g(1-p^{(t)})^{2-g}} {\sum_{h=0}^{2}L_i(h)\binom{2}{h}(p^{(t)})^h(1-p^{(t)})^{2-h}}. $$ Then update $$ p^{(t+1)}=\frac{1}{2n}\sum_{i=1}^{n}\sum_{g=0}^{2}g\,\tau_{ig}^{(t)}. $$ Here $\tau_{ig}^{(t)}$ is the posterior probability of genotype $g$ under the current frequency estimate.

Why the EM Update Has This Form

The complete data would include each observed data record $D_i$ and its unobserved genotype $G_i$. For fixed genotypes, the part of the complete-data log-likelihood that depends on $p$ is $\sum_i[G_i\log(p)+(2-G_i)\log(1-p)]$. The genotype likelihood terms and genotype combinatorial factors do not depend on $p$, so they do not affect the maximization.

In the expectation step, replace each unobserved count $G_i$ with its conditional expectation under the current estimate. Let $A=\sum_i\sum_g g\,\tau_{ig}^{(t)}$ be the resulting expected alternate-allele count. There are $2n$ allele copies in total, so the expected reference-allele count is $2n-A$. The part of the expected complete-data log-likelihood that depends on the candidate frequency $p$ is therefore

Q(p\mid p^{(t)})=A\log(p)+(2n-A)\log(1-p)+\text{constant}.

For an interior maximum, differentiate with respect to $p$, set the derivative to zero, and solve:

\frac{\partial Q}{\partial p} =\frac{A}{p}-\frac{2n-A}{1-p}=0 \quad\Longrightarrow\quad A(1-p)=p(2n-A) \quad\Longrightarrow\quad p=\frac{A}{2n}.

This is exactly the EM update. If $A=0$ or $A=2n$, the maximizer is at the corresponding boundary, $p=0$ or $p=1$. In ordinary applications, the updates are repeated until the change in $p$ is sufficiently small. EM does not guarantee that every possible likelihood has a unique maximum; the result can depend on initialization when the likelihood has multiple stationary points.

Worked Example: One EM Update From Genotype Likelihoods

Consider three hypothetical individuals with genotype likelihoods for $G=0,1,2$ equal to $(1,0.2,0.05)$, $(0.1,1,0.3)$, and $(0.05,0.4,1)$. Start the EM algorithm at $p^{(0)}=0.40$. These likelihood values need not sum to 1; multiplying all three likelihoods for one individual by the same positive constant would leave that individual’s posterior probabilities unchanged.

At $p^{(0)}=0.40$, the genotype probabilities under the Hardy–Weinberg model are $(0.36,0.48,0.16)$. For individual 1, the unnormalized weights are $(1(0.36),0.2(0.48),0.05(0.16))=(0.36,0.096,0.008)$, which sum to $0.464$. The normalized probabilities are approximately $(0.7759,0.2069,0.01724)$, so the expected count is $0(0.7759)+1(0.2069)+2(0.01724)=0.2414$.

For individual 2, the weights are $(0.1(0.36),1(0.48),0.3(0.16))=(0.036,0.48,0.048)$, with sum $0.564$. The normalized probabilities are approximately $(0.06383,0.8511,0.08511)$, giving expected count $0.8511+2(0.08511)=1.0213$.

For individual 3, the weights are $(0.05(0.36),0.4(0.48),1(0.16))=(0.018,0.192,0.16)$, with sum $0.370$. The normalized probabilities are approximately $(0.04865,0.5189,0.4324)$, giving expected count $0.5189+2(0.4324)=1.3838$.

The expected total is approximately $0.2414+1.0213+1.3838=2.6464$ alternate alleles. The first update is therefore

p^{(1)}=\frac{2.6464}{2(3)}=0.4411.

This is one iteration, not necessarily the final estimate. The next expectation step would use genotype probabilities calculated from $p^{(1)}=0.4411$, then another maximization step would update the frequency.

Why Posterior Probabilities Depend on the Frequency Model

The EM expectation step turns likelihoods into posterior genotype probabilities by combining them with the current population-frequency model. The same likelihoods can therefore produce different posterior probabilities at different values of $p$. This is a feature of the method, not a contradiction: the data evidence is unchanged, but its weighting by the genotype prior changes.

Worked Example: The Same Likelihoods Under Two Frequencies

For one hypothetical individual, suppose the genotype likelihoods are $(1,0.5,0.1)$. First use $p=0.20$, which gives genotype probabilities $(0.64,0.32,0.04)$. The unnormalized posterior weights are $(0.64,0.16,0.004)$, with total $0.804$. The expected alternate-allele count is

\frac{0(0.64)+1(0.16)+2(0.004)}{0.804} =\frac{0.168}{0.804} =0.2090.

Now use $p=0.50$, giving genotype probabilities $(0.25,0.50,0.25)$. The weights are $(0.25,0.25,0.025)$, with total $0.525$, so the expected count becomes

\frac{0(0.25)+1(0.25)+2(0.025)}{0.525} =\frac{0.30}{0.525} =0.5714.

The likelihoods are identical in both calculations; only the frequency-based genotype prior differs. This is why a supplied posterior probability should not be mistaken for a genotype likelihood.

Assumptions, Interpretation, and a Common Pitfall

The likelihood-based EM estimate is tied to its model. It assumes a shared allele frequency for the analyzed individuals, independent diploid genotypes, and Hardy–Weinberg genotype probabilities at the locus. It also relies on likelihoods that are calibrated conditional on genotype. If individuals are related, come from distinct populations with different frequencies, or do not follow the assumed genotype model, the update may not estimate the intended population frequency. In particular, Hardy–Weinberg assumptions should not be treated as a guarantee that the observed sample must have those genotype proportions.

A useful interpretation of the update is that each person contributes a fractional expected number of alternate alleles, rather than a forced count of 0, 1, or 2. The frequency is the average expected alternate-allele count per allele copy. It is not a statement that any individual has a fractional genotype.

Common mistake: “Choose the most probable genotype for each person, count the alternate alleles, and call the result the estimate that uses genotype probabilities.” This is incorrect because it discards the probabilities of the other genotypes. For example, probabilities $(0.20,0.50,0.30)$ imply an expected count of $0(0.20)+1(0.50)+2(0.30)=1.10$, even though the most probable genotype is $G=1$. The correction is to use the probability-weighted expected count when probabilities are available, and to use the frequency-weighted posterior probabilities in each EM iteration when starting from genotype likelihoods.

Finally, if genotype probabilities were calculated using a prior frequency from outside the current sample, averaging their expected counts gives an estimate conditional on those supplied probabilities. It is not automatically the maximum-likelihood estimate under the sample’s own frequency model. The EM procedure instead recalculates the genotype posterior probabilities using the current sample-frequency estimate at every iteration.

Key takeaway: With uncertain genotypes, estimate frequency by summing expected alternate-allele counts and dividing by $2n$. For genotype likelihoods, derive those expectations from likelihoods weighted by Hardy–Weinberg genotype probabilities at the current $p$, then iterate the resulting EM update.

Check Your Understanding

Answer each question using the genotype-count convention and model stated above.

For posterior probabilities $(0.10,0.70,0.20)$, what is the expected alternate-allele count for that individual?
Why are genotype likelihoods not themselves posterior genotype probabilities?
In the EM update, what quantities are calculated in the expectation step, and what is maximized in the maximization step?
For $n=5$ diploid individuals, why is the expected alternate-allele count divided by 10 to estimate $p$?
Name two assumptions of the likelihood-based EM procedure and explain how violating either could affect interpretation.