biohazard-0.6.6: doc/genotyping.tex
\documentclass{article}
\usepackage{a4}
\usepackage{amsmath}
\usepackage{todonotes}
\usepackage{xargs}
\newcommandx{\beware}[2][1=]{\todo[nolist,noline,inline,linecolor=red,backgroundcolor=red!25,bordercolor=red,#1]{#2}}
\newcommandx{\idea}[2][1=]{\todo[nolist,noline,inline,linecolor=blue,backgroundcolor=blue!25,bordercolor=blue,#1]{#2}}
\newcommandx{\oops}[2][1=]{\todo[nolist,noline,inline,linecolor=yellow,backgroundcolor=yellow!25,bordercolor=yellow,#1]{#2}}
\newcommandx{\result}[2][1=]{\todo[nolist,noline,inline,linecolor=green,backgroundcolor=green!25,bordercolor=green,#1]{#2}}
\newcommand{\argmax}{\operatornamewithlimits{arg\,max}}
\title{Deamination Aware Genotype Calling}
\author{Udo Stenzel}
\begin{document}
\maketitle
\section{``History'' Of Error Models}
I tried to track down the logic behind the \texttt{samtools} and
\texttt{maq} error models, which supposedly go back to \texttt{CAP3}.
Near as I can tell, there is absolutely no reasoning behind any of it.
\texttt{CAP3} may have originated the idea of setting the probability of
$k$ errors to $p^{f(k)}$ where $f$ is a function that grows more slowly
than the identity function. The paper cited by \texttt{samtools}
doesn't actually mention any of that, though.
\texttt{SOAPsnp}\cite{soapsnp} is the first(?) implementation of the idea. Bases are
dealt with in order of increasing(!) quality, the quality score of
observation $k$ is scaled by $\theta^k$, where $\theta$ is a constant
parameter.
\texttt{Maq}\cite{maq} follows the same idea, but attempts some combinatorial
simplification. The derivation is rather complicated: It starts out
with a simplification (counting similar bases), then proceeds to apply
approximations (equal error rates), then ends up being incomprehensible
(weird effective error rate derived from quality scores). By that time,
it is no longer obvious whether that derivation makes any sense, and in
some cases, according to publications about
\texttt{samtools}\cite{samtools}, it
doesn't and fails catastrophically instead.
\texttt{Samtools}\cite{samtools} improves upon the \texttt{maq} model, where the
claimed reason is that the \texttt{maq} model is ill-behaved at high
coverage and high error rate. Unfortunately, the fix in
\texttt{samtools} is only a different approximation in the last step of
an equally convoluted derivation. The chief difference seems to be that
\texttt{maq} computes a strange quantity based on a sort of average
error rate, while \texttt{samtools} deals with bases in order of
decreasing(!) quality. By that time, it's unclear that the
combinatorial contortions provide any benefit.
The take home message is that we model error dependency by having a more
slowly growing exponent, that errors happening on different strands are
independent(?) from each other, and that the combinatorial
contortions in both the \texttt{maq} and the \texttt{samtools} model do
not seem to be useful. The order in which we touch the bases is up for
grabs, since there are two cases of prior art. We'll go with a simple
implementation like \texttt{SOAPsnp} (or more recently
\texttt{BSNP}\cite{bsnp}), which needs to be generalized for ancient
DNA.
\section{Damage Model}
Ancient DNA is conceptually modelled as double stranded with heavily
deaminated single stranded ends, whose lengths are distributed
geometrically. We apply the same simplification as in \cite{mapdamage},
in that we don't try to estimate the overhang length for any molecule,
but compute an effective deamination probability for every position in a
read.
For classically prepared ``double strand'' libraries, we get C to T
transitions near the 5' end and G to A transitions near the 3' end, and
the overhang lengths are equal. For single strand libraries, we get C
to T transitions near both ends, and the expected overhang lengths may
be different. The same models work for libraries treated with UDG.
While we cannot see the actual overhangs anymore, the overall damage
looks the same, but with much shorter overhangs.
We can now construct a universal damage model that applies to every
library and estimate parameters for it. Let $\sigma_d$, $\delta_d$,
$\lambda_d$ be the single stranded deamination rate, the double
stranded deamination rate and the length parameter at the 5' end,
respectively, for a double stranded model. Let $\sigma_s$, $\delta_s$
and $\lambda_s$ be the corresponding parameters for a single stranded
model, and $\kappa_s$ be the overhang length parameter at the 3' end in
a single stranded model. The probabilities that a position $i$ in a
read of length $l$ is in an overhang in a single stranded model, in a 5'
overhang, or in a 3' overhang in a double stranded model, are then:
\begin{align*}
P_s(i|l) &= \lambda_s^{i+1} + \kappa_s^{l-i} - \lambda_s^{i+1} \kappa_s^{l-i} \\
P_5(i|l) &= \lambda_d^{i+1} \\
P_3(i|l) &= \lambda_d^{l-i}
\end{align*}
From these, we compute the rates of C to T transversions and G to A transversions:
\begin{align*}
P_{CT}(i|l) &= \sigma_s P_s(i|l) + \delta_s (1-P_s(i|l) + \sigma_d P_5(i|l) + \delta_d (1-P_5(i|l)) \\
P_{GA}(i|l) &= \sigma_d P_3(i|l) + \delta_d (1-P_3(i|l))
\end{align*}
This yields the complete damage matrix used in the follwing sections:
\begin{equation}
D_i = \left( \begin{array}{cccc}
1 & 0 & P_{GA}(i) & 0 \\
0 & 1-P_{CT}(i) & 0 & 0 \\
0 & 0 & 1-P_{GA}(i) & 0 \\
0 & P_{CT}(i) & 0 & 1
\end{array} \right)
\end{equation}
To fit parameters to real data, we assume simple errors and no divergence. We separately fit a pure single stranded model, a pure
double stranded model, and a model without damage. At the maximum likelihood parameters, we compute the probabilities of either the
single stranded or double stranded model being correct. (For brevity, all parameters not mentioned here are assumed to be zero.)
\begin{align*}
\hat{\sigma_s}, \hat{\delta_s}, \hat{\lambda_s}, \hat{\kappa_s} &:= \argmax_{\sigma_s, \delta_s, \lambda_s, \kappa_s}
P(D | \sigma_s, \delta_s, \lambda_s, \kappa_s ) \\
\hat{\sigma_d}, \hat{\delta_d}, \hat{\lambda_d} &:= \argmax_{\sigma_d, \delta_d, \lambda_d}
P(D | \sigma_d, \delta_d, \lambda_d ) \\
P_s(D) &:= \frac{
P(D | \hat{\sigma_s}, \hat{\delta_s}, \hat{\lambda_s}, \hat{\kappa_s}) }{
P(D | \hat{\sigma_s}, \hat{\delta_s}, \hat{\lambda_s}, \hat{\kappa_s}) +
P(D | \hat{\sigma_d}, \hat{\delta_d}, \hat{\lambda_d} ) +
P(D | \emptyset ) } \\
P_d(D) &:= \frac{
P(D | \hat{\sigma_d}, \hat{\delta_d}, \hat{\lambda_d} ) }{
P(D | \hat{\sigma_s}, \hat{\delta_s}, \hat{\lambda_s}, \hat{\kappa_s}) +
P(D | \hat{\sigma_d}, \hat{\delta_d}, \hat{\lambda_d} ) +
P(D | \emptyset ) }
\end{align*}
We might now simply select the most probable model, but this results in erratic behavior if the models cannot be distinguished
clearly. Instead, we use all of these models and roll their posterior probabilities into the damage probabilities. This yields the
final universal damage parameters
\begin{align*}
&\sigma_s = P_s(D) \hat{\sigma_s}, \quad
\delta_s = P_s(D) \hat{\delta_s}, \quad
\lambda_s = \hat{\lambda_s}, \quad
\kappa_s = \hat{\kappa_s}, \\
&\sigma_d = P_s(D) \hat{\sigma_d}, \quad
\delta_d = P_s(D) \hat{\delta_d}, \quad
\lambda_d = \hat{\lambda_d}
\end{align*}
\section{Genotype Calling with Simple Errors}
The following heavily borrows notation from the \texttt{BSNP}\cite{bsnp} paper. Consider a
genomic position $j$. Let $G_j$ be the true (unknown) genotype. For
convenience of notation, we write $G_j=\{1,0,0,0\}$ for \texttt{AA},
$G_j=\{\frac{1}{2},\frac{1}{2},0,0\}$ for \texttt{AC}, and so on; we'll
often drop the $j$ subscript. Let $X=(X_1, X_2, \ldots, X_n) \in
\{A,C,G,T\}^n$ be the base calls, $q=(q_1, q_2, \ldots, q_n)$ their
effective\footnote{We roll base quality, base alignment quality, and map
quality into one, see Appendix \ref{app_qualities}} quality scores, $Q=(Q_1, Q_2,
\ldots, Q_n)$ the corresponding error probabilities, $H=(H_1, H_2,
\ldots, H_n) \in \{A,C,G,T\}^n$ the (unobserved) haploid bases
sequenced in each read. The model is that the $H$ are obtained by
sampling from the $G$, possibly modified by chemical damage, then the
$X$ are obtained from the $H$ by application of sequencing error. In
general:
\begin{align*}
L(G) := P(X|G,Q) &= \prod_{i=1}^n P(X_i|G,Q_i,X_1,\ldots,X_{i-1}) \\
&= \prod_{i=1}^n \sum_{H_i} P(X_i|Q_i,H_i,X_1,\ldots,X_{i-1}) P(H_i|G) \\
&= \prod_{i=1}^n \sum_{H_i} P(X_i|Q_i,H_i,X_1,\ldots,X_{i-1}) (H_i \cdot D_i \cdot G_i)
\end{align*}
where $D_i$ is the damage matrix, which depends on the read,
specifically the position within the read. How to model damage is out
of scope here, but would probably follow \cite{mapdamage}.
For maximum likelihood fitting of parameters, we set $L_j = \sum_G L_j(G) P(G)$, for Bayesian \emph{maximum a posteriori} calling,
we set $P(G_j|X_j,Q_j) = L_j(G_j) P(G_j) / L_j$. The prior can be a complicated model, in which case the ML fit serves to
derive statistical parameters, but the simplest possible prior models only heterozygosity:
\begin{align*}
P(G=\{1,0,0,0\}) = P(G=\{0,1,0,0\}) = \cdots &= 1 - \frac{\pi}{4} \\
P(G=\{\frac{1}{2},\frac{1}{2},0,0\}) = P(G=\{\frac{1}{2},0,\frac{1}{2},0\}) = \cdots &= \frac{\pi}{6}
\end{align*}
The minimal error model has no dependency on other bases and directly applies the error rate from the quality score as $Q_i =
10^{-q_i/10}$, following \texttt{BSNP}\footnote{We assume a low quality base is random, as opposed to a random error. Both
are equivalent up to scaling of the error probability, see Appendix \ref{app_errprob}.}, we set:
\begin{align*}
P(X_i|H_i,Q_i,X_1,\ldots,X_{i-1}) := P(X_i|H_i,Q_i) = (1-Q_i) X_i H_i + \frac{Q_i}{4}
\end{align*}
A simple enhancement would be an error matrix modelling typical Illumina errors, another option would be three actual quality
scores from a suitable base caller.
\section{Genotype Calling w/ Dependent Errors}
Both the papers regarding \texttt{maq} and \texttt{BSNP} introduce dependency between errors by a ``dependency parameter'' $\theta$,
which could vary between 0 (totally independent) and 1 (totally dependent), and then raising quality scores to a power involving
$\theta$ and $k_i$, a counter of how many errors have happened so far:
\begin{align*}
Q_i := 10^{-0.1 q'} \qquad \mbox{and} \qquad q'_i := \theta^{k_i} q_i
\end{align*}
This is easy if we pretend that there is only one kind of error or that
we know which one happened (which is how \cite{samtools} gets away with
a very simple presentation). In the case of ancient DNA, we
do not necessarily know which error happened, since we cannot know what was the true base sequened. We have to generalize to
multiple kinds of errors, and count them fractionally. The above equation could be generalized by plugging in matrix exponentials,
but they are both expensive and unlikely to work in a predictable fashion. Instead we define a more general base likelihood:
\begin{align*}
P(X_i|H_i,Q_i) = \left\{ \begin{array}{ccc}
w_{X_i,H_i} 10^{-0.1 q_i \theta^{k_{i,X_i,H_i}}} & \mbox{if} & X_i \neq H_i \\
1 - \sum_{Y \neq X_i} P(Y|H_i,Q_i) & \mbox{if} & X_i = H_i
\end{array} \right.
\end{align*}
The $w_{X,H}$ allow for a substitution matrix or could all be set uniformly to
one quarter.\footnote{If we ever get four quality scores, this might produce
negative probabilities and may need to be modified. Such quality scores are
nowhere in sight, though.} For the $k_i$, we have to count errors
fractionally, so we set
\begin{align*}
k_{1,X_i,H_i} &= 0 \\
k_{i,X_i,H_i} &= k_{i-1,X_i,H_i} + P(H_i | X_i, Q_i, G) \\
&= k_{i-1,X_i,H_i} + \frac{P(X_i | Q_i, H_i) P(H_i | G)}{P(X_i|Q_i, G_i)} \\
&= k_{i-1,X_i,H_i} + \frac{P(X_i | Q_i, H_i) \left( H_i \cdot D_i \cdot G \right) }
{\sum_{H'} P(X_i| Q_i, H') \left( H' \cdot D_i \cdot G \right)}
\end{align*}
\oops{In an earlier version, this said $P(X_i,H_i|Q_i,G)$, which is
nonsense, since $X_i$ already happened. So we must condition on it and
write $P(H_i|X_i,Q_i,G)$.}
\section{Parameter Fitting for Single Sample}
We seek to fit the Allele Frequency Spectrum to a single sample, or in
other words, to estimate heterozygosity and divergence. Let $X,Y$ be
the genotype at some site and $R$ be the reference allele. Now set:
\begin{equation*}
P(x,y|r) = \left\{ \begin{array}{cl}
dh/3 & \mbox{if} \quad x \neq y \wedge (x=r \vee y=r) \\
d(h-1)/3 & \mbox{if} \quad x=y \wedge x \neq r \\
1 - d & \mbox{if} \quad x=y \wedge x=r \\
0 & \mbox{otherwise}
\end{array} \right.
\end{equation*}
Here we assume that all heterozygous mutations happen at the same
uniform rate, and all homozygous mutations at a different uniform rate.
We declare heterozygous genotypes with two alternative alleles to be
impossibru, getting us out of the need to fit a third parameter, which
would have little impact anyway. We get a divergence of $d$ and a
heterozygosity of $dh$. Labelling the genotype likelihoods as $G_{XY}$
and assuming the reference allele is $A$ for convenience, we can compute
the likelihood per site:
\begin{align*}
L &= G_{AA} \left( 1-d \right) + \frac{d(1-h)}{3} \left( G_{CC} + G_{GG} + G_{TT} \right)
+ \frac{dh}{3} \left( G_{AC} + G_{AG} + G_{AT} \right) \\
&= G_{AA} \left( 1-d + d(1-h) \frac{G_{CC} + G_{GG} + G_{TT}}{3 G_{AA}}
+ dh \frac{G_{AC} + G_{AG} + G_{AT}}{3 G_{AA}} \right)
\end{align*}
We now define three quantities of interest, which we sort by magnitude.
The smallest one is accumulated directly, the differences to the other
two are discretized and tabulated. Just six small tables are needed.
\begin{equation*}
R_i := \ln 3 G_{AA}, \quad
D_i := \ln G_{CC} + G_{GG} + G_{TT}, \quad
H_i := \ln G_{AC} + G_{AG} + G_{AT}
\end{equation*}
We further set $\delta = \ln \frac{d}{1-d}$ and $\eta = \ln
\frac{h}{1-h}$, because the log-odds-ratios pose fewer numerical
problems and provide cleaner confidence intervals. Thus we recover the
log-likelihood as follows:
\begin{align*}
L_l &= - \sum_{i} R_i - \sum_{i} \ln \left( \frac{1}{1+e^\delta} +
\frac{e^\delta}{1+e^\delta} \frac{1}{1+e^\eta} e^{D_i} +
\frac{e^\delta}{1+e^\delta} \frac{e^\eta}{1+e^\eta} e^{H_i} \right)
\end{align*}
\idea{Need to do something similar for indels.}
\result{Hand-crafting partial derivatives of this equation didn't seem to
yield anything that simplifies, so I applied Automatic
Differentiation and handed it over to Hager-Zhang.}
\section{Testing Method}
\subsection{Handcrafted Data}
To test for egregious bugs, we can write a couple of SAM or BAM files by
hand. This shouldn't really be called a test; it's ordinary, boring
debugging.
\result{Nothing to see here. Code runs.}
\subsection{Simulated Data}
For all of the simulations, the genome used does not matter at all. For
simplicity, the genome should be small (a megabase should be plenty) and
completely random. We start with a haploid reference genome, then we
apply divergence to get a diploid sample genome. Reads from the sample
genome are simulated along with their correct alignment, the reads are
then modified according to damage and error models.
\beware{We are \textbf{not} going to use the \textbf{human genome}
or some other monstrosity here. It's needlessly complicated and
provides no benefit.}
\beware{We specifically \textbf{do not} simulate reads, then
\textbf{align} them back. Doing so would introduce alignment problems
into the genotype calling, which makes testing harder while providing no
insight at all.}
\beware{It is not actually possible to test the precision of a genotype
caller by counting miscalls. In any such tests, a genotype callers that
doesn't try to call heterozygotes ``wins''. Since that's undesirable,
``simply counting errors'' is a terrible testing method.}
\paragraph{Simulated Modern Data}
Starting from a genome with known divergence and heterozygosity, we
simulate plain reads with some sequencing error and suitable quality
scores, then call genotypes.
Parameters (divergence, heterozygosity) should be fitted and compared to
their true values, particularly at low (roughly one or twofold)
coverage.
\beware{There is no point in simulating fancy sequencing error.
Here, we assume the simple model is correct and show that maximum
likelihood estimation of parameters works in this setting.}
\result{We accidentally skipped this part.}
\paragraph{Simulated Ancient Data}
This is the same idea, this time including damage with known parameters.
As before, genotype calls can be compared and parameters fitted. The
main goal is to get the estimates for heterozygosity, especially
heterozygous transitions, and damage, which interact, right.
Damage could be simulated either by chosing a position dependent damage
rate and damaging bases independently, or by chosing overhang lengths
according to a distribution, then damaging bases independently at
different constant rates for overhang vs. stem. The genotype caller
uses the former model, but the latter is more correct. We should test
both, both are expected to work reasonably well.
\beware{We will not use an empirical distribution for the damage
rates. The empirical distribution is no closer to reality than one
based on a formula, so there is nothing to be learned from it. We could
use an empirical distribution of overhang lengths, if that could be
obtained.}
\result{Since damage should not correlate with genotypes, estimating
damage in a separate first pass works and is a lot cheaper,
both conceptually and operationally, than co-estimating damage with
heterozygosity.}
\result{Homa simulated two diploid genomes with divergence and aligned,
damaged reads from these. I estimated damage from the reads, and the
estimate is nearly spot on. On a dataset with 20\% divergence, the
estimated divergence and heterozygosity are unaffected by damage,
probably because damage is a minor effect compared to divergence here.
On a 0.1\% divergent dataset, naive genotype calling overestimates
heterozygosity by a factor of 15, but gets homozygous divergence right
(estimated $D=1.56\%, H/D=0.93$, but should be $D=0.3\%, H/D=0.67$.
Running with the estimated damage parameters gives a practically perfect
result.}
\subsection{Real Sequencing Data}
\paragraph{Clean, high-coverage, modern, haploid data}
We need actual sequencing data from a haploid region at sensible
coverage. The goal is to test the two available error models in a
setting without confounding factors, especially heterozygosity. This
should be used to select the better error model and to fit the $\theta$
parameter by the maximum likelihood method, if the \texttt{Maq} model is
selected. It's a good idea to check whether $\theta$ varies between
sample sor sequencing runs.
The haploid region of choice might be the mitochondrion, which is
haploid, but the data will be somewhat contaminated with nuMT sequences.
Alternatively, a uniquely mappable region of the X or Y chromosomes of a
male specimen would work, here the difficulty is to find that unique
region.
\idea{A couple of variations on the error model need to be investigated:
handle in bases in ascending, descending, random or input order of
quality; treat errors on different strands as dependent or independent;
various values of $\theta$. For $\theta = 1$, all of these should match
the naïve error model.}
\paragraph{Clean, high-coverage, modern data, two mixed haploids}
Just like the previous test, but this time with two haploid samples
mixed in equal parts. Here, we simulate a diploid sample, but we know
from the individual calls which positions are heterozygous. Again, we
test which error model is better and assess the correctness of the
calls.
\beware{Counting miscalls suffers from the usual problem that it's
entirely unclear how to count errors. The plan is therefore to estimate
heterozygosity.}
\paragraph{Clean, high-coverage, diploid modern data}
We test the two error models and select the better one. This must be an
internal goodness-of-fit test, since we cannot know the true genotypes.
In principle, this is the only test strictly necessary to decide on an
error model.
\beware{In case the tests on haploid samples are inconclusive,
goodness-of-fit takes over. We don't mess with the particulars of the
above tests until we like the results.}
\paragraph{Clean, low-coverage, modern data}
Assuming we fixed the error model, assuming we can reliably estimate
difficult parameters like heterozygosity, here we investigate the
behaviour at low coverage. The sample can be a high coverage sample
suitably downsampled. In this case, we have a reasonable expectation
for the estimated parameters.
\paragraph{Ancient data, one mitochondrion}
To fit parameters to real data, using the best known model (at this
point, we should have selected an error model), without being confounded
by heterozygosity. Data should be clean, so we don't have to deal with
contamination.
\idea{FFPE samples could serve in place of ancient DNA, it
might come with less trouble due to bacterial contamination.}
\paragraph{Ancient data, two mixed mitochondria}
To investigate interaction of heterozygosity, deamination, error model
in a setting where true heterozygous genotypes are known. Data should
be clean and ideally from the same run (we don't want to deal with
additional contamination and different error profiles). The assumption
is that we can correctly call either sample on its own.
In principle, if we haven't encountered difficulties so far, this should
simply work. We can learn to which extent the parameter estimates are
confounded with each other in a simple setting.
\paragraph{Real Ancient Data}
Including autosomes, sex chromosomes, mitochondria of
varying coverage. To estimate parameters separately, to see possible
interactions. To demostrate the process makes sense.
\beware{I actually have no idea what to look for here. Any possible
result would have to be taken at face value.}
\section{Further Directions}
At this point, we should have a genotype caller that deals well with
deamination, heterozygosity, and recurrent errors. Thorough testing on
a real ancient diploid genome is not possible, because a validated data
set cannot be obtained. We also need a second calling step which
applies a prior and produces genotypes. After that, further ideas
include:
\begin{itemize}
\item Likelihoods involving contamination.
\item Priors involving covariance matrices.
\end{itemize}
\subsection{Dealing With Contamination}
A natural direction to go in is to model contamination\todo{Actual
equations}. To a first
approximation, a contaminated sample has four haplotypes, two of which
are endogenous and occur at the same high frequency, and two are from
the contaminant and occur at lower frequency\footnote{It's conceivable
that we could get away with modelling a haploid contaminant. However,
then allele frequencies are modelled incorrectly at contaminated
heterozygous sites, and the gains of this simplification are modest
anyway: instead of 100 distinct genotypes we'd have 40. I think the
full model is worthwhile here.}.
In a model with independent errors, our genotype likelihoods will now
depend on the contamination rate, which is a variable. However, the
dependency is linear\todo{Confirm}, so we can store each genotype
likelihood as two coefficients. This makes for 200 coefficients where
before we had 10 values. Bad, but no deal breaker.
If we have dependent errors, things get complicated. We can recover
simple formulas by assuming low contamination and then ignoring it where
we count repeated errors. This amounts to assuming that contamination
is low enough so that we never make the same error twice when sequencing
the contaminant. Even if that's not strictly true, the effect should be
minor.
Contamination is a property of a read, not really of a base. If we
tried to infer actual haplotypes, this could easily be taken into
account, but that seems too complicated to contemplate. Instead, when
considering a base in genotype calling, we could assign a probability of
coming from a contaminant to it, which is derived from length
distribution and deamination model. Strictly speaking, deamination
depends on genotype, so genotypes depend on other genotypes, which is
horrible. We recover a simple model by assuming everything a read
crosses matches the reference, except the base under consideration
\todo{Calculate and confirm}.
Fitting of a simple model on a single sample could recover deamination
parameters, length distributions and contamination rate. We're really
only interested in the contamination rate, which is to be treated as
preliminary. A later stage could try and fit both the sample and the
contaminant into a phylogeny, thereby learning a more precise
contamination rate\todo{Try and confirm}.
\subsection{Covariance-Matrix as Prior}
When co-calling individuals from multiple populations, the correct prior
for the genotypes would be based on a covariance matrix\footnote{We
restrict to site-by-site analysis and a couple more simplifying
assumptions.}. Estimating that matrix allows Treemix, Patterson's~D,
Pr\:ufer's Divergence, and PCA; possibly more. We also get a clean
method for imputation for free.
% And the replacement for TreeMix is "Blandskog" (== mixed forest)
Conceptually, it's easy: the covariance matrix serves as prior for the
allele frequencies in multiple populations, the allele frequency serves
as prior for the genotype. The goal is to maximize the likelihood with
respect to the covariance matrix. The likelihood function looks more or
less like this:
\begin{equation*}
L = \frac{1}{\sqrt{2^k \pi^k |\Sigma|}} \int_{0^k}^{1^k} f^T G f
\exp \left( -\frac{1}{2}(f-\mu)^T \Sigma^{-1} (f-\mu) \right) df
\end{equation*}
This appears to require integration over the space of possible
combinations of allele frequencies, which sounds impractical. If this
has to be done numerically (Monte-Carlo integration seems at least
feasible), it will be slow and horrible. It's not obvious how to do it
symbolically (it would only make sense if the genotypes can be
separated, otherwise marginalization over all \emph{combinations} of
genotypes remains), but even then, a few problems remain:
\begin{itemize}
\item The formula lacks a normalization factor arising from the
integration bounds.
\item PopGen models based on allele frequencies fail catastrophically if
the frequencies approach 1 or 0. It seems disingenious to integrate out
to these boundaries.
\item How to deal with $\mu$ is not obvious.
\end{itemize}
Instead, we sidestep these problems by taking inspiration from
SeqEm\cite{seqem}. We treat the allele frequencies $\bar{f}$ as
hidden parameters, marginalize over the genotypes(!), and maximize the
joint probability $P(\mathcal{D}, \bar{f} | \Sigma, \mu)$ with
respect to the variance-covariance matrix $\Sigma$ and the mean
distances $\mu$ using an EM algorithm\footnote{Effectively, we estimate
the allele frequency at every position for every sample. Literally,
this is of course impossible, but in aggregate makes sense for
populations.}.
The expectation step, which is finding estimates for the allele
frequencies, is much easier, as it requires only summation over the
possible genotypes of \emph{one} individual at a time and optimizes
\emph{one} allele frequency at a time, holding the others
constant\todo{Is this correct? At any rate, even if not, multivariate
maximization should do it.}. Maximization is finding a covariance
matrix from allele frequencies, and this is again easy. The obvious
moment based estimator should work.
Maximization (actually done by the moment based estimator for sample
mean and sample covariances):
\begin{equation*}
\hat{\mu}, \hat{\Sigma} := \argmax_{\mu,\Sigma}
P(\bar{f}_1,\ldots,\bar{f}_k| \mu,\Sigma)
\end{equation*}
Expectation (for each site):
\begin{align*}
\hat{f_1} \lell \hat{f_k} &= \argmax_{f_1,\ldots,f_k}
P(\mathcal{D} | f_1,\ldots,f_k) P(f_1,\ldots,f_k | \mu, \Sigma) \\
&= \argmax_{f_1,\ldots,f_k}
\prod_{i=1}^k \sum_{G_i} P(\mathcal{D} | G_i) P( G_i | f_i )
e^{ -\frac{1}{2} ((f_1,\ldots,f_k) - \mu)^T \Sigma^{-1} ((f_1,\ldots,f_k) - \mu) }
\end{align*}
Here, $P(\mathcal{D}|G_i)$ is a genotype likelihood, acting as a
constant, and $P( G_i | f_i)$ follows by assuming
Hardy-Weinberg-Equilibrium\footnote{Even if HWE is violated in practice,
it will be absorbed by the estimation of more drift along the affected
lineage.}:
\begin{equation*}
P(G=RR) = (1-f)^2 \quad P(G=RA) = 2f(1-f) \quad P(AA) = f^2
\end{equation*}
This \emph{looks} as if it can be optimized analytically. The summation
over the genotypes requires no nested sums, so it \emph{looks} tractable.
\idea{The practical implementation should probably not store the
aforementioned 600GB of likelihoods, but only 6GB or thereabout of
allele frequency data per sample. The likelihoods can be generated from
the smaller BAM files on the fly. That probably means 'pileup' needs to
get a lot faster. On second thought, maybe the frequencies don't need
to be stored at all.}
\idea{Dealing with contamination becomes obvious in this framework. For
a contaminated sample, every position has two allele frequencies (or two
sets of three allele frequencies). We need to assign a contaminant
probability to each read, which should derive from the alignment score
given the currently estimated allele frequencies.}
\subsection{Principal Component Analysis}
\beware{I personally think PCA is not an analysis, it barely(!) serves
as visualisation method. Nonetheless, it is frequently requested.}
``Modern'' PCA starts with a matrix where rows correspond to $m$
individuals and columns to $n$ markers, such as allele frequencies.
Means are subtracted from each column, then each column is divided by
its empirical standard deviation $\sqrt{f_j (1 - f_j)}$. For
multiallelic sites, one frequency is used per variant.
Since $m < n$, matrix $X = \frac{1}{n} M M^T$ is small and
\textsc{Lapack} can trivially perform PCA on it. It just so happens
that $X$ is the same covariance matrix we estimated above. (Before
David and Nick got involved, everybody PCA'd $\frac{1}{m} M^T M$
instead, which is much bigger. I don't know why.)
\appendix
\section{Random Base vs. Random Error}
\label{app_errprob}
In \texttt{BSNP}, likelihood of a base is calculated in a way that would be appropriate if the quality score described the
probability of any of the three possible errors. This is problematic, since a low quality score should imply that a base is
completely random. However:
\begin{align*}
L(X|H,Q) &= (1-Q)\delta_{H,X} + \frac{1}{3} Q (1-\delta_{H,X}) \\
&= \delta_{H,X} - Q\delta_{H,X} + \frac{1}{3}Q - \frac{1}{3}Q\delta_{H,X} \\
&= (1-\frac{4}{3}Q)\delta_{H,X} + \frac{1}{3}Q \\
&= (1-P) \delta_{H,X} + \frac{1}{4}P \qquad \mbox{with} \qquad P=\frac{4}{3}Q
\end{align*}
So if we replace $Q$ with $\frac{4}{3}Q$, we obtain the formula where the quality $P$ decribes the probability that an observation
is random. Therefore, the two approaches are equivalent up to scaling of the error probability. Both should be tested for any
given base caller, but no special coding is needed.
\section{Effective Quality}
\label{app_qualities}
We have two very different quality measures, base quality, which applies to bases, and map quality, which applies to reads.
\texttt{BSNP} tries to treat them separately, but since it treats different genomic location as independent, they become the same:
\begin{align*}
P(X|H,Q,Z=0) &= \frac{1}{4} \\
P(X|H,Q,Z=1) &= (1-Q)\delta_{X,H} + \frac{1}{4}Q \\
P(X|H,Q,M) &= P(Z=0|M) P(X|H,Q,Z=0) + P(Z=1|M) P(X|H,Q,Z=1) \\
&= \frac{M}{4} + (1-M)\left((1-Q)\delta_{X,H} + \frac{1}{4}Q \right) \\
&= (1-M)(1-Q)\delta_{X,H} + (1-M)\frac{Q}{4} + \frac{M}{4} \\
&= (1-M-Q+MQ) \delta_{X,H} + \frac{1}{4}(Q-MQ+M) \\
&= (1-Q_e) \delta_{X,H} + \frac{Q_e}{4} \qquad \mbox{with} \qquad Q_e = Q+M-MQ
\end{align*}
So we can handle map quality by defining an effective quality such that it describes at least of the two possible errors (sequencing
or mapping) happening, then computing everything with base qualities only. We can roll in base alignment quality (BAQ), which
doesn't even have a solid definition, too, and express it in quality scores: $q_{\operatorname{eff}} = \operatorname{softmin} \left[
q_{\operatorname{base}}, q_{\operatorname{map}}, q_{\operatorname{baq}} \right]$.
This treatment is not correct in that we try to model dependency between errors, which makes sense for base quality. Mapping
errors, however, are complex and probably even more dependent that sequencing errors. Considering that mapping quality is a crude
approximation anyway, this is not a major concern. In principle, PCR
error that happens before sequencing should also be modelled. However,
PCR error behaves similarly to mapping error, and mapping error is
always of at least the same magnitude. Therefore, we simply ignore PCR
error.
\listoftodos
\begin{thebibliography}{9}
\bibitem{bsnp}
Ilan Gronau et. al.,
\emph{Bayesian inference of ancient human demography from individual genome sequences}.
Nature Genetics 43, 1031---1034 (2011).
\bibitem{samtools}
Heng Li,
\emph{Mathematical Notes on SAMtools Algorithms}.
https://www.broadinstitute.org/gatk/media/docs/Samtools.pdf (2010).
\bibitem{soapsnp}
http://soap.genomics.org.cn/soapsnp.html
\bibitem{maq}
Heng Li, Jue Ruan, and Richard Durbin,
\emph{Mapping short DNA sequencing reads and calling variants using mapping quality scores}.
Genome Research 18, 1851--1858 (2008).
\bibitem{mapdamage}
Aurelien Ginolhac et. al.,
\emph{mapDamage: testing for damage patterns in ancient DNA sequences}.
Bioinformatics 27 (15), 2153--2155 (2011).
\bibitem{seqem}
E. R. Martin et. al.,
\emph{SeqEM: an adaptive genotype-calling approach for next-generation
sequencing studies}.
Bioinformatics 26 (22), 2803-2810 (2010).
\end{thebibliography}
\end{document}