Reference for the statistical methodology used in joint genotype calling

Question

I am searching for a good reference that explains the statistical methodology used in what is variously known as Joint Genotyping, Cohort Genotyping or Multi-sample Genotyping. One with a pseudo code implementation would also be ideal.

I have found several discussing the approach in very general terms and the relative benefits of the technique but nothing with actual equations or algorithms. Failing that, if there is an open source implementation that people would consider to be a reference implementation and that has excellent code documentation that would also be great.

Answer 1

To the best of my knowledge, the first joint calling model was published by Le & Durbin (2011) and refined in Li (2012). They come with detailed derivation. Both samtools and older GATK are using this model. GATK now uses a different model which only works with high-coverage data but is faster to compute. There are other alternatives (e.g. in freebayes).

Answer 2

Here is a recent example of the derivation of the genotype model from Octopus genotype calling paper:

All samples have known ploidy and copy number, so the likelihood function of reads $R$ given genotype $g$ is:

$p(R|g) = \mathop {\prod }\limits_{n = 1}^{|R|} \frac{1}{{|g|}}\mathop {\sum }\limits_{i = 1}^{|g|} p(r_n|g_i)$

where $|g|$ is the ploidy and $|R|$ is the number of reads. The joint genotype posterior for $S$ samples is therefore given by:

$p({\bf{g}}|{{R}},{\cal{M}}_{{g}}) \propto p({\bf{g}}|{\cal{M}}_g)\mathop {\prod }\limits_{s = 1}^S p({{R}}_s|{\bf{g}}_s)$

where the genotype prior model, $\cal{M}𝑔$, is either the uniform or HWE-coalescent prior. Unfortunately, the number of genotype combinations g grows exponentially in the number of samples S, so we cannot evaluate the full posterior distribution in general. Therefore, other than for trivial cases, we first approximate the sample marginal genotype posterior distribution under the HWE model (without mutations) $𝑝(𝐠𝑠|𝑅,\cal{M}𝑔)$ and use these marginal probabilities to select K genotype combinations $𝐠1,…,𝐠𝐾$ (K is user-defined) to evaluate under the full joint genotype model. The approximate posterior marginals are computed with expectation maximization.

The source code can be found here.