SNP calling: heuristics, likelihood and Bayes
The practice
A genome is sequenced, and you have to say what this individual carries at every position. You will do it three times — with thresholds you choose, with maximum likelihood, and with Bayes — and score all three on the same reads. Everything below is simulated, so the truth is known and every call can be marked.
The genome
A locus is polymorphic in the population with probability v, the variant rate. At a polymorphic locus the alternative allele has frequency q, and what this individual carries is drawn from
P(hom ref) = p² + F·p·q P(het) = 2pq(1 − F) P(hom alt) = q² + F·p·q
with p = 1 − q. F is Wright's inbreeding coefficient: at F = 0 the genotypes are in Hardy-Weinberg proportions, and as F rises heterozygotes are emptied into the two homozygous classes. Note what follows — a polymorphic locus can leave this individual homozygous for the reference, and at such a position there is nothing in the reads to find. Calls are scored against the individual's genotype, not against whether the locus is variable in the population.
The sequencing
Every position is sequenced, variable or not. It receives a Poisson number of reads with the mean depth you set, so some positions arrive with three reads and others with seventeen. Each read draws one of the two alleles the individual carries, with probability one half each, and carries a Phred quality Q of its own. The base is then misread with probability
e = 10−Q/10
becoming one of the three other bases. The quality is honest by construction: a base reported at Q20 really is wrong one time in a hundred. Whether a caller believes it is another matter.
What you have to build
Three callers, one in each tab, all scored on these same reads. Each reports the same two tables: whether a variant is there at all, and which genotype was called. Start at 10×, where every method looks good, then drop the depth to 4× and look again — the differences between the three live at low depth.
Heuristic
Count the reads carrying the alternative base, divide by the reads carrying either allele, and compare the ratio with two thresholds somebody chose: below the first it is homozygous reference, above the second homozygous alternative, in between a heterozygote. Positions with too few reads are left uncalled, and reads below the quality cutoff are thrown away rather than weighed — throwing information away is what a heuristic does with a quality value.
No probability is computed anywhere. The thresholds are the whole model, and finding out what they cost is the exercise.
Maximum likelihood
For each of the three genotypes, ask how probable these very reads would be if that genotype were the truth, and take the genotype that makes them most probable. One read showing base b from an individual carrying alleles A₁ and A₂ contributes
P(b | A) = 1 − e if b is A, else e/3 P(b | A₁A₂) = ½ P(b | A₁) + ½ P(b | A₂)
and the reads multiply, since each is independent. Every quality is used as a weight here rather than as a filter: a poor base argues for its genotype more weakly instead of not at all.
Bayes
The same likelihoods as the previous tab, multiplied by a prior over the three genotypes and normalised:
P(G | reads) ∝ P(G) · P(reads | G)
The prior is where the variant rate, the allele frequency and F belong. Most positions of a genome are not variable at all, so
P(hom ref) = (1 − v) + v(p² + F·p·q) P(het) = v·2pq(1 − F) P(hom alt) = v(q² + F·p·q)
These are the model's assumptions, not the genome's. They start at the values the simulation used; set them wrong on purpose and watch what a confident prior does to a weak piece of evidence. That is the lesson of the tab: a posterior is only as good as the prior it was built on.