Sequence comparison to determine how similar or different two sequences are is a fundamental operation of computational biology. Perhaps two laboratories want to compare their work on the same gene; or a biologist might want to determine similarities in gene sequences from different organisms, such as yeast and mouse. Programs that recreate a large DNA sequence from fragments use sequence comparison. As another example, some biologists are interested in searching databases for sequences that are locally similar to a given sequence in an attempt to discover gene function. (Refer to the "Genomic Data" module for background material for the current module.)

• +1 for a match
• -1 for a mismatch
• -2 for a space in one of the corresponding positions

For example, with indicated column scores, the alignment

Because we obtain the same column scores regardless of which sequence is in the first row, the order in which we write the sequences is irrelevant for determining an alignment score.  Thus, the similarity of two sequences does not depend on which we write first.  For example, the similarity of ATGAC and ACGC is the same as the similarity of ACGC and ATCAC.

• Diagonal: Starting with the best alignment of AT and A, prefixes of s and t, respectively, consider the impact of the match (+1 value), or in this case, the mismatch (-1 value) of G onto the prefix of s and C onto the prefix of t
 
• Above: Starting with the best alignment of AT and AC, prefixes of s and t, respectively, consider the impact (-2 value) of adding G to the prefix of s and blank to the prefix of t
 
• Left: Starting with the best alignment of ATG and A, prefixes of s and t, respectively, consider the impact (-2 value) of adding blank to the prefix of s and C to the prefix of t

To summarize, a[3][2] is the maximum of the following values:

• Diagonal value (a[2][1]) plus -1 for mismatch of G and C (or diagonal value plus +1 if the letters had matched)

• Above value (a[2][2]) plus (-2)

• Left value (a[3][1]) plus (-2)

• A diagonal line segment indicates a match or mismatch of the pair of bases that head the row and column. For example, the diagonal between a[4][3] and a[5][4] indicates that it is best to attach the Cs onto the prefixes aligned with similarity value a[4][3]. The diagonal between a[1][1] and a[2][2] says to align T and C onto the prefixes aligned with similarity value a[1][1].
• A vertical line segment matches the row heading in s with a space. For example, the horizontal line segment between a[3][3] and a[4][3] indicates to attach A to the prefix of s and a space to the prefix of t in the developing alignment.
• A horizontal line segment matches a space with the column heading in t.

Figure 5 displays the similarity matrix of Figure 4 with the path from a[5][4] to a[0][0] indicating the best alignment.  Using Figure 5, Figure 6 gives the optimum alignment of ATGAC and ACGC. Sometimes, more than one alignment yields the same similarity.

Figure 5. From Figure 4, path indicating best alignment of ATGAC and ACGC in bold

Figure 6. Optimal alignment of ATGAC and ACGC obtained from Figure 5

score(s, t, match Value, mismatch Value, spacePenalty, a)
	Procedure to determine array of similarity values for two sequences
		of bases
Pre:
	s and t are non-empty sequences of bases.
	match Value is the numeric value for two paired bases in an alignment
    		agreeing.
	mismatch Value is the numeric value for two paired bases in an alignment
        	disagreeing.
	spacePenalty is the numeric value for a base being paired with a space
        	in an alignment.
	match Value > mismatch Value
	match Value > spacePenalty
	a is an m-by-n array, where m is the length of s plus 1 and n is the length of t plus 1.

Post:
	a is an array of similarity values for prefixes of s and t. The value
        	in the bottom, right corner is the similarity of s and t.
Algorithm:
	m <-- length of s
	n <-- length of t
	for i from 0 through m, do the following:
		a[i][0] <-- i * spacePenalty
    	for j from 0 through n, do the following:
		a[0][j] <-- j * spacePenalty
    	for i from 1 through m, do the following:
		for j from 1 through n, do the following:
			if (s[i - 1] equals  t[j - 1])
				diagonalValue <-- match Value
			else
				diagonalValue <-- mismatch Value
			a[i][j] <-- max(a[i - 1][j - 1] + diagonalValue,
     				     a[i - 1][j]		+ spacePenalty,
     				     a[i][j - 1]		+ spacePenalty)

align(s, t, a, match Value, mismatch Value, spacePenalty, sAlignment, 
			tAlignment)
    Algorithm to determine optimal alignment of two sequences of bases

Pre:
    s and t are sequences of bases.
    a is the array of similarity values for prefixes of s and t, as completed by score.
    spacePenalty is the numeric value for a base being paired with a space in an alignment.
Post:
    sAlignment and tAlignment are sequences in the optimal alignment for
         sequences s and t, respectively, with a dash indicating a space.
Algorithm:
	sLength <-- length of s
	tLength <-- length of t
	recAlign(s, t, a, match Value, mismatch Value, spacePenalty, sAlignment,
        	tAlignment, sLength, tLength)

recAlign(s, t, a, match Value, mismatch Value, spacePenalty, sAlignment, 
       tAlignment, i, j)
    Algorithm to determine optimal alignment of two prefix sequences of
         bases
Pre:
    s and t are sequences of bases.
    a is the array of similarity values for prefixes of s and t, as completed by score.
    i is an index that is less than or equal to the length of s.
    j is an index that is less than or equal to the length of t.
   spacePenalty is the numeric value for a base being paired with a space in an alignment.

Post:
    sAlignment and tAlignment are sequences in the optimal alignment for
        prefixes s[0..i - 1] and t[0..j - 1], respectively, with a dash indicating 
        a space.
Algorithm:
    if (i equals 0 and j equals 0)
	sAlignment <--	null
	tAlignment <--	null
    else
    	if (i > 0 and a[i][j ] equals a[i - 1][j] + SpacePenalty	// above
    		recAlign(s, t, a, match Value, mismatch Value, spacePenalty, 
    			sAlignment, tAlignment, i - 1, j)
       		sAlignment <--	sAlignment + s[i - 1]
    		tAlignment <--	tAlignment + ‘-‘
    	else if (j > 0 and a[i][j ] equals a[i][j - 1] + spacePenalty)	// left
    		recAlign(s, t, a, match Value, mismatch Value, spacePenalty, 
    			sAlignment, tAlignment, i, j - 1)
    		sAlignment <--	sAlignment + ‘-‘
    		tAlignment <--	tAlignment + t[j - 1]
    	else					// diagonal
    		recAlign(s, t, a, MatchValue, MismatchValue, spacePenalty,
    			sAlignment, tAlignment, i - 1, j - 1)
    		sAlignment <--	sAlignment + s[i - 1]
    		tAlignment <--	tAlignment + t[j - 1]

1. Consider the two sequences ATGAC and ACGC.
    a. List every possible alignment in which ATGAC has no blanks.
    b. List every possible alignment in which ATGAC has at least one blank         and ACGC has only blank(s) at the end.

2. a. Trace through the score algorithm for the strings and similarity array in Quick Review Question 4 and Figure 4.
    b. Trace through align and recAlign algorithms for the strings and similarity array in Quick Review Question 5 and Figure 5.

3. a. Develop the similarity array for strings s = CGTA and t = GCATG using the scoring of +1, -1, -2 from the example in this module.
    b. What is the similarity score?
    c. Give the optimal alignment that align returns.
    d. Give all optimal alignments for these strings.

4. Repeat Exercise 1 for strings s = CGTA and t = GCATG.

5. Repeat Exercise 1 for strings s = ATGAC and t = ACGC using the scoring that a mismatch or space costs 1, while a match has score 2.

6. What is the complexity of score?

7. What is the complexity of recAlign?

8. Change the priority of the align algorithm so that when a choice exists, the algorithm gives priority to going vertically up, next to going on the diagonal, and finally to going to the left. Because priority is given to going up, the resulting alignment is called the "upmost alignment".  Give the additional parameters that align and recAlign need.

9. Change the priority of the recAalign algorithm so that when a choice exists, the algorithm gives priority to going to the left, next to going on the diagonal, and finally to going vertically up. Because priority is given to going to the left, the resulting alignment is called the "downmost alignment".  Give the additional parameters that align and recAlign need.