Marginal Likelihoods for Uncertain Genotypes

Why Marginalize Over Genotype?

Sequencing reads do not always identify a genotype unambiguously. At a site with two alleles, a diploid genotype can contain zero, one, or two copies of the alternate allele. A small or noisy set of reads may be compatible with more than one of these states. Replacing that uncertainty with a single genotype call discards information about how well each state fits the observed reads.

A likelihood for the observed reads should account for every possible genotype that could have produced them. It does this by combining the read likelihood conditional on each genotype with the probability assigned to that genotype, then summing across genotypes. The result is called a marginal likelihood because the unobserved genotype has been summed out.

Definition: Let $G\in\{0,1,2\}$ be the number of alternate alleles in a diploid genotype, and let $R$ denote the observed read data. Write $\pi_g=P(G=g)$ for the genotype probability and $L_g(R;\theta)=P(R\mid G=g;\theta)$ for the read likelihood under genotype $g$, where $\theta$ denotes parameters of the read model. The marginal likelihood of the observed reads is $$ L_{\mathrm{marg}}(\theta;R)=P(R;\theta)=\sum_{g=0}^{2}\pi_g L_g(R;\theta). $$ The genotype probabilities satisfy $\pi_g\geq 0$ and $\sum_{g=0}^{2}\pi_g=1$.

The read likelihoods describe how the reads would be distributed if a particular genotype were present. The genotype probabilities describe how much weight the model gives each possible genotype before seeing these reads. They may be specified directly or derived from a population model. For example, under Hardy–Weinberg equilibrium with alternate-allele frequency $p$, the probabilities for $G=0,1,2$ are $(1-p)^2$, $2p(1-p)$, and $p^2$. Whatever their source, the weights used in the marginal likelihood must be probabilities under the model being analyzed.

A Read-Count Model

A simple model for read counts illustrates the calculation. Suppose $n$ reads cover a site and $K=k$ of them show the alternate allele. Conditional on genotype $g$, assume the alternate-read count has a binomial distribution with alternate-read probability $q_g$:

K\mid G=g \sim \operatorname{Binomial}(n,q_g), \qquad P(K=k\mid G=g)=\binom{n}{k}q_g^k(1-q_g)^{n-k}.

This model assumes that, conditional on $g$, the reads are independent and have the same alternate-read probability $q_g$. In a simple symmetric-error model with per-read error probability $e$, the probabilities are $q_g=e$ for $g=0$, $q_g=1/2$ for $g=1$, and $q_g=1-e$ for $g=2$. Other read models may use different probabilities or additional parameters; the marginalization principle is unchanged.

When genotype $g$ is not observed, the probability of the read count is obtained by summing the three conditional binomial probabilities weighted by their genotype probabilities. This is the law of total probability applied to the unobserved genotype. The binomial coefficient is part of each read likelihood; it may be omitted only when comparing likelihoods for the same observed $n$ and $k$ and it is constant with respect to the parameter being compared.

Worked Example: Marginal Likelihood for One Read Count

Consider hypothetical data with $n=10$ reads and $k=6$ alternate reads. Suppose the conditional alternate-read probabilities for $G=0,1,2$ are $0.02, 0.50, 0.98$, respectively, and the genotype probabilities are $\pi_0=0.25$, $\pi_1=0.50$, and $\pi_2=0.25$. Calculate the marginal likelihood of the observed count.

For each genotype, the conditional read likelihood is $\binom{10}{6}q_g^6(1-q_g)^4$, with $\binom{10}{6}=210$. Substituting the genotype-specific probabilities gives:

\begin{aligned} L_0 &=210(0.02)^6(0.98)^4 \approx 1.240\times 10^{-8},\\ L_1 &=210(0.50)^6(0.50)^4 =0.2050781,\\ L_2 &=210(0.98)^6(0.02)^4 \approx 2.976\times 10^{-5}. \end{aligned}

Weight these likelihoods by their genotype probabilities and add:

L_{\mathrm{marg}} =0.25L_0+0.50L_1+0.25L_2 \approx 0.25(1.240\times10^{-8})+0.50(0.2050781)+0.25(2.976\times10^{-5}) \approx 0.1025465.

Thus, under the specified model, the probability of observing exactly six alternate reads among ten is approximately $0.1025$. This is a probability for the observed data under the model, not a probability that a particular genotype or parameter value is true.

Marginal Likelihood and Genotype Uncertainty

Marginalizing over genotype preserves the contribution from every possible state, including states with small but nonzero support. The posterior probability of a genotype is a different quantity: it describes the relative support for that genotype after combining the read evidence with the genotype probabilities. Bayes' rule gives

P(G=g\mid R;\theta) =\frac{\pi_g L_g(R;\theta)} {\sum_{h=0}^{2}\pi_h L_h(R;\theta)}.

The denominator is the marginal likelihood. The numerator is the weighted contribution from genotype $g$. Posterior probabilities sum to one across the three genotypes. They can summarize genotype uncertainty, but they are not themselves the marginal likelihood of the reads.

Worked Example: Posterior Genotype Probabilities

Use the hypothetical read-count model and data from the preceding example. The weighted contributions to the marginal likelihood are approximately $0.25L_0=3.099\times10^{-9}$, $0.50L_1=0.1025391$, and $0.25L_2=0.000007441$. Divide each contribution by the marginal likelihood, $0.1025465$, to find the posterior probabilities.

\begin{aligned} P(G=0\mid K=6)&\approx \frac{3.099\times10^{-9}}{0.1025465} \approx 0.0000000302,\\ P(G=1\mid K=6)&\approx \frac{0.1025391}{0.1025465} \approx 0.9999274,\\ P(G=2\mid K=6)&\approx \frac{0.000007441}{0.1025465} \approx 0.00007256. \end{aligned}

The heterozygous state has nearly all the posterior probability under these assumptions. Nonetheless, the marginal likelihood remains the sum of all three weighted contributions. A genotype call that reports only $G=1$ communicates the most supported state, but not the small remaining uncertainty.

What a Hard Call Changes

A hard genotype call selects one genotype, often the state with the largest posterior probability, and treats that state as if it were known. If the selected state is $\hat{g}$, the read likelihood conditional on that call is $L_{\hat{g}}(R;\theta)$. This is not generally the marginal likelihood $\sum_g\pi_g L_g(R;\theta)$: the hard-call quantity conditions on one state, whereas the marginal likelihood averages over states using their model probabilities.

In the first worked example, the most probable genotype is $G=1$. The read likelihood conditional on this hard call is $L_1=0.2050781$, while the marginal likelihood is $0.1025465$. These values differ because the marginal likelihood weights the read likelihoods for all three states by the genotype probabilities. Neither number is a posterior probability of the genotype.

The appropriate likelihood depends on what is observed and what model is being used. If the reads are observed but genotype is not, use the marginal likelihood. If genotype is genuinely known and the analysis conditions on it, a genotype-conditional read likelihood is appropriate. If a hard call was produced from the same reads, treating it as known is a data-reduction choice, not a way to retain all the read information. In particular, a genotype call should not be substituted into a model as though the uncertainty that produced the call did not exist.

Worked Example: The Effect of Genotype Probabilities

Suppose the likelihoods of a hypothetical read pattern under $G=0,1,2$ are $0.01,0.20,0.04$. Compare two models for genotype probabilities: Model A has $(\pi_0,\pi_1,\pi_2)=(0.80,0.19,0.01)$, and Model B has $(0.25,0.50,0.25)$. Calculate the marginal likelihood and posterior probability of each genotype under both models.

Under Model A, the marginal likelihood is:

L_{\mathrm{marg},A}=0.80(0.01)+0.19(0.20)+0.01(0.04) =0.008+0.038+0.0004=0.0464.

The posterior probabilities are the three weighted terms divided by $0.0464$: $P(G=0\mid R)=0.008/0.0464\approx0.1724$, $P(G=1\mid R)=0.038/0.0464\approx0.8190$, and $P(G=2\mid R)=0.0004/0.0464\approx0.008621$.

Under Model B, the marginal likelihood is:

L_{\mathrm{marg},B}=0.25(0.01)+0.50(0.20)+0.25(0.04) =0.0025+0.1000+0.0100=0.1125.

The corresponding posterior probabilities are $0.0025/0.1125\approx0.02222$, $0.1000/0.1125\approx0.8889$, and $0.0100/0.1125\approx0.08889$. Both models assign the highest posterior probability to $G=1$, but their marginal likelihoods and posterior distributions differ. The genotype probabilities are part of the model, so changing them can change both quantities even when the read likelihoods stay fixed.

Marginalizing Across Independent Individuals

For multiple individuals, each with an unobserved genotype and observed reads, marginalize each individual's genotype. If the individuals' genotypes are independent under the model and their read data are conditionally independent given their genotypes and the model parameters, the joint marginal likelihood is the product of the individual marginal likelihoods. Independence of the read data alone is not enough: the genotype probabilities must also be specified consistently with the model.

L_{\mathrm{marg}}(\theta;R_1,\ldots,R_m) =\prod_{i=1}^{m}\left\{\sum_{g=0}^{2}\pi_g L_g(R_i;\theta)\right\}.

This product can also be written as a sum over all possible combinations of the individuals' genotypes. For two individuals, independent genotype probabilities make the sum factor into a product of sums. The individual genotype states should not be collapsed into one shared state unless the model actually specifies a shared genotype.

Worked Example: Marginal Likelihood for Two Individuals

For two hypothetical individuals, suppose the conditional read likelihoods for $G=0,1,2$ are $(0.10,0.60,0.20)$ for individual A and $(0.40,0.20,0.05)$ for individual B. Both have genotype probabilities $(0.20,0.50,0.30)$, and assume their genotypes and read data are independent under the model. Individual A's marginal likelihood is

$$ m_A=0.20(0.10)+0.50(0.60)+0.30(0.20) =0.02+0.30+0.06=0.38. $$

Individual B's marginal likelihood is

$$ m_B=0.20(0.40)+0.50(0.20)+0.30(0.05) =0.08+0.10+0.015=0.195. $$

The joint marginal likelihood is therefore $m_A m_B=0.38(0.195)=0.0741$. To verify the sum-over-combinations interpretation, the term for genotype pair $(g_A,g_B)$ is $\pi_{g_A}\pi_{g_B}L_{g_A}(R_A)L_{g_B}(R_B)$. Summing these nine terms factors as

\sum_{g_A=0}^{2}\sum_{g_B=0}^{2} \pi_{g_A}\pi_{g_B}L_{g_A}(R_A)L_{g_B}(R_B) = \left(\sum_{g_A=0}^{2}\pi_{g_A}L_{g_A}(R_A)\right) \left(\sum_{g_B=0}^{2}\pi_{g_B}L_{g_B}(R_B)\right) =0.38(0.195)=0.0741.

The calculation depends on the stated independence assumptions. If genotype states are dependent, the joint genotype probabilities cannot in general be replaced by the product $\pi_{g_A}\pi_{g_B}$.

A Short Justification of the Marginalization Rule

The genotypes $G=0,1,2$ are mutually exclusive and exhaustive under the model. Therefore, the event $R$ can be partitioned into the three possibilities that the reads occurred with each genotype. The events “$R$ and $G=g$” are disjoint, so their probabilities add. Applying the definition of conditional probability to each term gives

\begin{aligned} P(R) &=\sum_{g=0}^{2}P(R,G=g)\\ &=\sum_{g=0}^{2}P(R\mid G=g)P(G=g)\\ &=\sum_{g=0}^{2}L_g(R;\theta)\pi_g. \end{aligned}

This argument requires that the genotype states cover all possibilities allowed by the model and that the genotype probabilities are defined for those states. It does not require the genotype to be assigned a single best value first. In fact, choosing a hard call before evaluating the likelihood replaces the sum with a single conditional contribution, which is a different calculation.

Common Mistake and Practical Checks

A common mistake is to compute the most probable genotype and then use its read likelihood as though it were the likelihood of the observed reads. The reasoning may be: “$G=1$ has the highest posterior probability, so $P(R\mid G=1)$ is the likelihood for these data.” The correction is that the posterior identifies a genotype supported by the data and genotype probabilities; it does not make that genotype observed. For uncertain genotypes, the observed-data likelihood still sums over states with their genotype-probability weights.

Before using a marginal likelihood, check that the conditional read model matches the data, the genotype probabilities sum to one, and the genotype states included in the sum are exhaustive. Also check which quantities are fixed and which are parameters: if a genotype probability depends on an unknown parameter, its contribution belongs in the marginal likelihood and must not be treated as a fixed number without justification. Finally, keep the target clear: a likelihood for observed reads is not a posterior probability for a genotype, and a hard-call likelihood is not generally the same quantity.

Key takeaway: When genotype is unobserved, the likelihood for the reads is the sum of the genotype-conditional read likelihoods weighted by genotype probabilities. A hard call selects one state and discards the uncertainty that marginalization retains.

Check Your Understanding

Use the definitions and assumptions in this tutorial to answer the following questions.

In a diploid two-allele model, what does $G=1$ represent, and what does $L_1(R;\theta)$ mean?
Why must the genotype probabilities sum to one before they can be used as weights in a marginal likelihood?
Suppose the largest posterior genotype probability is for $G=2$. Explain why the read likelihood conditional on $G=2$ is not automatically the marginal likelihood.
For two individuals, what independence assumptions justify multiplying their individual marginal likelihoods?
If two models have identical read likelihoods but different genotype probabilities, can their marginal likelihoods differ? Explain using the marginalization formula.