Linear time minimum segmentation enables scalable founder reconstruction

Background We study a preprocessing routine relevant in pan-genomic analyses: consider a set of aligned haplotype sequences of complete human chromosomes. Due to the enormous size of such data, one would like to represent this input set with a few founder sequences that retain as well as possible the contiguities of the original sequences. Such a smaller set gives a scalable way to exploit pan-genomic information in further analyses (e.g. read alignment and variant calling). Optimizing the founder set is an NP-hard problem, but there is a segmentation formulation that can be solved in polynomial time, defined as follows. Given a threshold L and a set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {R}} = \{R_1, \ldots , R_m\}$$\end{document}R={R1,…,Rm} of m strings (haplotype sequences), each having length n, the minimum segmentation problem for founder reconstruction is to partition [1, n] into set P of disjoint segments such that each segment \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[a,b] \in P$$\end{document}[a,b]∈P has length at least L and the number \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d(a,b)=|\{R_i[a,b] :1\le i \le m\}|$$\end{document}d(a,b)=|{Ri[a,b]:1≤i≤m}| of distinct substrings at segment [a, b] is minimized over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[a,b] \in P$$\end{document}[a,b]∈P. The distinct substrings in the segments represent founder blocks that can be concatenated to form \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\max \{ d(a,b) :[a,b] \in P \}$$\end{document}max{d(a,b):[a,b]∈P} founder sequences representing the original \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {R}}$$\end{document}R such that crossovers happen only at segment boundaries. Results We give an O(mn) time (i.e. linear time in the input size) algorithm to solve the minimum segmentation problem for founder reconstruction, improving over an earlier \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(mn^2)$$\end{document}O(mn2). Conclusions Our improvement enables to apply the formulation on an input of thousands of complete human chromosomes. We implemented the new algorithm and give experimental evidence on its practicality. The implementation is available in https://github.com/tsnorri/founder-sequences.


Background
A key problem in pan-genomics is to develop a sufficiently small, efficiently queriable, but still descriptive representation of the variation common to the subject under study [1]. For example, when studying human population, one would like to take all publicly available variation datasets (e.g. [2][3][4]) into account. Many approaches encode the variation as a graph [5][6][7][8][9][10] and then one can encode the different haplotypes as paths in this graph [11]. An alternative has been proposed [12] based on a compressed indexing scheme for a multiple alignment of all the haplotypes [13][14][15][16][17]. In either approach, scalability is hampered by the encoding of all the haplotypes.
We suggest to look for a smaller set of representative haplotype sequences to make the above pan-genomic representations scalable.
Finding such set of representative haplotype sequences that retain the original contiguities as well as possible, is known as the founder sequence reconstruction problem [18]. In this problem, one seeks a set of d founders such that the original m haplotypes can be mapped with minimum amount of crossovers to the founders. Here a crossover means a position where one needs to jump from one founder to another to continue matching the content of the haplotype in question. Unfortunately, this problem is NP-hard even to approximate within a constant factor [19].  14:12 For founder reconstruction to be scalable to the pangenomic setting, one would need an algorithm to be nearly linear to the input size. With this is mind, we study a relaxation of founder reconstruction that is known to be polynomial time solvable: Namely, when limiting all the crossovers to happen at the same locations, one obtains a minimum segmentation problem specific to founder reconstruction [18]. A dynamic programming algorithm solves this problem in O(n 2 m) time [18], where m is the number of haplotypes and n is the length of each of them.
In this paper, we improve the running time of solving the minimum segmentation problem of founder reconstruction to O(mn) (linear in the input size).
We also implement the new algorithm, as well as a further heuristic that aims to minimize crossovers over the segment boundaries (yielded by the optimal solution to the minimum segmentation problem). In our experiments, we show that the approach is practical on human genome scale setting. Namely, we apply the implementation on a multiple alignment representing 5009 haplotypes of human chromosome 6, and the result is 130 founder sequences with the average distance of two crossovers being 9624 bases. Preserving such long contiguities in just 2.5% of the original input space is promising for the accuracy and scalability of the short read alignment and variant calling motivating our study.
The main technique behind the improvement is the use of positional Burrows-Wheeler transform (pBWT) [20], and more specifically its extension to larger alphabets [21]. While the original dynamic programming solution uses O(nm) time to look for the best preceding segment boundary for each column of the input, we observe that at most m values in pBWT determine segment boundaries where the number of distinct founder substrings change. Minimums on the already computed dynamic programming values between each such interesting consecutive segment boundaries give the requested result. However, it turns out that we can maintain the minimums directly in pBWT internal structures (with some modifications) and have to store only the last L computed dynamic programming values, thus spending only O(m + L) additional space, where L is the input threshold on the length of each segment. The segmentation is then reconstructed by standard backtracking approach in O(n) time using an array of length n.
Preliminary version of this work appeared in WABI 2018 [22].

Notation and problem statement
For a string s = c 1 c 2 · · · c n , denote by |s| its length n. We write s [i] for the letter c i of s and s [i, j] for the substring c i c i+1 . . . c j . An analogous notation is used for arrays. For any numbers i and j, the set of integers {x ∈ Z : i ≤ x ≤ j} (possibly empty) is denoted by [i, j].
The input for our problem is the set R = {R 1 , . . . , R m } of strings of length n, called recombinants. A set F = {F 1 , . . . , F d } of strings of length n is called a founder set of R if for each string R i ∈ R , there exists a partition P i of the segment [1, n] into disjoint subsegments such that, for each The partition P i together with the mapping of the segments [a, b] ∈ P i to substrings F j [a, b] is called a parse of R i in terms of F , and a set of parses for all R i ∈ R is called a parse of R in terms of F . The integers a and b + 1 , for [a, b] ∈ P i , are called crossover points; thus, in particular, 1 and n + 1 are always crossover points.
It follows from the definition that, in practice, it makes sense to consider founder sets only for pre-aligned recombinants. Throughout the paper we implicitly assume that this is the case, although all our algorithms, clearly, work in the unaligned setting too but the produce results may hardly make any sense.
We consider the problem of finding a "good" founder set F and a "good" corresponding parse of R according to a reasonable measure of goodness. Ukkonen [18] pointed out that such measures may contradict each other: for instance, a minimum founder set obviously has size d = max j∈[1,n] |{R 1 [j], . . . , R m [j]}| , but parses corresponding to such set may have unnaturally many crossover points; conversely, R is a founder set of itself and the only crossover points of its trivial parse are 1 and n + 1 , but the size m of this founder set is in most cases unacceptably large. Following Ukkonen's approach, we consider compromise parameterized solutions. The minimum founder set problem is, given a bound L and a set of recombinants R , to find a smallest founder set F of R such that there exists a parse of R in terms of F in which the distance between any two crossover points is at least L (the crossover points may belong to parses of different recombinants, i.e., for [a, b] ∈ P i and [a ′ , b ′ ] ∈ P j , where P i and P j are parses of R i and R j , we have either a = a ′ or |a − a ′ | ≥ L).
It is convenient to reformulate the problem in terms of  [18] is, given a bound L and a set of recombinants R , to find a segmentation S of R such that max{|R[j, k]| : R[j, k] ∈ S} is minimized and the length of each segment from S is at least L; in other words, the problem is to compute where S L is the set of all segmentations in which all segments have length at least L.
The minimum founder set problem and the minimum segmentation problem are connected: any segmentation S with segments of length at least L induces in an obvious way a founder set of size max{|R[j, k]| : R[j, k] ∈ S} and a parse in which all crossover points are located at segment boundaries (and, hence, at distance at least L from each other); conversely, if F is a founder set of R and {j 1 , . . . , j p } is the sorted set of all crossover points in a parse of R such that j q − j q−1 ≥ L for q ∈ [2, p] , then S = {R[j q−1 , j q −1] : q ∈ [2, p]} is a segmentation of R with segments of length at least L and Our main result is an algorithm that solves the minimum segmentation problem in O(mn) time (linear in the input size). The solution normally does not uniquely define a founder set of R : for instance, if the built segmentation of R = {baaaa, baaab, babab} is S = {R [1,1], R [2,3], R [4,5]} , then the possible founder sets induced by S are F 1 = {baaab, babaa} and F 2 = {baaaa, babab} . In other words, to construct a founder set, one concatenates fragments of recombinants corresponding to the found segments in a certain order. We return to this ordering problem in the section describing experiments and now focus on the details of the segmentation problem.
Hereafter, we assume that the input alphabet is the set [0, | |−1] of size O(m), which is a natural assumption considering that the typical alphabet size is 4 in our problem. It is sometimes convenient to view the set R = {R 1 , . . . , R m } as a matrix with m rows and n columns. We say that an algorithm processing the recombinants R is streaming if it reads the input from left to right "columnwise", for each k from 1 to n, and outputs an answer for each set of recombinants The main result of the paper is the following theorem.

Minimum segmentation problem
Given a bound L and a set of recombinants R = {R 1 , . . . , R m } each of which has length n, Ukkonen [18] proposed a dynamic programming algorithm that solves the minimum segmentation problem in O(mn 2 ) time based on the following recurrence relation: It is obvious that M(n) is equal to the solution (1); the segmentation itself can be reconstructed by "backtracking" in a standard way [18]. We build on the same approach. For a given k ∈ [1, n] , denote by j k,1 , . . . , j k,r k the sequence of all positions j ∈ [1, k − L] in which the value of |R[j, k]| changes, i.e., 1 ≤ j k,1 < · · · < j k,r k ≤ k − L and We complement this sequence with j k,0 = 0 and j k,r k +1 = k − L + 1 , so that j k,0 , . . . , j k,r k +1 can be interpreted as a splitting of the range [0, k − L] into segments in which the value |R[j + 1, k]| stays the same: namely, for h ∈ [0, r k ] , one has , k]|, min j k,h ≤j<j k,h+1 M(j)} and, therefore, (2) can be rewritten as follows: Our crucial observation is that, for k ∈ [1, n] and j ∈ [1, k] , one has |R[j + 1, k]| ≤ |R[j, k]| ≤ m . Therefore, m ≥ |R[j k,1 , k]| > · · · > |R[j k,r k +1 , k]| ≥ 1 and r k < m . Hence, M(k) can be computed in O(m) time using (3), provided one has the following components: i. the sorted sequence j k,1 , . . . , j k,r k ii. the numbers |R[j k,h+1 , k]| , for h ∈ [0, r k ] iii. the values min{M(j) : j k,h ≤ j < j k,h+1 } , for h ∈ [0, r k ].
In the remaining part of the section, we describe a streaming algorithm that reads the strings {R 1 , . . . , R m } "columnwise" from left to right and computes the components (i), (ii), and (iii) immediately after reading each "column" , positions j on which these minima are attained (see below). Further details are straightforward and, thence, omitted.

Positional Burrows-Wheeler transform
Let us fix k ∈ [1, n] . Throughout this subsection, the . Given a set of recombinants R = {R 1 , . . . , R m } each of which has length n, a positional Burrows-Wheeler transform (pBWT), as defined by Durbin [20], is a pair of integer arrays a k [1, m] and d k [1, m]  Example 1 Consider the following example, where m = 6 , k = 7 , and = {a, c, t} . It is easy to see that the pBWT implicitly encodes the trie depicted in the right part of Fig. 1, and such interpretation drives the intuition behind this structure: The trie represents the reversed , read from right to left) in lexicographic order. Leaves (values a k ) store the corresponding input indices. The branches correspond to values d k (the distance from the root subtracted from k + 1 ). Our main algorithm in this paper makes implicitly a sweep-line over the trie stopping at the branching positions.
Durbin [20] showed that a k and d k can be computed from a k−1 and d k−1 in O(m) time on the binary alphabet. Mäkinen and Norri [21] further generalized the construction for integer alphabets of size O(m), as in our case. For the sake of completeness, we describe in this subsection the generalized solution [21] (see Algorithm 1), which serves then as a basis for our main algorithm. We also present a modification of this solution (see Algorithm 2), which, albeit seems to be slightly inferior in theory (we could prove only O(m log |�|) time upper bound), showed better performance in practice and thus, as we believe, is interesting by itself.
Proof Given a k−1 and d k−1 , we are to show that Algorithm 1 correctly computes a k and d k . Since, for lexicographically, it is easy to see that the array a k can be deduced from a k−1 by radix sorting the sequence of pairs Further, since, by definition of a k−1 , the second components of the pairs are already in a sorted order, it remains to sort the first components by the counting sort. Accordingly, in Algorithm 1, the first loop counts occurrences of letters in the sequence using an auxiliary array C[0, | |] ; as is standard in the counting sort, the second loop modifies the array C so that, for each letter , and performing the assignments The last three lines of the algorithm are responsible for computing d k . Denote the length of the longest common prefix of any strings s 1 and s 2 by LCP(s 1 , s 2 ) . The computation of d k relies on the following well-known fact: given a sequence of strings s 1 , . . . , s r such that s 1 ≤ · · · ≤ s r lexicographically, one has LCP(s 1 , Suppose that the last loop of the algorithm, which iterates through In order to calculate the maximums we build a range maximum query (RMQ) data structure on the array d k−1 [1, m] in O(m) time [23]. Therefore, the running time of Algorithm 1 is O(m). In practice the bottleneck of the algorithm is the RMQ data structure, which, although answers queries in O(1) time, has a sensible constant under the big-O in the construction time. We could naively compute the maximums by scanning the ranges d k−1 [P[b]+1, i] from left to right but such algorithm works in quadratic time since same ranges of d k−1 might be processed many times in the worst case. Our key idea is to store the work done by a simple scanning algorithm to reuse it in future queries. We store this information right in the arrays a k−1 and d k−1 rewriting them; in particular, since a k−1 is accessed sequentially from left to right in the last loop, the range a k−1 [1, i] is free to use after the ith iteration.
More precisely, after the ith iteration of the last loop, the subarrays a k−1 [1, i] and d k−1 [1, i] are modified so that the following invariant holds: for any j ∈ after this, we redirect the "jump pointers" in a k−1 to i + 1 and update the maximums in d k−1 accordingly. This idea is implemented in Algorithm 2. Notice the new line a k−1 [i] ← i + 1 in the main loop (it is commented), which erases a k−1 [i] and makes it a part of the "jump table". The correctness of the algorithm is clear. But it is not immediate even that the algorithm works in O(m log m) time. The next lemma states that the bound is actually even better, O(m log |�|).
The latter inequality follows from the fact that the query ranges correponding to the same symbol are non-overlapping.
We say that a position j is touched if the function maxd is called with its first argument equal to j. Since for each i the first call to maxd is with different j, it suffices to prove that the total number of touches is O(m logl) . While processing the query maxd(i−ℓ i , i) , we may have touched many positions. Denote the sequence of all such position, for the given i, by i 1 , . . . , i r ; in other words, at the time of the query maxd(i−ℓ i , i) , we have Fig. 2). We count separately the total number of scaling and nonscaling touches in all i.
For position j, denote by p(j) the number of non-scaling touches of j. We are to prove that P = m j=1 p(j) ≤ 2m logl . Let q h (j) denote the value of a k−1 [j] − j in the hth non-scaling touch of j, for h ∈ [1, p(j)] . Suppose that this hth touch happens during the processing of a query maxd(i − ℓ i , i) . By the definition, j + q h (j) follows j in the sequence of touched positions. Since the touch of j is non-scal-  whenever all p(j) are equal, i.e., m j=1 2 p(j) ≥ m j=1 2 P/m . Hence, once P > 2m logl , we obtain m j=1 p(j) h=1 q h (j) ≥ m j=1 2 P/m − m > ml 2 − m , which is larger than ml for l ≥ 2 (for the case l < 2 the claim follows directly), contradicting m j=1 p(j) h=1 q h (j) ≤ ml . Thus, P = m j=1 p(j) ≤ 2m logl. It remains to consider scaling touches. The definition implies that each query maxd(i−ℓ i , i) performs at most log ℓ i scaling touches. Thus, it suffices to upperbound m i=1 log ℓ i . Since the function log is concave, the sum m i=1 log ℓ i is maximized whenever all ℓ i are equal, i.e.,

Modification of the pBWT
We are to modify the basic pBWT construction algorithm in order to compute the sequence j k,1 , . . . , j k,r k of all positions j ∈ [1, k − L] in which |R[j, k]| � = |R[j + 1, k]| , and to calculate the numbers |R[j k,h+1 , k]| and min{M(j) : j k,h ≤ j < j k,h+1 } , for h ∈ [0, r k ] (assuming j k,0 = 0 and j k,r k +1 = k − L + 1 ); see the beginning of the section. As it follows from (3), these numbers are sufficient to calculate M(k), as defined in (2) and (3), in O(m) time. The following lemma reveals relations between the sequence j k,1 , . . . , j k,r k and the array d k .  Our modified algorithm does not store d k but stores the following four arrays (but we still often refer to d k for the sake of analysis):

Lemma 4 Consider recombinants
• s k [1, r] contains all distinct elements from d k [1, m] The arrays s k and e k together emulate d k . The array t k will be used to calculate some numbers |R[j, k]| required to compute M(k).   Suppose that ℓ > 0 . It suffices to show that By (4) and Lemma 5, one can calculate M(k) in O(m) time using the arrays t k and u k . It remains to describe how we maintain a k , e k , s k , t k , u k .   [1, m].

Lemma 6 Algorithm 3 computes the arrays
By definition, s k must contain only elements from d k , but this is not necessarily the case in line 14. In order to fix s k and t k , we simply have to remove all elements s k . (This is not completely correct since M ′ has changed as k was increased; namely, M ′ (k − L) was equal to +∞ but now is equal to M(k − L) ). As we join segments removing some elements from s k in the loop 16-23, the array u k must be fixed accordingly: If, as in our case, one does not need s k , t k , u k for all k, the arrays s k , t k , u k can be modified in-place, i.e., s k , t k , u k can be considered as aliases for s k−1 , t k−1 , u k−1 , and yet the algorithm remains correct. Thus, we really need only 7 arrays in total: a k , a k−1 , e k , e k−1 , s, t, u, where s, t, u serve as s k , t k , u k and the array tmp can be organized in place of a k−1 or e k−1 . It is easy to maintain along with each value u k [j] a corresponding position ℓ such that u k [j] = M ′ (ℓ) ; these positions can be used then to restore the found segmentation of R using backtracking (see the beginning of the section). To compute e k , instead of using an RMQ data structure, one can adapt in an obvious way Algorithm 2 rewriting the arrays a k−1 and e k−1 during the computation, which is faster in practice but theoretically takes O(m log |�|) time by Lemma 3. We do not discuss further details as they are straightforward.

From segmentation to founder set
Now we are given a segmentation S of R and we wish to produce a founder set F that obeys the segment boundaries. Recall that such founder set corresponds to a parse P of R with respect to segmentation S . We conjecture that finding an optimal parse/founder set that minimizes the number of crossovers at segment boundaries is an NP-hard problem, but unfortunately we have not been able to prove this claim. Therefore, we continue by proposing three natural strategies of which two latter have interesting theoretical properties. The first of the strategies is a naive baseline, second is a greedy strategy, and third one is based on maximum weight perfect matching in a bipartite graph analogous to one by Ukkonen [18]. This latter strategy provides an optimal solution for a special case, and greedy gives a 2-approximation for the same special case. We will present all the three strategies first for the special case and then describe how to turn the general case to this special case (however loosing all optimality guarantees while doing so). We compare the naive baseline with the perfect matching in our experiments. Assume (for our special case) that each segment in S induces exactly M(n) distinct substrings in R . Then the naive baseline strategy to produce a founder set is to concatenate the distinct substrings of segment 1 with the distinct substrings of segment 2 in random order, and continue this process form left to right until M(n) founder sequences of length n are produced. For the latter two strategies, the idea is that instead of a random permutation, we aim to find a permutation that gives a concatenation order that minimizes the number of crossovers at each segment boundary. For this purpose, it is sufficient to consider two consecutive segments  F[a, c] of some founder F, then this concatenation causes m − |A ∩ B| crossovers. Hence, to minimize crossovers, one seeks to maximize the intersection between two partitions, studied next.
Problem of maximum intersection between two partitions. Let a be an integer. Given two partitions E 1 and E 2 of {1, . . . , a} with |E 1 | = |E 2 | , the problem of Maximum Intersection Between two Partitions (MIBP) is to find the bijection f from E 1 to E 2 which maximizes By using the bipartite graph defined between the elements of E 1 and the elements of E 2 and such that for x ∈ E 1 and y ∈ E 2 , the weight of this edge is w(x, y) = |x ∩ y| , a maximum weight perfect matching of this graph gives an optimal solution of MIBP, and hence this problem can be solved in polynomial time.
We can define the greedy algorithm related to MIBP as the the greedy algorithm related to the problem of maximum weight perfect matching in the previous bipartite graph. As the greedy algorithm for maximum weight perfect matching is 1 2 -approximation [24], we have the same ratio of approximation for the greedy algorithm for MIBP.

Lemma 7
Let E 1 and E 2 be two partitions of {1, . . . , a} with |E 1 | = |E 2 | . We can compute the greedy algorithm for MIBP of E 1 and E 2 in O(a) time.
Proof Let E be a partition of {1, . . . , a} and ≺ be a total order on E, we denote by G E the array of elements of E of size a such that for all i, G E [i] = e i where i ∈ e i ∈ E . Let be x ∈ E 1 and y ∈ E 2 . We have w(x, y) = |x ∩ y| = |{i ∈ {1, . . . , a} | i ∈ x ∩ y}| = |{i ∈ {1, . . . , a} | G E 1 [i] = x and G E 2 [i] = y}| . It follows that the number of edges of no zero weight is at most a. By using Radix sort, we can compute in O(a) the sorted array of elements of {1, . . . , a} following the order where i < j iff With this array, as for all x ∈ E 1 and y ∈ E 2 w(x, y) ≤ a , we can compute (by further Radix sort and renaming steps) in O(a) the ordered list [(x 1 , y 1 ), . . . , (x q , y q )] such that w(x 1 , y 1 ) ≥ · · · ≥ w(x q , y q ) > 0 with q ≤ a . By taking the elements in the order of this list, we can compute in O(a) two arrays f and f −1 of size |E 1 | such that Optimal founder set for the special case. Now we can solve independently the MIBP problem for each pair of consecutive segments, resulting to the following theorems, where the first one follows directly also from earlier constructions [18], and the latter from Lemma 7. Theorem 8 ([18]) Given a segmentation S of R such that each segment induces exactly K distinct substrings in R , then we can construct an optimal parse P of R (and hence the corresponding set of founders) in polynomial time.
Theorem 9 Given a segmentation S of R such that each segment induces exactly K distinct substrings in R , then we can construct a greedy parse P of R (and hence the corresponding set of founders) that has at most twice as many crossovers than the optimal parse in O(|S| × m) time and O(|S| × m) space.
In the general case, there are segments inducing less than M(n) distinct substrings. We turn such segments to the special case by duplicating some of the substrings. The choices made have dependencies between segments, and this is the reason we believe this general case is NP-hard to solve optimally. Hence, we aim just to locally optimize the chances of minimizing crossovers by duplicating distinct substrings in proportion they cover R . That is, consider a segment inducing k < M(n) distinct substrings and the corresponding partitioning E of {1, . . . , m} . Consider the largest set x of E. We make k x = ⌈ |x| m (M(n) − k)⌉ copies of the corresponding distinct substring. We continue by decreasing cardinality, stop when the sum of the duplication counts is greater than or equal to M(n) and update the last one such that k + � i k i = M(n) (see Fig. 3). The fact that the corresponding partitioning is now a multi-partitioning (containing same set multiple times), does not affect the functioning of the greedy or perfect matching algorithms for the MIBP problem.

Results
We implemented the segmentation algorithm using Algorithm 2 to build the pBWT arrays and computed the minimum number of founders with the given value of L using the recursion in Eq. 3. This part of the implementation corresponds to Lemma 3, and thus the overall time complexity of the implemented approach is O(mn log |�|) . After computing the minimum number of founders, we use backtracking to determine the optimal segmentation. Since we use the pBWT arrays to determine the distinct substrings in each segment, as part of the first phase of building the arrays we also store samples and now update them to the segment boundary positions in parallel. We proceed to join adjacent segments from left to right until the number of distinct substrings in one segment would exceed the minimum number of founders, and finally we concatenate the substrings to generate founder sequences. The implementation outputs for each segment the distinct founder sequence fragments, and associates to each fragment the set of haplotypes containing that fragment as a substring at that location (these are easily deduced given the segmentation and the positional BWT structures). Our implementation uses integer vectors from the SDSL library [25].
As our goal is to produce reference sequences for aligning short reads, we wanted to find a good value of L to generate a segmentation suitable for this purpose. In particular, we wanted to have the length of most segments clearly above a typical read length, such that most reads could be aligned without hitting a recombination site.
We used the chromosome 6 variants from the phase 3 data of the 1000 Genomes Project [2] as the starting point. We converted the variant data to a multiple sequence alignment with vcf2multialign, 1  c a c g c g c g a t g ...

Set of Founders
Optimal perfect matching Fig. 3 The duplication of the fragments and the link between optimal solution of perfect matching and the concatenation of the fragments to obtain the set of founder sequences tool, we discarded columns of identical characters as they would not affect the number of recombination sites. This reduced each sequence to approximately 5.38 million characters.
We used an increasing number of the generated sequences as an input to our tool with the value of L fixed to 10 to verify the usability of the tool in terms of running time and memory consumption. The tests were run on a Ubuntu Linux 16.04 server. The server had 96 Intel Xeon E7-4830 v3 CPUs running at 2.10GHz and 1.4 TB of memory. In addition to our own RMQ data structure, we tested with a general-purpose RMQ from the SDSL library. As seen in Fig. 4, our special-purpose RMQ data structure performed somewhat better in terms of speed compared to the general-purpose library implementation. From this experiment it is conceivable that processing of thousands of complete human genomes takes only few CPU days. As we did not optimize the memory usage of our tool, the maximum resident set size with 5009 inputs was around 257 GB which corresponds to approximately 10.25 bytes per input character. We expect that the memory consumption may be reduced without much affecting the performance.
Our second experiment was to see the effect of the minimum length L on the number of founders as well as the length of the segments. The results have been summarized in Table 1. We tested with a number of values of L ranging from 10 to 80. After generating the founders, we mapped the segment co-ordinates back to the original sequences to determine the segment lengths. The results are shown in Figs. 5 and 6. We note that while the average segment length of 2395 bases with L = 10 is fitting our purpose, there is a peak of short segments of approximately 250 bases. The peak is magnified in Fig. 7. We also tested smaller values of L to conclude that decreasing L further rapidly makes the situation more difficult. On the other hand, setting L = 10 resulted in only 130 founders, which makes aligning reads much faster than using all of the haplotypes for indexing.
We proceeded with two tests in which we measured the number of recombinations needed to express each of the original sequences with the generated founder sequences depending on the method of concatenating the fragments into the set of founder sequences. Using the method given earlier, we began by duplicating some fragments so that each segment had exactly the same amount of fragments. For these tests, we implemented the three concatenation strategies: a Random matching which corresponds to concatenating the consecutive fragments in random order, a Perfect matching which Fig. 4 The running time of our implementation plotted against the number of input sequences with L = 10 and using either our RMQ data structure or rmq_succinct_sct from SDSL. The data points have been fitted with a least-squares linear model, and the grey band shows the 95% confidence interval

Table 1 Summarized results with 5009 input sequences
We measured the average segment length from the segmentation, median number of recombinations from mapping the input sequences to the founder sequences, and average distance between recombinations by dividing the length of the original sequences by the average number of recombinations. The last three columns report the results for the perfect matching approach takes an optimal solution of the maximum weight perfect matching problem as the order for the concatenation of the fragments, and a Greedy matching which solves the matching problem greedily. For evaluating the different concatenation strategies, we mapped each one of the original sequences to the founders, using a simple greedy algorithm that is also optimal [19]. In the first test, we fixed the value of L to 10 and mapped an increasing number of input sequences to a set of founder sequences generated with the same input sequences. In the second one, we used all of the 5009 input sequences and varied the value of L. The results are shown in Figs. 8 and 9. Considering the 17768 and 43333 recombinations achieved with perfect and random matching, respectively, given 5009 input sequences and L = 10 (see Table 1), we conclude that the heuristic part of optimizing the concatenation of founder blocks yields an improvement of around 2.44 compared to a random concatenation of segments with  duplications. Greedy approach works even slighly better than perfect matching in our experiments: the number of recombinations on the same setting is 17268. As the numbers are very close, we refer to perfect matching numbers in the sequel.
The results look promising, as using 130 founders instead of 5009 haplotypes as the input to our pangenome indexing approach [12] will result into significant saving of resources; this solves the space bottleneck, and the preprocessing of founder reconstruction also saves time in the heavy indexing steps.
Our intention was to compare our tool to an implementation of Ukkonen's algorithm [19]. However, initial testing with four input sequences showed that the latter implementation is not practical with a data set of this size.

Conclusions
As our experiments indicate that one can reduce 5009 haplotypes down to 130 founders with the average distance of two crossovers being 9624 bases, one can expect short read alignment and variant calling to become practical on such pan-genomic setting. We are investigating this on our tool PanVC [12], where one can simply replace its input multiple alignment with the one made of the founder sequences. With graph-based approaches, slightly more effort is required: Input variations are encoded with respect to the reference, so one first needs to convert variants into a multiple alignment, apply the founder reconstruction algorithm, and finally convert the multiple alignment of founder sequences into a directed acyclic graph. PanVC toolbox provides the required conversions. Alternatively, one can construct the pan-genome graph using other methods, and map the founder sequences afterwards to the paths of the graph: If original haplotype sequences are already spelled as paths, each founder sequence is a concatenation of existing subpaths, and can hence be mapped to a continuous path without alignment (possibly requiring adding a few missing edges).
Finally, it will be interesting to see how much the contiguity of the founder sequences can still be improved with different formulations of the segmentation problem. We are investigating a variant with the number of founder sequenced fixed.