The EM Algorithm for Haplotype-Frequency Estimation

Why Estimate Haplotypes from Unphased Genotypes?

At two loci, a diploid individual carries two haplotypes, one on each chromosome. Genotyping can identify the two alleles at each locus without identifying which alleles lie together on the same chromosome. That missing chromosome pairing is phase ambiguity. For most genotype classes there is only one compatible pair of haplotypes, but a double heterozygote has two possibilities. The expectation-maximization (EM) algorithm uses the observed genotypes and current frequency estimates to allocate these ambiguous individuals probabilistically, then updates the frequencies.

The goal here is to estimate the population frequencies of the four possible haplotypes at two biallelic loci. The algorithm does not determine the phase of each person with certainty. Rather, it estimates the population frequencies by accounting for uncertainty about phase under a specified model.

Definition: Let the alleles at locus 1 be \(A\) and \(a\), and those at locus 2 be \(B\) and \(b\). The four haplotypes are \(AB\), \(Ab\), \(aB\), and \(ab\), with frequencies \(f_{AB}\), \(f_{Ab}\), \(f_{aB}\), and \(f_{ab}\), respectively. They are nonnegative and sum to 1. A double heterozygote has genotype \(Aa,Bb\) and can have either coupling phase \(AB/ab\) or repulsion phase \(Ab/aB\).

Model and Likelihood for the Ambiguous Genotype

The standard two-locus EM model assumes a randomly mating population in Hardy–Weinberg equilibrium with respect to haplotypes, independent sampled individuals, and correctly observed genotypes at both loci. Under this model, an individual's two haplotypes are independent draws from the four haplotype frequencies. The genotype likelihood is the probability of the observed genotype under those draws. The method also assumes the four haplotype categories are the relevant possibilities; it does not model genotype error, missing calls, or population mixture.

For a genotype that identifies its haplotype pair uniquely, the compatible pair contributes its probability under the model. For example, an \(AB/AB\) individual has probability \(f_{AB}^2\), and an \(AB/Ab\) individual has probability \(2f_{AB}f_{Ab}\), because either chromosome can carry either haplotype. For a double heterozygote, both phase arrangements are compatible. Its probability is the sum of the probabilities of those two arrangements.

P(Aa,Bb)=2f_{AB}f_{ab}+2f_{Ab}f_{aB}.

The factor 2 in each term accounts for the two chromosome orders within a diploid individual. Given the genotype and current frequencies, the conditional probability that a double heterozygote is in coupling phase is therefore

q = P(AB/ab\mid Aa,Bb) = \frac{f_{AB}f_{ab}} {f_{AB}f_{ab}+f_{Ab}f_{aB}}.

The repulsion-phase probability is \(1-q\). These are posterior probabilities conditional on the genotype and the current frequency estimates; they are not fixed properties of a person independent of the model.

Deriving the EM Update

If phase were known for every individual, the frequency estimate for a haplotype would be its chromosome count divided by the total number of chromosomes. With unphased data, the double heterozygotes contribute expected rather than known chromosome counts. Let \(U_{AB},U_{Ab},U_{aB},U_{ab}\) denote the observed counts of chromosomes from genotypes whose haplotype pair is unambiguous. Let \(m\) be the number of double heterozygotes. At current frequencies, each has probability \(q\) of coupling phase and \(1-q\) of repulsion phase.

A coupling individual contributes one \(AB\) and one \(ab\) chromosome. A repulsion individual contributes one \(Ab\) and one \(aB\) chromosome. Thus the expected counts contributed by the ambiguous individuals are \(mq\) to each coupling haplotype and \(m(1-q)\) to each repulsion haplotype. If \(n\) individuals were genotyped, the total number of chromosomes is \(2n\).

EM update: At iteration \(t\), calculate the coupling posterior \(q^{(t)}\) using the current frequencies. Then replace each frequency by its expected chromosome count divided by \(2n\):

\begin{aligned} q^{(t)} &= \frac{f_{AB}^{(t)}f_{ab}^{(t)}} {f_{AB}^{(t)}f_{ab}^{(t)}+f_{Ab}^{(t)}f_{aB}^{(t)}},\\ f_{AB}^{(t+1)}&=\frac{U_{AB}+m q^{(t)}}{2n}, & f_{ab}^{(t+1)}&=\frac{U_{ab}+m q^{(t)}}{2n},\\ f_{Ab}^{(t+1)}&=\frac{U_{Ab}+m(1-q^{(t)})}{2n}, & f_{aB}^{(t+1)}&=\frac{U_{aB}+m(1-q^{(t)})}{2n}. \end{aligned}

These equations can also be justified from the complete-data likelihood. If the phase assignments were observed, the log likelihood, apart from a constant, would be the sum of the haplotype chromosome counts multiplied by the log of their respective frequencies. In the E-step, unknown phase indicators are replaced by their conditional expected counts. In the M-step, maximizing that expected log likelihood subject to the frequencies summing to 1 gives each frequency as its expected count divided by the total count.

\ell_{\mathrm{complete}} = C_{AB}\log f_{AB} +C_{Ab}\log f_{Ab} +C_{aB}\log f_{aB} +C_{ab}\log f_{ab}.

To see the maximization, let \(E_h\) be the expected count for haplotype \(h\). Add a Lagrange multiplier term \(\lambda(1-\sum_h f_h)\) to \(\sum_h E_h\log f_h\). Differentiating with respect to \(f_h\) gives \(E_h/f_h-\lambda=0\), so \(f_h=E_h/\lambda\). Since the frequencies sum to 1, \(\lambda=\sum_h E_h=2n\), giving \(f_h=E_h/(2n)\). This completes the update justification.

Worked Iterations

Worked Example: One E-Step and M-Step

Suppose a hypothetical sample has \(n=15\) individuals. Across the 10 individuals whose phases are unambiguous, the haplotype chromosome counts are \(U_{AB}=10\), \(U_{ab}=4\), \(U_{Ab}=3\), and \(U_{aB}=3\). The other \(m=5\) individuals are double heterozygotes. Start with frequencies \(f_{AB}=0.40\), \(f_{ab}=0.10\), \(f_{Ab}=0.25\), and \(f_{aB}=0.25\). These sum to 1, as required.

E-step. The coupling product is \(0.40(0.10)=0.0400\), and the repulsion product is \(0.25(0.25)=0.0625\). Therefore,

q=\frac{0.0400}{0.0400+0.0625}=0.3902, \qquad 1-q=0.6098.

The five double heterozygotes contribute an expected \(5(0.3902)=1.951\) chromosomes to each of \(AB\) and \(ab\), and \(5(0.6098)=3.049\) to each of \(Ab\) and \(aB\).

M-step. There are \(2n=30\) chromosomes. The updated frequencies are

\begin{aligned} f_{AB}^{(1)}&=\frac{10+1.951}{30}=0.3984,& f_{ab}^{(1)}&=\frac{4+1.951}{30}=0.1984,\\ f_{Ab}^{(1)}&=\frac{3+3.049}{30}=0.2016,& f_{aB}^{(1)}&=\frac{3+3.049}{30}=0.2016. \end{aligned}

The updated values sum to 1 (up to rounding). They are expected-count estimates, not the result of assigning each double heterozygote to whichever phase is more common.

Worked Example: Posterior Allocation at a Later Iteration

Consider a hypothetical iteration whose current frequency estimates are \(f_{AB}=81/140\), \(f_{ab}=25/140\), and \(f_{Ab}=f_{aB}=17/140\). They are nonnegative and sum to 1. For one double heterozygote, calculate the phase probabilities. The coupling product is \((81/140)(25/140)=2025/19600\). The repulsion product is \((17/140)^2=289/19600\). Hence

q = \frac{2025/19600}{2025/19600+289/19600} = \frac{2025}{2314} \approx 0.8751.

The posterior probability of repulsion phase is \(1-q=289/2314\approx0.1249\). If a hypothetical sample contains eight double heterozygotes, the expected number in coupling phase is \(8(0.8751)\approx7.001\), and the expected number in repulsion phase is \(8(0.1249)\approx0.9991\). The corresponding expected chromosome contributions are about 7.001 to each of \(AB\) and \(ab\), and 0.9991 to each of \(Ab\) and \(aB\). These are fractional expected counts; they do not imply fractional chromosomes in an individual.

Worked Example: Updating Again and Checking the Counts

Continue the first example from its updated frequencies, approximately \(0.3984,0.1984,0.2016,0.2016\) in the order \(AB,ab,Ab,aB\). For a clear check, retain the unrounded first-step values: \(f_{AB}=10.39024/30\), \(f_{ab}=5.39024/30\), and \(f_{Ab}=f_{aB}=6.10976/30\). The next coupling probability is

q^{(1)} = \frac{(10.39024)(5.39024)} {(10.39024)(5.39024)+(6.10976)^2} \approx \frac{56.014}{93.344} \approx0.6000.

(The final digits depend on rounding the displayed intermediate counts; using the exact first-step values gives approximately \(0.6000\).) The next expected counts are \(10+5(0.6000)=13.000\) for \(AB\), \(4+5(0.6000)=7.000\) for \(ab\), and \(3+5(0.4000)=5.000\) for each of \(Ab\) and \(aB\). Dividing by 30 gives the next frequencies \(0.4333,0.2333,0.1667,0.1667\), which sum to 1. This illustrates how an updated estimate changes the posterior allocation at the following iteration.

Running the Algorithm in Practice

For a real analysis, the EM steps are repeated until the frequency estimates change by less than a chosen numerical tolerance. The starting frequencies must be valid: they should be nonnegative and sum to 1. The update itself preserves these properties when the observed chromosome counts are valid, because the expected contributions from the double heterozygotes total \(2m\), and the unambiguous counts total \(2(n-m)\).

The algorithm's result can depend on initialization in difficult or weakly informative data, so a practical analysis may use several starting values and compare converged solutions. This is a computational safeguard, not a substitute for checking the model assumptions. EM iterations provide a way to find a likelihood-based estimate; they do not by themselves establish that the sampled population satisfies Hardy–Weinberg equilibrium or that the genotypes are error-free.

Key takeaway: Each iteration first uses current haplotype frequencies to calculate the posterior phase probability for every double heterozygote, then updates frequencies from expected chromosome counts. Unambiguous genotypes contribute their observed counts; ambiguous genotypes contribute expected counts.

Common Mistakes and Why the Method Matters

Assigning every double heterozygote to one phase. Choosing the more probable phase discards uncertainty and can produce biased chromosome counts. Use the posterior probability as a fractional expected allocation.
Forgetting the two chromosomes per individual. The denominator for frequency updates is \(2n\), not \(n\). The expected haplotype counts must sum to \(2n\).
Using allele counts in place of haplotype counts. Allele frequencies at each locus do not uniquely determine the four joint haplotype frequencies. The EM update tracks the four haplotypes separately.
Misreading a posterior phase probability. A value such as \(q=0.8751\) is conditional on the observed double-heterozygote genotype, the current frequency estimates, and the model. It is not a statement that 87.51% of the sample is in coupling phase unless applied to the relevant count.
Ignoring the model. Population mixture, nonrandom mating, genotyping error, or departures from the assumed sampling model can make the estimates misleading. EM handles unobserved phase under its model; it does not repair a misspecified model.

The practical value of the method is that it extracts information about linkage between loci even when phase is unavailable for some individuals. The result supplies haplotype-frequency estimates needed to describe two-locus variation. The central discipline is to preserve the uncertainty in phase rather than silently turning an ambiguous genotype into a definite chromosome pair.

Check Your Understanding

Try these questions using the definitions and update equations above.

Why does the probability of a double heterozygote include two terms rather than just the coupling-phase term?
If \(f_{AB}f_{ab}=0.06\) and \(f_{Ab}f_{aB}=0.02\), what is the posterior probability of coupling phase?
A sample has 40 individuals. What is the denominator for the haplotype-frequency update, and why?
Explain how one double heterozygote contributes expected counts to the four haplotypes when its coupling probability is 0.7.
Why is assigning every ambiguous individual to the more probable phase not the EM update?