Introduction to Bioinformatics Lecture V Alignment Counting and

Introduction to Bioinformatics: Lecture V Alignment Counting and Alignment Algorithms Jarek Meller Division of Biomedical Informatics, Children’s Hospital Research Foundation & Department of Biomedical Engineering, UC JM - http: //folding. chmcc. org 1

Outline of the lecture Complexity of inexact string matching: exercises in counting alignments with gaps l The dynamic programming algorithm for sequence alignment: how it works l The dynamic programming algorithm for sequence alignment: why it works l Limitations and faster heuristic approaches: BLAST l JM - http: //folding. chmcc. org 2

Web watch: “Genes and Disease” and other NCBI resources Genes (proteins) work in herds. Being co-localized may imply co-expression and interactions. Always check the neighbors of your favorite gene! http: //www. ncbi. nlm. nih. gov/books/bv. fcgi? call=bv. View. . Show. TOC&rid=gnd. TOC&depth=2 Additional reading materials regarding sequence alignment: check out the web site for the course … JM - http: //folding. chmcc. org 3

How many alignments with gaps are there? All the possible alignments (with gaps, however without the unnecessary alignment of two gaps against each other) may be represented in the form of a grid with only three steps (South, East, Southeast) allowed, i. e. , there is a bijection between the set of walks on such a grid and the set of alignments. a 1 a 2 1 1 3 5 7 1 5 13 25 1 7 25 63 0 b 1 b 2 b 3 a 3 321 1683 JM - http: //folding. chmcc. org 4

How many alignments with gaps are there? JM - http: //folding. chmcc. org 5

How many alignments with gaps are there? JM - http: //folding. chmcc. org 6

How many alignments with gaps are there? JM - http: //folding. chmcc. org 7

Time for the main idea of the algorithm … Suppose we knew best alignments (and their scores) up to the nodes which delineate the yellow part of the DP graph. What would be the best extension, given that we have three choices: 1) Align the last characters in each string (diagonal extension) and add the score for this pair, e. g. , s(a 3, b 3)= -5 2) Align a 3 to a gap (horizontal extension), s(a 3, -)= -8 3) Align b 3 to a gap (vertical extension), s(-, b 3)= -8 a 1 Note well, the score for an extension does not depend on the alignment up to this point. a 2 a 3 b 1 b 2 5 12 b 3 -3 JM - http: //folding. chmcc. org 4 8

Tracing back optimal local extensions given the best alignment up to the previous node in the graph … a 1 a 2 a 3 b 1 b 2 b 3 -3 s(a 2, b 1)= 10 s(a 3, b 2)= 10 s(a 1, b 1)= -5 2 -6 5 12 b 1 a 2 a 3 -8 -16 -5 2 -6 5 12 4 b 2 b 3 -3 JM - http: //folding. chmcc. org 4 9

This conceptual step may be reversed to obtain the best score and alignment up to a given point The score is a sum of independent piecewise scores, in particular, the score up to a certain point is the best score up to the point one step before plus the incremental score of the new step: l F(i-1, j-1) F(i-1, j) F(i, j-1) (i, j) Global alignment (Needleman-Wunsch): F(0, 0) = 0; F(k, 0) = F(0, k) = - k d; F(i, j) = max { F(i-1, j-1)+s(ai, bj) ; F(i-1, j)-d ; F(i, j-1)-d } l Local alignment (Smith-Waterman): F(0, 0) = 0; F(k, 0) = F(0, k) = 0; F(i, j) = max { 0 ; F(i-1, j-1)+s(ai, bj) ; F(i-1, j)-d ; F(i, j-1)-d } JM - http: //folding. chmcc. org 10

The general scheme of the NW algorithm l Use the recurrence relations, starting from the left upper corner (convention). l Find the highest score in the DP table (last, bottom right cell in the global alignment by definition) l Trace back the alignment using the pointers in the DP graph that show the best local steps led to the best overall alignment. JM - http: //folding. chmcc. org 11

Examples of pairwise scores from the Blosum 50 matrix H E A G A W G H E P -2 -1 -1 -2 -1 -4 -2 -2 -1 A -2 -1 5 0 5 -3 0 -2 -1 W -3 -3 -3 15 -3 -3 -3 H 10 0 -2 -2 -2 -3 -2 10 0 E 0 6 -1 -3 -3 0 6 JM - http: //folding. chmcc. org 12

An example of DP table for global alignment H E A G A W G H E 0 <-8 <-16 <-24 <-32 <-40 <-48 <-56 <-64 <-72 P ^-8 *-2 *-9 *-17 <-25 *-33 <-41 <-49 <-57 *-65 A ^-16 ^-10 *-3 *-4 <-12 *-20 <-28 <-36 <-44 <-52 W ^-24 ^-18 ^-11 *-6 *-7 *-15 *-5 <-13 <-21 <-29 H ^-32 *-14 *-18 *-13 *-8 *-9 ^-13 *-7 *-3 <-11 E ^-40 ^-22 *-8 <-16 ^-16 *-9 *-12 ^-15 *-7 *3 HEAGAWGHE --P-AW-HE JM - http: //folding. chmcc. org 13

An example of DP table for local alignment H E A G A W G H E 0 0 0 0 0 P 0 0 0 0 0 A 0 0 0 *5 0 0 W 0 0 *20 <12 <4 0 H 0 *10 <2 0 0 0 ^12 *18 *22 <14 E 0 ^2 *16 <8 0 0 ^4 ^10 *18 *28 AWGHE AW-HE JM - http: //folding. chmcc. org 14

Why does it work? l All the possible alignments (with gaps) are represented in the DP table (graph) l The score is a sum of independent piecewise scores, in particular, the score up to a certain point is the best score up to the point one step before plus the incremental score of the new step l Once the best score in the DP table is found the trace back procedure generates the alignment since only the best “past” leading to the present score is represented by the pointers between the cells JM - http: //folding. chmcc. org 15

Why does it work? Formally, one needs to show that the walk (alignment) found using the NW DP recurrence relations and the traceback procedure is indeed optimal, i. e. , it maximizes the alignment score. An argument instead of a proof In case of global alignment each path starts at cell (0, 0) and must end at cell (n, m). Consider the latter cell and the immediate past that led to it through one (most favorable together with the cost of the incremental step) of the 3 neighboring cells. Changing the last step (e. g. from initially chosen, optimal diagonal step) to an alternative one does not affect the scores at the preceding cells that represent the best trajectory up to this point. Clearly we get suboptimal solution if we assume that optimal solutions have been found in the previous steps. Hence, formalizing this argument we get proof by induction with reductio ad absurdum. Problem How to modify the NW algorithm for suffix-prefix matches? Problem What is the meaning of the cut off threshold for the SW algorithm? JM - http: //folding. chmcc. org 16

Approximate, heuristic solutions may be nearly as good and much faster: BLAST algorithm l BLAST approach: gapless seeds (High Scoring Pairs with well defined confidence measures), DP extensions: JM - http: //folding. chmcc. org 17