biohazard-0.6.1: 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}}
\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} &:= \arg \max_{\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} &:= \arg \max_{\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{Testing Method}
\subsection{Handcrafted Data}
To test for egregious bugs, we can write a couple of SMA or BAM files by
hand. This shouldn't really be called a test; it's ordinary, boring
debugging.
\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.}
\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 genotype call.
Called genotypes can be compared to the correct genome, but more
importantly, 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.}
\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.}
\idea{Since damage should not correlate with genotypes, estimating
damage in a separate first pass might work and would be lot cheaper,
both conceptually and operationally, than co-estimating damage with
heterozygosity. This is a good time to try it.}
\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 if the \texttt{Maq} model is selected.
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.
\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.
\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, he we investigate the bahviour
at low coverage. The sample can be a high coverage sample suitably
downsampled. In this case, we have an 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 assumtion
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. Estimating
that matrix allows Treemix, Patterson's~D and Pruefer's Divergence.
Conceptually, it's easy: the covariance matrix serves as prior for the
allele frequencies in multiple populations, the allele frequency
(together with a small term for new mutations) serves as prior for the
genotypes\todo{Equation!}. Maximizing the covariance matrix is
straight-forward, but it would require integrating over the space of
possible combinations of allele frequencies, which sounds impractical,
and becomes more impractical the more samples are considered.
Instead, we can estimate the joint probability (genotype(s), allele
frequency) and maximize that, which is much easier\footnote{Effectively,
we estimate the allele frequency at every position for every sample.
Which is impossible, but the aggregate makes sense for populations.}.
Only summation over the possible genotypes is necesssary, which is just
10 per individual, and individuals are independent; allele frequencies
and covariance matrix are co-estimated using something resembling the EM
algorithm. (One idea would be to not store the aforementioned 600GB of
likelihoods, but only 6GB or thereabout of allele frequency data. The
likelihoods can be generated from the smaller BAM files on the fly.)
\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).
\end{thebibliography}
\end{document}