Transforming match bonus into cost

Tricks with match bonus or how to fool Dijkstra’s limitations Link to heading

The reader is assumed to have basic knowledge about pairwise alignment and graph theory.

In this post we will show how to transform traditional alignment score models that include a match bonus into an equivalent cost model using only non-negative costs. We generalize the relationship in Eizenga and Paten (2022) and show that it applies to both affine and asymmetric cost models. Unlike their method, our method does not increase the base penalties by a factor $2$, leading to a lower cost alignment, and hence a potentially faster runtime of WFA.

Tl;dr: Negative match costs can be transformed into an equivalent non-negative cost model. This allows WFA to solve such instances of global alignment, and may provide up to 2x speedup over an earlier transformation.

Edit graph Link to heading

The default global pairwise alignment task can be considered as the task of finding a shortest path in the edit graph, as shown in Figure 1.

Edges in this graph correspond to the cost you need to ‘pay’ for an insertion (cost $I$), deletion (cost $D$), or mismatch (cost $X$) of a letter. If two letters are the same, you can go through a ‘match’ edge for free (cost $M=0$). Since matches are free, it is allowed to greedily match characters whenever possible without considering the insertion and deletion edges. This principle is the basis of the WFA algorithm, which can extend diagonals efficiently because it does not check other operations.

In some cases, a negative match cost $M<0$ is used, where using a match decreases the cost of the alignment, as in Figure 2. We call the positive number $B:=-M$ the match bonus.

Algorithms Link to heading

Let’s briefly recall what algorithms can be used to compute the shortest path in an edit graph. (Note: $n$ and $m$ are sequence lengths here, and $s$ is the distance between them.)

• Needleman-Wunsch: $O(nm)$ approach based on dynamic programming that works for any cost model.
• Exponential band: $O(ns)$ by doubling the searched band on each iteration. Only for non-negative costs.
• Dijkstra: $O(ns)$ when used with a bucket-queue. Only works for non-negative costs.
• WFA: A more efficient variant of Dijkstra, $O(ns)$ worst case and $O(s^2)$ expected case. Only works for non-negative edges.

As you can see, only the Needleman-Wunsch DP can handle the negative cost edges created by introducing a match bonus.¹¹

That’s because NW visits all vertices in topologically sorted order, and hence considers all possible path to each vertex. 

All other algorithms break.

Potentials Link to heading

A classic variant of Dijkstra’s algorithm that can handle negative-cost edges is Johnson’s algorithm (Johnson 1977). The idea is to give a potential $\varphi (u)$ to each vertex $u$, and change the cost of an edge from $u$ to $v$ from $E(u,v)$ to²²

Note that most sources use $\varphi$ with an opposite sign. We find our choice to be more natural though: think of the potential as potential energy (height). When going up, you pay extra for this, and can use this energy/cost reduction later when going down. 

$$E^{\prime }(u,v)=E(u,v)-\varphi (u)+\varphi (v).$$ Regardless of how the potential function $\varphi$ is chosen, the length of each path from any vertex $s$ to any vertex $t$ will change by exactly $\varphi (t)-\varphi (s)$. This means that the transformation preserves shortest paths.

It turns out $\varphi$ can always be chosen in such a way that all new edge costs $E^{\prime }$ are non-negative when the graph does not contain negative cycles. In particular, in Jonson’s algorithm $-\varphi (u)=d(s,u)$ is the length of the shortest path from a fixed vertex $s$ to $u$. This somewhat defeats our purpose though: computing shortest paths is what we are trying to do in the first place.

In our case, the edit graph is weighted, directed, without cycles, and “all roads lead to Rome” in it – any way you can take in the edit graph will end in the target point to which we are finding a shortest path. To get to this terminus we need to cross the graph both vertically and horizontally as we go from the upper left corner to the bottom right one. All paths take the same number of steps down, and the same number of steps to the right, where an insertion is one step down, a deletion is one step right, and a match and mismatch are counted as one step down and one step right.

We can increase the cost of each horizontal step by ${\mathrm{\Delta }}_D$ and the cost of each vertical step by ${\mathrm{\Delta }}_I$ (Figure 3), so that (mis)matches become ${\mathrm{\Delta }}_I+{\mathrm{\Delta }}_D$ more expensive. This way, the total cost of the alignment goes up by $n\cdot {\mathrm{\Delta }}_D+m\cdot {\mathrm{\Delta }}_I$.

Now, a simple way to get rid of negative match edges is to simply require that ${\mathrm{\Delta }}_I+{\mathrm{\Delta }}_D\ge B$.

Let $u=⟨i,j⟩$ denote the vertex in column $i$ and row $j$ (Figure 4). Then, the above is equivalent to a potential $\varphi$ given by $$\varphi (u)=\varphi ⟨i,j⟩:=i\cdot {\mathrm{\Delta }}_D+j\cdot {\mathrm{\Delta }}_I.$$

Indeed, the cost of a deletion goes up by ${\mathrm{\Delta }}_D$ $$D^{\prime }:=D-\varphi ⟨i,j⟩+\varphi ⟨i+1,j⟩=D+{\mathrm{\Delta }}_D,$$ and the cost of an insertion edge goes up by ${\mathrm{\Delta }}_I$: $$I^{\prime }:=I-\varphi ⟨i,j⟩+\varphi ⟨i,j+1⟩=I+{\mathrm{\Delta }}_I.$$ Similarly, for matches and mismatches, the cost goes up by ${\mathrm{\Delta }}_I+{\mathrm{\Delta }}_D$,

$$\begin{array}{}\text{(1)}& X^{\prime }& :=X+{\mathrm{\Delta }}_D+{\mathrm{\Delta }}_I,\text{(2)}& M^{\prime }& :=M+{\mathrm{\Delta }}_D+{\mathrm{\Delta }}_I.\end{array}$$

We would like all edges to have a non-negative cost, so our choice of ${\mathrm{\Delta }}_D$ and ${\mathrm{\Delta }}_I$ must satisfy the following constraints:

$$\begin{array}{}\text{(3)}& {\mathrm{\Delta }}_D& \ge -D,\text{(4)}& {\mathrm{\Delta }}_I& \ge -I,\text{(5)}& {\mathrm{\Delta }}_D+{\mathrm{\Delta }}_I& \ge -M=B\ge 0.\end{array}$$

Multiple variants Link to heading

The above equations give us some flexibility in choosing ${\mathrm{\Delta }}_I$ and ${\mathrm{\Delta }}_D$. We will restrict ourselves to the case where $M^{\prime }=M+{\mathrm{\Delta }}_D+{\mathrm{\Delta }}_I=0$, i.e. where exact matches are free.

There are a few natural choices of ${\mathrm{\Delta }}_D$ and ${\mathrm{\Delta }}_I$ that we cover in the table below.

The symmetric option in the first row is the most natural choice, and roughly corresponds to the transformation suggested in Eizenga and Paten (2022). It differs in that all costs are divided by $2$ and the half-integer costs in the odd $B$ case are distributed between $I^{\prime }$ and $D^{\prime }$.

The bottom two rows are even applicable when matches are already free ($B=0$), and transfer the cost of horizontal steps to vertical steps or vice-versa by increasing $D$ and decreasing $I$ by the same amount.

Some notes on algorithms Link to heading

WFA Link to heading

For WFA this cost transformation is a life-saver because with $M^{\prime }=0$ it allows for “greedy matching” again—the core of the algorithm that enables fast extending diagonals. This is important, because one of the main WFA data structures is a wavefront, which covers all states that can be reached with specific cost. If the cost could go down, it could disrupt previous wavefronts, which would break the logic of the algorithm. The drawback of the cost transformation is that it increases the cost of the optimal alignment, and hence increases the number of wavefronts, leading to a slower WFA execution.

Thus, we expect that the lower the potential of the target state, the faster the WFA algorithm runs.

A* Link to heading

The A* algorithm by itself can not handle negative costs. However, a heuristic function $h$ can be modified to achieve exactly the same result as the potentials introduced above by using

$$\begin{array}{}\text{(6)}& h^{\prime }(u):=h(u)+\varphi (u).\end{array}$$

This works in our case because it ensures that $f(u)=g(u)+h^{\prime }(u)$ can never decrease when taking an edge.

Extending to different cost models Link to heading

Affine costs Link to heading

The potentials defined above naturally extend to affine costs. Each state in an affine layer naturally corresponds to a state $⟨i,j⟩$ in the main layer, and can use the corresponding potential.

This means that the delete-extend cost increases with ${\mathrm{\Delta }}_D$, and similarly the insert-extend cost increases with ${\mathrm{\Delta }}_I$, while gap-open costs remain the same.

Substitution matrices Link to heading

The BLOSUM matrix (Figure 5) specifies a match score for each pair of amino acids, with some entries being positive (indicating similarity) and some being negative, so we are maximizing the score. This can be transformed into a cost model by simply negating all scores, which allows us the previous techniques.

In general, let $B$ be the maximum score assigned to any pair of letters in the substitution matrix $S$,³³

I’m using $S$ instead of $M$ to indicate that this is a score instead of a cost. 

i.e. the maximum entry in the matrix. Again we choose ${\mathrm{\Delta }}_D$ and ${\mathrm{\Delta }}_I$ such that ${\mathrm{\Delta }}_D+{\mathrm{\Delta }}_I\ge B$. Then, we simply replace each matrix element $S_{xy}$ by $S_{xy}^{\prime }=S_{xy}-{\mathrm{\Delta }}_D-{\mathrm{\Delta }}_I$, ensuring that all the scores are at most $0$.

But not local alignment Link to heading

We will stress here that this idea does not work for local alignment. The reason it works for global alignment is that the cost of each path is increased by the same amount. Two local alignments can have completely different lengths and thus span a different number of rows/columns of the table. That means that the cost increase $\varphi (end)-\varphi (start)$ of these two local alignments is different, and hence their final scores are not directly comparable.

Evaluations Link to heading

Now let’s see how WFA performs when using different transformation variants.

We will use the following cost model:

• $I=1$ - insertion cost
• $D=1$ - deletion cost
• $X=1$ - substitution cost
• $B=2$ - match bonus

Unequal string length Link to heading

Firstly, let’s make an experiment on sequences where the first sequence is longer than the second, i.e. $|A|=n>m=|B|$.

Given the match bonus, we can take different values for ${\mathrm{\Delta }}_I$ and ${\mathrm{\Delta }}_D$. According to our hypothesis of keeping the cost increase minimal, the optimal option should be ${\mathrm{\Delta }}_I=2$ and ${\mathrm{\Delta }}_D=0$, since ${\mathrm{\Delta }}_I$ is multiplied by $m=min(m,n)$. Let’s check that!⁴⁴

The author’s laptop was made in times of the Second Punic War, so he decided not to take very long sequences to save his working station. 

^, 55

Also, the code is not optimized in the first place. Look at relative timings only. 

[TODO: Also compare to the transformation in Eizenga and Paten (2022).]

Note: the last row in the table depictes experiment on cost model with match bonus introduced in Eizenga and Paten (2022), so ${\mathrm{\Delta }}_D+{\mathrm{\Delta }}_I$ is not required to be equal to $b$ in this particular case.

The numbers are eloquent enough, so let’s look at the pixels. As you can see the second sequence has a suffix of the first one deleted.

We can clearly see that increasing ${\mathrm{\Delta }}_I$ (the vertical cost) instead of ${\mathrm{\Delta }}_D$ (the horizontal cost) yields less downward expansion, and fewer visited states in total.

Equal string lengths Link to heading

Let’s see how things work out in case of equal length strings, where the total score will be the same independent of how we choose ${\mathrm{\Delta }}_D$ and ${\mathrm{\Delta }}_I$.

Note: the last row in the table depictes experiment on cost model with match bonus introduced in Eizenga and Paten (2022), so ${\mathrm{\Delta }}_D+{\mathrm{\Delta }}_I$ is not required to be equal to $b$ in this particular case.

In this case, we can see in the experiments that it is preferred to split the match bonus equally between ${\mathrm{\Delta }}_D$ and ${\mathrm{\Delta }}_I$. Pictures for this case are so:

The middle picture clearly visits fewer states, and seems to benefit from being pushed towards the end, without needless exploration to the sides.

Conclusion Link to heading

We have shown that alignment score/cost models that include a match bonus can be transformed into an equivalent cost model using only non-negative costs when doing global alignment. This transformation even works for affine and asymmetric costs. It avoids doubling costs (to preserve integral values) by distributing an odd match bonus unevenly over the horizontal and vertical steps in the edit graph.

To summarize: when the cost model is given by match cost $M=-B$, mismatch cost $X$, insertion cost $I$, and deletion cost $D$, the equivalent cost model for global alignment with non-negative costs is given by:

$$\begin{array}{}\text{(7)}& M^{\prime }& =0,\text{(8)}& X^{\prime }& =X+B,\text{(9)}& D^{\prime }& =D+\lfloor B/2\rfloor ,\text{(10)}& I^{\prime }& =I+\lceil B/2\rceil .\end{array}$$

References Link to heading