An average-case sublinear forward algorithm for the haploid Li and Stephens model

Background Hidden Markov models of haplotype inheritance such as the Li and Stephens model allow for computationally tractable probability calculations using the forward algorithm as long as the representative reference panel used in the model is sufficiently small. Specifically, the monoploid Li and Stephens model and its variants are linear in reference panel size unless heuristic approximations are used. However, sequencing projects numbering in the thousands to hundreds of thousands of individuals are underway, and others numbering in the millions are anticipated. Results To make the forward algorithm for the haploid Li and Stephens model computationally tractable for these datasets, we have created a numerically exact version of the algorithm with observed average case sublinear runtime with respect to reference panel size k when tested against the 1000 Genomes dataset. Conclusions We show a forward algorithm which avoids any tradeoff between runtime and model complexity. Our algorithm makes use of two general strategies which might be applicable to improving the time complexity of other future sequence analysis algorithms: sparse dynamic programming matrices and lazy evaluation.


Background
Probabilistic models of haplotypes describe how variation is shared in a population. One application of these models is to calculate the probability P(o|H), defined as the probability of a haplotype o being observed, given the assumption that it is a member of a population represented by a reference panel of haplotypes H. This computation has been used in estimating recombination rates [1], a problem of interest in genetics and in medicine. It may also be used to detect errors in genotype calls.
Early approaches to haplotype modeling used coalescent [2] models which were accurate but computationally complex, especially when including recombination. Li and Stephens wrote the foundational computationally tractable haplotype model [1] with recombination. Under their model, the probability P(o|H) can be calculated using the forward algorithm for hidden Markov models (HMMs) and posterior sampling of genotype probabilities can be achieved using the forward-backward algorithm. Generalizations of their model have been used for haplotype phasing and genotype imputation [3][4][5][6][7].

The Li and Stephens model
Consider a reference panel H of k haplotypes sampled from some population. Each haplotype h j ∈ H is a sequence (h j,1 , . . . , h j,n ) of alleles at a contiguous sequence 1, . . . , n of genetic sites. Classically [1], the sites are biallelic, but the model extends to multiallelic sites [8].
Consider an observed sequence of alleles o = (o 1 , . . . , o n ) representing another haplotype. The monoploid Li and Stephens model (LS) [1] specifies a probability that o is descended from the population represented by H. LS can be written as a hidden Markov model wherein the haplotype o is assembled by copying (with possible error) consecutive contiguous subsequences of haplotypes h j ∈ H. Rosen and Paten Algorithms Mol Biol (2019) 14:11 Definition 1 (Li and Stephens HMM) Define x j,i as the event that the allele o i at site i of the haplotype o was copied from the allele h j,i of haplotype h j ∈ H . Take parameters and from them define the transition and recombination probabilities We will write µ i (j) as shorthand for p(o i |x j,i ) . We will also define the values of the initial probabilities p(x j,1 , o 1 |H ) = µ 1 (j) k , which can be derived by noting that if all haplotypes have equal probabilities 1 k of randomly being selected, and that this probability is then modified by the appropriate emission probability.
Let P(o|H) be the probability that haplotype o was produced from population H. The forward algorithm for hidden Markov models allows calculation of this probability in O(nk 2 ) time using an n × k dynamic programming matrix of forward states The probability P(o|H) will be equal to the sum j p n [j] of all entries in the final column of the dynamic programming matrix. In practice, the Li and Stephens forward algorithm is O(nk) (see "Efficient dynamic programming" section).

Li and Stephens like algorithms for large populations
The O(nk) time complexity of the forward algorithm is intractable for reference panels with large size k. The UK Biobank has amassed k = 500, 000 array samples. Whole genome sequencing projects, with a denser distribution of sites, are catching up. Major sequencing projects with k = 100, 000 or more samples are nearing completion. Others numbering k in the millions have been announced. These large population datasets have significant potential benefits: They are statistically likely to more accurately represent population frequencies and probability of any recombination between sites i − 1 and i (2) µ i probability of a mutation from one allele to another at site i those employing genome sequencing can provide phasing information for rare variants.
In order to handle datasets with size k even fractions of these sizes, modern haplotype inference algorithms depend on models which are simpler than the Li and Stephens model or which sample subsets of the data. For example, the common tools Eagle-2, Beagle, HAPI-UR and Shapeit-2 and -3 [3][4][5][6][7] either restrict where recombination can occur, fail to model mutation, model longrange phasing approximately or sample subsets of the reference panel.
Lunter's "fastLS" algorithm [8] demonstrated that haplotypes models which include all k reference panel haplotype could find the Viterbi maximum likelihood path in time sublinear in k, using preprocessing to reduce redundant information in the algorithm's input. However, his techniques do not extend to the forward and forwardbackward algorithms.

Our contributions
We have developed an arithmetically exact forward algorithm whose expected time complexity is a function of the expected allele distribution of the reference panel. This expected time complexity proves to be significantly sublinear in reference panel size. We have also developed a technique for succinctly representing large panels of haplotypes whose size also scales as a sublinear function of the expected allele distribution.
Our forward algorithm contains three optimizations, all of which might be generalized to other bioinformatics algorithms. In "Sparse representation of haplotypes" section, we rewrite the reference panel as a sparse matrix containing the minimum information necessary to directly infer all allele values. In "Efficient dynamic programming" section, we define recurrence relations which are numerically equivalent to the forward algorithm but use minimal arithmetic operations. In "Lazy evaluation of dynamic programming rows", we delay computation of forward states using a lazy evaluation algorithm which benefits from blocks of common sequence composed of runs of major alleles. Our methods apply to other models which share certain redundancy properties with the monoploid Li and Stephens model.

Sparse representation of haplotypes
The forward algorithm to calculate the probability P(o|H) takes as input a length n vector o and a k × n matrix of haplotypes H. In general, any algorithm which is sublinear in its input inherently requires some sort of preprocessing to identify and reduce redundancies in the data. However, the algorithm will truly become effectively sublinear if this preprocessing can be amortized over many iterations. In this case, we are able to preprocess H into a sparse representation which will on average contain better than O(nk) data points. This is the first component of our strategy. We use a variant of column-sparse-row matrix encoding to allow fast traversal of our haplotype matrix H. This encoding has the dual benefit of also allowing reversible size compression of our data. We propose that this is one good general data representation on which to build other computational work using very large genotype or haplotype data. Indeed, extrapolating from our single-chromosome results, the 1000 Genomes Phase 3 haplotypes across all chromosomes should simultaneously fit uncompressed in 11 GB of memory.
We will show that we can evaluate the Li and Stephens forward algorithm without needing to uncompress this sparse matrix.

Sparse column representation of haplotype alleles
Consider a biallelic genetic site i with alleles {A, B} . Consider the vector h 1,i , h 2,i , . . . , h k,i ∈ {A, B} k of alleles of haplotypes j at site i. Label the allele A, B which occurs more frequently in this vector as the major allele 0, and the one which occurs less frequently as the minor allele 1. We then encode this vector by storing the value A or B of the major allele 0, and the indices j 1 , j 2 , . . . of the haplotypes which take on allele value 1 at this site.
We will write φ i for the subvector h j 1 ,i , h j 2 ,i , . . . of alleles of haplotypes consisting of those haplotypes which possess the minor allele 1 at site i. We will write |φ i | for the multiplicity of the minor allele. We call this vector φ i the information content of the haplotype cohort H at the site i.

Relation to the allele frequency spectrum
Our sparse representation of the haplotype reference panel benefits from the recent finding [9] that the distribution over sites of minor allele frequencies is biased towards low frequencies. 1 Clearly, the distribution of |φ i | is precisely the allele frequency spectrum. More formally,

Lemma 1
Let E[f ](k) be the expected mean minor allele frequency for k genotypes. Then

Dynamic reference panels
Adding or rewriting a haplotype is constant time per site per haplotype unless this edit changes which allele is the most frequent. It can be achieved by addition or removal or single entries from the row-sparse-column representation, wherein, since our implementation does not require that the column indices be stored in order, these operations can be made O(1) . This allows our algorithm to extend to uses of the Li and Stephens model where one might want to dynamically edit the reference panel. The exception occurs when φ i = k 2 -here it is not absolutely necessary to keep the formalism that the indices stored actually be the minor allele.

Implementation
For biallelic sites, we store our φ i 's using a length-n vector of length |φ i | vectors containing the indices j of the haplotypes h j ∈ φ i and a length-n vector listing the major allele at each site (see Fig. 1 panel iii) Random access by key i to iterators to the first elements of sets φ i is O(1) and iteration across these φ i is linear in the size of φ i . For multiallelic sites, the data structure uses slightly more space but has the same speed guarantees.
Generating these data structures takes O(nk) time but is embarrassingly parallel in n. Our "*.slls" data structure doubles as a succinct haplotype index which could be distributed instead of a large vcf record (though genotype likelihood compression is not accounted for). A vcf → slls conversion tool is found in our github repository.

Efficient dynamic programming
We begin with the recurrence relation of the classic forward algorithm applied to the Li and Stephens model [1]. To establish our notation, recall that we k . For i > 1 , we may then write: We will reduce the number of summands in (8) and reduce the number indices j for which (7) is evaluated. This will use the information content defined in "Sparse column representation of haplotype alleles" section. Rosen  Proof Suppose first that µ i (j) = µ i for all j. Then Now suppose that µ i (j) = 1 − µ i for some set of j. We must then correct for these j. This gives us The same argument holds when we reverse the roles of µ i and 1 − µ i . Therefore we can choose which calculation to perform based on which has fewer summands. This gives us the following formula: We note another redundancy in our calculations. For the proper choices of µ ′ i , µ ′′ i among µ i , 1 − µ i , the recurrence relations (7) are linear maps R → R of which there are precisely two unique maps, f i corresponding to the recurrence relations for those x j such that j ∈ φ i , and F i to those such that j / ∈ φ i .
Proof Equation (12) lets us calculate S i−1 without knowing any p i−1 [j] for any j / ∈ φ i−1 . From S i−1 we also have f i and F i . Therefore, we can calculate provided that j ′ � = j . This then shows us that we can calculate p i [j ′ ] for all j ′ ∈ φ i without knowing any j such that j / ∈ φ i and j / ∈ φ i−1 . Finally, the first statement follows from another application of (12) (Fig. 2). (8) and the minimum set of recurrences (7) needed to compute (8)

Corollary 2 The recurrences
We address the assumption on prior calculation of the necessary p i−1 [j] 's in "Lazy evaluation of dynamic programming rows" section.

Time complexity
Recall that we defined E[f ](k) as the expected mean minor allele frequency in a sample of size k. Suppose that it is comparatively trivial to calculate the missing p i−1 [j] (15)

Lazy evaluation of dynamic programming rows
Corollary 2 was conditioned on the assumption that specific forward probabilities had already been evaluated. We will describe a second algorithm which performs this task efficiently by avoiding performing any arithmetic which will prove unnecessary at future steps. 2

Then the dynamic programming matrix entries
need not be calculated in order to calculate S ℓ , S ℓ+1 , . . . , S i−1 .
Proof By repeated application of Lemma (3).

Corollary 3 Under the same assumption on j,
need not be calculated in order to calculate F ℓ+1 , . . . , F i . This is easily seen by definition of F i .

Lemma 5 Suppose that
can be calculated in the time which it takes to calculate It is immediately clear that calculating the p i [j] lends well to lazy evaluation. Specifically, the x j / ∈ φ i are data which need not be evaluated yet at step i. Therefore, if we can aggregate the work of calculating these data at a later iteration of the algorithm, and only if needed then, we can potentially save a considerable amount of computation.

Definition 2 (Longest major allele suffix classes) Define
That is, let E ℓ→i−1 be the class of all haplotypes whose sequence up to site i − 1 shares the suffix from ℓ to i − 1 inclusive consisting only of major alleles, but lacks any longer suffix composed only of major alleles.
Remark 1 E ℓ→i−1 is the set of all h j where p ℓ−1 [j] was needed to calculate S ℓ−1 but no p (·) [j] has been needed to calculate any S (·) since.
Note that for each i, the equivalence classes E ℓ→i−1 form a disjoint cover of the set of all haplotypes h j ∈ H.

The lazy evaluation algorithm
Our algorithm will aim to: To accomplish these goals, at each iteration i, we maintain the following auxiliary data. The meaning of these are clarified by reference to Figs. 3, 4 and 5.  Stephens HMM recurrence relations. ii Our procedure specified in Eq. (12). Black lines correspond to arithmetic operations; operations which cannot be parallelized over j are colored yellow Rosen and Paten Algorithms Mol Biol (2019) 14:11 These are used to rapidly extend prefixes F ℓ→m into prefixes F ℓ→i−1 .
Finally, we will need the following ordering on tuples T ℓ to describe our algorithm: Definition 4 Impose a partial ordering < on the Fig. 4.
We are now ready to describe our lazy evaluation algorithm which evaluates p i [j] = f i • F ℓ→i−1 (p ℓ−1 [j]) justin-time while fulfilling the aims listed at the top of this section, by using the auxiliary state data specified above.  The algorithm is simple but requires keeping track of a number of intermediate indices. We suggest referring to the Figs. 3, 4 and 5 as a visual aid. We state it in six steps as follows.
Step 1: Identifying the tuples containing φ-O(φ i ) time complexity Identify the subset U (φ) of the tuples T ℓ for which there exists some h j ∈ φ i such that h j ∈ E ℓ→i−1 .
Step 2: Identifying the preparatory map suffix calculations to be performed-O(φ i ) time complexity Find the maximum depth d of any T ℓ ∈ U (φ) with respect to the partial ordering above. Equivalently, find the minimum m such that T ℓ = (E ℓ→i−1 , F ℓ→m , m) ∈ U (φ) . See  In this case, we wish to calculate p 6 [2] . This is a member of the equivalence class E 2→5 , since it hasn't needed to be calculated since time 1. In step 4 of the algorithm, we therefore must update the whole tuple T 2 by post-composing the partially completed prefix F 2→4 of the map F 2→5 which we need using our already-calculated suffix map F 5 . In step 5, we use F 2→5 to compute p 6 [2] = f 6 • F 2→5 (p 1 [j]) . In step 6, we update the tuple T 2 to reflect its loss of h 2 , which is now a member of E 6→6 Step 4: Performing the deferred calculations for the tuples containing h j ∈ φ i -O(φ i ) time complexity If not already done in Step 3.2, for every T ℓ ∈ U (φ) , extend its map element from (E ℓ→i−1 , F ℓ→m , m) to The following two trivial lemmas allow us to bound d by k such that the aggregate time complexity of the lazy evaluation algorithm cannot exceed O(nk) . Due to the irregularity of the recursion pattern used by the algorithm, is likely not possible to calculate a closed-form tight bound on i d , however, empirically it is asymptotically dominated by i φ i as shown in the results which follow.

Lemma 6
The number of nonempty equivalence classes E ℓ→i−1 in existence at any iteration i of the algorithm is bounded by the number of haplotypes k.
Proof Trivial but worth noting.

Lemma 7
The number of unique indices m in existence at any iteration i of the algorithm is bounded by the number of nonempty equivalence classes E ℓ→i−1 .

Implementation
Our algorithm was implemented as a C++ library located at https ://githu b.com/yohei rosen /subli near-Li-Steph ens. Details of the lazy evaluation algorithm will be found there. We also implemented the linear time forward algorithm for the haploid Li and Stephens model in C++ as to evaluate it on identical footing. Profiling was performed using a single Intel Xeon X7560 core running at 2.3 GHz on a shared memory machine. Our reference panels H were the phased haplotypes from the 1000 Genomes [10] phase 3 vcf records for chromosome 22 and subsamples thereof. Haplotypes o were randomly generated simulated descendants.

Minor allele frequency distribution for the 1000 Genomes dataset
We found it informative to determine the allele frequency spectrum for the 1000 Genomes dataset which we will use in our performance analyses. We simulated haplotypes o of 1,000,000 bp length on chromosome 22 and recorded the sizes of the sets φ i (o i ) for k = 5008 . These data produced a mean |φ i (o i )| of 59.9, which is 1.2% of the size of k. We have plotted the distribution of |φ i (o i )| which we observed from this experiment in (Fig. 6). It is skewed toward low frequencies; the minor allele is unique at 71% of sites, and it is below 1% frequency at 92% of sites.

Comparison of our algorithm with the linear time forward algorithm
In order to compare the dependence of our algorithm's runtime on haplotype panel size k against that of the standard linear LS forward algorithm, we measured the CPU time per genetic site of both across a range of haplotype panel sizes from 30 to 5008. This analysis was achieved as briefly described above. Haplotype panels spanning the range of sizes from 30 to 5008 haplotypes were subsampled from the 1000 Genomes phase 3 vcf records and loaded into memory in both uncompressed and our column-sparse-row format. Random sequences were sampled using a copying model with mutation and recombination, and the performance of the classical forward algorithm was run back to back with our algorithm for the same random sequence and same subsampled haplotype panel. Each set of runs was performed in triplicate to reduce stochastic error. Figure 7 shows this comparison. Observed time complexity of our algorithm was O(k 0.35 ) as calculated from the slope of the line of best fit to a log-log plot of time per site versus haplotype panel size.
For data points where we used all 1000 Genomes project haplotypes ( k = 5008 ), on average, time per site is 37 μs for our algorithm and 1308 μs for the linear LS algorithm. For the forthcoming 100,000 Genomes Project, these numbers can be extrapolated to 251 μs for our algorithm and 260,760 μs for the linear LS algorithm.

Lazy evaluation of dynamic programming rows
We also measured the time which our algorithm spent within the d-dependent portion of the lazy evaluation subalgorithm. In the average case, the time complexity of our lazy evaluation subalgorithm does not contribute to the overall algebraic time complexity of the algorithm (Fig. 8, right). The lazy evaluation runtime also contributes minimally to the total actual runtime of our algorithm (Fig. 8, left).

Sparse haplotype encoding Generating our sparse vectors
We generated the haplotype panel data structures from "Sparse representation of haplotypes" section using the vcf-encoding tool vcf2slls which we provide. We built indices with multiallelic sites, which increases their time and memory profile relative to the results in "Minor allele frequency distribution for the 1000 Genomes dataset" section but allows direct comparison to vcf records. Encoding of chromosome 22 was completed in 38 min on a single CPU core. Use of M CPU cores will reduce runtime proportional to M.

Size of sparse haplotype index
In uncompressed form, our whole genome *.slls index for chromosome 22 of the 1000 genomes dataset was 285 MB in size versus 11 GB for the vcf record using uint16_t's to encode haplotype ranks. When compressed with gzip, the same index was 67 MB in size versus 205 MB for the vcf record.
In the interest of speed (both for our algorithm and the O(nk) algorithm) our experiments loaded entire chromosome sparse matrices into memory and stored haplotype indices as uint64_t's. This requires on the order of 1 GB memory for chromosome 22. For long chromosomes or larger reference panels on low memory machines, the algorithm can operate by streaming sequential chunks of the reference panel.

Discussions and Conclusion
To the best of our knowledge, ours is the first forward algorithm for any haplotype model to attain sublinear time complexity with respect to reference panel size. Our algorithms could be incorporated into haplotype inference strategies by interfacing with our C++ library. This opens the potential for tools which are tractable on haplotype reference panels at the scale of current 100,000 to 1,000,000+ sample sequencing projects.

Applications which use individual forward probabilities
Our algorithm attains its runtime specifically for the problem of calculating the single overall probability P(o|H , ρ, µ) and does not compute all nk forward probabilities. We can prove that if m many specific forward probabilities are also required as output, and if the time complexity of our algorithm is O( i |φ i |) , then the time complexity of the algorithm which also returns the m forward probabilities is O( i |φ i | + m).
In general, haplotype phasing or genotype imputation tools use stochastic traceback or other similar sampling algorithms. The standard algorithm for stochastic traceback samples states from the full posterior distribution and therefore requires all forward probabilities. The algorithm output and lower bound of its speed is therefore O(nk) . The same is true for many applications of the forward-backward algorithm.
There are two possible approaches which might allow runtime sublinear in k for these applications. Using stochastic traceback as an example, first is to devise an O(f (m)) sampling algorithm which uses m = g(k) forward probabilities such that O(f • g(k)) < O(k) . The second is to succinctly represent forward probabilities such that nested sums of the nk forward probabilities can be queried from O(φ) < O(nk) data. This should be possible, perhaps using the positional Burrows-Wheeler transform [11] as in [8], since we have already devised a forward algorithm with this property for a different model in [12].

Generalizability of algorithm
The optimizations which we have made are not strictly specific to the monoploid Li and Stephens algorithm. Necessary conditions for our reduction in the time complexity of the recurrence relations are   The reduction in time complexity of the recurrence relations depends on the Markov property, however we hypothesize that the delayed evaluation needs only the semi-Markov property.

Other haplotype forward algorithms
Our optimizations are of immediate interest for other haplotype copying models. The following related algorithms have been explored without implementation.

Example 1 (Diploid Li and Stephens)
We have yet to implement this model but expect average runtime at least subquadratic in reference panel size k. We build on the statement of the model and its optimizations in [13]. We have found the following recurrences which we believe will work when combined with a system of lazy evaluation algorithms:

Lemma 8 The diploid Li and Stephens HMM may be expressed using recurrences of the form which use on the intermediate sums defined as
where α (·) , β (·) , γ (·) depend only on the diploid genotype o i . Implementing and verifying the runtime of this extension of our algorithm will be among our next steps.
Example 2 (Multipopulation Li and Stephens) [14] We maintain separate sparse haplotype panel representations φ A i (o i ) and φ B i (o i ) and separate lazy evaluation mechanisms for the two populations A and B. Expected runtime guarantees are similar.
This model, and versions for > 2 populations, will be important in large sequencing cohorts (such as NHLBI TOPMed) where assuming a single related population is unrealistic.
(17) p i [j 1 , j 2 ] = α p p i−1 [j 1 , j 2 ] + β p (S i−1 (j 1 ) + S i−1 (j 2 )) + γ p S i−1 Example 3 (More detailed mutation model) It may also be desirable to model distinct mutation probabilities for different pairs of alleles at multiallelic sites. Runtime is worse than the biallelic model but remains average case sublinear.
Example 4 (Sequence graph Li and Stephens analogue) In [12] we described a hidden Markov model for a haplotype-copying with recombination but not mutation in the context of sequence graphs. Assuming we can decompose our graph into nested sites then we can achieve a fast forward algorithm with mutation. An analogue of our row-sparse-column matrix compression for sequence graphs is being actively developed within our research group.
While a haplotype HMM forward algorithm alone might have niche applications in bioinformatics, we expect that our techniques are generalizable to speeding up other forward algorithm-type sequence analysis algorithms.