Pairwise Sequence Alignment and Database Searching

1. Pairwise Sequence Alignment and Database Searching

3. Objectives • What is the function of this gene? • Do other genes have this functional motif? • Can I predict the higher order structure of this protein? • Is this gene a member of a known gene family? • Do other organisms have this gene?

4. Objectives • What is the function of this gene? • Do other genes have this functional motif? • Can I predict the higher order structure of this protein? • Is this gene a member of a known gene family? • Do other organisms have this gene?

5. Intuition • “Similar” sequences should have (long) regions of similar/identical residues. • Why? • Evolution: descent from a common ancestral sequence • Functional/structural convergence

6. A sequence similarity search is an application of pairwise sequence alignment!

7. Dotplots G A G C T T A G G G C C T T T G G G A A A

8. General Alignment/Search Issues • Search using amino acid sequence if possible • Why? Protein evolution is slower than DNA sequence evolution • Ask the program to translate your query sequence in all 6 possible reading frames. • Statistical theory is based on unrealistic assumptions; consider searches as exploratory analyses.

9. Substitutions Insertions Deletions A B Alignment Jargon ancestor Evolutionarily related sequences differ from one other because of several processes: Observed sequences

10. ACG GCG ACG Alignment Jargon GCG || ACG Substitution A-G • 1 mismatch • 2 matches

11. ACG ATCG ACG Alignment Jargon ATCG | || A-CG Insertion T • 0 mismatches • 3 matches • 1 gap

12. Deletion ATCG T ATCG ACG Alignment Jargon ATCG | || A-CG • 0 mismatches • 3 matches • 1 gap

13. Alignment Jargon Results of insertion and deletion events can be indistinguishable. Indel: INsertion or DELetion

14. Sequence Alignment Sequence alignment is simply the “optimal” assignment of substitution and indel events to a pair of sequences. • Global alignment: align entire sequences • Local alignment: find best matching regions of sequences

15. Dotplots G A G C T T A G G G C C T T T G G G A A A

16. Measuring Alignment Quality Good alignments should have … • “many” exact matches • “few” mismatches • “many” of the mismatches should be similar residues • “few” gaps

17. Measuring Alignment Quality Begin with... QTRPQNVLNPP ||| STRQNVINPWAAQ Longest Exact Match S = 3a S=alignment score a=match score

18. Measuring Alignment Quality … allow some mismatches QTRPQNVLNPP ||| || STRQNVINPWAAQ S=alignment score a=match score b=mismatch penalty S = 5a - 1b

19. Measuring Alignment Quality QTRPQNVLNPP || ||| || STR-QNVINPWAAQ …and finally, introduce some gaps S=alignment score a=match score b=mismatch penalty c=gap penalty S = 7a - 1b -1c

20. Scoring Issues • Relative costs of matches, mismatches, and gaps should depend on their probabilities (rare events receive higher penalties) • In practice, the appropriate costs are rarely known. • A variety of scoring matrices are available.

21. Scoring Matrices • Scoring matrix specifies a score, sij, for aligning aa i with aa j. • Choice of matrix depends (ideally) on the divergence level of desired/expected hits. • Examples: PAM, BLOSUM • Both can be modified for different divergence levels (eg, BLOSUM40, BLOSUM62) • Advice: try several matrices when possible.

22. Dayhoff Family of Matrices • Based on an empirically-derived, probabilistic model of protein evolution. • Using closely-related sequences, count the frequencies of different types of amino acid substitutions, and use these frequencies to construct the scoring matrix. • Extrapolate to higher levels of divergence

23. Dayhoff Family of Matrices • Dayhoff model measures sequence evolution in units of “PAMs” • 2 sequences separated by k PAMs are expected to have experienced k substitutions per 100 amino acid sites since they diverged from a common ancestral sequence • Mutability of an aa is its relative rate of change (amino acids with high mutabilities are more likely to change) • Mutability of alanine was defined to be 100.

24. Dayhoff Family of Matrices Problems with the original Dayhoff scheme • It does not consider the genetic code. • Not all amino acid substitutions can occur by a single nucleotide substitution event. • Parameters were estimated from a small sample of closely related proteins. • Evolution at the “average site” of the “average protein” is being modeled.

25. BLOSUM Scoring Matrices • Construct a database of “blocks”: ungapped, aligned conserved regions of proteins) • Cluster sequences within a block that are more similar than a chosen threshold (eg, 62% for BLOSUM62) • Represent each cluster of sequences using their “average” sequence

26. BLOSUM Scoring Matrices • After the averaging procedure, each modified block consists of average sequence(s) and sequences that did not cluster with the others • Examine all pairs of sequences in the modified blocks. Tabulate pij the probability of observing aa i in one sequence and aa j in a second sequence.

27. BLOSUM Scoring Matrices • Similarly, tabulate i, the probability that a particular position in a block will be aa i. • The ratio of the probability of observing the residue pair ij in in two related sequences to the same probability for unrelated sequences is where

28. BLOSUM Scoring Matrices • The logarithm of this ratio serves as a score for pairing aa i with aa j. • Why 62%? • Henikoff and Henikoff found that use of this threshold in conjunction with BLAST was more effective than other values for correctly identifying which sequences belonged in which groups

29. Scoring Indels GCCTATTG GCCTATTG | ||| | ||| AC--TTTG A-C-TTTG Using a scoring scheme of match (a), mismatch (-b), gap (-c), each of these alignments would score: 4a - 2b -2c

30. A more biologically realistic scoring scheme is to have a relatively large penalty for opening, or initiating, a gap, and a smaller penalty for each position the gap is extended. GCCTATTG GCCTATTG | ||| | ||| A-C-TTTG AC--TTTG Using a scoring scheme of match (a), mismatch (-b), gap open (-i), gap extend (-e), the alignment scores become: 4a - 2b - 2i and 4a - 2b - i - e

31. BLAST (www.ncbi.nlm.nih.gov/BLAST)Basic Local Alignment Search Tool Main idea: Good alignments are likely to have two or more short (3+ residues) high-scoring words. Observation: Similarities of interest are usually longer than a single word, so look for multiple hits on the same diagonal, separated by a short distance

32. BLAST: basic strategy • 1. Find pairs of high scoring words (>T), separated by no more than A positions, on the same diagonal (ie, no gaps). • 2. Ungapped extension (fast) • 3. If extended HSP scores >Sg, then • 4. Perform gapped extension (slow) • 5. Report if “E-val” is lower than c • T, A, c, and Sg are adjustable constants

33. Why do BLOSUM matrices seem to outperform Dayhoff matrices? • Good guess: Dayhoff matrices are based on changes seen among closely related sequences. • Searches tend to target ancient homologies. • Different sites in proteins evolve with different rates and patterns. • The remaining sequence similarities after long periods of time are likely to be in functionally constrained regons

34. Why do BLOSUM matrices seem to outperform Dayhoff matrices? • Changes counted by Dayhoff methods may tend to be in unconstrained, quickly evolving regions

35. Statistical Significance: E-values Let Pi be the frequency of aa i. For unrelated sequences, an alignment of i with j has probability Pij. Given Pij and sij, we can calculate normalized scores (“bit scores”) from the raw score, S: K is a function of the database size,  is a function of the scoring matrix

36. Statistical Significance: E-values When 2 random sequences of length m and n are aligned, the expected number of HSPs with normalized scores greater than S´ is approximately

37. Statistical Significance: E-values

38. Statistical Significance: E-values • The preceding theory is strictly correct for ungapped alignments only. • Empirical observations suggest that it may apply approximately to gapped alignments as well.

39. Ramble about deriving distribution of BLAST scores if time allows.

40. How many alignments? Consider aligning 2 sequences of length n and m. , where [Perspective: approximately 1080 particles in the universe]

41. Dynamic Programming Alignment of 2 Sequences Consider the global alignment of sequences a and b. Define d(x,y) = cost of aligning an x with a y. d(A,A) = 5 match score d(A,T) = -2 mismatch penalty d(A,–) = -3 gap penalty

42. We can think of alignment as a series of decisions. At each step, we decide among three possible alternatives: • Add the next residues from both sequences (match or mismatch) • Add the next residue from a (gap in b) • Add the next residue from b (gap in a)

43. A–CC ATCC A–CC ATCC The “best” alignment is the series of decisions that maximizes the total score:

44. Let Dij be the score of aligning the first i residues of a with the first j residues of b. Then,

45. Match:+5 Mismatch: -2 Gap: -3 0 A G C T T A 0 -6 -9 -12 -15 -18 G 2 -1 -4 -7 -10 C -6 -1 7 4 -1 -4 T -9 -4 4 12 9 6 G -12 -11 -3 1 9 10 7 A -15 -7 -6 -2 6 7 15 0 -3 -3 -2 -5 -8

46. The optimal global alignments is the following: A: AGCTTA ||| | B: -GCTGA It has a score of 15: 4 matches, 1 mismatches, and 1 gap. (20-2-3=15)

47. Align a=ACTCG with b=ACCTG. Use match score +4, mismatch penalty -1, and gap penalty -2. 0 A C C T G 0 A C T C G