We showed how to calculate Edit distance

“Distance” means that smaller numbers are better

If the distance is 0, then the two sequences are identical

If the distance is small, the two sequences are similar

If the distance is big, the sequences are different

The strategy is called Dynamical Programming

Here “Programming” means “finding the best option”

• divides a large problem into a series of smaller problems
• uses the smaller solutions to build a solution to the larger problem

We will talk about the Needleman–Wunsch Algorithm

Published in 1970 by Saul B. Needleman and Christian D. Wunsch

They used dynamic programming to find a score for the alignment

The philosophy is the same, but we look for big numbers instead of small numbers

Initially there is no difference in the resulting alignment, but later it will be important.

Let’s say we have two sequences q and s

They have length m and n, respectively

• q="TTAATAATCAGTTTTCAT"
• s="TTAACATTAGTATTTCA"

q[1:3] is "TTA". s[5:7] is "CAT"

We see that q[1:m] is the same as q, and s[1:n]=s

We want to evaluate the score between q and s

We call it Score(q, s). We can write Score(q[1:m], s[1:n])

Since we divide the problem into smaller ones, we need to calculate Score(q[1:i], s[1:j]) for different values of i and j

To simplify, we write M[i, j] instead of Score(q[1:i], s[1:j])

Therefore the edit distance will be

Score(q, s) = Score(q[1:m], s[1:n]) = M[m, n]

If we pair two letters x and y, we need to know how they contribute to the score value. For example,

• if x == y the score increases by 1
• if x != y the score decreases by 1
  • (increases by -1)

In general, we have a substitution matrix C[x, y] for each x and y

If one of the letters is - (a gap), the score reduces by a value G

Since Score(q, s) = M[m, n] we just need to evaluate M[m, n]

Let’s assume that we can solve the problem just before the last position of the alignment

We want to align the last two letters. We have three options

• pairing q[m] with s[n], match or mismatch
• pairing q[m] with a gap -, i.e. insertion
• pairing a gap - with s[n], i.e. deletion

We see that to calculate the optimal alignment score we need to calculate all the values of the matrix M[i,j]

In every step we choose the option that produces the highest score

M[i, j] = max(
              M[i-1, j-1] + C[q[i], s[j]],
              M[i-1, j] - G,
              M[i, j-1] - G
             )

To get M[m,n] we need M[m-1,n-1] and others

To get M[m-1,n-1] we need M[m-2,n-2] and others

What happens when we get to M[0,0], M[0,1] or M[1,0]?

That is, what happens when we look before the beginning of the sequences?

• If one of the sequences is larger than the other
  • or if one is contained in the other
• then we can put gaps outside the sequence

----------AACCTACCTCAAAGAGGGGGATA------------------
ACACGTGGGCAACCTGCCTCAAAGTGGGGGATAGCCTTCCGAAAGGAAGAT

To represent that, we fix the values of M[i,0] and M[0,j]

If we follow the idea that both sequences must be equal, then each gap outside will contribute with -G to the score

M[0,0] = 0, M[1,0] = -G, M[2,0] = -2*G, M[i,0] = -i*G

M[0,0] = 0, M[0,1] = -G, M[0,2] = -2*G, M[0,j] = -j*G

This strategy is called Global Alignment

If instead we assume that one of the sequences is smaller than the other, then each gap outside the sequence must not change the score

M[0,0] = 0, M[1,0] = 0, M[2,0] = 0, M[i,0] = 0

M[0,0] = 0, M[0,1] = 0, M[0,2] = 0, M[0,j] = 0

This strategy is called Semi-Global Alignment

In semi-global alignment it is not true that Score(q, s) = M[m, n]

Instead, we have to look at all the values in the last row

Score_semiglobal(q, s) = max(M[m, 1:n])

The scoring matrix (also known as a substitution matrix) has one row and one column for each letter in our alphabet

For example 4 rows and 4 columns for DNA sequences

In DNA it is common to have +1 for match and -2 for mismatch

Sometimes we can use (+1,-3), (+1,-1) or even (4,-5)

For proteins the alphabet has 20 letters, so the scoring matrix is 20x20

The scores of these matrices are based on the relative frequency of each mutation in nature

Obviously, we only care about non-lethal mutations

More specifically, we care about the Point Mutations that are Accepted

These are called PAM matrices (Margaret Dayhoff, 1978)

We prefer this alignment

GGGTAACCTACCTC
GGGCAACCTGCCTC

instead of this other alignment

GGGT-AACCTA-CCTC
GGG-CAACCT-GCCTC

Thus, gap penalty must be greater than mismatch penalty

We prefer this alignment

TCAAAGAG---GATA
TCA--GAGGGGGATA

instead of this other alignment

TCAAAGA-G-G-ATA
TC-A-GAGGGGGATA

We want few long gaps instead of many short gaps

Gap values must reflect how real insertions and deletions occur in nature

We observe that, once an indel event starts, it can easily grow

If the polymerase jumps, then it can jump a long distance

To represent this, we use affine gaps

So far we considered only linear gaps

The penalty of n consecutive gaps is n*G

Now we consider affine gaps, where the first gap is expensive, but the consecutive are cheap

The penalty of n consecutive gaps is G + n*E

G is the initial gap penalty, E is the gap extension penalty

Modify the Excel spreadsheet that we did last class to handle border gaps

• in global alignment
• in semi-global alignment