Pangenome-based genome inference using integer programming

Overview

Our method, pangenome-based haplotype inference (PHI), takes as input a pangenome graph reference and short-read sequencing data from a target haploid genome. We assume that the given pangenome graph has been constructed using high-quality, haplotype-resolved genome assemblies. This graph is modeled as a DAG, in which each node is labeled with a DNA sequence, and edges represent the adjacency of these sequences along known haplotypes. Any path through this graph can represent an existing haplotype or a recombinant mosaic of segments from multiple haplotypes.

In the Li–Stephens haplotype model (Li and Stephens 2003), a new haplotype is formed by copying segments from reference haplotypes, occasionally switching between them. Inspired by this model, we aim to identify a path through the graph that best explains the observed sequencing reads by maximizing the number of shared k-mers between the reads and the sequence encoded by the chosen path. To account for recombination, our method allows switching between haplotypes in the graph but imposes a fixed penalty for each recombination event. We refer to this computational task as the path inference problem (illustrated in Fig. 1).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 1.

A schematic to illustrate the path inference problem. At the top, we show a haplotype-resolved pangenome graph containing two haplotype paths, h₁ and h₂. These paths are shown in pink and blue, respectively. An inferred path containing a single recombination is illustrated using a thick dashed line. At the bottom, we show a set of k-mers observed in the sequencing reads assuming k = 3. A k-mer is marked with a green tick if it is a substring of the sequence spelled by the inferred path AGGTTACTGAAGTT. Two k-mers (TTC and ATC) are not present in this sequence. Our optimization framework identifies an optimal inferred path with minimum cost given user-defined penalties for recombinations and for missing k-mers.

We show that solving this optimization task exactly is computationally intractable (NP-hard). To address this, we develop two efficient integer programming algorithms that leverage structural properties of the problem. A key component in our method is the construction of an expanded graph, which is a modified version of the pangenome graph in which each reference haplotype is represented as a separate path, and all haplotype paths start and end at common source and sink vertices, respectively. Within this expanded graph, the path inference problem can be reinterpreted as a constrained network flow problem: identifying a unit flow from source to sink that minimizes the combined penalties for recombinations and unmatched k-mers. In the following subsections, we describe our algorithms formally.

Notations and problem formulation

Let denote a DAG representing a haplotype-resolved pangenome reference. Function σ assigns a string label over alphabet Σ = {A, C, G, T} to each vertex. A path (u₁, u₂, …, u_n) in G spells string σ(u₁) ○ σ(u₂) ○ … ○ σ(u_n), where denotes the concatenation of strings s₁ and s₂. denotes a set of paths in G such that each of these paths spells a reference haplotype sequence used in the pangenome reference. We refer to these paths as haplotype paths. We assume that each haplotype path is described by an array; that is, h_i[1] is the first vertex in h_i, h_i[2] is the second vertex in h_i, etc. The length of a haplotype path h_i, that is, the count of vertices in h_i, is denoted as |h_i|. The set of haplotype paths covering vertex v ∈ V is denoted as haps(v). We assume that, for each edge (u, v) ∈ E, there exists a haplotype path such that u and v are consecutive vertices in h_i. In other words, each edge is supported by at least one haplotype path. Functions N⁺(v) and N⁻(v) denote the sets of out-neighbors and in-neighbors of vertex v, respectively.

Genotyping a haploid genome sample can be framed as finding a path in the pangenome graph that contains the sample's variants (Paten et al. 2017). Such a path can be modeled as a mosaic path formed by recombinations between the haplotype paths of the graph. Based on this intuition, we define an inferred path P of length n as an ordered set (a₁, a₂, …, a_n), where each a_i is a two-tuple (u, h), . For a tuple a_i, we use a_i.u to denote the respective vertex and a_i.h to denote the respective haplotype path. By using this notation, we keep track of the haplotype path being used alongside the vertex indices. In an inferred path, (a_i.u, a_i+1.u) ∈ E for all i ∈ [1, n). If a_i.h = a_i+1.h, then we also require that a_i.u and a_i+1.u must be consecutive vertices in haplotype path a_i.h. The string spelled by path P is σ(a₁.u) ○ σ(a₂.u) ○ … ○ σ(a_n.u). We denote this string as , with a slight abuse of notation.

We say a recombination, or a haplotype switch, occurs between two consecutive vertices a_i.u and a_i+1.u in P if a_i.h ≠ a_i+1.h. We use to denote the count of recombinations in P. Because an inferred path must be a path in the pangenome graph by our definition, this setup models homologous recombination events during meiosis. Our goal is to determine the best path in the pangenome graph that is consistent with the set of strings (e.g., k-mers) present in the sequencing reads. We define the problem as follows.

The intuition behind our formulation is to maximize the number of string matches along the inferred path while minimizing the number of recombinations. This approach yields an inferred path that incorporates the majority of strings from S as a substring with a finite number of recombinations, constrained by a recombination penalty c. Set S in the above formulation can be a set of either k-mers or minimizers observed in the sequencing reads. The optimization is robust to sequencing errors. Errors in reads typically introduce novel k-mers that are not present anywhere in the graph. Because such erroneous k-mers cannot be matched to any path in the graph, the associated cost is unavoidable. As a result, they do not affect the optimization outcome.

We assume a constant cost per recombination. Ideally, the model should use variable, locus-specific recombination costs based on known recombination rates (Petes 2001). However, recombination maps analogous to those used on linear reference genomes are not yet available for pangenome graphs.

Next, we consider the complexity of path inference problem. We show that the problem is NP-hard, even when restricted to binary alphabets.

We refer the reader to Supplemental Note S1 for a proof (see also Supplemental Fig. S1). This result suffices to establish the theoretical hardness of finding an optimal inferred path under the given cost function. The hardness result holds specifically for our formulation (Problem 1) and does not apply to other possible formulations for genome inference on pangenome graphs. Note that there exist many versions of Problem 1 that could be practically useful and whose computational complexity remains open, for example, when a small (e.g., polylogarithmic) bound is given on the number of allowed recombinations in the inferred path or when the set of haplotype paths in the pangenome graph is small. We leave this investigation for future research.

Construction of an expanded graph

Before developing our integer programming solutions to Problem 1, it is first helpful to define an additional graph representation, which we call as expanded graph. In pangenome graphs, multiple haplotype paths share vertices if the sequences are conserved, whereas in the expanded graph, we will split all haplotypes into separate paths (Fig. 2A,B). The expanded graph enables us to model Problem 1 as a sort of network flow problem. In particular, the inferred path will be reconstructed from a flow of value one in the expanded graph. We will assign weights to edges to account for recombination penalty. Additional constraints will be used to capture how many strings in S occur in the resulting inferred path.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 2.

Illustration of an expanded graph. (A) A pangenome graph with three haplotype paths: h₁, h₂, and h₃. (B) The corresponding expanded graph, which includes three disjoint paths, one for each haplotype path. The recombination edges are shown in purple; these edges have a weight of c. We consider only the useful recombinations (Lemma 1). The edges that are not recombination edges in the expanded graph have a weight of 0. (C) The corresponding optimized expanded graph.

In an inferred path we call a recombination a_i.h ≠ a_i+1.h useful except when a_i.u and a_i+1.u are consecutive vertices in haplotype path a_i.h. Lemma 1 allows us to only consider useful recombinations when finding an optimal solution to Problem 1.

Next, we present a definition of the expanded graph in which we will consider only the useful recombinations. For technical reasons, we preprocess each edge in E, splitting it and adding a new vertex labeled with the empty string ε. Each added vertex inherits the haplotype paths that supported the edge it was formed from. This added step is to prevent recombinations from a haplotype to itself when we build our expanded graph.

We use G_E = (V_E, E_E, σ_E) to denote the expanded graph. In G_E, vertices are string-labeled and edges are weighted. Vertex set V_E is defined as (1) where h_j[i] is a vertex in V (which may appear within multiple paths in H), and (h_j[i])^j is a newly created vertex in V_E. We refer to the ordered vertex set (h_j[1])^j · · · (h_j[|h_j|])^j as a haplotype path for h_j in G_E. Note that the superscript is used to indicate to which haplotype path the vertex is designated. This results in disjoint vertices for each haplotype path in H. The vertex set also contains a source and sink vertex, s and t, respectively (Fig. 2A,B).

We denote weighted edges in E_E as tuples of the form (start, end, weight). The weighted edge set is (2) (3) (4) (5)

Next, we give some intuition for each line (2)–(5) in the above construction of E_E.

• (2) Weight zero edges are created from s to the start of each haplotype path in G_E.

• (3) Weight zero edges are created from the end of each haplotype path in G_E to t.

• (4) Weight zero edges are created between adjacent vertices in each haplotype path. That is, in the path for h_j, an edge is created from (h_j[i])^j to (h_j[i + 1])^j.

• (5) Weight c edges are used to represent the useful recombinations described in Lemma 1. We call these recombination edges.

  We use ε to denote the empty string. The vertex labels are defined as follows: (6) (7)

• (6) The vertices in a haplotype path are labeled according to the corresponding vertex label in G. These labels will be used to identify matches.

• (7) The source and sink do not require vertex labels and are hence labeled with the empty string ϵ.

Integer linear programming solution

We assume that the maximum length of any vertex label is upper bounded by the length of any string in S; that is, . This condition can be easily enforced in the input graph by adjusting the lengths of vertex labels, for example, by splitting a vertex with a long label into two, while ensuring that the graph's topology is preserved. We assume .

The basis for our solution is to find an st-flow with a flow of one through the expanded graph G_E. Our integer programs will utilize binary decision variable x_uv for each edge. The variable x_uv will take the value one if edge (u, v) ∈ E_E is part of the solution flow and zero otherwise. Because these are binary variables, the flow will always be a path. From the solution path in G_E, it is straight forward to recover the corresponding inferred path P. We use binary decision variable z_r for each string such that z_r will take the value one if the solution flow includes a subpath from hits(r). We also use variable z_rω for each .

Letting weight(u, v) denote the weight of an edge (u, v) ∈ E_E, our integer linear programming (ILP) formulation is as follows: (10) subject to (11) (12) (13) (14) In the ILP formulation, the Equation 10 models . The summation over weight(u, v) · x_uv imposes penalty c for each recombination. This is because of the two c/2 weighted recombination edges that must traversed when the path switches between haplotype paths in G_E (Fig. 2C). In the second summation, the term (1 − z_r) adds a penalty of one to the objective for every , where . Constraint 11 enforces flow conservation, allowing a unit flow from the source vertex s to the sink vertex t, ensuring that the ILP formulation selects a single path in the expanded graph.

To explain the function of Constraint 12, termed as linear string-hit constraint, and Constraint 13, observe that in an optimal solution, whenever possible the variable z_r is set to one. This is because the term (1 − z_r) in the objective function adds a penalty of zero whenever z_r = 1. However, this is only possible when z_rω is equal to one for some ω ∈ hits(r). This, in turn, is only possible if , meaning r occurs as a substring in the inferred path. Also note that at most one z_rω variable can equal one in Constraint 13. Other variables, where and , can have a value of zero, even if , justifying the use of equality in Constraint 13.

A weakness of the proposed ILP formulation is that the number of string-hit constraints equals the total number of string matches, that is, . We design another formulation with quadratic constraints in which fewer constraints are needed.

Integer quadratic programming solution

In our integer quadratic programming (IQP) formulation, Equation 10, and Constraints 11, 13, and 14 remain unchanged from the ILP formulation. Constraints in Equation 12 are replaced by quadratic constraints defined as (15) which is the quadratic string-hit constraint. Again, because of Constraint 13 at most one z_rω variable can be one. The expression sums to one when the subpath ω is contained in the flow. In this case z_r will take the value one, and no penalty is paid in the objective. Conversely, if some of the edges for ω are not in the flow, the expression will sum to zero or less. If this is the case for each ω ∈ hits(r), then Constraint 15 can only be satisfied by setting z_r = 0 and z_rω = 0 for each ω ∈ hits(r). Because z_r = 0, a penalty is paid in the objective. The total number of quadratic string-hit constraints is . In our experiments, we observe that IQP formulation solves the problem faster, albeit while requiring more memory.

As a further improvement, we relax the variables x_uv for all (u, v) ∈ E_E to continuous values x_uv ∈ [0, 1] in Constraint 14, following Lemma 2.

If z_r = 0 for all in ϕ_cont, then ϕ can be trivially obtained as a single haplotype path in G_E without recombination penalties. In such a case, all edge variables are assigned either zero or one.

For the remaining cases, we introduce the following terms:

• – ω ∈ hits(r) is a used hit-subpath if z_rω = 1.

• – A flow between vertices u and v can be decomposed into uv-paths each assigned some positive flow and called flow subpaths.

• – ω is the first used hit-subpath if there is a flow subpath from vertex s to the first vertex of ω without passing through another used hit-subpath.

• – ω is the last used hit-subpath if there is a flow subpath from the last vertex of ω to vertex t without passing through another used hit-subpath.

• – ω and are consecutive used hit-subpaths if there is a flow subpath between them without passing through a third used hit-subpath, where and .

Now, if z_r = 1 in ϕ_cont for some , there exists a used hit-subpath. We obtain ϕ as following. The flow used to reach the first hit-subpath avoids recombination penalties by following a single haplotype path. Similarly, the flow from the end vertex on the last used hit-subpath to t avoids recombinations penalties by staying on a single haplotype path. Next, consider two consecutive used hit-subpaths ω and , with u and v as their respective end and start vertices. If u and v are on different haplotype paths, any flow subpaths between u and v must minimize the recombination penalty. The same minimum recombination cost can be achieved by replacing the potentially multiple fractional flow subpaths with a single path that incurs the same recombination penalty. We can select any flow subpath from u to v and assign its edge variables to one. Edge variables on edges used on the flow from u to v and not on this selected path are set to zero.