String Alignment, Dynamic Programming, & DNA

Most modern runtime environments provide an implementation of some basic functionality for working with strings. We’re used to leveraging simple tasks like splicing together strings or checking whether one string contains another, but what if we wanted to do something more complex and determine whether two strings have a shared etymology or are likely to have a similar meaning? The problem is interesting for a few reasons. First, there are some very interesting applications for predicting how closely related strings are. In terms of spoken and written language, being able to glean information about the evolutionary origin of words can serve as one building block for more complex problems like both natural language processing and translation. There are also very important uses for this kind of algorithm in related fields like DNA sequencing, and even in comparing the source code of programs. Second, the algorithms at the core of the solution to the problem provide a neat example of Dynamic Programming, which is a handy tool for any developer to have in their toolbox.

To get started, let’s consider a simple real world example. Suppose that we need to programmatically decide whether the English word “dictionary” and the Spanish word “diccionario” are related in order to draw some conclusions about the meaning of each word in absence of a… well, dictionary. We begin by noting that in practice, words tend to evolve over time through the insertion, deletion, and mutation of characters. To capitalize on this observation, we need to start by trying to find the most likely alignment between the two words by inserting gaps to make similar characters match up. The first step in finding that alignment is to establish a scoring system between characters in each word. The scoring system should be based on data and should reward similarity while penalizing differences between the strings. For our example, let’s start with a very simple scoring system where letters that match exactly are worth +2, vowels that don’t match are still worth +1 (and we’ll include the character ‘y’ with vowels), and letters that are mismatched (or matched with gaps) are worth -3. This scoring system takes advantage of a somewhat anecdotal observation that vowels tend to turn into other vowels more frequently than other mutations.

Given these words and the scoring system, our brains can actually scan the two words rather intuitively and posit a valid alignment. We can represent the alignment by showing the two words with dashes inserted for gaps, and a third row where letters indicate a direct match and plus signs indicate a non-match that was still assigned a positive score:

DICTIONAR-Y
DIC IONAR +
DICCIONARIO

Our alignment has 8 letters that match perfectly and a score of 11, and at a glance it seems like a very strong match. It also happens to be a global alignment, which I’ll touch on later. The gap in “dictionary” is also somewhat arbitrary in that the score would be identical if the gap followed the ‘y’ instead of preceding it, which means that in this case there are multiple alignments that are equally valid.

So we’ve produced an alignment, but how could we automate the alignment finding process? To take things a step further, how can we calculate the probability that the alignment is meaningful and not just coincidence? Let’s deal with the former problem first. At first blush, it looks like we could just implement a recursive top down algorithm that finds the optimal alignment by drilling down to individual characters and then building back up by adding together the maximum scoring alignments of each possible combination of substrings and gaps. If we want to align “dictionary”, we could first find the possible alignments and scores of ‘d’ and ‘i’ and then find the alignments and scores for “di” by picking the from individual character alignments and gaps, and so on. The problem with this approach is that we would end up considering each character in the first string with each character in the second string, and then we would also consider each combination of characters in the first string with each combination in the second string. The result is a less than desirable runtime complexity of O(n^2m) for strings with length n and m, respectively. Unfortunately that won’t scale well with bigger strings (like DNA, which we’ll touch on later) when we have to run a large number alignments to try to calculate the probability that the alignment is significant (which we’ll also touch on later).

An alternative is to use Dynamic Programming (DP): breaking the problem up into smaller sub-problems and solving from the bottom up. The beauty of DP is that we can optimize so that we don’t have to redo some of the simple calculations that get repeated with recursion, and we can also take advantage of nuances of the problem in the way that we setup the algorithm. For string alignment the particular nuance that we can capitalize on is the fact that subsequent characters of the string must occur in order, so at each character we can only either insert the next character or insert a gap in one string or the other.

I’ll pause briefly before explaining the DP algorithm to clarify that there are two flavors to string alignment. Global alignment requires that we use each string in it’s entirety. Local alignment requires that we find only the most aligned substring between the two strings. The most widely used global alignment algorithm is called Needleman-Wunsch, while the local equivalent is an algorithm called Smith-Waterman. Both algorithms are examples of DP, and begin by building up a two dimensional matrix of integers with dimensions of the size of each respective string plus one (because we can start each aligned string with either a character or a gap). I’m going to use the local alignment algorithm in my example, but global alignment is very similar to implement in practice.

The first step in finding a local alignment is building up the matrix of integers, where each location in the matrix represents a score as we traverse the alignment. We start by initializing the values along the top and left of the matrix to zero, since local alignment can start at any point in the string. If we were finding a global alignment we would have initialized the edges so that at each subsequent location we imposed the mismatch or gap penalty, because we would need to include the entire string in the result. Next we start with the top left corner that remains and work either right or down (it doesn’t matter as long as we stay consistent), and at each location we enter the maximum of four values: the value to the left plus the score of the top character and a gap, the value to the top plus the score of the left character and a gap, the value up and left plus the score of the characters to the top and left, or zero. If that explanation isn’t clear or the intuition isn’t there, consider what each decision is implying with respect to the alignment that we’re scoring. The first two options imply inserting a gap in one of the strings, the third implies aligning two characters, and the forth implies ending a local alignment. For our contrived example and scoring system, the matrix looks like the following:

  – D I C C I O N A R I O
– 0 0 0 0 0 0 0 0 0 0 0 0
D 0 2 0 0 0 0 0 0 0 0 0 0
I 0 0 4 1 0 2 0 0 0 0 2 0
C 0 0 1 6 3 0 0 0 0 0 0 0
T 0 0 0 3 3 0 0 0 0 0 0 0
I 0 0 2 0 0 5 2 0 0 0 2 0
O 0 0 0 0 0 2 7 4 1 0 0 4
N 0 0 0 0 0 0 4 9 6 3 0 1
A 0 0 0 0 0 0 1 611 8 5 2
R 0 0 0 0 0 0 0 3 81310 7
Y 0 0 0 0 0 0 0 0 51010 7

Apologies if my somewhat crude tables don’t format well on your device, I was lazy and copied/pasted them from program output without applying a lot of formatting. You can see that building up the matrix takes O(nm) time for strings with length n and m, so we’ve reduced our runtime complexity by an order of magnitude.

The next step is retrieving the alignment from the matrix we start with the location with the maximum value, which we can store as we build up the matrix. From that location we trace backwards and figure out which action resulted in that value by considering the locations to the top, left, and top left of the current location and calculating which value was a valid part of the calculation for the current value. In some cases there will be multiple valid alignments, as indicated by either multiple starting locations with identical scores or multiple valid locations that could have contributed to the calculation for the current location. In that case it’s important to consider the context and decide whether it makes sense to return multiple alignments, pick a single alignment arbitrarily, or create a criteria to pick between alignments (for example penalizing longer alignments). Since we’re finding a local alignment, we trace until we reach a location where the score is equal to or less than zero. If we were doing global alignment we would use the same logic, but we would start in the bottom right location and go until we reach the top left location. The same matrix below shows our trace back with the path that we’ve selected in bold:

  – D I C C I O N A R I O
– 0 0 0 0 0 0 0 0 0 0 0 0
D 0 2 0 0 0 0 0 0 0 0 0 0
I 0 0 4 1 0 2 0 0 0 0 2 0
C 0 0 1 6 3 0 0 0 0 0 0 0
T 0 0 0 3 3 0 0 0 0 0 0 0
I 0 0 2 0 0 5 2 0 0 0 2 0
O 0 0 0 0 0 2 7 4 1 0 0 4
N 0 0 0 0 0 0 4 9 6 3 0 1
A 0 0 0 0 0 0 1 611 8 5 2
R 0 0 0 0 0 0 0 3 81310 7
Y 0 0 0 0 0 0 0 0 51010 7

At this point we’ve got an alignment and a score of 13, which is pretty rad. But how can we tell whether the score is statistically significant? Clearly we can’t make hard and fast rules, because the number that the score produces is relative to a lot of factors including the size of the matrix and the alphabet that was used. One option is to create a bunch (in practice a bunch usually means a minimum of 10^4) of random permutations of one of the strings that reuses all the characters and is the same length, align them with the other string, and then see how many alignments score equal to or better than the initial match. The resulting “P-Value” can be calculated by taking the number of better matches plus one divided by the number of attempts plus one (because we include the initial match in our set). Interpretation of the value depends on the context, but in our case we may decide that a probability value of .5 may indicate that the relationship is random while a value of 10^-4 may suggest a more compelling relationship between the strings. For the given example I tried running 10^5 different alignments while randomizing a string and didn’t see a single alignment with a score that was equal or better, which isn’t a surprise because the alignment that we produced is almost a direct match. The obvious conclusion is that the words almost certainly did evolve from the same origin at some point in human history, and thus likely share meaning.

Now that we’re armed with our nifty algorithm, where else can we put it to use? Some of the more interesting applications for string alignment occur when trying to sequence and compare DNA. At the end of the day, DNA is just a long string that consists of a 4 letter alphabet of nucleotides. DNA produces RNA, which in turn also produces proteins that are made up of a chain of amino acids. Each protein is also just a string that consists of a 20 letter alphabet of amino acids. DNA (and as a result, the proteins that it produces) has evolved through a series of insertions, deletions, and mutations just like written or spoken language. DNA evolution is actually much more pure than that of written and spoken language, because it’s much more difficult for an entire section of DNA to instantly mutate than it is for a brand new word to become adopted in a culture. A gene is the DNA equivalent of a word: it’s a section of DNA that is ultimately responsible for the production of a protein that may contribute to a physical trait. By aligning either DNA or protein sequences, we can make very educated guesses about whether genes are evolutionarily linked, and whether the proteins that they produce fold and function similarly. By looking at highly conserved regions of related proteins over time and making statistical observations about the likelihood of certain amino acids and common mutations, researchers have produced scoring systems like BLOSUM that are very effective when comparing proteins.

Another interesting application for our algorithm is in examining iterations to the source code of a program, or comparing multiple implementations of a program. A program is really just a string where the alphabet is generally defined by the context free grammar of a programming language. Suppose that we want to try to figure out whether a student in a computer science class copied his assignment from a peer, or did the work on his own. In this case we can simply align each two programs, look at the probabilities that the programs are related, and apply some manual investigation and common sense investigation to any results that appear closely related.

That’s a lot of information to pack into a single post, but I hope you find it helpful and are able to think of other creative applications for the algorithms and concepts. As always, I’ll look forward to reading your comments!

Leave a Reply