Sequence Alignment I Lecture #2

1. Sequence Alignment ILecture #2 Background Readings: The second chapter (pages 12-45) in the text book, Biological Sequence Analysis, Durbin et al., 2001. This class has been edited from Nir Friedman’s lecture which is available at www.cs.huji.ac.il/~nir. Changes made by Dan Geiger. Thanks to Shlomo Moran for several improvements. .

2. Sequence Comparison Much of bioinformatics involves sequences • DNA sequences • RNA sequences • Protein sequences We can think of these sequences as strings of letters • DNA & RNA: alphabet of 4 letters • Protein: alphabet  of 20 letters

3. g1 g2 Sequence Comparison (cont) • Finding similarity between sequences is important for many biological questions For example: • Find similar proteins • Allows to predict function & structure • Locate similar subsequences in DNA • Allows to identify (e.g) regulatory elements • Locate DNA sequences that might overlap • Helps in sequence assembly

4. Sequence Alignment Input: two sequences over the same alphabet Output: an alignment of the two sequences Example: • GCGCATGGATTGAGCGA • TGCGCCATTGATGACCA A possible alignment: -GCGC-ATGGATTGAGCGA TGCGCCATTGAT-GACC-A

5. Alignments -GCGC-ATGGATTGAGCGA TGCGCCATTGAT-GACC-A Three elements: • Perfect matches • Mismatches • Insertions & deletions (indel)

6. Choosing Alignments There are many possible alignments For example, compare: -GCGC-ATGGATTGAGCGA TGCGCCATTGAT-GACC-A to ------GCGCATGGATTGAGCGA TGCGCC----ATTGATGACCA-- Which one is better?

7. Scoring Alignments Intuition: • Similar sequences evolved from a common ancestor • Evolution changed the sequences from this ancestral sequence by mutations: • Replacements: one letter replaced by another • Deletion: deletion of a letter • Insertion: insertion of a letter • Scoring of sequence similarity should examine how many and which operations took place

8. Simple Scoring Rule Score each position independently: • Match: +1 • Mismatch : -1 • Indel -2 Score of an alignment is sum of position scores

9. Example Example: -GCGC-ATGGATTGAGCGA TGCGCCATTGAT-GACC-A Score: (+1x13) + (-1x2) + (-2x4) = 3 ------GCGCATGGATTGAGCGA TGCGCC----ATTGATGACCA-- Score: (+1x5) + (-1x6) + (-2x11) = -23

10. More General Scores • The choice of +1,-1, and -2 scores is quite arbitrary • Depending on the context, some changes are moreplausible than others • Exchange of an amino-acid by one with similar properties (size, charge, etc.) vs. • Exchange of an amino-acid by one with opposite properties • Probabilistic interpretation: (e.g.) How likely is one alignment versus another ?

11. Additive Scoring Rules • We define a scoring function by specifying a function • (x,y) is the score of replacing x by y • (x,-) is the score of deleting x • (-,x) is the score of inserting x • The score of an alignment is the sum of position scores

12. The Optimal Score • The optimal (maximal) score between two sequences is the maximal score of all alignments of these sequences, namely, • Computing the maximal score or actually finding an alignment that yields the maximal score are closely related tasks with similar algorithms. • We now address these problems.

13. Computing Optimal Score • How can we compute the optimal score ? • If |s| = n and |t| = m, there are many alignments ! • Exercise : Show that the number of legal alignments denoted by A(m,n), where n > m, satisfies (legal means not to a “-” versus another “-”) : • The additive form of the score allows us to perform dynamic programming to compute optimal score efficiently.

14. Recursive Argument • Suppose we have two sequences:s[1..n+1] and t[1..m+1] The best alignment must be one of three cases: 1. Last match is (s[n+1],t[m +1] ) 2. Last match is (s[n +1],-) 3. Last match is (-, t[m +1] )

15. Recursive Argument • Suppose we have two sequences:s[1..n+1] and t[1..m+1] The best alignment must be one of three cases: 1. Last match is (s[n+1],t[m +1] ) 2. Last match is (s[n +1],-) 3. Last match is (-, t[m +1] )

16. Recursive Argument • Suppose we have two sequences:s[1..n+1] and t[1..m+1] The best alignment must be one of three cases: 1. Last match is (s[n+1],t[m +1] ) 2. Last match is (s[n +1],-) 3. Last match is (-, t[m +1] )

17. Recursive Argument Define the notation: • Using our recursive argument, we get the following recurrence for V:

18. T S Recursive Argument • Of course, we also need to handle the base cases in the recursion: AA - - versus We fill the matrix using the recurrence rule:

19. T S Dynamic Programming Algorithm We continue to fill the matrix using the recurrence rule

20. T S +1 -2 -A A- -2 (A- versus -A) Dynamic Programming Algorithm versus

21. T S Dynamic Programming Algorithm

22. T S Conclusion: d(AAAC,AGC) = -1 Dynamic Programming Algorithm

23. T S Reconstructing the Best Alignment • To reconstruct the best alignment, we record which case(s) in the recursive rule maximized the score

24. T S AAAC AG-C Reconstructing the Best Alignment • We now trace back a path that corresponds to the best alignment

25. T S AAAC -AGC AAAC AG-C Reconstructing the Best Alignment • Sometimes, more than one alignment has the best score AAAC A-GC

26. T S Time Complexity Space: O(mn) Time: O(mn) • Filling the matrix O(mn) • Backtrace O(m+n)

27. Space Complexity • In real-life applications, n and m can be very large • The space requirements of O(mn) can be too demanding • If m = n = 1000, we need 1MB space • If m = n = 10000, we need 100MB space • We can afford to perform extra computation to save space • Looping over million operations takes less than seconds on modern workstations • Can we trade space with time?

28. A G C 1 2 3 0 0 0 -2 -4 -6 1 -2 1 -1 -3 A 2 -4 -1 0 -2 A 3 -6 -3 -2 -1 A 4 -8 -5 -4 -1 C Why Do We Need So Much Space? To compute V[n,m]=d(s[1..n],t[1..m]), we need only O(min(n,m)) space: • Compute V(i,j), column by column, storing only two columns in memory (or line by line ). Note however that • This “trick” fails when we need to reconstruct the optimizing sequence. • Trace back information requires O(mn) memory bytes.

29. t • Construct two alignments • A= s[1,n/2] vs t[1,j] • B= s[n/2+1,n] vs t[j+1,m] s • Return the concatenated alignment AB Space Efficient Version: Outline Input: Sequences s[1,n] and t[1,m] to be aligned. Idea: perform divide and conquer • Find position (n/2, j) at which the best alignment crosses a midpoint

30. We need to find j that maximizes this score. Such j determines a point (n/2,j) through which the best alignment passes. t • Thus, we need to compute these two quantities for all values of j s Finding the Midpoint • The score of the best alignment that goes through j equals: d(s[1,n/2],t[1,j]) + d(s[n/2+1,n],t[j+1,m])

31. t s Finding the Midpoint (Algorithm) Define • V[i,j] = d(s[1..i],t[1..j]) (As before) • B[i,j] = d(s[i+1..n],t[j+1..m]) (Symmetrically) • F[i,j] + B[i,j] = score of best alignment through (i,j)

32. Computing V[i,j] V[i,j] = d(s[1..i],t[1..j]) As before: Requires linear space complexity

33. Computing B[i,j] B[i,j] = d(s[i+1..n],t[j+1..m]) Symmetrically (replacing i with i+1, and j with j+1): Requires linear space complexity

34. t s Time Complexity Analysis • Time to find a mid-point: cnm (c - a constant) • Size of recursive sub-problems is (n/2,j) and (n/2,m-j-1), hence T(n,m) = cnm + T(n/2,j) + T(n/2,m-j-1) Lemma: T(n,m)  2cnm Proof (by induction): T(n,m)  cnm + 2c(n/2)j + 2c(n/2)(m-j-1)  2cnm. Thus, time complexity is linear in size of the problem. At worse, twice the cost of the regular solution.

35. Local Alignment Consider now a different question: • Can we find similar substrings of s and t • Formally, given s[1..n] and t[1..m] find i,j,k, and l such that d(s[i..j],t[k..l]) is maximal

36. Local Alignment • As before, we use dynamic programming • We now want to setV[i,j] to record the best alignment of a suffix of s[1..i] and a suffix of t[1..j] • How should we change the recurrence rule? • Same as before but with an option to start afresh • The result is called the Smith-Waterman algorithm

37. Local Alignment New option: • We can start a new match instead of extending a previous alignment Alignment of empty suffixes

38. T S Local Alignment Example s =TAATA t =TACTAA

39. T S Local Alignment Example s =TAATA t =TACTAA

40. T S Local Alignment Example s =TAATA t =TACTAA

41. T S Local Alignment Example s =TAATA t =TACTAA A maximal alignment starts at a maximum entry and follows backwards till a zero entry.

42. T S Local Alignment Example s =TAATA t =TACTAA

43. Variants of Sequence Alignment We have seen two variants of sequence alignment: • Global alignment • Local alignment Other variants in the book and in tutorial time: • Finding best overlap • Using an affine cost d(g) = -d –(g-1)e for gaps of length g. The –d is for introducing a gap and –e for continuing the gap. We used d=e=2. We could use smaller e. These variants are based on the same basic idea of dynamic programming.

44. Remark: Edit Distance Instead of speaking about the score of an alignment, one often talks about an edit distance between two sequences, defined to be the “cost” of the “cheapest” set of edit operations needed to transform one sequence into the other. • Cheapest operation is “no change” • Next cheapest operation is “replace” • The most expensive operation is “add space”. Our goal is now to minimize the cost of operations, which is equivalent to maximizing the corresponding score (e.g., change scores to edit costs as follows: 1 -1,-1 1,-22).