1 / 66

Sequence Alignment

Understand the concept of sequence alignment and its importance in comparing DNA, RNA, and protein sequences. Discover the scoring and additive scoring rules used in aligning sequences and the dynamic programming algorithm for finding optimal alignments.

ajake
Download Presentation

Sequence Alignment

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Sequence Alignment AGGCTATCACCTGACCTCCAGGCCGATGCCC TAGCTATCACGACCGCGGTCGATTTGCCCGAC -AGGCTATCACCTGACCTCCAGGCCGA--TGCCC--- TAG-CTATCAC--GACCGC--GGTCGATTTGCCCGAC Definition Given two strings x = x1x2...xM, y = y1y2…yN, an alignment is an assignment of gaps to positions 0,…, N in x, and 0,…, N in y, so as to line up each letter in one sequence with either a letter, or a gap in the other sequence

  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. Sequence Comparison • Finding similarity between sequences is important for many biological questions • During evolution biological evolves (mutation, deletion, duplication, addition, move of subsequences…) • Homologous (share a common ancestor) sequences are (relatively) similar • Algorithms try to detect similar sequence that possibly share a common function

  4. g1 g2 Sequence Comparison (cont) 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

  5. 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

  6. Alignments -GCGC-ATGGATTGAGCGA TGCGCCATTGAT-GACC-A Three elements: • Matches • Mismatches • Insertions & deletions (indel)

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

  8. Scoring Alignments Intuition: • Similar sequences evolved from a common ancestor • Evolution changed the sequences from this ancestral sequence by mutations: • Substitution: 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

  9. Simple Scoring Rule Score each position independently: • Match m: +1 • Mismatch s: -1 • Indel d: -2 Score of an alignment is sum of position scores Scoring Function: Match: m m≥0 Mismatch: s s≤0 Gap: d s≤0 Score F = (#matches)m + (#mismatches)s + (#gaps)d

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

  11. 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.) • Exchange of an amino-acid by one with opposite properties • Probabilistic interpretation: (e.g.) How likely is one alignment versus another ?

  12. 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

  13. How do we compute the best alignment? AGTGCCCTGGAACCCTGACGGTGGGTCACAAAACTTCTGGA • Alignment is a path from cell (0,0) to cell (m,n) • Too many possible alignments: • O( 2M+N) AGTGACCTGGGAAGACCCTGACCCTGGGTCACAAAACTC

  14. The Optimal Score • The optimal alignment score between two sequences is the maximal score over all alignments of these sequences: • 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 two problems.

  15. Alignment is additive Observation: The score of aligning x1……xM y1……yN is additive Say that x1…xi xi+1…xM aligns to y1…yj yj+1…yN The two scores add up: V(x[1:M], y[1:N]) = V(x[1:i], y[1:j]) + V(x[i+1:M], y[j+1:N])

  16. Dynamic Programming • We will now describe a dynamic programming algorithm Suppose we wish to align x1……xM y1……yN Let V(i,j) = optimal score of aligning x1……xi y1……yj

  17. Dynamic Programming Notice three possible cases: • xi aligns to yj x1……xi-1 xi y1……yj-1 yj 2. xi aligns to a gap x1……xi-1 xi y1……yj - • yj aligns to a gap x1……xi - y1……yj-1 yj m, if xi = yj V(i,j) = V(i-1, j-1) + s, if not V(i,j) = V(i-1, j) + d V(i,j) = V(i, j-1) + d

  18. Dynamic Programming • How do we know which case is correct? Inductive assumption: V(i, j-1), V(i-1, j), V(i-1, j-1) are optimal Then, V(i-1, j-1) + σ(xi, yj) V(i, j) = max V(i-1, j) + d V( i, j-1) + d Where σ(xi, yj) = m, if xi = yj; s, if not

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

  20. 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:

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

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

  23. T S Dynamic Programming Algorithm

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

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

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

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

  28. The Needleman-Wunsch Matrix x1 ……………………………… xM Every nondecreasing path from (0,0) to (M, N) corresponds to an alignment of the two sequences y1 ……………………………… yN An optimal alignment is composed of optimal subalignments

  29. The Needleman-Wunsch AlgorithmGlobal Alignment Algorithm • Initialization. • F(0, 0) = 0 • F(0, j) = j  d • F(i, 0) = i  d • Main Iteration. Filling-in partial alignments • For each i = 1……M For each j = 1……N F(i-1,j-1) + s(xi, yj) [case 1] F(i, j) = max F(i-1, j) + d [case 2] F(i, j-1) + d [case 3] DIAG, if [case 1] Ptr(i,j) = LEFT, if [case 2] UP, if [case 3] • Termination. F(M, N) is the optimal score, and from Ptr(M, N) can trace back optimal alignment

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

  31. 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?

  32. 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 if lines are shorter). Note however that • This “trick” fails when we need to reconstruct the optimizing sequence. • Trace back information requires O(mn) memory bytes.

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

  34. s t 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]) • Thus, we need to compute these two quantities for all values of j

  35. Finding the Midpoint (Algorithm) Define • F[i,j] = d(s[1,i],t[1,j]) • B[i,j] = d(s[i+1,n],t[j+1,m]) • F[i,j] + B[i,j] = score of best alignment through (i,j) • We compute F[i,j] as we did before • We compute B[i,j] in exactly the same manner, going “backward” from B[n,m] • Requires linear space complexity because there is no need to keep trace back information in this step

  36. 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 worst, twice the cost of the regular solution.

  37. A variant of the NW algorithm: • Maybe it is OK to have an unlimited # of gaps in the beginning and end: ----------CTATCACCTGACCTCCAGGCCGATGCCCCTTCCGGC GCGAGTTCATCTATCAC--GACCGC--GGTCG-------------- • Then, we don’t want to penalize gaps in the ends

  38. Different types of overlaps

  39. The Overlap Detection variant Changes: • Initialization For all i, j, V(i, 0) = 0 V(0, j) = 0 • Termination maxi V(i, N) VOPT = max maxj V(M, j) x1 ……………………………… xM y1 ……………………………… yN

  40. Overlap Alignment Example s =PAWHEAE t =HEAGAWGHEE • Scoring system: • Match: +4 • Mismatch: -1 • Indel: -5

  41. Overlap Alignment • Initialization:V[i,0]=0,V[0,j]=0 • Recurrence: as in global alignment • Score: maximum value at the bottom line and rightmost line in the matrix

  42. Overlap Alignment Example s =PAWHEAE t =HEAGAWGHEE • Scoring system: • Match: +4 • Mismatch: -1 • Indel: -5

  43. Overlap Alignment Example s =PAWHEAE t =HEAGAWGHEE • Scoring system: • Match: +4 • Mismatch: -1 • Indel: -5

  44. Overlap Alignment Example The best overlap is: PAWHEAE------ ---HEAGAWGHEE Pay attention! A different scoring system could yield a different result, such as: ---PAW-HEAE HEAGAWGHEE-

  45. The local alignment problem Given two strings x = x1……xM, y = y1……yN Find substrings x’, y’ whose similarity (optimal global alignment value) is maximum e.g. x = aaaacccccgggg y = cccgggaaccaacc

  46. Why local alignment • Genes are shuffled between genomes • Portions of proteins (domains) are often conserved

  47. Cross-species genome similarity • 98% of genes are conserved between any two mammals • >70% average similarity in protein sequence hum_a : GTTGACAATAGAGGGTCTGGCAGAGGCTC--------------------- @ 57331/400001 mus_a : GCTGACAATAGAGGGGCTGGCAGAGGCTC--------------------- @ 78560/400001 rat_a : GCTGACAATAGAGGGGCTGGCAGAGACTC--------------------- @ 112658/369938 fug_a : TTTGTTGATGGGGAGCGTGCATTAATTTCAGGCTATTGTTAACAGGCTCG @ 36008/68174 hum_a : CTGGCCGCGGTGCGGAGCGTCTGGAGCGGAGCACGCGCTGTCAGCTGGTG @ 57381/400001 mus_a : CTGGCCCCGGTGCGGAGCGTCTGGAGCGGAGCACGCGCTGTCAGCTGGTG @ 78610/400001 rat_a : CTGGCCCCGGTGCGGAGCGTCTGGAGCGGAGCACGCGCTGTCAGCTGGTG @ 112708/369938 fug_a : TGGGCCGAGGTGTTGGATGGCCTGAGTGAAGCACGCGCTGTCAGCTGGCG @ 36058/68174 hum_a : AGCGCACTCTCCTTTCAGGCAGCTCCCCGGGGAGCTGTGCGGCCACATTT @ 57431/400001 mus_a : AGCGCACTCG-CTTTCAGGCCGCTCCCCGGGGAGCTGAGCGGCCACATTT @ 78659/400001 rat_a : AGCGCACTCG-CTTTCAGGCCGCTCCCCGGGGAGCTGCGCGGCCACATTT @ 112757/369938 fug_a : AGCGCTCGCG------------------------AGTCCCTGCCGTGTCC @ 36084/68174 hum_a : AACACCATCATCACCCCTCCCCGGCCTCCTCAACCTCGGCCTCCTCCTCG @ 57481/400001 mus_a : AACACCGTCGTCA-CCCTCCCCGGCCTCCTCAACCTCGGCCTCCTCCTCG @ 78708/400001 rat_a : AACACCGTCGTCA-CCCTCCCCGGCCTCCTCAACCTCGGCCTCCTCCTCG @ 112806/369938 fug_a : CCGAGGACCCTGA------------------------------------- @ 36097/68174 “atoh” enhancer in human, mouse, rat, fugu fish

  48. 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

  49. The Smith-Waterman algorithm Idea: Ignore badly aligning regions Modifications to Needleman-Wunsch: Initialization: V(0, j) = V(i, 0) = 0 0 Iteration: V(i, j) = max V(i – 1, j) + d V(i, j – 1) + d V(i – 1, j – 1) + σ(xi, yj)

  50. The Smith-Waterman algorithm Termination: • If we want the best local alignment… VOPT = maxi,j V(i, j) • If we want all local alignments scoring > t ?? For all i, j find V(i, j) > t, and trace back Complicated by overlapping local alignments

More Related