Genetics and population analysis Haplotype matching in large cohorts using the Li and Stephens model

Motivation: The Li and Stephens model, which approximates the coalescent describing the pattern of variation in a population, underpins a range of key tools and results in genetics. Although highly efﬁcient compared to the coalescent, standard implementations of this model still cannot deal with the very large reference cohorts that are starting to become available, and practical implementations use heuristics to achieve reasonable runtimes. Results: Here I describe a new, exact algorithm (‘fastLS’) that implements the Li and Stephens model and achieves runtimes independent of the size of the reference cohort. Key to achieving this runtime is the use of the Burrows-Wheeler transform, allowing the algorithm to efﬁciently identify partial haplotype matches across a cohort. I show that the proposed data structure is very similar to, and generalizes, Durbin’s positional Burrows-Wheeler transform.


Introduction
The genetic variation in a population of interbreeding individuals is highly structured. Kingman (1982) introduced the canonical model that describes this structure mathematically, known as Kingman's coalescent, later extended by Hudson (1983) and Griffiths and Marjoram (1997) to include recombination. Although mathematically elegant, it is challenging to use these models directly for statistical inference. Li and Stephens (2003) introduced a model (LS) that is both a good approximation to the coalescent with recombination, and computationally tractable. As a result, LS now underpins a large range of key tools and scientific findings (Beaumont, 2010;Howie et al., 2009;The International HapMap Consortium, 2005;The Wellcome Trust Case Control Consortium, 2007). Depending on whether the input sequence is haploid or diploid, LS in its straightforward implementation as a hidden Markov model (HMM) runs in linear or quadratic time in the number of reference haplotypes. While this is orders of magnitude more efficient than algorithms based directly on Kingman's coalescent or the ARG, the recent availability of affordable DNA sequencing technology has provided access to very large reference sets, on which even the LS model is intractable in its standard implementation, so that current implementations of LS use heuristics to cope with datasets encountered in practice (Howie et al., 2009).
A very different algorithm that is making an impact in genomics was introduced by Burrows and Wheeler (1994). Known as the Burrows-Wheeler transform (BWT), it permutes an arbitrary text in such a way that the original text can be recovered, while at the same time improving the compressibility of the transformed text by increasing simple repetitions. In addition, the transformed text, even in compressed form, serves as an index that allows rapid searching in the original text. In genomics this idea has so far been used mainly for fast alignment of short reads against a large and relatively repetitive reference genome (Langmead et al., 2009;Li and Durbin, 2009). More recently, Durbin (2014) introduced a variant of the BWT, termed the Positional Burrows-Wheeler Transform (PBWT), that exploits the additional structure that exists in a set of haplotypes in a population sample. These data, which are usually encoded as a series of 0 and 1 s representing the absence or presence in a sample of particular genetic variants along a reference sequence, have a natural representation as a matrix, where rows represent samples and columns represent the particular positions in a reference. Local matches between samples are only relevant at matching positions, and exploiting this restriction leads to improvements over a standard application of the BWT. The resulting data structure again allows for fast haplotype searches against a database, and achieves very high compression ratios.

Approach
There are two main results in this paper. First, I establish a formal connection between the standard and positional BWT, showing how the PBWT as introduced in Durbin (2014) is a special case of the BWT. This connection also shows how the PBWT can be slightly generalized to cope with the multiallelic case. Besides providing an additional perspective on the positional BWT algorithms, which helps to better understand them, it also provides a mechanical way to 'lift' existing algorithms operating on the BWT data structure to their positional equivalent, allowing the large literature on BWT algorithms to be applied to the current data structure. I show how this works by deriving the haplotype search algorithm from the equivalent BWT algorithm.
The second contribution consist of algorithms that implement the LS model on top of the BWT. More precisely, I present algorithms that compute maximum-likelihood ('Viterbi') paths through the LS hidden Markov model, providing a parsimonious description of a given sequence as an imperfect mosaic of reference haplotypes. The ability to efficiently identify matches in the database of reference haplotypes result in considerable improvements in runtime over the standard implementation, reducing the linear and quadratic asymptotic runtime to empirical constant time, independent of the number of reference haplotypes. More precisely, for H samples of n loci each, the standard implementation runs in OðHnÞ time for a haploid input sequence, and OðH 2 nÞ for a diploid input sequence, while the proposed algorithms run in empirical O(n) time in both cases. This allows the Li and Stephens model to be used on very large reference panels, without recourse to approximations.

Materials and methods
3.1 Haplotype matching using the BWT Let x 0 ; . . . ; x HÀ1 be H haplotype sequences, each consisting of n symbols from the alphabet A representing the possible allelic states at a locus; for simplicity I will often use A ¼ f0; 1g in this paper. A straightforward way of identifying haplotype matches would be to use the BWT on the concatenation x 0 x 1 Á Á Á x HÀ1 of haplotype sequences. It turns out that a more efficient algorithm is obtained, in terms of time and memory use, by embedding this sequence of Hn characters into a sequence of 2Hn characters taken from a much larger alphabet. The increase in sequence length and alphabet size is offset by the additional structure in the BWT that results from the chosen embedding. This in turn translates into better compression and a streamlined search algorithm.
I will write x½j for the jth symbol in the sequence x, and x½j; kÞ for the subsequence starting at position j and ending at k -1. I will also use ½i; jÞ to denote the half-open interval fi; i þ 1; . . . ; j À 1g, and if M ij is a matrix, M k ½i; jÞ is the subsequence M k;i ; M k;iþ1 ; . . . ; M k;jÀ1 of the kth row of the matrix. Throughout this paper, all indices start at 0.
Let p 0 ; . . . ; p nÀ1 be n additional symbols in the alphabet, ordered such that p 0 < Á Á Á < p nÀ1 < 0 < 1. Introduce a new sequence X of length 2Hn by inserting a symbol p j after each symbol x i ½j and concatenating the resulting sequences into a single sequence of the form (To impose a particular initial ordering I will later on replace the last symbol p nÀ1 by H symbols p 0 nÀ1 < Á Á Á < p HÀ1 nÀ1 , but to avoid cluttering the notation I ignore this detail for now.) Consider all cyclic shifts Let M be the matrix obtained by writing X k on the kth row of a square matrix, and sorting the resulting rows lexicographically. Let p be the permutation that sorts the rows, so that X pð0Þ < X pð1Þ < Á Á Á < X pð2HnÀ1Þ , and M ij ¼ X pðiÞ ½j. The Burrows-Wheeler transform of X is the last column of this matrix: BWTðXÞ½i ¼ X pðiÞ ½2Hn À 1. Note that this is almost the traditional BWT of the sequence X, except that there is no special 'end' character. This character is used to identify the start of the sequence; here, the special structure of X is sufficient to navigate BWT(X). Now consider how the matrix M may be constructed. The position symbols p i determine the coarse structure of M, which is independent of the data x i apart from the haplotype frequencies f 0 i and f 1 i (see Fig. 1). The fine-scale structure of M within each 'block' of H rows is determined by the data. More precisely, rows in the block starting at index iH are those cyclic shifts of X that start with symbol p i and end with x k ½i for some k 2 ½0; HÞ, such that these rows are ordered lexicographically within the block. Let j7 !a i j denote the permutation of ½0; HÞ that describes this order within block i, so that row iH þ j ends with symbol x a i j ½i. Determining M therefore boils down to determining the n permutations a i j for i 2 ½0; nÞ, since these determine the top half of M, and those in turn determine the remaining rows (see Fig. 1 and the explanation).
The permutations a i j are determined recursively, working from i ¼ n À 1 backwards. Because we imposed the special ordering p 0 nÀ1 < Á Á Á < p HÀ1 nÀ1 on the final position symbols, the permutation for block n -1 is given by the identity permutation j7 !a nÀ1 j ¼ j. Now suppose the permutation a i j for block i has been determined. The sequences in block i -1 are formed from those in block i by moving two characters from the end to the front. The first character in any sequence of this new block is p iÀ1 , which does not influence the ordering within the block. The second character is an allele marker x a i j ½i. To sort the sequences in block i -1 in lexicographic Algorithm 1 Calculating BWT(X) Input: sequences x 0 ; . . . ; x HÀ1 , each of length n; alphabet A Output: Block permutations j7 !a i j ; i ¼ 0; . . . ; n À 1. 1: i n; a nÀ1 j j for j 2 ½0; HÞ 2: While i > 0: 3: order, it is therefore sufficient to list those sequences that start a 0 symbol first, followed by those starting with a 1 symbol (followed by other symbols if the locus is multiallelic), and otherwise leave the original order undisturbed. Doing this results in Algorithm 1.
To show that the proposed construction is equivalent to the positional Burrows-Wheeler transform, Algorithm 1 is given both for general alphabets A and specialized for the case A ¼ f0; 1g, since that in that case the inner loop is precisely Algorithm 1 in Durbin (2014) (except that the proposed algorithm runs back-to-front, as is usual for BWT algorithms). As in the PBWT algorithm, the permutations a i j play the role of the suffix array in the ordinary BWT algorithm. Note that the output includes a permutation a À1 j , which encodes how the very first characters x j ½0 influence the permutation of the cyclic shifts X k ; this permutation is used in Algorithm 5. Following Durbin (2014) I now define the PBWT of x 0 ; . . . ; x HÀ1 as the first half of BWT(X), which is availably implicitly as Figure 1 shows that the second half of BWT(X) is determined by the allele frequencies f c i ; i 2 ½0; nÞ, which can be computed easily from the relevant block in the first half of BWT(X), so that the PBWT of x 0 ; . . . ; x HÀ1 is in fact equivalent to BWT(X).

Substring searching
Algorithm 1 calculates BWT(X) in linear time by exploiting the special structure of X, and is not a specialization of an existing, general algorithm to calculate the BWT. By contrast, Algorithm 2, which performs a substring search, can be derived directly from its analogous algorithm for a general BWT.
To describe the algorithm, let M be the sorted matrix of cyclic shifts of an arbitrary sequence X of length n, so that BWTðXÞ½i ¼ M i ½n À 1, and let R a ðiÞ (the 'a-rank' for row i) be the number of times that a appears in BWTðXÞ½0; iÞ. This function can be calculated efficiently from BWT(X), particular if the data is stored in compressed form. Finally, let C(a) (the cumulative symbol frequency) be the number of symbols in X that are less than a. This notation makes it possible to write down Algorithm 2, for substring searching. (The symbol . is used throughout to mark comments and invariants in the algorithms.) To understand the algorithm, consider all rows of M that end with a symbol a. If these rows are cyclically shifted rightward, so that the last symbol becomes the first and all others are moved one position to the right, all rows will now start with a, and the relative order in which they appear in M (which they must as M contains all cyclic shifts of X) is the same as before the shift since they were ordered lexicographically to start with. Suppose that M k is a row that ends with a, and that after right-shifting it ends up as row M k 0 ; then the above observation means that the rank R a ðkÞ of the symbol a in M k in the last column of M, is the same as the rank in the first column of M of the symbol a in M k 0 . Because M is sorted lexicographically, the rows that start with a form a contiguous block in M, so that the first-column rank of the symbol a in row M k 0 is k 0 À CðaÞ, so that R a ðkÞ ¼ k 0 À CðaÞ or LFðk; aÞ : The function k7 !LFðk; aÞ, mapping row k to the row corresponding to its right-shifted counterpart k 0 , is called the last-to-first mapping because it maps the last (rightmost) symbol of M k to the corresponding symbol in the first (leftmost) position of M k 0 . It is repeatedly used to identify the interval of rows corresponding to sequences that match one additional character of w.
Note that the mapping is well-defined whether or not M k ½n À 1 ¼ a. This makes it possible to think of k as representing a possible location between two entries (k and k -1) in M where a sequence (or sequence prefix) x not necessarily represented in M would be inserted; this is the view taken in the search algorithm. Alternatively, when k is thought of as a particular row in M, that row's initial character a can be obtained from the CðÁÞ function, and since the mapping (2) is invertible when restricted to the set of rows k ending in a, this makes the mapping k7 !LFðk; M k ½n À 1Þ invertible for all k. The existence of this inverse mapping also follows is the rightmost column of M, while the positional BWT of the sequences x0; . . . ; x HÀ1 is the upper half of the same column (see text). The column indices are determined by f a i , the allele frequency of symbol a at locus i, and F a i :¼ P i j¼0 f a i , the cumulative frequency of symbol a across loci 0; . . . ; i. Note that ordering of rows ðn À 1ÞH to Hn À 1 is determined by the special position symbols p 0 nÀ1 < Á Á Á < p HÀ1 nÀ1 , but to avoid cluttering the notation these are all written as p nÀ1 Algorithm 2 General subsequence search Input: Sequence w½0; jÞ; BWTðXÞ of sequence X½0; nÞ Output: Indices s, e such that M k ½0; jÞ ¼ w for k 2 ½s; eÞ 1: s 0; e n; i j 2: While s < e and i > 0: " w½i; jÞ matches M k ½0; j À iÞ for k 2 ½s; eÞ 3: i i À 1 4: s Cðw½iÞ þ R w½i ðsÞ 5: e Cðw½iÞ þ R w½i ðeÞ directly from the observation that it corresponds to rotating the sequence one position leftward; it could be called the first-to-last mapping, k7 !FLðkÞ, and is used in Algorithm 5.
To derive the corresponding algorithm for matching a sequence in the PBWT data structure, it is enough to track the bounding variables for two steps through the standard BWT algorithm acting on the 'lifted' sequence X, matching a haplotype character and a position character. The first step identifies the new range depending on the haplotype character to be matched, and points these variables to the second half of the matrix. The next step moves the bounding variable back into the first half by moving a position character in front. Because of the regular form of BWT(X) (see Fig. 1), these two steps can be followed algebraically and combined into a single update step. The derivation, which is straightforward but requires additional notation, is presented in the Appendix. The resulting combined update step is given by a modified last-to-first mapping function, which now additionally depends on the current position i: LFðk; a; iÞ : or for an arbitrary alphabet, LFðk; a; iÞ ¼ r a i ðkÞ þ P c < a f c i . Here r a i ðkÞ is the positional analogue of R a ðiÞ, and counts how often a appears in the first k rows of the ith block of PBWTðx 0 ; . . . ; x HÀ1 Þ, or equivalently, in BWTðXÞ½Hi; Hi þ kÞ ¼ x a i 0 ½i; . . . ; x a i kÀ1 ½i, and f a i is the (haplotype) frequency of a at position i. This leads to Algorithm 3.

Haploid Li and Stephens
The Li and Stephens (2003) model approximates the coalescent model describing the relationship between DNA sequences in a population, by generating a new sequence as a mosaic of imperfect copies of existing sequences The popularity of the model stems from the fact that it is both a good approximation to the full coalescent model with recombination, as well as fast to compute in its natural implementation as a hidden Markov model, running in OðHnÞ time for H sequences of length n. However, for very large population samples this is still too slow in practice.
Here I describe an algorithm to compute the maximum likelihood path through the LS hidden Markov model (HMM) in empirical O(n) time. The approach is not to consider single sequences to copy from, but groups of sequences that share a common subsequence. Like the Viterbi algorithm for HMMs, the proposed algorithm traverses the sequence to be explained, but rather than using a dynamic programming approach, it uses a branch-and-bound approach considering (groups of) potential path prefixes to a maximum likelihood path. Where at each iteration the Viterbi algorithm must consider all possible sequences that a potential path prefix Algorithm 3 PBWT subsequence search Input: Sequence w½0; jÞ, PBWT of x 0 ; . . . ; x HÀ1 Output: Indices s, e such that x a 0 k ½0; jÞ ¼ w for k 2 ½s; eÞ 1: s 0; e H; i j 2: While s < e and i > 0: " w½i; jÞ matches x a i k ½i; jÞ for k 2 ½s; eÞ 3: i i À 1 4: s LFðs; w½i; iÞ " see equation (3)  a; gm idx FLðgm idx; iÞ 12: If a 6 ¼ x½i: gm gm À l 13: i i þ 1 14: If gm ¼ t score: 15: gm idx t idx; gm gm À q; path:appendðði; a iÀ1 gm idx ÞÞ 16: Return path could end with, the proposed algorithm in principle considers all extensions of the current potential path prefixes (the 'branch' part), but ignores prefixes that cannot be part of an optimal path (the 'bound' part). For instance, if a prefix can be extended with a matching nucleotide, a recombination does not have to be considered, since the recombination can be postponed at no cost. Below I will show this more formally. This formal approach is perhaps not necessary (or even helpful) for the haploid case, but becomes useful when I introduce the diploid Li and Stephens algorithm.
First some definitions. A placed character is a character c at a sequence position i; it is equivalent to a pair cp i where p i is the position symbol introduced before. Two placed characters are contiguous if they occupy neighbouring positions; subsequences of placed characters are contiguous if every pair of neighbouring characters is; and two or more subsequences are contiguous if their concatenation is. A path p of m parts through a set of sequences X ¼ fx 0 ; . . . ; x HÀ1 g is a contiguous sequence of m subsequences s 0 ; . . . ; s mÀ1 such that each s i is a subsequence of some x j . I will write a path as where c i is a character placed at position i, and k 0 ; k 1 ; . . . ; k mÀ2 are the recombination breakpoints identified by the symbol R (which is not part of the alphabet), and l is the length of the path. The (sequence) group associated with p is the set GðpÞ of all sequences x 2 X for which the subsequences x½k mÀ2 ; lÞ agree with the suffix c kmÀ2 Á Á Á c lÀ1 that follows the last recombination in p. The extension pc l (of length l þ 1) is the path ðc 0 Á Á Á Rc kmÀ2 Á Á Á c lÀ1 c l Þ, if it exists; since by definition all subsequences that make up a path are subsequences of some x j , existence of an extension implies that its group is nonempty. The extension pR (of length l) is defined as ðc 0 Á Á Á Rc kmÀ2 Á Á Á c lÀ1 RÞ, and always exists; its group is X. Finally, the path prefix p½0; tÞ is the path ðc 0 Á Á Á c tÀ1 Þ including any R symbols for recombinations between positions 0 and t -1; a path prefix never ends with an R symbol.
For a given sequence x and a path p, the Li and Stephens model assigns a joint likelihood to the event that p occurred and gave rise to sequence x. If p has m parts and has k mismatches to x, this likelihood is where p q is the probability of recombining into a particular other sequence, and p l is the probability of a mutation to one of the three other nucleotides. The negative log likelihood takes a particularly simple form, where C is a constant, q ¼ Àlog ðp q =ð1 À np q ÞÞ and l ¼ Àlog ðp l = ð1 À 3p l ÞÞ. This motivates defining the path score as s x ðpÞ ¼ mqþ kl, where m and k are defined as above. I drop the subscript x from s x ðpÞ when this is possible without creating confusion. Suppose we want to calculate a path p that minimizes sðpÞ. This can be done by iteratively constructing path prefixes p 0 , so that at each step one of them is a prefix of a full path p that minimizes sðpÞ. Note that the minimum score achievable by a path p that has p 0 as its prefix depends on the prefix score sðp 0 Þ and the prefix group Gðp 0 Þ, but not on the rest of the prefix. This is because Gðp 0 Þ is the set of sequences the Li and Stephens model could be copying from at the end of p 0 , and the Markov property of the model implies that the minimum score only depends on the sequence being copied from (and the prefix score). This justifies the definition of state of a path (prefix) p 0 to be the pair ðGðp 0 Þ; sðp 0 ÞÞ.
The key observation for the algorithm is that some states (G, s) can be ignored, because any of their extensions give rise to paths and scores that are also achievable via other states. To make this precise I need one more definition. A set S of path prefixes, all of length l, is a full prefix set for x½0; lÞ if for any sequence x 0 whose prefix x 0 ½0; lÞ agrees with x½0; lÞ, there exists a path p that achieves the minimum score (i.e. s x 0 ðpÞ ¼ min p 0 s x 0 ðp 0 Þ) and whose prefix p½0; lÞ is in S. If we can somehow find a way to iteratively construct full prefix sets of increasing length, the problem of finding a minimum-score path is solved, because the required path will be an element of the full prefix set for the full-length sequence x. The following theorem shows how to do this: THEOREM 1. Suppose S is a full prefix set for x½0; lÞ; S 0 a set of prefixes of length l þ 1, and let s min ¼ min p2S sðpÞ and s 0 min ¼ min p2S 0 sðpÞ. Then S 0 is a full prefix set for x½0; l þ 1Þ if the following conditions hold: a. a For all p 2 S and all a 2 f0; 1g so that pa is an extension and sðpaÞ < s 0 min þ q we have pa 2 S 0 ; and b. b If there is no p 2 S so that sðpÞ ¼ s min and px½l is an extension, then S 0 contains a path of the form pRx½l with sðpÞ ¼ s min .
In other words, certain extensions are not required to be in S 0 : extensions pa whose score exceed the minimum plus q can be left out (since a recombination from the minimum-scoring prefix would give a path that is at least as good), and recombinations can be ignored altogether as long as any current lowest-scoring path has a matching extension (since otherwise postponing the recombination would again be at least as good) -and if not, only a single recombination from a lowest-scoring path needs to be considered.
Algorithm 4 implements these ideas. It does not actually construct prefix sets of paths, but sets of states of paths in prefix sets. This is sufficient since the state determines how paths can be extended. By using the PBWT, these states can be represented efficiently, using just the score and a pair of indices into the PBWT that correspond to a set of subsequence matches to sequences in X, similar to how the variables s and e in Algorithm 3 represent the interval ½s; eÞ corresponding to a set of subsequence matches. Another difference with the description above is that the algorithm scans the sequence back-to-front, extending partial matches leftward, so that the invariant refers to the full suffix set, rather than the full prefix set.
The algorithm computes gm ¼ s min , and keeps a running minimum score gm 0 that bounds s 0 min , ignoring states whose new score are not less than gm 0 þ q. At the end of an iteration, states whose score are not lower than the now updated gm 0 plus q are not immediately removed, but are instead ignored in the next iteration. The algorithm implicitly considers both score bounds implied by gm and gm 0 , but in each situation uses only the tighter bound of the two to decide which states to ignore.
It is possible for different paths to result in overlapping or identical states, resulting in duplicate or otherwise redundant entries in the st array. Although redundant entries do not impact the correctness of the algorithm, they can dramatically reduce efficiency. A practical implementation therefore includes a step that occasionally removes redundant states.
The algorithm can be generalized a little by allowing the mutation score l ! 0 to depend on the position. The path score is then defined as sðpÞ ¼ mq þ P i:x½i6 ¼p½i l i . Theorem 1 continues to hold, and so does Algorithm 4, with the obvious changes. The current approach does not lend itself easily to generalize to a positiondependent recombination probability, as the proof of Theorem 1 relies on delaying the recombination without changing the score, which is only possible if q is constant along the sequence.
Note that the algorithm can be simplified when l i ! 2q, because a mismatch can always be circumvented by two recombinations (before and after the offending locus), so that only exact matches need to be considered. In human genetics polymorphisms are sparse, and recombinations can only be localized to within hundreds or thousands of positions. Even when a maximum likelihood path is sought it is natural to marginalize over these positions, and this makes the probability of a recombination between two polymorphic sites at least an order of magnitude higher than the probability of a mutation, so that l ) q. However, in the presence of phasing errors the probability of a mismatch can be much higher than that of a mutation, so that the regime l < 2q is of practical importance.
Algorithm 4 only computes the optimal score, and to obtain an optimal-scoring path p itself a backtracking step is needed (Algorithm 5). Here it is useful that Algorithm 4 works in the backward direction, so that the result of the backtracking is oriented in the natural direction. To track an optimal path along a sequence, the PBWT index corresponding to that sequence can be tracked using the 'first-to-last' mapping, inverting the steps in lines 6 and 12 in Algorithm 4, and the minimum score of the remaining suffix is updated whenever a difference between this sequence and x is found.
Recombinations are followed greedily, as it is always correct to follow a feasible recombination, and it is never clear whether a particular recombination is the last feasible one for a particular sequence. Algorithm 4 collects information about recombinations in the traceback list, and when a recombination and score is identified that forms a feasible suffix to the path so far, it is followed.
The naive implementation of Algorithm 5 is somewhat slower than the haploid Li and Stephens algorithm itself, due to the FL function which takes Oðlog HÞ time in the implementation shown. In practice the PBWT will be stored in compressed form using runlength encoding, which allows a faster implementation of FL.

Diploid Li and Stephens
Where the haploid Li and Stephens algorithm computes a single haplotype path maximizing the probability of a given haploid sequence, the diploid Li and Stephens algorithm aims to find a pair of haplotype paths that maximizes the probability of a sequence of diploid genotypes under the same model. The approach used to derive the haploid algorithm also works in this case, but the details are more involved.
Let x be a sequence of genotypes, encoded as values 0, 1 or 2 at each position representing homozygous ancestral, heterozygous and homozygous derived genotypes. The aim is to compute a pair of paths a, b that minimizes a score. As before this score contains terms for recombinations and mismatches, but the mismatch term now considers genotypes rather than haplotypes. More precisely, the score associated to the pair fa; bg is defined as sða; bÞ ¼ qmðaÞ þ qmðbÞ þ lkða; bÞ, where mðaÞ represents the number of parts of path a, as before, and k ¼ P i ja½i þ b½i À x½ij counts the number of mismatches of the paths a and b to the genotype sequence x.
The approach of the algorithm is similar to the haploid case, again sequentially building full prefix sets for ever longer sequence prefixes until a minimum path pair is found. To describe the
The sequence group associated to an unordered pair of paths fa; bg is defined as Gða; bÞ ¼ ffx; ygjx 2 GðaÞ; y 2 GðbÞg. Similarly, using the same justification as before, the state of an (unordered) path pair fa; bg is defined to be the pair ðGða; bÞ; sða; bÞÞ. A full prefix set S for x½0; lÞ is defined as a set of (unordered) pairs of path prefixes such that for any sequence x 0 that extends x½0; lÞ, there exists a path pair fa; bg that achieves the minimum score s x 0 ða; bÞ ¼ min a 0 ;b 0 s x 0 ða 0 ; b 0 Þ and whose prefix pair fa½0; lÞ; b½0; lÞg is in S. Finally, to formulate the theorem it is handy to introduce the notation S to denote the set of 'haplotype' paths in S, or formally S ¼ fajfa; bg 2 Sg. Algorithm 6 implements these ideas. The core of the algorithm is formed by lines 11 and 14 that consider regular extensions with a pair of characters (a 1 , a 2 ); lines 22-23 and 27 that consider single recombinations; and line 34 that considers simultaneous recombinations in both haplotypes. The remainder of the algorithm is concerned with implementing the conditions of Theorem 2 to ensure that redundant extensions are ignored. The variables gm and gm 0 keep track of the current and next global minimum score s min and s 0 min , while the associative arrays lm½ and lm 0 ½ keep track of s min ðaÞ and s 0 min ðaÞ respectively. The associative array extended½ keeps track which paths a have a partner b 0 that achieves the minimum score s min ðaÞ, and for which both a and b 0 have extensions required in conditions b and c; whether the extension is appropriate is computed by the function consider recomb. Finally, the variable double recomb is used to ensure that at most one double recombination is considered at every iteration.
The traceback algorithm for diploid Li and Stephens is similar to the haploid algorithm. Again, the traceback list contains records describing the recombinations that have been considered. These records now additionally contain a pair s x , e x that represent the range of PBWT indices corresponding to the sequence that does not undergo a recombination. As with the haploid algorithm, the traceback algorithm follows a recombination only if the path scores agree, but now also ensures that the index of the non-recombining path is contained in the range ½s x ; e x Þ. Double recombinations are encoded by setting e x ¼ À1, and for such recombinations only the scores need to agree. A pseudocode implementation is given as Algorithm 7.

Performance
For testing the fastLS algorithms were implemented in Cþþ, with all tables stored in uncompressed form in memory. To validate the implementations and to compare runtimes, standard Viterbi Algorithm 7 Diploid traceback Input: Sequence x½0; nÞ, PBWT of x 0 ; . . . ; x HÀ1 , scores l ! 0; q ! 0, minimum score gm, corresponding indices gm idx1; gm idx2, traceback list traceback. Output: Representation of a minimum-scoring diploid path 1: Function FL(k, i): " "First-to-last" mapping 2: lo 0; hi H; a 0 if k < f 0 i else 1 3: While lo < hi: " LFðj; a; iÞ k 8j < lo and LFðj; a; iÞ > k 8j ! hi 4: mid bðlo þ hiÞ=2c 5: If LFðmid; a; iÞ k: lo mid þ 1 6: Else hi mid 7: Return a; lo À 1 8: i 0; path1 ½ði; a iÀ1 gm idx1 Þ; path2 ½ði; a iÀ1 gm idx2 Þ 9: For ðt locus; t start; t end; t idx; t scoreÞ in reverseðtr 0 backÞ: 10: While i t locus: 11: a1; gm idx1 FLðgm idx1; iÞ 12: a2; gm idx2 FLðgm idx2; iÞ 13: gm gm À lja1 þ a2 À x½ij; i i þ 1 14: If gm ¼ t score: 15: If t end ¼ À1: " Double recombination 16: gm idx1 t start; gm idx2 t idx 17: path1:appendðði; a iÀ1 gm idx1 ÞÞ; path2:appendðði; a iÀ1 gm idx2 ÞÞ 18: gm gm À 2q 19: Else If t start gm idx1 < t end:" Single rec. in path 2 20: gm idx2 t idx; path2:appendðði; a iÀ1 gm idx2 ÞÞ 21: gm gm À q 22: Else if t start gm idx2 < t end: " Single rec. in path 1 23: gm idx1 t idx; path1:appendðði; a iÀ1 gm idx1 ÞÞ 24: gm gm À q 25: Return path1; path2 (a) ( b) Fig. 2. Running time for inferring inheritance patterns under the haploid (dashes) and diploid Li and Stephens model over a simulated reference set of n (horizontal axis) haploid sequences, using the Viterbi (red) and fastLS (green) algorithms, using q=l ¼ 2. Dots represent measurements, curves show quadratic fits. (a) Results for a simulated reference population of n samples. (b) Results for a fixed simulated reference population of 100 000, subsampled to n samples algorithms for the haploid and diploid LS model were also implemented. Traceback was included in the fastLS algorithms, but was excluded from the Viterbi implementations because of memory constraints. Two sets of simulations were performed. For the first, 30 Mb of sequence in populations of size 100 to 10 000 were simulated by scrm (Staab et al., 2015) using the 'standard simulation' model of Li and Durbin (2011) which roughly resembles the demography of the European population. For each population I simulated an additional 50 samples to serve as input sequences. This resulted in a number of segregating sites ranging from 129 945 for the 150sample case, to 436 361 for 10 050 samples. For the second set, I simulated a single population of 100 000 samples under the same model (resulting in 621 156 segregating sites) and sub-sampled reference populations of 100 to 10 000 samples from these (Fig. 2). The run-times of the Viterbi algorithms show the expected linear and quadratic dependence on H. The fastLS algorithms show a weak dependence on H. In the case of the sub-sampled population, which have a fixed number of loci (not all of which segregate in the sample), the dependence on H is weakest, and in fact the diploid algorithm becomes faster for larger populations, probably because longer haplotype matches can be found in larger populations, resulting in more efficient pruning of the prefix sets.