Write a program to compute the optimal sequence alignment of two DNA strings. This program will introduce you to the emerging field of computational biology in which computers are used to do research on biological systems. Further, you will be introduced to a powerful algorithmic design paradigm known as dynamic programming.

Biology review.  A genetic sequence is a string formed from a four-letter alphabet {Adenosine (A), Thymidine (T), Guanosine (G), Cytidine (C)} of biological macromolecules referred to together as the DNA bases. A gene is a genetic sequence that contains the information needed to construct a protein. All of your genes taken together are referred to as the human genome, a blueprint for the parts needed to construct the proteins that form your cells and, by extension, your body. Each new cell produced by your body receives a copy of the genome. This copying process, as well as natural wear and tear, introduces a small number of changes into the sequences of many genes. Among the most common changes are the substitution of one base for another and the deletion of a substring of bases; such changes are generally referred to as point mutations. As a result of these point mutations, the same gene sequenced from closely related organisms will have slight differences.

The problem.  Through your research you have found the following sequence of a gene in a previously unstudied organism.

A A C A G T T A C C

What is the function of the protein that this gene encodes? You could immediately begin a series of uninformed experiments in the lab to determine what role this gene plays. However, there is a good chance that it is a variant of a known gene in a previously studied organism. Since biologists and computer scientists have laboriously determined (and published) the genetic sequence of many organisms (including humans), you would like to leverage this information to your advantage. We'll compare the above genetic sequence with one which has already been sequenced and whose function is well understood.

T A A G G T C A

Edit-distance.  We measure the similarity of two genetic sequences by using a very popular method known as the edit distance, a concept which is also widely used in spell checking, speech recognition, plagiarism detection, file revision, and computational linguistics. We align the two sequences, but we are permitted to insert gaps in either sequence (e.g., to make them have the same length). We pay a penalty for each gap that we insert and also for each pair of characters that mismatch in the final alignment. Intuitively, these penalties model the relative likeliness of point mutations arising from deletion/insertion and substitution. We produce a numerical score according to the following simple rule, which is widely used in biological applications:

Penalty per gap	2
Penalty per mismatch	1
Penalty per match	0

As an example, two possible alignments for the genetic sequences discussed above are:

Sequence 1	A	A	C	A	G	T	T	A	C	C
Sequence 2	T	A	A	G	G	T	C	A	-	-
Penalty	1	0	1	1	0	0	1	0	2	2

Sequence 1	A	A	C	A	G	T	T	A	C	C
Sequence 2	T	A	-	A	G	G	T	-	C	A
Penalty	1	0	2	0	0	1	0	2	0	1

The first alignment has a score of 8, while the second one has a score of 7. The edit-distance is the score of the best possible alignment between the two genetic sequences over all possible alignments. In this example, the second alignment is in fact optimal, so the edit-distance between the two strings is 7. Computing the edit-distance between two strings is a nontrivial computational problem because we must find the best alignment among exponentially many possibilities. For example, if both strings are 100 characters long, then there are more than 10^75 possible alignments.

A recursive solution.  Your job is to write a program to compute the edit-distance and the optimal alignment of two genetic sequences. We can compute the edit-distance recursively by breaking up the sequence alignment problem on the two original strings x and y into many alignment problems on the suffixes of the two strings. We use the notation x[i..M] to refer to the suffix of x consisting of the characters x[i], x[i+1], ..., x[M]. Note that x[M] is the string termination character '\0'. For example, consider the two strings x = "AACAGTTACC" and y = "TAAGGTCA" of length M = 10 and N = 8, respectively. Then, x[2..M] is "CAGTTACC" and y[8..N] is the empty string.

Now, we return to our recursive solution technique. Consider the zeroth column in an optimal alignment of x[0..M] with y[0..N]. There are three possibilities:

1. The optimal alignment matches x[0] = 'A' up with y[0] = 'T'. In this case, we pay a penalty of 1 for a mismatch and still need to align the string x[1..M] with y[1..N]. What is the best way to do this? This subproblem is exactly the same as the original sequence alignment problem, except that the two inputs are each suffixes of the original inputs. We can solve this subproblem recursively.

2. The optimal alignment matches the x[0] = 'A' up with a gap. In this case, we pay a penalty of 2 for a gap and still need to align the string x[1..M] with y[0..N]. This subproblem is identical to the original sequence alignment problem, except that the first input is a suffix of the original input.

3. The optimal alignment matches the y[0] = 'T' up with a gap. In this case, we pay a penalty of 2 for a gap and still need to align the string x[0..M] with y[1..N]. This subproblem is identical to the original sequence alignment problem, except that the second input is a suffix of the original input.

A dynamic programming approach.  A direct implementation of the above recursive scheme will work, but it is spectacularly inefficient. If both input strings have N characters, then the number of recursive calls will exceed 2^N.

To overcome this performance bug, we will use dynamic programming. (Read Sedgewick 5.3 for an introduction to this technique.) Dynamic programming is a powerful algorithmic paradigm that forms the core computational engine of many programs, including BLAST (the sequence alignment program almost universally used by molecular biologist in their experimental work). The key idea of dynamic programming is to break a large computational problem up into smaller subproblems, store the answers to those smaller subproblems, and, eventually, using the stored answers to solve the original problem. This avoids recomputing the same quantity over and over again. Specifically, we maintain an array, say knownOpt[i][j], that contains a table (or cache) of the solutions to all of the subproblems that the recursive function has already solved. The recursive function computes and fills in the appropriate table entry when it is presented with a new subproblem. However, when the function is presented with a previously solved subproblem, it simply looks up the appropriate value in the table and returns it.

Finding the alignment itself.  The above procedure above indicates how to compute the value of the optimal alignment. We now describe how to find the optimal alignment itself. In order to reconstruct the optimal alignment, we maintain a character matrix, say sol[i][j], to keep track of where the minimum value for aligning x[i..M] with y[j..N] came from. For example, if the minimum came from aligning x[i] with y[j] (in which case we solved the subproblem opt(i+1, j+1), then we can record this fact by drawing an arrow from (i, j) to (i+1, j+1). We can obtain a crude ASCII picture of such an arrow by storing one of the three characters '\', '-', or '|' into sol[i][j]. We interpret the three symbols as arrows emanating from (i, j) and terminating at (i+1, j+1), (i+1, j), and (i, j+1), respectively. For the example above, we get the the following solution matrix:

       |  0  1  2  3  4  5  6  7  8
   x\y |  T  A  A  G  G  T  C  A  -
-----------------------------------
  0  A |  \  \  \  \  |  \  |  \  | 
  1  A |  \  \  \  \  |  \  |  \  | 
  2  C |  \  \  |  \  |  |  \  |  | 
  3  A |  -  \  \  \  |  |  \  \  | 
  4  G |  \  \  \  \  \  |  \  |  | 
  5  T |  \  \  \  \  \  \  \  |  | 
  6  T |  \  \  \  \  \  \  \  |  | 
  7  A |  -  \  \  \  \  \  |  \  | 
  8  C |  \  \  \  \  \  \  \  \  | 
  9  C |  -  -  -  -  -  -  \  \  | 
 10  - |  -  -  -  -  -  -  -  -  . 

Input and output.  Read the two input strings from standard input, one per line. Print the edit-distance to standard output. If both strings have length at most 40, also print out a table of optimal values, optimal choices, and an optimal alignment. Test your program using the example provided above, as well as the genomic data sets from GenBank.

Analysis.  Estimate the running time and memory usage of your program as a function of the lengths of the two input strings M and N.