Sequence Alignment II Lecture 3 Background Readings The

Sequence Alignment II Lecture #3 Background Readings: The second chapter (pages 12 -45) in the text book, Biological Sequence Analysis, Durbin et al. , 2001. Chapter 3. 5. 1 in Introduction to Computational Molecular Biology, Setubal and Meidanis, 1997. 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. .

Reminder u Last class we discussed dynamic programming algorithms for l global alignment l local alignment u All of these assumed a scoring rule: that determines the quality of perfect matches, substitutions, insertions, and deletions. 2

Alignment in Real Life u One of the major uses of alignments is to find sequences in a “database. ” u The current protein database contains about 108 residues ! So searching a 103 long target sequence requires to evaluate about 1011 matrix cells which will take about three hours in the rate of 10 millions evaluations per second. u Quite annoying when, say, one thousand target sequences need to be searched because it will take about four months to run. 3

Heuristic Search u Instead, most searches rely on heuristic procedures u These are not guaranteed to find the best match u Sometimes, they will completely miss a high-scoring match We now describe the main ideas used by the best known of these heuristic procedures. 4

Basic Intuition u Almost all heuristic search procedures are based on the observation that real-life matches often contain long strings with gap-less matches. u These heuristic try to find significant gap-less matches and then extend them. 5

Banded DP that we have two strings s[1. . n] and t[1. . m] such that n m u If the optimal alignment of s and t has few gaps, then path of the alignment will be close to diagonal u Suppose t s 6

Banded DP u To find such a path, it suffices to search in a diagonal band of the matrix. u If the diagonal band consists of k diagonals (width k), then dynamic programming takes O(kn). u Much faster than O(n 2) of standard DP. t s V[i, i+k/2] Out of range V[i+1, i+k/2] V[i+1, i+k/2 +1] Note that for diagonals, i-j = constant. k 7

Banded DP for local alignment Problem: Where is the banded diagonal ? It need not be the main diagonal when looking for a good local alignment. t How do we select which subsequences to align using banded DP? k s We heuristically find potential diagonals and evaluate them using Banded DP. This is the main idea of FASTA. 8

Overview of FASTA Input: strings s and t, and a parameter ktup Output: A highly scored local alignment. 1. 2. 3. Find pairs of matching substrings s[i. . . i+ktup]=t[j. . . j+ktup] Extend to ungapped diagonals Extend to gapped matches using banded DP 9

Finding Potential Diagonals Suppose there exists a relatively long gap-less match S=****AGCGCCATGGATTGAGCGA* T=**TGCGACATTGATCGACCTA** u Each such sequence defines a potential diagonal as follows. If the first sequence starts at location i (e. g. , 5 above) and the second starts at location j (e. g. , 3 above), then the potential diagonal starts at location (i, j). j t u Can we identify potential diagonals quickly? i s u Such diagonals can then be evaluated using Banded DP. 10

Identifying Potential Diagonals Assumption: High scoring gap-less alignments contain several “seeds” of perfect matches S=****AGCGCCATGGATTGAGCGA* T=**TGCGACATTGATCGACCTA** Since this is a gap-less alignment, all perfect match regions reside on the same diagonal (defined by i-j). j t i s How do we find seeds efficiently ? 11

Formalizing the task Let ktup be a parameter denoting the seed length of interest. Task at hand (Identifying seeds): Find all pairs (i, j) such that s[i. . . i +ktup] = t[j. . . j+ktup] 12

Finding Seeds Efficiently Prepare an index table of the database sequence S such that for any sequence of length ktup, one gets the list of its positions in S. u March on the query sequence T while using the index table to list all matches with the database sequence S. u Index Table (ktup =2) AA AC AG 5, 19 AT 11, 15 CA 10 CC 9, 21 CG 7 … TT 16 S=****AGCGCCATGGATTGAGCGA* 5 10 15 20 7 89 T=**TGCGACATTGATCGACCTA** (-, 7) No match (10, 8) One match (11, 9), (15, 9) Two matches These steps take linear time: O(|s|+|t|). 13

Comments The maximal size of the index table is | |ktup where is the alphabet size (4 or 20). For small ktup, the entire table is stored. For large ktup values, one should keep only entries for tuples actually found in the database, so the index table size is indeed linear. In this case, hashing is needed. Typical values of ktup are 1 -2 for Proteins and 4 -6 for DNA. Tradeoffs of these values to be discussed. The index table is prepared for each database sequence ahead of users’ matching requests, at compilation time. So matching time is O(t). Index table AA AC AG AT CA CC CG … TT 5, 19 11, 15 10 9, 21 7 16 14

Identifying Potential Diagonals u. Input: Sets of pairs. E. g, (6, 4), (10, 8), (14, 12), (15, 10), (20, 4) … u. Task: Locate sets of pairs that are on the same u Method: Sort according to the difference i-j. 6 10 14 4 8 12 diagonal. 20 S=***AGCGCCATGGATTGAGCGA* T=**TGCGACATTGATCGACCTA** i-j = 2; 6 -4 ; 10 -8; 14 -12 i-j = 20 -4=16 j t i s 15

Processing Potential Diagonals For high i-j offset frequency, namely, diagonals with many pieces, combine the pieces into regions by extending pieces greedily along the diagonal as long as the score improves (and never below some score value). t s 16

FASTA’s Final steps: using banded DP l l List the highest scoring diagonal matches Run banded DP on regions containing a high scoring diagonal (say with width 12). t 2 s 3 1 Hence, the algorithm may combine some diagonals into gapped matches. In the example above it could combine diagonals 2 and 3). 17

FASTA- practical choices Most applications of FASTA use very small ktup (1 -2 for proteins, and 4 -6 for DNA). Higher values yield less potential diagonals. Hence to search around potential diagonals (DP) is faster. But the chance to miss an optimal local alignment is increased. t s Some implementation choices /tricks have not been explicated herein. 18

BLAST Overview Based on similar ideas described earlier (High scoring pairs rather than exact k tuples as seeds). Uses an established statistical framework to determine thresholds. The new PSI-BLAST (Position Specific Iterated – BLAST ) is the state of the art sequence comparison software. Iterative Procedure l Perform BLAST on a database l Uses significant alignments to construct “position specific” score matrix. l This matrix is used in the next round of database searching until no new significant alignments are found. Can sometime detect remote homologs. 19

Where do scoring rules come from ? We have defined an additive scoring function by specifying a function ( , ) such that l (x, y) is the score of replacing x by y l (x, -) is the score of deleting x l (-, x) is the score of inserting x But how do we come up with the “correct” score ? Answer: By encoding experience of what are similar sequences for the task at hand. Similarity depends on time, evolution trends, and sequence types. 20

Why use probability to define and/or interpret a scoring function ? • Similarity is probabilistic in nature because biological changes like mutation, recombination, and selection, are not deterministic. • We could answer questions such as: • How probable two sequences are similar? • Is the similarity found significant or random? • How to change a similarity score when, say, mutation rate of a specific area on the chromosome becomes known ? 21

A Probabilistic Model u For now, we will focus on alignment without indels. u For now, we assume each position (nucleotide /amino-acid) is independent of other positions. u We consider two options: M: the sequences are Matched (related) R: the sequences are Random (unrelated) 22

Unrelated Sequences u Our random model of unrelated sequences is simple l Each position is sampled independently from a distribution q( ) over the alphabet u Then: 23

Related Sequences u We assume that each pair of aligned positions (s[i], t[i]) evolved from a common ancestor. u Let p(a, b) be the probability that some ancestral letter evolved into this particular pair of letters. 24

Odd-Ratio Test for Alignment If Q > 1, then the two strings s and t are more likely to be related (M) than unrelated (R). If Q < 1, then the two strings s and t are more likely to be unrelated (R) than related (M). 25

Log Odd-Ratio Test for Alignment Taking logarithm of Q yields Score(s[i], t[i]) If log Q > 0, then s and t are more likely to be related. If log Q < 0, then they are more likely to be unrelated. How can we relate this quantity to a score function ? 26

Probabilistic Interpretation of Scores u We define the scoring function via u Then, the score of an alignment is the log-ratio between the two models: l Score > 0 l Score < 0 Model is more likely Random is more likely 27

Constructing Scoring Rules The formula suggests how to construct a scoring rule: u Estimate p(·, ·) and q(·) from the data u Compute (a, b) based on the estimated p(·, ·) and q(·) u How to estimate these parameters is the subject matter of parameter estimation in Statistics. 28

Modeling Assumptions u It is important to note that this interpretation depends on our modeling assumption!! u For example, if we assume that the letter in each position depends on the letter in the preceding position, then the likelihood ratio will have a different form. u Next week we will start to discus languages that allow us to employ more complex models (FSA, and Bayesian networks). 29

Estimating p(·) u Suppose we are given a long string s[1. . n] of letters from (say, the concatenation of all sequences in a database). u We want to estimate the distribution q(·) Likelihood function: ML parameters MAP parameters Number of times a appears in s 30

Estimating p(·, ·) Intuition: u Find pair of aligned sequences s[1. . n], t[1. . n], u Estimate probability of pairs: Number of times a is aligned with b in (s, t) u Again, s and t can be the concatenation of many aligned pairs from a database 31

Estimating p(·, ·) for proteins Generate a large diverse collection of accepted mutations. An accepted mutation is a mutation due to an alignment of closely related protein sequences. For example, Hemoglobin alpha chain in humans and other organisms (homologous proteins). Let pa = na/n where na is the number of occurrences of letter a and n is the total number of letters in the collection, so n = ana. Mutation counts be the number of mutations a b, be the total number of mutations that involve a, be the total number of amino acids involved in a mutation. Note that f is twice the number of mutations. 32

PAM-1 matrices Define Mab to be the symmetric probability matrix for switching between a and b. We set, Maa = 1 – ma, so that ma is the probability that a is involved in a change. We define Mab, such that only 1% of amino acids change according to this matrix or 99% don’t. Hence the name, 1 -Percent Accepted Mutation (PAM). In other words, 33

PAM-1 matrices We wish that ma will be proportional to the relative mutability of letter a compared to other letters. where K is a proportional constant. We select K to satisfy the PAM-1 definition: So K=100 for PAM-1 matrices. Note that K=50 yields 2% change, etc. 34

Evolutionary distance The choice that 1% of amino acids change (and that K =100) is quite arbitrary. It could fit specific set of proteins whose evolutionary distance is such that indeed 1% of the letters have mutated. This is a unit of evolutionary change, not time because evolution acts differently on distinct sequence types. What is the substitution matrix for k units of evolutionary time ? 35

Model of Evolution We make some assumptions: 1. Each position changes independently of the rest 2. The probability of mutations is the same in each position 3. Evolution does not “remember” T A T C C C t t+ t+2 t+3 G G t+4 Time 36

Model of Evolution u How do we model such a process? u This process is called a Markov Chain A chain is defined by the transition probability u P(Xt+ =b|Xt=a) - the probability that the next state is b given that the current state is a u We often describe these probabilities by a matrix: M[ ]ab = P(Xt+ =b|Xt=a) 37

Multi-Step Changes on Mab, we can compute the probabilities of changes over two time periods u Based Using Conditional independence (No memory) u Thus u By M[2 ] = M[ ] induction (HMW exercise): M[n ] = M[ ] n 38

Longer Term Changes M[ ] = M (PAM-1 matrices) u Use M[n ] = Mn (PAM-n matrices) u Define u Estimate u Use this quantity to define the score for your application of interest. 39

Comments regarding PAM Historically researchers use PAM-250. (The only one published in the original paper. ) u Original PAM matrices were based on small number of proteins (circa 1978). Later versions use many more examples. u Used to be the most popular scoring rule, but there are some problems with PAM matrices. u 40

Degrees of freedom in PAM definition With K=100 the 1 -PAM matrix is given by With K=50 the basic matrix is different, namely: Thus we have two different ways to estimate the matrix M[4 ] : Use the 1 -PAM matrix to the fourth power: M[4 ] = M[ ] 4 Or Use the K=50 matrix to the second power: M[4 ] = M[2 ] 2 41

Problems in building distance matrices u How do we find pairs of aligned sequences? u How far is the ancestor ? l earlier divergence low sequence similarity l later divergence high sequence similarity E. g. , M[250 ] is known not reflect well long period changes. u Does one letter mutate to the other or are they both mutations of a third letter ? 42

BLOSUM Outline Idea: use aligned ungapped regions of protein families. These are assumed to have a common ancestor. Similar ideas but better statistics and modeling. u Procedure: l Cluster together sequences in a family whenever more than L% identical residues are shared. l Count number of substitutions across different clusters (in the same family). l Estimate frequencies using the counts. u Practice: Blosum 50 and Blosum 62 are wildly used. (See page 43 -44 in the text book). u Considered state of the art nowadays. 43