Haplotype-aware sequence alignment to pangenome graphs

Modern pangenome graphs are built using haplotype-resolved genome assemblies. During read mapping to a pangenome graph, prioritizing alignments that are consistent with the known haplotypes has been shown to improve genotyping accuracy. However, the existing rigorous formulations for sequence-to-graph co-linear chaining and alignment problems do not consider the haplotype paths in a pangenome graph. This often leads to spurious read alignments to those paths that are unlikely recombinations of the known haplotypes. In this paper, we develop novel formulations and algorithms for haplotype-aware sequence alignment to an acyclic pangenome graph. We consider both sequence-to-graph chaining and sequence-to-graph alignment problems. Drawing inspiration from the commonly used models for genotype imputation, we assume that a query sequence is an imperfect mosaic of the reference haplotypes. Accordingly, we extend previous chaining and alignment formulations by introducing a recombination penalty for a haplotype switch. First, we solve haplotype-aware sequence-to-graph alignment in O(|Q| | E| |ℋ|) time, where Q is the query sequence, E is the set of edges, and ℋ is the set of haplotypes represented in the graph. To complement our solution, we prove that an algorithm significantly faster than O(|Q| | E| |ℋ|) is impossible under the Strong Exponential Time Hypothesis (SETH). Second, we propose a haplotype-aware chaining algorithm that runs in O(|ℋ| N log |ℋ|N) time after graph preprocessing, where N is the count of input anchors. We then establish that a chaining algorithm significantly faster than O(|ℋ|N) is impossible under SETH. As a proof-of-concept of our algorithmic solutions, we implemented the chaining algorithm in the Minichain aligner (https://github.com/at-cg/minichain). We demonstrate the advantage of the algorithm by aligning sequences sampled from human major histocompatibility complex (MHC) to a pangenome graph of 60 MHC haplotypes. The proposed algorithm offers better consistency with ground-truth recombinations when compared to a haplotype-agnostic algorithm.


Introduction
Pangenome graphs represent the maximum possible diversity that exists in the population [55].A variety of methods have been developed that use pangenome graphs for common applications including genotyping, variant calling, haplotype reconstruction, etc. (see [10] for a review).Efforts toward using pangenome graphs have been further catalyzed by the advent of long and accurate reads.The recent progress in long-read technologies and assembly algorithms enables high-quality haplotype-phased assembly of human genomes [7,41,48,62].Building a pangenome graph directly from phased assemblies instead of variant calls (VCF files) allows a more comprehensive representation of the variation [61].Accordingly, the latest methods for pangenome graph construction use phased assemblies to construct a graph [9,19,15,31,28].The input haplotypes are stored as paths in the graph.
Pangenome graphs raise fundamental questions about the trade-offs between complexity and usability.An arbitrary path in a pangenome graph corresponds to either the original haplotypes or their recombinations.The number of sequences spelled by a graph increases combinatorially with the number of variants.This issue has been addressed previously by using different techniques, e.g., by limiting the amount of variation in the graph [23,44,58], artificially simplifying complex regions [16], or restricting the alignment to either one [38,57,56] or two haplotype paths [4].A more principled approach to tackle this problem may be to leverage the correlations between two or more genetic variants, i.e.where individuals tend to have same genotype [9,18].For example, PanGenie is an alignment-free genotyping algorithm that leverages long-range haplotype information inherent in the phased genomes [9].
A pangenome graph is represented either as a directed cyclic graph or a directed acyclic graph (DAG), where each vertex is labeled by a sequence.The primary formulation for sequence-to-graph alignment problem seeks a path in the graph which spells a sequence with minimum edit distance from the query sequence [2,14,24,25,40,36,46].O(|V |+|Q||E|)-time algorithms exist for both exact and approximate pattern matching problems for graphs, where Q denotes the query sequence, E denotes the set of edges, and V denotes the set of vertices.These formulations do not consider the associations between genetic variants and may lead to alignments with spurious recombinations in variant-dense regions [44].Co-linear chaining is another common technique used in modern aligners.It is used to identify a coherent subset of anchors (short exact matches) that can be joined together to produce an alignment.The existing formulations for chaining on graphs share the same limitation of not considering the associations between genetic variants [6,32,35,45,49].Some of these chaining algorithms run in O(KN log KN ) time after graph preprocessing [6,32], where K denotes the minimum number of paths required to cover all the vertices and N denotes the count of input anchors.
In this paper, we address the above limitations by introducing 'haplotype-aware' formulations for (i) sequence-to-DAG alignment and (ii) sequence-to-DAG chaining problems.Our formulations use the haplotype path information available in modern pangenome graphs.The formulations are inspired from the classic Li-Stephens haplotype copying model [30].The Li-Stephens model is a probabilistic generative model which assumes that a sampled haplotype is an imperfect mosaic of known haplotypes.Similarly, our formulation for haplotype-aware sequence-to-DAG alignment problem optimizes the number of edits and haplotype switches simultaneously (Figure 1).We give O(|Q||E||H|) time algorithm for solving the problem, where H denotes the set of haplotypes represented in the pangenome DAG.We further prove that the existence of a significantly faster algorithm than O(|Q||E||H|) is not possible under the strong exponential time hypothesis (SETH) [59].We also formulate the haplotype-aware co-linear chaining problem.We solve it in O(|H|N log |H|N ) time, assuming that a one-time preprocessing of the DAG is done in O(|E||H|) time.To provide evidence of the near optimality of this algorithm, we show that it is impossible to solve this problem significantly faster than O(|H|N ) under SETH.Fig. 1: An example of an acyclic pangenome graph built using three haplotypes.The left figure shows the optimal alignment of the query sequence to the DAG with minimum edit distance.Accordingly, the edit distance is zero and the count of recombinations is two.Next, suppose that we use recombination penalty −5 (formally defined in Section 2).The right figure shows the new optimal alignment, where there is no recombination because the alignment path is consistent with Haplotype 2.
We implemented the proposed chaining algorithm in Minichain and evaluated it using simulated and real sequences.We built a pangenome DAG using 60 publicly-available complete major histocompatibility complex (MHC) sequences [27].We simulated query MHC haplotypes as mosaics of the 60 haplotypes.We demonstrate that introducing recombination penalty in the formulation leads to a much better consistency between the observed recombination events and the ground truth.We achieved Pearson correlation up to 0.94 between the observed recombination count and the ground-truth recombination count, whereas the correlation remained below 0.32 if recombinations are not penalized.We also tested Minichain using publicly-available PacBio and Oxford Nanopore long reads from CHM13 human genome [41,60].Here we demonstrate a better consistency between the ground-truth haplotype (CHM13) and the selected haplotype in the output.Our experiments suggest that haplotype-aware pattern matching to pangenome graphs can be useful to improve read mapping and variant calling accuracy.
2 Methods for haplotype-aware sequence-to-graph alignment

Technical background and notations
Haplotype-aware pangenome DAG.Let G(V, E, σ, H) be a DAG representing a pangenome reference.We assume that the DAG is connected and |E| ≥ |V | − 1.Let Σ denote an alphabet.Function σ : V → Σ assigns a character label to each vertex.H = {H 1 , H 2 , . . ., H |H| } denotes the set of haplotype sequences represented in the DAG.Each haplotype sequence can be identified using a path in the DAG.We define haps(v) as {i | H i includes vertex v}.We assume that haps(v) ̸ = ϕ for all v ∈ V and haps(u) ∩ haps(v) ̸ = ϕ for all edges (u, v) ∈ E. In other words, each vertex and edge must be supported by at least one haplotype.Given a query sequence Q ∈ Σ + , we will find its exact and approximate matches in the DAG.For brevity, we use [m] to denote the set {1, 2, . . ., m}, m ∈ N. Let an alignment path of length l in the DAG be denoted as an ordered sequence (w 1 , w 2 , . . ., w l ), where each w i is a pair of the form (v, h) such that v ∈ V , h ∈ haps(v), and . Thus, in our haplotype-aware setting, an alignment path specifies a path in the DAG and the indices of the selected haplotypes along the path.We say that a recombination has occurred in between w i and w i+1 if w i .h̸ = w i+1 .h.Given an alignment path P in the DAG, we use functions s(P) and r(P) to denote the sequence spelled by the path and the count of recombinations in the path, respectively.
Hardness result for the Edit Distance problem.Our hardness result (Lemma 3) will build on top of the seminal work of Backurs and Indyk [5] .The Orthogonal Vector (OV) Problem asks: given two sets of d-dimensional binary vectors A and B, |A| = |B| =: n, decide if there is a pair a ∈ A, b ∈ B such that a • b = 0.If strong exponential time hypothesis (SETH) [59] holds and d = ω(log n), there is no ϵ > 0 for which an O(n 2−ϵ poly(d)) algorithm can solve the OV problem [63].This result is frequently used to study the hardness of other problems [5,13,14,17,20].For any two sequences P 1 and P 2 , let EDIT(P 1 , P 2 ) be the edit distance between them.Let us define another distance measure between two sequences as following: x is a substring of P2 EDIT(P 1 , x) Backurs and Indyk [5] proved the hardness of computing edit distance by designing two reductions: (i) OV to PATTERN, and (ii) PATTERN to EDIT.We recall the properties of their first reduction gadget.
Lemma 1 (ref.[5]).Given an OV instance with sets A = {a 1 , . . ., a n } and B = {b 1 , . . ., b n } comprising d-dimensional binary vectors, there exist transformations of sets A and B into two sequences S 1 and S 2 respectively over the alphabet {0, 1, 2} such that the following properties hold: -Computing PATTERN(S 1 , S 2 ) determines the existence of orthogonal vectors.There exist constants We denote the injective function that generates the gadget sequence S 1 from set A as f 1 .Similarly, injective function f 2 generates sequence S 2 from set B. The exact definitions of f 1 and f 2 are available in [5].

Problem definitions
The input to each of the following problem is a DAG G(V, E, σ, H), a query Q ∈ Σ + , and a user-specified threshold k ∈ N. In each problem, we seek a match of the query in the DAG while controlling the number of recombinations.
-Problem 1 (Exact matching, limited recombinations): Determine if there is an alignment path P in DAG G such that EDIT (s(P), Q) = 0 and r(P) < k.
-Problem 2 (Approximate matching, no recombination): Determine if there is an alignment path P in DAG G such that EDIT (s(P), Q) < k and r(P) = 0.
-Problem 3 (Approximate matching, limited recombinations): Let k, c 1 , c 2 be user-specified thresholds, where k ∈ N and c 1 , c 2 ∈ R ≥0 .Here c 1 indicates cost per edit and c 2 indicates cost per recombination.Determine if there is an alignment path P in DAG G such that c 1 • EDIT (s(P), Q) + c 2 • r(P) < k.
In Problem 3, we assume a constant cost per recombination but one can also consider using different cost per locus based on known rates of recombination [43].In the above set of problems, we have omitted the case with exact matching and no recombination.One way to solve this case is by doing exact sequence matching against each haplotype independently.
Proof.We show how to solve Problem 3 using dynamic programming.The same can be extended to Problems 1 and 2 as well.Our approach is similar to the standard sequence-to-DAG alignment algorithm [25] except that we also track recombinations.We add a dummy source vertex v ϵ with an empty label.It has zero incoming edges and |V | outgoing edges to all v ∈ V .Assume haps(v ϵ ) = [|H|].Let D(i, v, h) denote the minimum cost for aligning the (possibly empty) prefix of sequence Q ending at index i against an alignment path that ends at pair (v, h), v ∈ V ∪ {v ϵ }, h ∈ haps(v).We also maintain a second table D ′′ such that D ′′ (i, v) will store min h∈haps(v) D(i, v, h).Initialize D(0, v, h) = D ′′ (0, v) = 0 for all v ∈ V ∪{v ϵ }, h ∈ haps(v).Initialize the remaining cells in table D ′′ to ∞.For i ≥ 1, v ∈ V ∪ {v ϵ }, h ∈ haps(v), update D(i, v, h) and D ′′ (i, v) using the following recursion: The recursion computes the optimal value of D(i, v, h) by considering the possibility of character deletion, insertion, match, and mismatch.The product It equals c 1 in case of a mismatch.The other expressions involve addition of c 1 to consider the possibility of an insertion or a deletion.The expressions involving the addition of c 2 evaluate the possibility of a recombination.After updating D(i, v, h) for all h ∈ haps(v), we update D ′′ (i, v) using the above equation.We compute values in table D in the increasing order of i, and then the increasing topological order of v ∈ V ∪ {v ϵ } for a fixed i.The output, i.e., the minimum alignment cost is solves the OV problem for instance (A, B).Next, we describe the construction of our gadget comprising query sequences Q 1 , . . ., Q √ n and haplotypeaware pangenome DAG G(V, E, σ, H).Assume alphabet {0, 1, 2} for defining the query sequences and the vertex labels in the DAG.Define Using Lemma 1, we know that all our haplotype sequences have uniform length, i.e., |H Next, we construct a DAG where these haplotype sequences can be represented.We build a DAG with |H i | 'layers'.Each layer comprises three vertices labelled with characters '0', '1' and '2' respectively (see Figure 2a).Each layer (except the last) is fully connected to its next adjacent layer using 3 2 = 9 directed edges.Subsequently, we identify the unique path that spells each of the √ n haplotype sequences in this DAG.Finally, we discard the vertices and edges that are not used by any haplotype sequence (Figure 2b).
Next, the above gadget will allow us to solve the OV problem by invoking the haplotype-aware sequence to graph alignment algorithm √ n times, once for each of the √ n query sequences.An alignment path with no recombination in the DAG spells a substring of a haplotype.Using Lemma 1, we make the following claim.There exist two orthogonal vectors a ∈ A and b ∈ B if and only if we have at least one query sequence q ∈ {Q 1 , . . ., Q √ n } for which an alignment path P exists such that r(P) = 0 and EDIT (s(P), q) < β 1

Methods for haplotype-aware chaining on graphs
Most alignment tools use seed-chain-extend heuristic to compute the alignments quickly [53,54].Given a set of seed matches (anchors) as input, co-linear chaining is a rigorous optimization technique to identify promising alignment regions in a reference.It identifies a coherent subset of anchors such that their coordinates are ordered on the query and the reference.The selected anchors are subsequently combined together to form an alignment.Several versions of co-linear chaining problems have been studied for aligning two sequences [1,22,34,39,11,12,42].Recent works have further studied the extension of chaining on acyclic [35,32,6,49] and cyclic pangenome graphs [3,45] but these formulations do not consider the haplotype paths.

Technical background and notations
From hereon we assume that vertices of a haplotype-aware pangenome DAG are labeled with sequences.
Sequence labeled vertices permit a more compact graph representation which will be useful in the context of chaining.It is trivial to transform a graph from sequence-labeled form to character-labeled form, and vice versa.Let function σ : .y] in the DAG (Figure 3).The weight function assigns a user-specified weight to each anchor.Anchor weights are typically set proportional to the length of the matching substrings in practice.We use R − (v) ⊆ V to denote the subset of vertices that can reach v using a path in the DAG.Vertex v is also included in R − (v).Next, we define the precedence relationship between a pair of anchors.
Definition 1 (Precedence).Given two anchors Definition 2 (Chain).A chain is an ordered sequence (s 1 , s 2 , . . ., s p ), where each s j is a pair of the form In our definition, a chain specifies the indices of the selected haplotypes alongside the anchors (Figure 3).We say that a recombination has occurred in between s j and s j+1 if s j .h̸ = s j+1 .h.Given a chain S, we use function r(S) to denote the count of recombinations in S. Let γ ∈ R ≥0 be a user-specified parameter that specifies cost per recombination.
Similar to [35], we will use path cover to facilitate efficient sparse dynamic programming on the DAG.A path cover is a set of paths in the DAG such that every vertex in V is included in at least one path.Previous chaining algorithms, e.g., in [6,32,35], identify a path cover with the minimum number of paths because the runtime of sparse dynamic programming increases with the size of the path cover [35].In the haplotype-aware setting, we will directly use the paths corresponding to the haplotype sequences as path cover.Let P 1 , P 2 , . . ., P |H| be the paths corresponding to the haplotype sequences H 1 , H 2 , . . ., H |H| respectively.{P 1 , P 2 , . . ., P |H| } is a path cover of the DAG because each vertex in the DAG is included in at least one haplotype (Section 2.1).Suppose last2reach(v, h), v ∈ V, h ∈ [|H|] denotes the last vertex in path P h that belongs to R − (v).last2reach(v, h) does not exist if no vertex in path P h can reach vertex v.We precompute last2reach(v, h) for all v ∈ V and h ∈ [|H|] in O(|E||H|) time by using the algorithm from Makinen et al. [35].This is a one-time preprocessing step that will be useful to solve the chaining problems efficiently.Next, we also use the following standard data structure in our chaining algorithm.
Lemma 4 (ref.[33]).Let n be the number of leaves in a tree, each storing a (key, value) pair.The data structure uses a balanced binary search tree.It supports update and RMQ (range maximum query) operations, defined below, in O(log n) time: update(k, val): For the leaf w with key = k, value(w) ← − max(value(w), val).
-RMQ(l, r): Return max{value(w) | l < key(w) < r}.Given n (key, value) pairs, the tree can be constructed in O(n log n) time and O(n) space.

Problem definitions
-Problem 4 (Limited recombinations): Determine an optimal chain S = s 1 , s 2 , . . ., s p that maximises the chaining score, defined as The second term γ • r(S) in the above scoring function corresponds to the penalty for haplotype switches in a chain.While scoring a chain, it is also essential to penalise the gap corresponding to the distance (in basepairs) between adjacent pair of anchors.Accordingly, we formulate a gap-sensitive version of the above problem.We assume the same definition of function gap(M -Problem 5 (Limited recombinations, gap-sensitive): Determine an optimal chain S = s 1 , s 2 , . . ., s p that maximises the chaining score, defined as

Proposed algorithms and complexity analysis
First, it is useful to discuss a simple dynamic programming solution for Problem 4. Let C(i, h) denote the optimal score of a chain ending at pair (i, h), where i ∈ [N ], h ∈ haps(M [i].v).We use C ′′ (i) to store max h∈haps(M [i].v)C(i, h).A single anchor is a valid chain of size one; therefore, we initialize C(i, h) We update C(i, h) to its optimal value using the following recursion: , thus the inequality holds.Next, consider the case when h ′ ̸ = h.Let S 1 = s 1 , s 2 , . . ., s p−1 , s p be an optimal chain that ends at pair (i, h ′ ).We consider another chain S 2 = s 1 , s 2 , . . ., s p−1 , s p ′ such that s p ′ = (i, h).Suppose the score of chain S 2 is x.Both chains S 1 and S 2 differ only by the haplotype used in their last pairs.Therefore, the difference in the chaining scores of S 1 and S 2 is at most γ (one recombination).Thus, x ≥ C ′′ (i) − γ.Also, Proof.See Algorithm 1 for an outline of our approach.We extend the sparse dynamic programming framework of Makinen et al. [35] to our haplotype-aware setting.Our path cover {P The score C(k, h ′ ) would have already been recorded in search tree T h ′ (Line 26).
All the anchors that lack a preceding anchor are optimally processed at the initialization stage (Line 2).Suppose an optimal chain ending at (i, h) is s 1 , s 2 , . . ., s p−1 , s p such that s p−1 = (k, h ′ ) and s p = (i, h).Assume that C(k, h ′ ) is already computed optimally.While calculating C(i, h), the algorithm handles the following three cases: -(Case 1) h = h ′ : In this case, we optimally compute C(i, h) by executing a range query in search tree In Lines 10-12, we add tuples in array Z to ensure that C l2r (i) is updated by using a range query in search tree T h ′ (Line 22).We also put a penalty γ due to recombination in this case.This update to C l2r (i) happens after all anchors in vertex last2reach(M [i].v, h ′ ) have been processed and search tree T h ′ has been updated because of our definition and ordering of the tuples in array Z (Lines 12, 15).Later when we update C(i, h) in Line 19, it gets its optimal value using C l2r (i)., (i, h) is an optimal chain ending at (i, h), then s 1 , s 2 , . . ., (k, h ′ ), (i, h ′ ) must be an optimal chain that ends at (i, h ′ ).Accordingly, Accordingly, C(i, h) is updated to its optimal value in Line 24.The total runtime of the algorithm is dominated by sorting of array Z. ⊓ ⊔ In the above algorithm, the optimal chain can be obtained easily by adding backtracking pointers.We next show that any algorithm to solve Problem 4 has a time complexity that is not polynomially faster than O(|H|N ) under SETH.Lemma 7.For any constant ϵ > 0, there is no algorithm that solves Problem 4 in O(|H| 1−ϵ N + |H|N 1−ϵ ) time, unless SETH is false.
13 end 14 end 15 for z ∈ Z in lexicographically ascending order based on the key (rank(v), pos, task) do Proof.The reduction from the OV problem is as follows.Suppose an OV instance is given as A = {a 1 , . . ., a n } and B = {b 1 , . . ., b n } and each set is comprised of d-dimensional binary vectors.We first define component gadgets: For two sequences s 1 and s 2 , we use s 1 ⃝ s 2 or s 1 • s 2 to denote concatenation of s 1 and s 2 .For a symbol c and an integer i, we use c i to denote symbol c repeated i times.Using the above component gadgets, we define our vector gadgets: ).The query sequence is defined as To create the vertices of graph G we start with: a vertex denoted v lef t with label 0 nd ; 2d vertices denoted v x,0 , v x,1 for 1 ≤ x ≤ d where each v x,0 has label CG(0) and each v x,1 has label CG(1); and a vertex denoted v right with label 0 nd .For edges, we add edges from v lef t to v 1,0 and from v lef t to v 1,1 ; for 1 ≤ x < d, we add edges from v x,0 to v x+1,0 and v x+1,1 , and from v x,1 to v x+1,0 and v x+1,1 ; we then add edges from v d,0 to v right and from v d,1 to v right .Next, we embed each haplotype in G by following its only possible path from v lef t to v right and deleting any vertices and edges not supported by some haplotype.See Figure 4. We call the subgraph of G excluding only v lef t and v right the center of G.
. This part of the chain has a value of 2d.Lastly, for j ′ > j, 1 ≤ x ≤ d, take anchors M right j ′ ,x , this part of the chain has value (n − j)d.This forms a valid chain and requires no recombinations.Thus, the total score is If a chain exists with score 2d + (n − 1)d, then in this chain, no recombinations are used, and anchors starting at all the anchor starting positions in Q are used.
Proof.For a given x, only one anchor to either v x,0 or v x,1 can be used while maintaining the precedence condition.Since each anchor to the center of G has weight 2, the total contribution of the center of G to the score is at most 2d.The remaining anchors in the chain have a weight of 1, and there are at most (n − 1)d of them after d anchors are used in the center.Hence, the maximum possible score is 2d + (n − 1)d.
The claim that there are no recombinations follows, as the maximum possible score using recombination becomes 2d + (n − 1)d − 1.The second claim follows since any chain using fewer than dn anchors can not achieve the maximum score of 2d + (n − 1)d.
⊓ ⊔ Lemma 10.If a chain exists with the score 2d + (n − 1)d, then in this chain, all anchors for some vector gadget in Q are to the center of G.
Proof.If anchors from two different vector gadgets in Q are to the center of G, then since the center is contributing 2d to the score (the only way a score of 2d + (n − 1)d is achieved), there must be anchors M and M ′ from two different vector gadgets that occur in adjacent vertices in the center, say with M ≺ M ′ .However, this implies that there are at least d − 1 anchor positions in Q between M and M ′ are not used in the solution, contradicting Lemma 9. ⊓ ⊔ Lemma 11.If there exists a chain with score 2d + (n − 1)d, then there exists orthogonal vectors a i and b j .
Proof.By Lemma 10, some vector gadget V G(b j ) has all of its anchors going to the center of G.By Lemma 9, no recombinations are used, and these must lie on some path for a haplotype constructed from a vector a i .By design of the component gadgets, for 1 ≤ x ≤ d if a i [x] = 1 we get anchor between v x,1 on the haplotype path for h i and V G(b j ) if and only if b j [x] = 0. We conclude that a i and b j are orthogonal.

⊓ ⊔
Observe that the total number of M lef t anchors is dn, the total number of M right anchors is dn, and the total number of M center anchors is at most 2dn, thus N = O(dn).The number of haplotypes is |H| = n.Problem 5 can be solved by extending Algorithm 1 (Supplementary Document).The extension requires a few additional terms corresponding to gap cost between adjacent anchors.This algorithm also runs in O(|H|N log |H|N ) time.The reduction presented in Lemma 7 can be easily extended for Problem 5 as well.

Implementation and benchmarking details
Implementation details.We implemented Algorithm 1 and the modifications to include gap penalty in Minichain v1.3 (https://github.com/at-cg/minichain).Our new chaining implementation replaces the haplotype-agnostic chaining implementation proposed in our previous work [6].Minichain supports an endto-end seed-chain-extend workflow to align long reads and contigs to acyclic pangenome graphs.Similar to the previous version of the software, we use the seeding and base-to-base alignment routines from Minigraph [26].We leave the implementation of our O(|Q||E||H|) time base-to-base alignment algorithm (Section 2.3) as future work because it requires further practical optimizations to reduce the runtime, e.g., using banded alignment or wavefront heuristics [37,47,64].
In Minichain, we parse the set of vertices, edges, and haplotype paths from a graphical fragment assembly (GFA) file [21].A read may be sampled from the same or the opposite orientation with respect to the graph.Accordingly, for each connected component in the input graph, we create another component of the same size with reverse-complemented vertex labels and reversed edges.We use the same path cover for the second component as the original component but the direction of each path is reversed.We compute anchors by using Minigraph's seeding method, which finds (w, k)-minimizer matches [51] between the query sequence Q and the vertex string labels σ(v) for all v ∈ V .The default values of w and k are 11 and 17 respectively.We set the weight of each anchor as 200 • k by default based on our experimental observations.Thus for k = 17, each anchor in a chain contributes 3400 to the total chain score.Minichain can report multiple best scoring chains and assign a mapping quality [29] to each chain; the criteria for this remains same as before [6].Pangenome graph construction.We used 60 publicly-available complete major histocompatibility complex (MHC) haplotype sequences [27].This set of 60 MHC haplotype sequences comprises the CHM13 MHC haplotype [41] and 59 other haplotypes from the Human Pangenome Reference Consortium [31].The length of these sequences varies from 4.8 Mbp to 5.1 Mbp.The MHC region is known to have significant polymorphism, inter-gene sequence similarity and long-range haplotype structures [8].We used Minigraph (v0.20-r559) [28] to build an acyclic pangenome graph using the 60 MHC haplotypes.The graph construction algorithm in Minigraph augments structural variants (SVs) of size ≥ 50 bp in the graph.Our MHC pangenome graph contained 984 vertices, 1, 387 edges, and 210 SVs.After computing the graph, we realigned each haplotype to the graph to obtain the corresponding haplotype paths.Simulation of query sequences.We simulated 135 MHC haplotype sequences.Each sequence was simulated as an imperfect mosaic of the 60 reference haplotypes (Figure 5).We used the following procedure for simulation: (i) Select a random haplotype path in the graph, (ii) Read the first l bases from the chosen haplotype path in the graph, (iii) If the (l + 1) th base is in vertex v and |haps(v)| ≥ 2, then switch to another randomly chosen haplotype path at vertex v, (iv) Read the next l bases, and repeat the procedure until we hit the last base of a haplotype path.After generating a sequence, we added substitutions at randomly chosen loci in the sequence.These substitutions approximately mimic true mutations and sequencing errors in real data.For each substitution rate of 0.1%, 1% and 5%, we simulated 45 sequences.In each case, three sets of 15 sequences were simulated using l = 1 Mbp, l = 2 Mbp and l = 3 Mbp respectively.Using this procedure, the number of recombination events in all the simulated sequences ranged from 1 to 4. Our pangenome graph and the query sequences are available on Zenodo (DOI: 10.5281/zenodo.10663934).

Experimental results
Comparison with haplotype-agnostic chaining using simulated data.We demonstrate that integrating recombination penalty in the sequence-to-graph chaining problem formulation leads to an improved concordance of the haplotypes used in our chains and the ground-truth.A positive finite recombination penalty in our problem formulation allows the user to limit the recombinations in an optimal chain.If we set recombination penalty γ = 0, then our algorithm is equivalent to the haplotype-agnostic chaining algorithm in [6].Recombination penalty γ = ∞ corresponds to haplotype-restricted chaining, where a chain can use only one of the reference haplotypes.We evaluated Minichain by using simulated MHC haplotype sequences and different values of recombination penalty γ = 0, 10 3 , 10 4 , 10 5 , 10 6 , ∞. Minichain computes an optimal chain as a sequence of pairs of the selected anchors and the haplotype indices.Thus, we know the sequence of haplotype indices and the number of recombinations in each chain.
We show the Pearson correlation coefficients between the observed number of recombinations and the true number of recombinations using different values of substitution rates and recombination penalties (Figure 6).The results suggest that γ = 10 4 gives the best correlation across different substitution rates.The correlation is weak when γ = 0 or γ = 10 6 .The correlation is not defined when γ = ∞ because the observed number of recombinations remained zero; thus the standard deviation of the observed values remained zero.γ = 10 3 is too low for setting recombination penalty because a single anchor in a chain contributes 3400 to the chaining score (Section 4).The algorithm may still favor an additional anchor even if it involves a haplotype switch.We evaluated the performance by using different substitution rates and recombination penalties.
Next, we compared the list of query sequences' source haplotypes to the selected haplotypes in the output chains.Suppose a query sequence was simulated using haplotypes in the order H i1 , H i2 , . . ., H in , then we recorded the "true haplotype recombination pairs" as a set of pairs {(start, H i1 ), (H i1 , H i2 ), . . ., (H in−1 , H in ) , (H in , end)}.Similarly, we recorded the sets corresponding to the "observed haplotype recombination pairs" from Minichain's output chains.We compared these two sets for each query and calculated the F1-scores (Figure 7 and Supplementary Table 1).We achieved higher F1-scores with γ = 10 4 , 10 5 compared to the haplotype-agnostic (γ = 0) and haplotype-restricted (γ = ∞) modes.These results suggest that haplotypeaware chaining and alignment algorithms can be useful to prevent spurious recombinations.Our MHC pangenome graph comprises SVs only.We expect further improvements in the accuracy of our algorithm if substitution and short indel variation are also augmented in a pangenome graph [19].This is because more variants will help to distinguish between near-identical or closely-related haplotypes.Doing this will also require improvements in Minichain's seeding and chaining implementation to allow the use of anchors that cover small bubbles in a graph.This is possible by considering a more flexible definition of anchors that can span multiple vertices [50,49,52].
The runtime of Minichain ranged from 5 to 11 minutes for aligning a simulated MHC sequence to the pangenome graph.The worst-case time complexity of the proposed haplotype-aware chaining algorithm is higher compared to the existing haplotype-agnostic chaining algorithms.This is because we require haplotype paths as the path cover, whereas the haplotype-agnostic chaining algorithms use a path cover of minimum size [6,32,35].For example, our graph with 60 haplotype paths has a minimum path cover of size 5. Thus, we observed that our haplotype-aware chaining algorithm in Minichain was about 10× slower compared to the haplotype-agnostic chaining algorithm [6].Space requirements of our algorithm are also proportional to the count of haplotype paths (e.g., to store the score table C).As a result, we observed that the haplotype-aware algorithm required about 15× more memory compared to the haplotype-agnostic algorithm [6].Evaluation using real data.We evaluated Minichain by using real long-read sequencing datasets from CHM13 human genome.Our MHC pangenome graph includes CHM13 as one of the 60 haplotypes.Since the CHM13 genome is effectively haploid, we expect all chains to be consistent with the CHM13 haplotype path in the graph.We used two PacBio HiFi datasets (SRX5633451, SRR11292121) and one Oxford Nanopore dataset (SRR23365080).We filtered the reads sampled from MHC region by mapping the reads to T2T-CHM13v2.0genome reference (GCF_009914755.1)using minimap2 v2.26 [26].Next, we discarded reads of length ≤ 1 kbp.After these filtering steps, the N50 read lengths of the three datasets were 11 kbp, 20 kbp, and 83 kbp, respectively.We considered the chaining output of a read as correct if all the anchors in the highest scoring chain overlapped with the CHM13 haplotype path in the pangenome graph.The fraction of correctly chained reads using the haplotype-agnostic objective (γ = 0) on the three datasets were 97.3%, 96.1%, and 91.1%, respectively.Using γ = 10 5 as recombination penalty, the fraction improved to 97.8%, 97.2%, and 94.2%, respectively.These results, although preliminary, suggest that the proposed haplotype-aware sequence-to-graph alignment algorithms can be useful for genotyping and variant calling using pangenome references.

Discussion
Pangenome graph has emerged as a powerful alternative representation as a reference during resequencing of a genome.These graphs can be represented in either a haplotype-agnostic or haplotype-aware manner.The former excludes the information of the reference haplotypes, whereas the latter incorporates it.Leveraging the haplotype path information during read-to-graph or contig-to-graph mapping can be useful to compute accurate alignments and avoid unlikely haplotype recombinations.This will become even more important in near future as more reference-quality human genomes are produced and included in a human pangenome graph.Prior work on haplotype-aware mapping has focused on restricting the recombination count to either zero or one [38,57,56].
In this paper, we showed how the routine chaining and alignment tasks can be rigorously formulated and solved in a haplotype-aware manner.Our formulations focus on optimizing the standard metrics, e.g., edit distance for alignment and chaining score for co-linear chaining, while considering the switches used between the reference haplotypes.Thus, the proposed algorithms will help address the combinatorial challenges as more variants are included in a graph.The other key theoretical contribution is that we proved that our alignment and chaining algorithms are near-optimal assuming SETH.Subsequently, we demonstrated the practical value of our haplotype-aware chaining algorithm.We developed a proof-of-concept implementation in Minichain (v1.3).Our experiments show that the algorithm produces alignment with fewer false recombinations and better selection of the correct haplotypes.
More work is needed, especially to establish that the proposed haplotype-aware mapping algorithms can lead to improvements in routine genomic applications.This requires a more careful engineering of Minichain to make it useful for graphs that comprise SNPs and indels.The current implementation of Minichain is restricted to acyclic pangenome graphs that contain structural variation only.It will also be interesting to investigate whether the genotype imputation accuracy can be improved using the proposed algorithms on low-pass whole-genome sequencing data.and the Intel India Research Fellowship.We used computing resources at the C-DAC National PARAM Supercomputing Facility, India and the National Energy Research Scientific Computing Center, USA.
Table 1: We show the true recombinations between haplotypes and the recombinations observed in co-linear chains using recombination penalty γ = 10 4 .We show this data for 45 query sequences that were generated using 0.1% substitution rate.'S' indicates the start and 'E' indicates the end.
is a function of n and d.Similarly, |S 2 | is a function of n and d.Both |S 1 | and |S 2 | are in O(n • poly(d)).-Time to compute sequences S 1 and S 2 from sets A and B is in O(n • poly(d)).

√ n . Fig. 2 :
Fig. 2: An illustration of the pangenome DAG used in the reduction proof of Lemma 3.
Sorting an array of size O(|H|N ) takes O(|H|N log |H|N ) time.Initializing |H| search trees takes O(|H|N log N ) time.The algorithm uses O(|H|N ) number of O(log N )-time RMQ and update operations on the search trees (Lemma 4).
Furthermore, |Q| = O(dn) and the graph G has O(d) vertices and edges, and the sum of vertex label lengths is O(dn).Hence, the reduction takes O(dn) time.Combined with Lemmas 8 and 11, it follows that an algorithm for Problem 4 running in time O(|H| 1−ϵ N + |H|N 1−ϵ ), solves the OV problem in time O(n 2−ϵ • poly(d)).This completes the proof of Lemma 7.

Fig. 5 :
Fig.5: An example to illustrate simulation of a query sequence as a mosaic of reference haplotypes.In this example, l = 19.

Fig. 6 :
Fig.6: Pearson correlation between the number of recombinations in Minichain's output chain and the true count.We evaluated the performance by using different substitution rates and recombination penalties.

Fig. 7 :
Fig.7: Box plots showing the levels of consistency between the haplotype recombination pairs in Minichain's output chain and the ground-truth using three different sets of simulated MHC sequences with substitution rates (a) 0.1%, (b) 1%, and (c) 5%.We tested using different recombination penalties.Each red dot in the plots corresponds to a query sequence.The median values are highlighted in green.
The runtime of the algorithm is dominated by the recursion.Solving the recursion takes O |Q|Σ v∈V ∪{vϵ} (|haps(v)| • (δ in (v) + 1)) time, where δ in (v) denotes the in-degree of vertex v.As a result, the time complexity of the algorithm is O(|Q||E||H|).⊓ ⊔ Proof.Problem 3 is generalized version of Problem 2. Therefore, it suffices to prove the hardness of Problem 2. We derive the reduction from OV problem.Suppose we have an OV instance with equal sized sets A and B comprising d-dimensional binary vectors.Define n = |A| = |B|.We assume w.l.o.g. that n is a perfect square.We partition set A into √ n subsets A 1 , A 2 , . . ., A √ n , each containing √ n vectors.Similarly, we partition set B into √ n subsets B 1 , B 2 , . . ., B √ n , each containing √ n vectors.Observe that solving the OV problem for all and path P h covers both the vertices M [k].v and M [i].v, then the chaining score obtained using C(k, h) is considered without recombination penalty.The third term in the expression considers the score with recombination penalty γ.This is needed when (k, h ′ ), h ′ ̸ = h precedes (i, h).After updating C(i, h) for all h ∈ haps(M [i].v), we also update C ′′ (i) using the above equation.Let function rank(v) specify the topological ordering of vertex v in the DAG.The C(i, h) values should be computed in lexicographically ascending order based on the key (rank(M[i].v),M[i].x).For a fixed i, C(i, h) values may be computed in any order.The above algorithm fills up to |H|N values in table C. Computing each value using a linear scan would require considering all N number of anchors and filtering those which satisfy the precedence criteria.Even if we ignore the time for filtering and checking precedence, the worst-case asymptotic runtime of this algorithm grows at least as fast as |H|N 2 .This makes it too slow when N is in millions, which typically happens when aligning megabase-long de novo assembled sequences to human pangenome DAGs.To address this, we propose a faster O(|H|N log |H|N ) time algorithm for solving Problem 4. First, we note the following inequality.Proof.Consider any value of i ∈ [N ] and h ∈ haps(M [i].v).Let h ′ = argmax k∈haps(M [i].v)C(i, k), breaking ties arbitrarily if necessary.Then, C ′′ Lemma 5. C(i, h) ≥ C ′′ (i) − γ for all i ∈ [N ], h ∈ haps(M [i].v).
After sorting of Z (Line 15), these tuples ensure that the search trees are updated and queried in a proper order.At the time of processing C 1 , P 2 , . . ., P |H| } contains |H| paths corresponding to the |H| haplotype sequences.We use |H| search trees T 1 , T 2 , . . ., T |H| , one per path.Each search tree T h will record C(i, h) scores specific to that haplotype.We initialise each search tree with keys {M [i].d |1 ≤ i ≤ N } and values −∞ (Line 1).In Lines 4-13, we create an array Z with O(|H|N ) tuples.