Genotype Likelihoods for Rare Sequencing Calls

Why a Rare-Variant Call Can Be Uncertain

At a rare variant site, a small number of alternate-allele reads may be compatible with more than one diploid genotype. One alternate read, for example, might come from a true heterozygote, or it might occur while the person is homozygous for the reference allele. A sequencing pipeline must summarize that ambiguity, not just count reads and announce a genotype.

A genotype likelihood records how compatible the observed reads are with each possible genotype. A posterior genotype probability combines that read evidence with prior probabilities for the genotypes. These are different quantities: the likelihood describes evidence from the reads, while the posterior is a probability distribution over genotypes conditional on the reads and the chosen prior.

Definition: For a diploid, biallelic site, let \(G\in\{0,1,2\}\) be the number of alternate alleles, and let \(R\) denote the observed sequencing reads. The genotype likelihood for genotype \(g\) is \(L_g=P(R\mid G=g)\). Given genotype prior probabilities \(\pi_g=P(G=g)\), the posterior genotype probability is \(P(G=g\mid R)\). The three posterior probabilities are nonnegative and sum to 1.

A Read Model and Genotype Priors

To make the calculations explicit, consider a simple model for a single diploid, biallelic site. Suppose \(D\) reads cover the site and \(K\) of them show the alternate allele. Let \(q_g\) be the probability that a read shows the alternate allele if the true genotype is \(g\). For the examples, assume \(q_0=0.01\), \(q_1=0.50\), and \(q_2=0.99\). This symmetric model allows a small probability of observing the opposite allele from either homozygote and assumes a heterozygote produces alternate reads half the time.

The model assumes that, conditional on the genotype, reads are independent and have the same alternate-read probability. It treats the read depth and the values \(q_g\) as known, and does not model allele-specific mapping or other sources of read dependence. Those assumptions make the calculations transparent; actual likelihoods may use more detailed read-level information.

Read-count likelihood: Under this binomial read model, the likelihood of observing \(K\) alternate reads among \(D\) reads under genotype \(g\) is

L_g=P(K\mid G=g)=\binom{D}{K}q_g^K(1-q_g)^{D-K}.

The binomial coefficient is included here for clarity. Because it is the same for all three genotypes at a given site, it can be omitted when calculating posterior probabilities: it cancels in the normalization. Genotype likelihoods are sometimes reported after rescaling, for example by dividing each likelihood by the largest one. Such a common rescaling preserves their relative evidence, but does not turn the likelihoods into genotype probabilities.

A prior is needed to turn likelihoods into posterior probabilities. If the alternate-allele frequency in the relevant population is \(p\), and genotypes follow Hardy–Weinberg proportions, the genotype priors are \(P(G=0)=(1-p)^2\), \(P(G=1)=2p(1-p)\), and \(P(G=2)=p^2\). These priors rely on the stated frequency and population assumptions; they are not automatically appropriate in every sample or ancestry group.

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

The denominator is the total probability of the observed read data after averaging over the possible genotypes using their priors. It ensures that the three posterior probabilities sum to one. In particular, a genotype can have a relatively large likelihood but a modest posterior if its prior probability is small.

Worked Calculations

Worked Example: One Alternate Read at a Rare Site

Suppose a hypothetical sample has \(D=8\) reads and \(K=1\) alternate read. Take \(p=0.01\), use the Hardy–Weinberg priors, and use the read model above. Calculate the posterior probability of a heterozygous genotype.

The genotype priors are \(\pi_0=(0.99)^2=0.9801\), \(\pi_1=2(0.01)(0.99)=0.0198\), and \(\pi_2=(0.01)^2=0.0001\). The likelihoods are \(L_0=\binom{8}{1}(0.01)(0.99)^7\), \(L_1=\binom{8}{1}(0.50)^8\), and \(L_2=\binom{8}{1}(0.99)(0.01)^7\).

Multiplying each likelihood by its prior gives the unnormalized weights: \(L_0\pi_0\approx0.07308138\), \(L_1\pi_1=8(0.50)^8(0.0198)=0.00061875\), and \(L_2\pi_2\approx7.92\times10^{-18}\). Their sum is approximately \(0.07370013\). Therefore,

P(G=1\mid K=1)= \frac{0.00061875}{0.07370013} \approx 0.008396.

Under this model and prior, the posterior probability of a heterozygote is about 0.84%. The single alternate read is not enough to overcome the much larger prior probability of the reference homozygote. This is a conditional probability under the specified model, not a general probability that a rare variant call is correct.

Worked Example: Evidence Strongly Favoring a Heterozygote

Now suppose a hypothetical sample has \(D=8\) and \(K=5\), with the same read model and \(p=0.01\). The priors remain \(0.9801\), \(0.0198\), and \(0.0001\). Calculate the posterior genotype probabilities.

For \(G=0\), the weighted likelihood is \(\binom{8}{5}(0.01)^5(0.99)^3(0.9801)\approx5.326\times10^{-9}\). For \(G=1\), it is \(\binom{8}{5}(0.50)^8(0.0198)=0.00433125\). For \(G=2\), it is \(\binom{8}{5}(0.99)^5(0.01)^3(0.0001)\approx5.326\times10^{-9}\).

The total weight is approximately \(0.004331260651\). Dividing each weight by that total gives \(P(G=0\mid K=5)\approx0.000001230\), \(P(G=1\mid K=5)\approx0.99999754\), and \(P(G=2\mid K=5)\approx0.000001230\). The read evidence strongly favors the heterozygous genotype, even though the prior probability of a heterozygote is only 0.0198.

The likelihoods for the two homozygotes are equal in this example because the read model is symmetric and five of eight reads show the alternate allele. Their posterior probabilities are also equal because their prior weights, \(0.9801\) and \(0.0001\), are not equal? In fact, the weighted likelihoods shown are equal under these specific priors and the symmetric read model: the likelihood ratio between \(G=0\) and \(G=2\) is offset by their prior ratio. This equality is a feature of these chosen values, not a general rule.

From a Posterior Distribution to a Genotype Call

One way to turn the posterior distribution into a single genotype is the maximum a posteriori (MAP) call: choose the genotype with the largest posterior probability. This rule identifies the most probable genotype under the model and prior, but does not indicate how much probability remains assigned to alternatives. A MAP call is therefore not the same thing as certainty.

A different useful summary is the posterior expected genotype dosage, the mean number of alternate alleles under the posterior distribution. It preserves some information about uncertainty: probabilities spread across genotypes contribute proportionally rather than being replaced by a single category.

Posterior expected dosage: If \(P(G=g\mid R)\) is available for \(g=0,1,2\), the posterior expected alternate-allele count is

E(G\mid R)=\sum_{g=0}^{2}g\,P(G=g\mid R) =P(G=1\mid R)+2P(G=2\mid R).

Worked Example: A Hard Call Can Hide Uncertainty

For a hypothetical site with two reads, one alternate read, and the same read model, the likelihoods are \(L_0=2(0.01)(0.99)=0.0198\), \(L_1=2(0.50)(0.50)=0.50\), and \(L_2=2(0.99)(0.01)=0.0198\). With \(p=0.01\), the unnormalized posterior weights are \(0.0198(0.9801)=0.01940598\), \(0.50(0.0198)=0.0099\), and \(0.0198(0.0001)=0.00000198\). Their sum is \(0.02930796\).

After normalization, the posterior probabilities are approximately \(P(G=0\mid R)=0.66214\), \(P(G=1\mid R)=0.33779\), and \(P(G=2\mid R)=0.0000676\). The MAP call is \(G=0\), since its posterior probability is largest. But the posterior expected dosage is \(0(0.66214)+1(0.33779)+2(0.0000676)\approx0.33793\) alternate alleles.

The hard call says only “reference homozygote” and discards the substantial posterior probability assigned to a heterozygote. The expected dosage records that the evidence is uncertain, although it is not a replacement for the full three-probability posterior when that full distribution is needed.

Why Retain Genotype Uncertainty?

A hard call is convenient for tables and some downstream analyses, but converting every posterior distribution to its most likely genotype can cause information loss. Individuals with posterior distributions such as \((0.51,0.48,0.01)\) and \((0.99,0.01,0)\) would both receive a reference-homozygote MAP call, even though the evidence for that call differs greatly. Retaining the likelihoods or posterior probabilities lets later analyses distinguish these cases.

When posterior probabilities are used, their interpretation depends on the prior as well as the read model. Genotype likelihoods can be retained when a later analysis will specify its own model for genotype frequencies; posterior probabilities already incorporate the prior used to calculate them. A downstream method must account for that distinction rather than treating a posterior as if it were a prior-free property of the reads.

Expected dosage is a compact uncertainty-aware summary, but it also has limits. Two individuals can have the same expected dosage and different posterior distributions. For example, a distribution concentrated on \(G=1\) and a distribution split between \(G=0\) and \(G=2\) can have the same mean. If the analysis depends on more than the mean genotype count, retaining the full probability vector is more informative.

Deriving the Posterior Odds

The way prior information changes the comparison between two genotypes follows directly from Bayes’ rule. For any two possible genotypes \(g\) and \(h\), write the posterior probability for each as its likelihood times its prior divided by the same normalizing denominator. Taking their ratio cancels that denominator:

\frac{P(G=g\mid R)}{P(G=h\mid R)} = \frac{P(R\mid G=g)\pi_g}{P(R\mid G=h)\pi_h} = \frac{P(R\mid G=g)}{P(R\mid G=h)} \frac{\pi_g}{\pi_h}.

Thus, posterior odds equal the likelihood ratio multiplied by the prior odds. Read evidence can favor a heterozygote, while a rare-allele prior can favor the reference homozygote; the posterior reflects both. This identity also explains why the posterior in the one-read example was small: the read evidence did not supply enough relative support to offset the prior odds against a heterozygote.

Common Mistakes

Calling a likelihood a probability of the genotype. \(P(R\mid G=g)\) describes the probability of the observed reads if genotype \(g\) were true. It is not \(P(G=g\mid R)\).
Ignoring the prior in a posterior calculation. A likelihood can be multiplied by a genotype prior, but the resulting weights must also be divided by their sum to obtain posterior probabilities.
Treating the MAP call as certain. The largest posterior probability identifies the most probable genotype under the model. It can still leave considerable probability on other genotypes.
Assuming an expected dosage is a full genotype distribution. The mean summarizes the posterior but cannot generally reconstruct its three probabilities.
Interpreting a posterior without its assumptions. The posterior depends on the read model, the genotype prior, and the population-frequency assumptions used to obtain that prior.

Key takeaway: Genotype likelihoods quantify how well each genotype explains the reads; posterior probabilities combine those likelihoods with genotype priors. A hard call keeps only the most probable genotype, whereas retaining posterior probabilities—or a suitable dosage summary—preserves information about uncertainty.

Check Your Understanding

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

In your own words, distinguish \(P(R\mid G=g)\) from \(P(G=g\mid R)\).
Under Hardy–Weinberg proportions, what are the prior probabilities of \(G=0,1,2\) when the alternate-allele frequency is \(p\)?
Why must the prior-weighted likelihoods be normalized before they can be interpreted as posterior genotype probabilities?
If the posterior probabilities for \(G=0,1,2\) are \(0.60, 0.35, 0.05\), what is the MAP call and the posterior expected alternate-allele dosage?
Give one reason to retain the full posterior distribution rather than keeping only a MAP genotype or expected dosage.