Dynamic
Programming Algorithms
Overview
Dynamic programming can be employed as a means to
align two sequences. The premise behind
the algorithm is that an optimal alignment can be constructed by adding two
aligned characters to a previously obtained optimal alignment. This insight means that it is not necessary
to search all possible alignments in order to obtain the optimal one (given two
N-lengh sequences, this amounts to time proportional to N4N);
rather, you must sequentially add characters to an optimal alignment that you
have already constructed.
Optimal
alignments are made using a sum matrix where any given element is given a score
based on whether or not it represents a match, plus the score of the optimal
alignment to which it could be added, minus a gap penalty incurred for adding
spaces in the sequence. At the end, a
traceback is performed where you start with the highest scoring element in the
final column or row, and follow a path back to the end of the sequence such
that every element contributing to the highest score is aligned. For two N-length sequences, this amounts to
searching an average of N values for each of the N2 spots in the
matrix. This implies that the time will
be proportional to N3 – a significant improvement over searching all
possible alignments. However, given any
cell, you have already calculated the maximum score for two cells adjacent to
it, so only two comparisons need to be made for each matrix element, reducing
the necessary time to O(N2).
Implementation
This project includes to
programs: one which uses an algorithm that should work in N3
(“align”) time and one that should work in N2 (“alignfast”) time by
caching previous values.
Consider two N length strings (for simplicity), A and
B. Both implementations create an N by
N matrix of structures M, where M(i,j) is an alignment between the ith element
of A and the jth element of B.
Each structure contains the score, its i and j coordinates (for easy
reference during traceback), and a pointer to the highest scoring element in
M(i+1, j+1) to M(i+1, N) or M(i+1, j+1) to M(N, j+1). In the event of a tie, the element that introduces fewer gaps is
chosen. If there are two equivalent
paths, the “high road” (M(i+1, j+1) to M(i+1, N) ) is taken.
After the sum matrix is constructed, the maximum value
in the top row or column closest to the origin (again, picking the top row in a
tie) is chosen as the starting point.
To complete the traceback, the pointer in the current high scoring structure
is taken as the next alignment.
The only difference in the two
implementations is that align searches each element along M(i+1, j+1) to M(i+1,
N) or M(i+1, j+1) to M(N, j+1) for every matrix element; alignfast stores the
previous high row and column scores, and then compares the stored scores with
M(i+1,j+1) to determine the final score of the element. The high scores are stored in another matrix
of structures of length N-1 by N-1 (since you don’t need to compare anything
for edges) as the sum matrix is being computed.
Documentation
Included in the tar file are six c files. Utility.h and Utilityfast.h are header files
for align and align fast, respectively, which contain structure definitions and
function prototypes. Utility.c and
utilityfast.c contain definitions for the functions. Align.c and alignfast.c contain the main source code. The tar also contains a Makefile so that
“make”, invoked with no arguments produces the executables “align” and “alignfast”.
Both programs work identically
with the exception of time. Two strings
are taken from stdin, with a single carriage return separating the first from
the second, and either a carriage return or “end-of-file” sign (which can be
simulated by typing Ctrl-D at stdin) after the second. Previously stored strings in a text file can
be supplied to the programs by an input redirect with the following
syntax: align < your-text-file. (alignfast could also be
used)
Output will consist of the two aligned strings (one
above the other), with “|” linking aligned characters that are identical, and
“*” linking aligned characters that don’t match. Spaces in the sequences are represents by “-“. After the aligned strings, the total score
(representing the total number of matches with no gap penalties) is printed.
The programs can be timed using: time align <
your-text-file. (type man time for details)
Outcome
In order to determine the length dependence of the
programs, we can look at the ratio of the time it takes to run a program with
strings of length N to the time for the program with strings of length 2N. If the program runs in O(N3),
this ratio should be 23 = 8, since proportionality constants cancel
in the ratio. For a program that takes
O(N2) time, the ratio should be 22 = 4.
I created several test files (containing and equal
number of “a”s in each string) to examine the different algorithms. The results were as follows:
The ratio for the align implementation is close to 8
as expected. Any discrepancy is likely
due to quadratic and linear factors in the time function. In the limit as N approaches infinity, we
expect these other factors to disappear (since N3+N2+N is
approximately equal to N3 for large N); as N increases in the table
above the factor does indeed get closer to 8.
The ratio for the alignfast
implementation is close to 4, also as expected. Again, as N increases the ratio seems to get closer to 4; the
anomalous ratio for the two smallest strings is likely a problem of precision
since the time function only reports to hundredths of seconds.
Finally, since these programs are
being run on a remote linux server, the percentage of CPU time being given to
your particular terminal will alter the time they take to execute. For all of the above measurements, the time
function returned a cpu percentage of about 100% +- .5%, so this consideration
is probably unimportant.
2
sample runs:
cobra:~/dynamic% make
gcc -c align.c
gcc -c utility.c
gcc -o align align.o utility.o
gcc -c alignfast.c
gcc -c utilityfast.c
gcc -o alignfast alignfast.o utilityfast.o
cobra:~/dynamic% align
abcnyrqclcrpm
aycynrckcrbp
abcny-rqclcr-pm
|*| | | |*|| |
ayc-ynr-ckcrbp-
Score: 8
cobra:~/dynamic% time alignfast < test
aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa
aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa
Score: 158
0.010u 0.010s 0:00.03 66.6% 0+0k 0+0io 87pf+0w
cobra:~/dynamic%