Sequence Alignment II

1. Sequence Alignment II K-tuple methodsStatistics of alignments

2. Database searches • What is the problem? • Large number of sequences to search your query sequence against. • Various indexing schemes and heuristics are used, one of which is BLAST. • heuristic is a technique to solve a problem that ignores whether the solution can be proven to be correct, but usually produces a good solution, are intended to gain computational performance or conceptual simplicity potentially at the cost of accuracy or precision. http://en.wikipedia.org/wiki/Heuristics#Computer_science

3. K-tuple methods http://creativecommons.org/licenses/by-sa/2.0/

4. Concepts of Sequence Similarity Searching • The premise: • The sequence itself is not informative; it must be analyzed by comparative methods against existing databases to develop hypothesis concerning relatives and function.

5. Important Terms for Sequence Similarity Searching with very different meanings • Similarity • The extent to which nucleotide or protein sequences are related. In BLAST similarity refers to a positive matrix score. • Identity • The extent to which two (nucleotide or amino acid) sequences are invariant. • Homology • Similarity attributed to descent from a common ancestor.

6. Sequence Similarity Searching: The Approach • Sequence similarity searching involves the use of a set of algorithms (such as the BLAST programs) to compare a query sequence to all the sequences in a specified database. • Comparisons are made in a pairwise fashion. Each comparison is given a score reflecting the degree of similarity between the query and the sequence being compared.

7. Blast QUERY sequence(s) BLAST results BLAST program BLAST database

8. Topics: BLAST program • There are different blast programs • Understanding the BLAST algorithm • Word size • HSPs (High Scoring Pairs) • Understanding BLAST statistics • The alignment score (S) • Scoring Matrices • Dealing with gaps in an alignment • The expectation value (E)

9. The BLAST algorithm • The BLAST programs (Basic Local Alignment Search Tools) are a set of sequence comparison algorithms introduced in 1990 for optimal local alignments to a query. • Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ (1990) “Basic local alignment search tool.” J. Mol. Biol. 215:403-410. • Altschul SF, Madden TL, Schaeffer AA, Zhang J, Zhang Z, Miller W, Lipman DJ (1997) “Gapped BLAST and PSI-BLAST: a new generation of protein database search programs.” NAR 25:3389-3402.

10. blastp blastx tblastn tblastx http://www.ncbi.nlm.nih.gov/BLAST blastn

11. Other BLAST programs • BLAST 2 Sequences (bl2seq) • Aligns two sequences of your choice • Gives dot-plot like output

12. More BLAST programs • BLAST against genomes • Many available • BLAST parameters pre-optimized • Handy for mapping query to genome • Search for short exact matches • BLAST parameters pre-optimized • Great for checking probes and primers

13. How Does BLAST Work? • The BLAST programs improved the overall speed of searches while retaining good sensitivity (important as databases continue to grow) by breaking the query and database sequences into fragments ("words"), and initially seeking matches between fragments. • Word hits are then extended in either direction in an attempt to generate an alignment with a score exceeding the threshold of “T".

14. Picture used with permission from Chapter 11 of “Bioinformatics: A Practical Guide to the Analysis of Genes and Proteins”

15. Each BLAST “hit” generates an alignment that can contain one or more high scoring pairs (HSPs)

16. Each BLAST “hit” generates an alignment that can contain one or more high scoring pairs (HSPs)

17. Where does the score (S) come from? • The quality of each pair-wise alignment is represented as a score and the scores are ranked. • Scoring matrices are used to calculate the score of the alignment base by base (DNA) or amino acid by amino acid (protein). • The alignment score will be the sum of the scores for each position.

18. What’s a scoring matrix? • Substitution matrices are used for amino acid alignments. These are matrices in which each possible residue substitution is given a score reflecting the probability that it is related to the corresponding residue in the query.

19. PAM vs. BLOSUM scoring matrices • BLOSUM 62 is the default matrix in BLAST 2.0. Though it is tailored for comparisons of moderately distant proteins, it performs well in detecting closer relationships. A search for distant relatives may be more sensitive with a different matrix.

20. The PAM Family PAM matrices are based on global alignments of closely related proteins. The PAM1 is the matrix calculated from comparisons of sequences with no more than 1% divergence. Other PAM matrices are extrapolated from PAM1. The BLOSUM family BLOSUM matrices are based on local alignments. BLOSUM 62 is a matrix calculated from comparisons of sequences with no less than 62% divergence. All BLOSUM matrices are based on observed alignments; they are not extrapolated from comparisons of closely related proteins. PAM vs BLOSUM scoring matrices

21. What happens if you have a gap in the alignment? • A gap is a position in the alignment at which a letter is paired with a null • Gap scores are negative. Since a single mutational event may cause the insertion or deletion of more than one residue, the presence of a gap is frequently ascribed more significance than the length of the gap. • Hence the gap is penalized heavily, whereas a lesser penalty is assigned to each subsequent residue in the gap.

22. Percent Sequence Identity • The extent to which two nucleotide or amino acid sequences are invariant A C C T G A G – A G A C G T G – G C A G mismatch indel 70% identical

23. BLAST algorithm • Keyword search of all words of length w in the query of default length n in database of length m with score above threshold • w = 11 for nucleotide queries, 3 for proteins • Do local alignment extension for each hit of keyword search • Extend result until longest match above threshold is achieved and output

24. BLAST algorithm (cont’d) keyword Query: KRHRKVLRDNIQGITKPAIRRLARRGGVKRISGLIYEETRGVLKIFLENVIRD GVK 18 GAK 16 GIK 16 GGK 14 GLK 13 GNK 12 GRK 11 GEK 11 GDK 11 Neighborhood words neighborhood score threshold (T = 13) extension Query: 22 VLRDNIQGITKPAIRRLARRGGVKRISGLIYEETRGVLK 60 +++DN +G + IR L G+K I+ L+ E+ RG++K Sbjct: 226 IIKDNGRGFSGKQIRNLNYGIGLKVIADLV-EKHRGIIK 263 High-scoring Pair (HSP)

25. Local alignment • Find the best local alignment between two strings, over the recurrence:

26. Local alignment (cont’d) • Input: strings v and w and scoring matrix d • Output: substrings of v and w whose global alignment as defined by d, is maximal among all global alignments of all substrings of v and w

27. Original BLAST • Dictionary • All words of length w • Alignment • Ungapped extensions until score falls below statistical threshold T • Output • All local alignments with score > statistical threshold

28. Original BLAST: Example A C G A A G T A A G G T C C A G T • w = 4, T = 4 • Exact keyword match of GGTC • Extend diagonals with mismatches until score falls below a threshold • Output result • GTAAGGTCC • GTTAGGTCC C T G A T C C T G G A T T G C G A From lectures by Serafim Batzoglou (Stanford)

29. Gapped BLAST: Example A C G A A G T A A G G T C C A G T • Original BLAST exact keyword search, THEN: • Extend with gaps in a zone around ends of exact match • Output result GTAAGGTCCAGT GTTAGGTC-AGT C T G A T C C T G G A T T G C G A From lectures by Serafim Batzoglou (Stanford)

30. Gapped BLAST : Example (cont’d) A C G A A G T A A G G T C C A G T • Original BLAST exact keyword search, THEN: • Extend with gaps around ends of exact matchuntil score <T, then merge nearby alignments • Output result GTAAGGTCCAGT GTTAGGTC-AGT C T G A T C C T G G A T T G C G A From lectures by Serafim Batzoglou (Stanford)

31. Topics: BLAST databases • The different blast databases provided by the NCBI • Protein databases • Nucleotide databases • Genomic databases • Considerations for choosing a BLAST database • Custom databases for BLAST

32. blastp db BLAST protein databases available at through blastp web interface @ NCBI

33. Considerations for choosing a BLAST database • First consider your research question: • Are you looking for an ortholog in a particular species? • BLAST against the genome of that species. • Are you looking for additional members of a protein family across all species? • BLAST against nr, if you can’t find hits check wgs, htgs, and the trace archives. • Are you looking to annotate genes in your species of interest? • BLAST against known genes (RefSeq) and/or ESTs from a closely related species.

34. When choosing a database for BLAST… • It is important to know your reagents. • Changing your choice of database is changing your search space completely • Database size affects the BLAST statistics • record BLAST parameters, database choice, database size in your bioinformatics lab book, just as you would for your wet-bench experiments. • Databases change rapidly and are updated frequently • It may be necessary to repeat your analyses

35. Topics: BLAST results • Choosing the right BLAST program • Running a blastp search • BLAST parameters and options to consider • Viewing BLAST results • Look at your alignments • Using the BLAST taxonomy report

36. BLAST parameters and options to consider: conserved domains Entrez query E-value cutoff Word size

37. More BLAST parameters and options to consider: filtering gap penalities matrix

38. Run your BLAST search: BLAST

39. The BLAST Queue: click for more info Note your RID

40. Formatting and Retrieving your BLAST results: Results options

41. A graphical view of your BLAST results:

42. The BLAST “hit” list: Score E-Value GenBank alignment EntrezGene

43. The BLAST pairwise alignments Identity Similarity

44. Sample BLAST output • Blast of human beta globin protein against zebra fish Score E Sequences producing significant alignments: (bits) Value gi|18858329|ref|NP_571095.1| ba1 globin [Danio rerio] >gi|147757... 171 3e-44 gi|18858331|ref|NP_571096.1| ba2 globin; SI:dZ118J2.3 [Danio rer... 170 7e-44 gi|37606100|emb|CAE48992.1| SI:bY187G17.6 (novel beta globin) [D... 170 7e-44 gi|31419195|gb|AAH53176.1| Ba1 protein [Danio rerio] 168 3e-43 ALIGNMENTS >gi|18858329|ref|NP_571095.1| ba1 globin [Danio rerio] Length = 148 Score = 171 bits (434), Expect = 3e-44 Identities = 76/148 (51%), Positives = 106/148 (71%), Gaps = 1/148 (0%) Query: 1 MVHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESFGDLSTPDAVMGNPK 60 MV T E++A+ LWGK+N+DE+G +AL R L+VYPWTQR+F +FG+LS+P A+MGNPK Sbjct: 1 MVEWTDAERTAILGLWGKLNIDEIGPQALSRCLIVYPWTQRYFATFGNLSSPAAIMGNPK 60 Query: 61 VKAHGKKVLGAFSDGLAHLDNLKGTFATLSELHCDKLHVDPENFRLLGNVLVCVLAHHFG 120 V AHG+ V+G + ++DN+K T+A LS +H +KLHVDP+NFRLL + + A FG Sbjct: 61 VAAHGRTVMGGLERAIKNMDNVKNTYAALSVMHSEKLHVDPDNFRLLADCITVCAAMKFG 120 Query: 121 KE-FTPPVQAAYQKVVAGVANALAHKYH 147 + F VQ A+QK +A V +AL +YH Sbjct: 121 QAGFNADVQEAWQKFLAVVVSALCRQYH 148

45. Sample BLAST output (cont’d) • Blast of human beta globin DNA against human DNA Score E Sequences producing significant alignments: (bits) Value gi|19849266|gb|AF487523.1| Homo sapiens gamma A hemoglobin (HBG1... 289 1e-75 gi|183868|gb|M11427.1|HUMHBG3E Human gamma-globin mRNA, 3' end 289 1e-75 gi|44887617|gb|AY534688.1| Homo sapiens A-gamma globin (HBG1) ge... 280 1e-72 gi|31726|emb|V00512.1|HSGGL1 Human messenger RNA for gamma-globin 260 1e-66 gi|38683401|ref|NR_001589.1| Homo sapiens hemoglobin, beta pseud... 151 7e-34 gi|18462073|gb|AF339400.1| Homo sapiens haplotype PB26 beta-glob... 149 3e-33 ALIGNMENTS >gi|28380636|ref|NG_000007.3| Homo sapiens beta globin region (HBB@) on chromosome 11 Length = 81706 Score = 149 bits (75), Expect = 3e-33 Identities = 183/219 (83%) Strand = Plus / Plus Query: 267 ttgggagatgccacaaagcacctggatgatctcaagggcacctttgcccagctgagtgaa 326 || ||| | || | || | |||||| ||||| ||||||||||| |||||||| Sbjct: 54409 ttcggaaaagctgttatgctcacggatgacctcaaaggcacctttgctacactgagtgac 54468 Query: 327 ctgcactgtgacaagctgcatgtggatcctgagaacttc 365 ||||||||| |||||||||| ||||| |||||||||||| Sbjct: 54469 ctgcactgtaacaagctgcacgtggaccctgagaacttc 54507

46. What do the Score and the e-value really mean? • The quality of the alignment is represented by the Score. • Score (S) • The score of an alignment is calculated as the sum of substitution and gap scores. Substitution scores are given by a look-up table (PAM, BLOSUM) whereas gap scores are assigned empirically . • The significance of each alignment is computed as an E value. • E value (E) • Expectation value. The number of different alignments with scores equivalent to or better than S that are expected to occur in a database search by chance. The lower the E value, the more significant the score.

47. E value • E value (E) • Expectation value. The number of different alignments with scores equivalent to or better than S expected to occur in a database search by chance. The lower the E value, the more significant the score.

48. Assessing sequence homology • Need to know how strong an alignment can be expected from chance alone • “Chance” is the comparison of • Real but non-homologous sequences • Real sequences that are shuffled to preserve compositional properties • Sequences that are generated randomly based upon a DNA or protein sequence model (favored)

49. High Scoring Pairs (HSPs) • All segment pairs whose scores can not be improved by extension or trimming • Need to model a random sequence to analyze how high the score is in relation to chance

50. Expected number of HSPs • Expected number of HSPs with score >S • E-value E for the score S: • E = Kmne-lS • Given: • Two sequences, length n and m • The statistics of HSP scores are characterized by two parameters K and λ • K: scale for the search space size • λ: scale for the scoring system