From Sequencing Reads to a Genotype Call
A sequencing-based genotype call is an inference from reads that cover a genomic position. At a diploid site with two alleles, a person’s underlying genotype has zero, one, or two copies of the alternate allele. Sequencing does not reveal that genotype without uncertainty: a read can contain an error, and the number of reads supporting each allele can vary. The central question is therefore not simply which allele appears most often, but how compatible the observed reads are with each possible genotype.
Read depth is the number of usable sequencing reads covering the position. The allele-specific read counts record how many of those reads support each allele. For a biallelic site, write \(n_R\) for the reference-supporting reads and \(n_A\) for the alternate-supporting reads, so the total depth is \(n=n_R+n_A\). These counts provide evidence for a genotype; they are not themselves the genotype.
The observed alternate fraction, \(n_A/n\), is a useful description of the reads, but it is not a genotype probability. A fraction near one-half often supports a heterozygous call, while a fraction near zero or one often supports a homozygous call. How strong that support is depends on depth, sequencing error, and whether the model’s assumptions are reasonable. For example, one alternate read among two total reads and 50 alternate reads among 100 total reads have the same fraction but do not provide the same amount of evidence.
Genotype Likelihoods Under a Simple Model
To compare genotypes, specify how reads are expected to arise under each one. Consider a deliberately simple model for a diploid, biallelic site. Assume reads are independent, the probability of a sequencing error is a known constant \(\epsilon\), errors are equally likely to change either allele to the other, and there is no allele-specific mapping or sequencing bias. Under these assumptions, the probability that a read supports the alternate allele is \(q_0=\epsilon\) for \(g=0\), \(q_1=1/2\) for \(g=1\), and \(q_2=1-\epsilon\) for \(g=2\).
The heterozygote has expected alternate probability one-half because each read is assumed to sample either chromosome with equal probability; symmetric errors do not change that probability. In real analyses, callers may use more detailed models for read quality, mapping, or bias. The simple model is useful for seeing what a genotype likelihood means and what assumptions are involved.
This is the binomial likelihood for \(n_A\) alternate-supporting reads out of \(n\), with alternate-read probability \(q_g\). The full binomial probability also contains the factor \(\binom{n}{n_A}\); since it is the same for every candidate genotype, it cancels when comparing or normalizing their likelihoods. A larger likelihood means the observed counts are more compatible with that genotype under the chosen model. It is not, by itself, the probability that the genotype is true.
Worked Example: Comparing Genotypes at Low Depth
Suppose a hypothetical site has 9 reference-supporting reads and 1 alternate-supporting read, for depth \(n=10\). Use the simple model with \(\epsilon=0.01\). Compare the three genotypes and, for this calculation, assign them equal prior probabilities.
For \(g=0\), \(q_0=0.01\), so \(L(0)=0.01^1(0.99)^9=0.009135\). For \(g=1\), \(q_1=0.5\), so \(L(1)=0.5^1(0.5)^9=0.0009766\). For \(g=2\), \(q_2=0.99\), so \(L(2)=0.99^1(0.01)^9=9.900\times10^{-19}\).
The largest likelihood is for \(g=0\), so the maximum-likelihood call is reference homozygous. With equal priors, the normalized likelihoods are also the posterior genotype probabilities for this calculation. The probability assigned to \(g=0\) is \(0.009135/(0.009135+0.0009766+9.900\times10^{-19})\approx0.9034\). Thus the probability assigned to the other genotypes combined is about \(0.0966\). The call favors reference homozygosity, but the one alternate read leaves appreciable uncertainty.
Why Depth Changes the Strength of Evidence
At the same allele fraction, more reads usually provide stronger evidence because the count pattern is less likely to arise by chance under a competing genotype. This conclusion depends on the reads behaving approximately as the model assumes. If reads are correlated, errors are more common than assumed, or one allele is systematically underrepresented, simply counting more reads can overstate the evidence.
Worked Example: Same Fraction, Different Depth
Compare two hypothetical samples at sites modeled with \(\epsilon=0.01\): sample A has 9 reference and 1 alternate read, while sample B has 36 reference and 4 alternate reads. Both have a 10% alternate fraction. For each, compare reference homozygosity with heterozygosity, using equal priors and the same independent-read model.
For sample A, the likelihoods are \(L_A(0)=0.01(0.99)^9=0.009135\) and \(L_A(1)=0.5^{10}=0.0009766\). Their ratio is \(L_A(1)/L_A(0)=0.0009766/0.009135\approx0.1069\). The heterozygote is less supported, but its likelihood is not negligible relative to the reference homozygote likelihood.
For sample B, \(L_B(0)=0.01^4(0.99)^{36}\approx6.964\times10^{-9}\), while \(L_B(1)=0.5^{40}\approx9.095\times10^{-13}\). Thus \(L_B(1)/L_B(0)\approx0.0001306\). Under equal priors, and with the alternate homozygote likelihood negligible in both cases, the probability assigned to a non-reference-homozygous genotype is about 0.0966 for sample A and 0.0001306 for sample B. The corresponding Phred-scaled error values are about 10.15 and 38.84. In this model, the greater depth gives much stronger support for the reference-homozygous call, despite the matching allele fractions.
A likelihood ratio is a comparison of evidence, not a direct measure of posterior probability unless prior probabilities are also specified. To obtain posterior probabilities, multiply each genotype likelihood by its prior probability and normalize the resulting values to sum to one. The priors represent genotype probabilities before considering these reads. They may affect the final probabilities, particularly when the read evidence is weak.
Worked Example: Likelihoods and Priors Are Different
Suppose a hypothetical site has 8 reference-supporting reads and 2 alternate-supporting reads, with \(\epsilon=0.01\). The likelihoods are \(L(0)=0.01^2(0.99)^8\approx0.00009227\), \(L(1)=0.5^{10}=0.0009766\), and \(L(2)=0.99^2(0.01)^8\approx9.801\times10^{-17}\). The heterozygote has the largest likelihood.
Now suppose, purely for illustration, the prior probabilities for \(g=0,1,2\) are 0.90, 0.09, and 0.01. Multiplying likelihood by prior gives weights of approximately \(0.00008304\), \(0.00008789\), and \(9.801\times10^{-19}\). After dividing each weight by their sum, the posterior probability for \(g=1\) is about \(0.514\). It is still the most probable genotype, but its posterior probability is only about 51.4%, not the likelihood value \(0.0009766\). The other genotypes together have probability about 48.6%.
This example shows why a likelihood cannot be read as the probability that a genotype is correct. The likelihood is conditional on a genotype; a posterior probability also depends on the prior probabilities and the set of genotypes being considered.
Genotype Quality as a Phred-Scaled Error Probability
A genotype call is a selected genotype, while genotype quality describes uncertainty in that call. A common convention expresses the probability that the selected call is wrong on the Phred scale. If \(P_{\mathrm{error}}\) is the probability that the call is incorrect under the caller’s model, genotype quality is
For example, an error probability of \(0.01\) gives \(GQ=20\), since \(-10\log_{10}(0.01)=20\). An error probability of \(0.001\) gives \(GQ=30\). Each increase of 10 Phred units corresponds to an error probability ten times smaller. These values are measures conditional on the model and the information used by the caller; they do not guarantee that the call is correct.
When posterior genotype probabilities are available, a natural choice for \(P_{\mathrm{error}}\) is one minus the posterior probability of the selected genotype. If likelihoods alone are used, an implementation must define how it turns them into an error probability; likelihoods should not be silently treated as posterior probabilities. Details can vary with the calling method, so a quality score should be interpreted according to how it was produced.
Deriving the Read-Count Likelihood
The binomial form follows directly from the independent-read assumption. Under a proposed genotype \(g\), each read independently supports the alternate allele with probability \(q_g\), and supports the reference allele with probability \(1-q_g\). For one particular ordering of the reads, the probability of observing \(n_A\) alternate and \(n_R\) reference reads is \(q_g^{n_A}(1-q_g)^{n_R}\). There are \(\binom{n}{n_A}\) possible orderings with those same counts, so the probability of the counts is
The combinatorial factor depends on the observed counts but not on \(g\). Therefore it does not affect which genotype has the greatest likelihood, nor does it affect normalized probabilities across genotypes when the same data and model are used. Removing that common factor yields \(L(g)\) above. This justifies using the simpler expression for genotype comparisons.
Assumptions and Common Pitfalls
The simple calculations here rely on a biallelic diploid site, correctly assigned reference and alternate reads, independent reads, a specified error rate, and no allele-specific bias. Real read data can depart from these assumptions because of mapping ambiguity, local sequence context, duplicated reads, contamination, or unequal representation of alleles. Such departures can make a call appear more certain than it should. The model is an aid to interpreting evidence, not a substitute for checking whether the data and assumptions are credible.
- Mistake: treating the most frequent allele as the genotype. A heterozygote can produce more reads from one allele than the other by chance. Compare the observed counts with the expected read distributions for all candidate genotypes.
- Mistake: saying a likelihood is the probability the genotype is true. A likelihood is \(P(\text{reads}\mid\text{genotype})\), not \(P(\text{genotype}\mid\text{reads})\). The latter also depends on priors.
- Mistake: interpreting a high quality score as certainty. The score reflects estimated error under a model. A misspecified error rate or biased reads can undermine that estimate.
- Mistake: comparing allele fractions without considering depth. Fractions alone discard information about how many reads support the pattern.
Check Your Understanding
Use the definitions and assumptions above to answer the following questions.
- A site has 24 reference-supporting reads and 6 alternate-supporting reads. What is its depth and alternate-read fraction, and why are these not themselves a genotype call?
- Under the simple model with \(\epsilon=0.01\), what alternate-read probabilities \(q_g\) correspond to \(g=0,1,2\)?
- Explain why two samples with the same alternate-read fraction can provide different strengths of evidence for a genotype.
- In one sentence, distinguish a genotype likelihood from a posterior genotype probability.
- If a call has estimated error probability \(0.001\), calculate its Phred-scaled genotype quality and interpret that value.