Getting an intuition about dynamic programming in general (and alignment in particular)

In this module, you have seen several examples of problems that are solved through dynamic programming: making change, finding a path through a graph, finding the longest common subsequence, and aligning (DNA or protein) strings.

A commonality that all of these algorithms share is that they can be broken up into sub-problems that are related to each other in a systematic way.

For example, when making change for $10, a possible sub-problem is making change for a smaller dollar amount, say $5. If you know how to make change for $5, then you just need to figure out all the ways of getting from $5 to $10, making the problem easier. Similarly, in the graph problem, where you are trying to get from A to B following the shortest or longest path, the shortest/longest path from A to some node C between A and B constitutes one of the sub-problems you can consider.

The key question is coming up with sub-problems that have a meaningful relationship between each other, in the sense that they can be combined to give you the answer to the problem you started with. In the graph case, the shortest path between A and B that passes through C must follow the shortest path between A and C. If there are multiple choices for C, you can just try each of them separately and select the one that gives you the best path between A and B (the shortest path between A and B goes through some node C, so you just iterate through all of them).

This may still be a bit confusing, but the key point is that setting up a problem to be solved as a dynamic programming algorithm requires you to think about what the sub-problems are, and how they are related to each other. For any sub-problem, you essentially have to think about all the other sub-problems that can be used to solve it in a simple way. Following the chain of sub-problems is what gives you, in the end, the solution to the problem you started with.

To be concrete,let's again start with the question of the best way of making change for some dollar amount, say $48.73 . While the answer is not obvious, you can focus on the simplest transition between sub-problems where you focus on just one coin or bill. The best way to make change for $48.73 could be to use a penny and the best way of making change for $48.72; or a dime and the $48.63 sub-problem, or a $20 bill and the $28.73 sub-problem. Then you can proceed in the same way for each of the sub-problems identified, building a virtual directed acyclic graph (DAG) that chains together increasingly simpler problems into a solution for the problem you started with.

Of course, you have to stop somewhere, so you need to define a base case - a problem or set of problems for which the answer is obvious (say, making change from a penny).

How does all of this lead into the alignment question? More precisely, assume we are given strings S1 and S2 and we are asked to find the maximum scoring global alignment between them.

The sub-problem we typically use asks how the ith and jth prefixes of the two strings align to each other, i.e., solving the alignment problem for S1[0..i] and S2[0..j]. The graph that relates different problems to each other, or its matrix representation, contain a node for each of the possible choices of i and j, or MN nodes if len(S1) = M and len(S2) = N.

So far, we've effectively defined the sub-problem structure we will work with. Now, how about the relationship between the different sub-problems? Let's think for a moment about how the optimal alignment for the two strings may look like. It could end with the last two characters of the string being aligned to each other:

....X
....Y

or one or the other of the character aligned to a gap character:

...X  ...-
...-  ...Y

There aren't any other choices, as it makes no sense to have two gap characters aligned to each other.

What does this mean for the sub-problem structure?

If X and Y are aligned to each other, then we can say that the best alignment between S1 and S2 can be derived from the best alignment of S1[0..N-1] and S2[0..M-1] plus the score/penalty of aligning X and Y to each other (which you get from the BLOSUM matrix in the case of the programming assignment).

What if X is aligned to a gap character? Then the best alignment of S1 and S2 is derived from the best alignment of S1[0..N-1] to the full S2, plus the cost of aligning X to a gap character (or deleting X).

Making this more abstract, what all of this is saying is that sub-problem (M, N) can be solved by looking at sub-problems (M-1, N-1), (M-1, N), and (M, N-1). In the graph model, the node (M,N) will be connected by edges to nodes (M-1, N-1), (M-1, N), and (M, N-1), and this ultimately leads to the graphs presented in the book. The scores of the edges, or of the transitions between nodes in the matrix, is defined by the score of aligning two characters to each other or of deleting one or the other character. The final score is simply the maximum over the three possible choices.

We are just left to define one more item - the base case. While you could work out the alignment of one letter against one, two, three, etc. letters, it is a lot easier to think of deleting one or more letters from a string, essentially assuming that one string is entirely empty and you align increasingly longer prefixes of the other string to this empty string. In the matrix formulation, this corresponds to the first row and column of the matrix, or nodes (i, 0) and (0, j) in the graph formulation. As each horizontal or vertical edge has the cost of -sigma in the graph, the value of the matrix cells (i, 0) and (0, j) is exactly -i sigma and -j sigma, respectively.

To summarize - defining the dynamic programming algorithm requires 3 parts:

Defining the sub-problem (what does a cell in the matrix or a node in the graph represent)
Defining the recurrence rules - how you determine the solution for a sub-problem by using previously-completed sub-problems. This also defines the connections between nodes in the graph.
The base case - the sub-problems for which the answer is obvious.

To compute the score of the optimal solution, the algorithm simply starts with the base-case then iteratively computes the values of more sub-problems using the recurrence rules.

Ok, there's actually a fourth requirement - and that's the starting point for the backtracking process, or the sub-problem that represents the solution to the problem you've started with. For global alignment, that is the sub-problem (M, N), or the bottom-right corner of the graph/matrix.

Variants sequence alignment

Now that the basics are in place, we can discuss how small changes to the 4 elements described above can alter the nature of the solutions.

Let's take the 4th element - the location of the answer. If we look for a global alignment, then the answer is the alignment of the two full strings, hence the bottom right corner of the matrix.

If we look for a local alignment, we are looking for substrings within the two strings that align to each other optimally. In other words, we are no longer looking for the full two strings, but we are looking for some sub-problem (I, J), or the alignment of the Ith and Jth prefixes of the two strings. Which particular sub-problem? Since we want the highest scoring alignment, we want the I and J that refer to the maximum value in the matrix.

What if we want to find an alignment of string 2 within string 1? In other words, we have to use the whole length of string 2 but it can align somewhere within string 1? Then we are looking for a sub-problem (I, M), where I is again selected so that the choice of (I, M) gives the highest score.

And so on.

Similarly, by changing the recurrence equation, we can also change the nature of the alignment. Assume we allow a "free ride" in addition to the three choices we had before, i.e., we allow each node in the graph to be connected to the (0, 0) node with an edge of weight 0. In that case, the alignment does not need to extend to the beginning of each string, rather you can stop in the middle, yielding a local alignment.

Note that this effectively sets all base cases to 0, as you can always get a better score by jumping from (0, 0) to (0, i) or (i, 0) for a cost of 0 rather than following a series of edges with penalty -sigma.

Finally, by playing with the base case we can allow alignments to start in the middle of one or the other string - if the first row is all zeroes, that means we can start an alignment at any point within the first string without penalty. Keeping the global alignment penalties (-sigma, -2 sigma, -3 sigma, etc.) on the first column forces the alignment to "use up" the other string. By coupling this trick with the idea of finding the answer somewhere else than the bottom right corner allows you to come up with a range of different alignments:

Global alignment

First column and row have gap penalties
Recurrence equation disallows free rides
Answer found in bottom right corner

Local alignment

First row and column have 0s
Recurrence equation allows free rides
Answer found in the middle of matrix

Fitting alignment (one string matches within the other)

Assume the longer string is at the top

The first row has 0s, the first column has gap penalties
Recurrence equation disallows free rides
Answer found at highest value on the bottom row

And so on.

The trick with backtracking

The reason backtracking is necessary is because the "forward" part of the dynamic programming approach effectively explores a large number of possible solutions and you don't know which solution you will select until you get the end of the calculation. The best way to understand this is in the context of local alignment. There can be multiple local alignments between two strings, and you have to fill out the full matrix before finding out which one has the highest score. Once you do, you have to retrace the steps that took you to that highest value, essentially reverse the calculation you did as part of the recurrence step.

The best way to think about the recurrence is to ask at each position in the matrix "where did the value at this position come from?". You may have ties (e.g., the score is the same if you delete a character in either string), but it doesn't matter which one you choose - they will all yield optimal alignments. Note that there may be an exponential number of alignments with the same score, but that's a story for another day.

Dealing with gaps

Assume for now that we want to score a string of k gaps in a row with a function g(k), where the function could really be anything (e.g., k²+2k-5). How would you structure the global alignment algorithm to handle such a scoring function? More precisely, which one of the 4 elements of the algorithm will you have to change?

The sub-problems are the same - alignments of prefixes of the two strings. The initialization is the same - you have to penalize gaps so that the strings align end to end. And the answer is found in the bottom right corner. Thus, the only thing that may need to be changed is the recurrence equation.

How do you change it? The question is what are the possible ways to get to a particular solution. Before you had:

diagonal - align to characters to each other
vertical - align one character to a gap
horizontal - align the other character to a gap

But now, we want to penalize different numbers of gaps differently. So we could have multiple vertical or horizontal "moves" that jump over 2, 3, ... k gaps. Each of the jumps would be penalized by g(2), g(3), ... g(k).

Now the question is how to pick the correct link. Just like the other example, you simply take the max over all possible choices, yielding something like:

max(diagonal, vertical(1), vertical(2), vertical(3),... vertical(l), horizontal(1), horizontal(2), ..., horizontal(m)).

adding an O(n) loop inside the two nested loops. That's what yields an O(n^3) algorithm that can compute the alignment of two strings using an arbitrary gap penalty function.

For the special case where g(k) = ak + b you have the affine gap penalties and you can solve those in O(n^2) time by decomposing the dynamic programming table into 3 - two of which model consecutive gaps on the horizontal or vertical directions.

Finding the middle in the linear space alignment

First, I hope it is clear that if all you want to compute is the score of the alignment, it is trivial to run the global alignment algorithm in linear space. At any point in the algorithm you only need information from two rows (or columns, or diagonals, depending on how you traverse the matrix), and you use the information in one of the rows to fill in the information in the next. Thus, you can simply alternate between the rows until you reach the bottom right corner of the matrix.

Doing the backtracking, however, requires you to have the full matrix as you cannot "back-compute" the table starting from the bottom right corner. That's where the recursive trick described in the book comes in.

There are two key elements to the effectiveness of this approach. The first is that at each iteration you do half the amount of work that you did in the prior iteration. Since $LaTeX: \sum_{i=0..N}\frac{1}{2^i}\:<\:2$ , the overall runtime remains O(N²), just the constant is twice as large.

The second element is the construction of the alignment itself. The way the algorithm works is to break the problem into two sub-problems, and construct a small piece of the alignment in between the two sub-problems. Here is where finding a middle column or a middle row comes in, or, more precisely, finding the place where the optimal alignment crosses this middle column or row.

To better understand how this works, you need to notice that the optimal alignment must cross this middle column at some point. The question is whether one can figure out where this crossing point is, without actually computing the alignment.

First, note that the score of the alignment is decomposable - if I break up the optimal alignment in the middle, the sum of the scores for either side of the alignment add up to the full score. This is because, in the end, the score is simply the sum of the scores of the aligned characters.

Since the optimal alignment has the best score, then the sum of the scores on either side of the middle column must also be the highest, and this is the key to finding the place where the optimal alignment crosses the middle column.

On the left side of the middle column, you simply have the traditional dynamic programming alignment, but now performed between half of one string and the entirety of the other string. Note, however, that we only care about a prefix of the second string aligning to the half of the first, as the suffix will align to the second half. Thus, we want to find the best score along the middle column, as we've discussed above. Since we can compute the score using just linear space as just described, this is trivial - we fill in column by column until we get the scores for the middle column.

Now on the right of the middle column things are a bit more complicated. It appears that we cannot run the dynamic programming going forward as we need the information from the left side of the alignment. But, if you think about it - we are simply aligning the second half of one string to the whole of the second string. The score and the alignment won't change if we flip both strings around (essentially reading both from the end to the front), and that effectively translates into filling the dynamic programming table from the bottom right corner moving towards the left. Again, we can use this approach to compute the scores along the middle column, now coming from the right.

Now, on the middle column, we'll have a set of scores coming from the left - the partial alignments of the first half of the first string - and a set of scores coming from the right - the partial alignments of the second half of the first string. The optimal alignment is one that maximizes the sum of these scores, meaning that the optimal alignment crosses the middle column at the place where the sum of the two scores is the highest.

There's a bit more that we need to take care of - specifically accounting for the fact that the optimal alignment may entail more than one cell on the middle column - that's something you can figure out by following the backtracking on the left and the right until you fall off the middle column.

Then you simply repeat the process on the left and on the right, but focusing on just the top or bottom of the corresponding side of the matrix, as defined by the point where the alignment crossed the column. Iteration by iteration you will then fill in additional bits of the optimal alignment.

Since the alignment itself only requires linear space (it's two strings with possible gaps inserted in them), and the score computation can be done in linear space, the total space usage is linear. As we already discussed, the many iterations do not blow up run time too much since at each iteration we halve the amount of work that needs to be done.