1 / 41

Lecture 2: Sequence Alignment BMI/IBGP 730

Lecture 2: Sequence Alignment BMI/IBGP 730. Kun Huang Department of Biomedical Informatics Ohio State University. Major issues in genomics Homology Format Search Alignment as an optimization problem Dynamical programming Hashing and BLAST Tools and examples. Charles Darwin

quanda
Télécharger la présentation

Lecture 2: Sequence Alignment BMI/IBGP 730

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. Lecture 2: Sequence Alignment BMI/IBGP 730 Kun Huang Department of Biomedical Informatics Ohio State University

  2. Major issues in genomics • Homology • Format • Search • Alignment as an optimization problem • Dynamical programming • Hashing and BLAST • Tools and examples

  3. Charles Darwin (1809-1882) “I think …”

  4. Homology A Working Definition: Sequences or structures which share a common ancestor

  5. Homology Defined "The same organ in different animals under a variety of form and function." Sir Richard Owen, Lectures on the Comparative Anatomy and Physiology of the Invertebrate Animals, 1843. "The mechanism of homology is heredity." Allan Boyden, Homology and Analogy: A century after the definitions of "homologue" and "analogue" of Richard Owen, 1943. "Homology is a relation bearing on recency of common ancestry.“ Olivier Rieppel, Homology and logical fallacy, 1993.

  6. Sequence Homology • Genes in separate species derived from the same ancestral genes in the last common ancestor of those two species are orthologs. • Related genes resulted from a gene duplication event within a single genome--and are likely to have diverged in their function--are paralogs. • Both orthologs and paralogs are homologs, a general term to cover both types of relationships

  7. Recognizing Sequence Homology • Relies primarily on understanding random sequence similarity • Only by knowing what random similarity looks like can we tell when two sequences are significantly similar • Understanding mutational regularity and sequence evolution increases the significance • Closely-related: Transitions/transversions • Distantly-related: PAM mutation probabilities • Even distantly-related sequences can be recognized • "Significant Similarity" is not a definition of homology.

  8. Functional Homology Biochemical Function What the protein does on a biochemical level, almost always conserved between homologues. Example: Serine Kinase Physiological Function What role the protein plays within the organism, only conserved for orthologous proteins. Example: Glycogen Synthase Kinase

  9. Databases • GenBank • EMBL • DDBJ • SWISSPROT • …

  10. Data Format FASTA (.fasta file) >gi|33469954|ref|NM_000240.2| Homo sapiens monoamine oxidase A (MAOA), nuclear gene encoding mitochondrial protein, mRNAGGGCGCTCCCGGAGTATCAGCAAAAGGGTTCGCCCCGCCCACAGTGCCCGGCTCCCCCCGGGTATCAAAA GAAGGATCGGCTCCGCCCCCGGGCTCCCCGGGGGAGTTGATAGAAGGGTCCTTCCCACCCTTTGCCGTCC CCACTCCTGTGCCTACGACCCAGGAGCGTGTCAGCCAAAGCATGGAGAATCAAGAGAAGGCGAGTATCGC GGGCCACATGTTCGACGTAGTCGTGATCGGAGGTGGCATTTCAGGACTATCTGCTGCCAAACTCTTGACT GAATATGGCGTTAGTGTTTTGGTTTTAGAAGCTCGGGACAGGGTTGGAGGAAGAACATATACTATAAGGA ATGAGCATGTTGATTACGTAGATGTTGGTGGAGCTTATGTGGGACCAACCCAAAACAGAATCTTACGCTT GTCTAAGGAGCTGGGCATAGAGACTTACAAAGTGAATGTCAGTGAGCGTCTCGTTCAATATGTCAAGGGG AAAACATATCCATTTCGGGGCGCCTTTCCACCAGTATGGAATCCCATTGCATATTTGGATTACAATAATC TGTGGAGGACAATAGATAACATGGGGAAGGAGATTCCAACTGATGCACCCTGGGAGGCTCAACATGCTGA CAAATGGGACAAAATGACCATGAAAGAGCTCATTGACAAAATCTGCTGGACAAAGACTGCTAGGCGGTTT GCTTATCTTTTTGTGAATATCAATGTGACCTCTGAGCCTCACGAAGTGTCTGCCCTGTGGTTCTTGTGGT ATGTGAAGCAGTGCGGGGGCACCACTCGGATATTCTCTGTCACCAATGGTGGCCAGGAACGGAAGTTTGT … >gi|4557735|ref|NP_000231.1| monoamine oxidase A [Homo sapiens] MENQEKASIAGHMFDVVVIGGGISGLSAAKLLTEYGVSVLVLEARDRVGGRTYTIRNEHVDYVDVGGAYV GPTQNRILRLSKELGIETYKVNVSERLVQYVKGKTYPFRGAFPPVWNPIAYLDYNNLWRTIDNMGKEIPT DAPWEAQHADKWDKMTMKELIDKICWTKTARRFAYLFVNINVTSEPHEVSALWFLWYVKQCGGTTRIFSV TNGGQERKFVGGSGQVSERIMDLLGDQVKLNHPVTHVDQSSDNIIIETLNHEHYECKYVINAIPPTLTAK IHFRPELPAERNQLIQRLPMGAVIKCMMYYKEAFWKKKDYCGCMIIEDEDAPISITLDDTKPDGSLPAIM GFILARKADRLAKLHKEIRKKKICELYAKVLGSQEALHPVHYEEKNWCEEQYSGGCYTAYFPPGIMTQYG RVIRQPVGRIFFAGTETATKWSGYMEGAVEAGERAAREVLNGLGKVTEKDIWVQEPESKDVPAVEITHTF WERNLPSVSGLLKIIGFSTSVTALGFVLYKYKLLPRS

  11. Data Format • Other formats • NBRF/PIR (.pir file) • Begin with “>P1;” for protein sequence and “>N1;” for nucleotide. • GDE (.gde file) • Similar to FASTA file, begin with “%” instead of “>”.

  12. Major issues in genomics • Homology • Format • Search • Alignment as an optimization problem • Dynamical programming • BLAST • Tools and examples

  13. Aligning Text Strings T C A T G C A T T G 2 matches, 0 gaps T C A T G | | C A T T G 4 matches, 1 inserted gaps T C A - T G | | | | C A T T G 3 matches, 2 end gaps T C A T G | | | C A T T G 4 matches, 1 inserted gaps T C A T - G | | | | C A T T G Optimal solution: with respect to what criteria / cost function?

  14. Alignment as An Optimization Problem • Optimization criteria / cost function • Parameters to be adjusted • Search algorithm / process • Exhaustive testing • Suboptimal solutions • Computational cost / complexity • Statistical significance

  15. Alignment as An Optimization Problem • Optimization criteria / cost function • What sort of alignment should be considered • Scoring system • Additive model • Based on probability compared with random sequence (PAM, BLOSUM) • Assumption of independence • More complicated cases • Gap penalty – linear (s = -gd) • affine (s = -d – (g-1)e)

  16. Alignment as An Optimization Problem • Optimization criteria / cost function • What sort of alignment should be considered • Scoring system • Additive model • Based on probability compared with random sequence (PAM, BLOSUM) • Assumption of independence • More complicated cases • Gap penalty – linear (s = -gd) • affine (s = -d – (g-1)e)

  17. Alignment as An Optimization Problem • Parameters to be adjusted • Shift • Number of gaps • Position of gaps 3 matches, 2 end gaps T C A T G | | | C A T T G 4 matches, 1 inserted gaps T C A - T G | | | | C A T T G

  18. Alignment as An Optimization Problem • Search algorithm / process • Exhaustive testing • Try all possible configuration of parameters. • E.g., sequence a with length m, sequence b with length n. Try all m+n shifts. 0 matches, 0 gaps T C A T G C A T T G 3 matches, 2 end gaps T C A T G | | | C A T T G 2 matches, 0 gaps T C A T G | | C A T T G

  19. Alignment as An Optimization Problem • Search algorithm / process • Computational cost / complexity • Big O notation – running time (number of basic computations) and memory usage • O(m+n) running time for the above example. It means the running time is asymptotically proportional to m+n. 0 matches, 0 gaps T C A T G C A T T G 3 matches, 2 end gaps T C A T G | | | C A T T G 2 matches, 0 gaps T C A T G | | C A T T G

  20. Alignment as An Optimization Problem • Search algorithm / process • Computational cost / complexity • What if we allow gaps? 0 matches, 0 gaps T C A T G C A T T G 3 matches, 2 end gaps T C A T G | | | C A T T G 2 matches, 0 gaps T C A T G | | C A T T G

  21. Many possible alignments to consider • Without gaps, there are are n+m possible alignments between sequences of length n and m • Once we start allowing gaps, there are many possible arrangements to consider:abcbcd abcbcd abcbcd | | | | | | | | || || abc--d a--bcd ab--cd • This becomes a very large number when we allow mismatches, since we then need to look at every possible pairing between elements: there are roughly nm possible alignments.

  22. Exponential computations get big fast • If n=m=100, there are 100100 = 10200 = 100,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000 different alignments. • And 100 amino acids is a small protein!

  23. Alignment as An Optimization Problem • Statistical significance • Not only are there many possible gapped alignments, but introducing too many gaps makes nonsense alignments possible:s--e-----qu---en--ce sometimesquipsentice • Need to distinguish between alignments that occur due to homology, and those that could be expected to be seen just by chance. • Define a score function that accounts for both element mismatches and a gap penalty

  24. Major issues in genomics • Homology • Format • Search • Alignment as an optimization problem • Dynamical programming • BLAST • Tools and examples

  25. Dot matrix sequence comparison • Write one sequence across top of matrix, the other across left side, then put a dot where character on line i equals one in column j • Examples below: DNA and amino acid sequences of the phage  cI (vertical axis) and phage P22 c2 (horizontal axis) repressors

  26. Dynamic programming • The name comes from an operations research task, and has nothing to do with writing programs. Programming – use tabular structure for computing. • The key idea is to start aligning the sequences left to right; once a prefix is optimally aligned, nothing about the remainder of the alignment changes the alignment of the prefix. • We construct a matrix of possible alignment scores (nxm2 calculations worst case) and then "traceback" to find the optimal alignment. • Called Needleman-Wunch (for global matching) or Smith-Waterman (for local matching)

  27. 18 3 20 5 25 B 2 A 11 Dynamic programming • The name comes from an operations research task, and has nothing to do with writing programs. Programming – use tabular structure for computing. 17

  28. Dynamic programming matrix • Each cell has the score for the best aligned sequenceprefix up to that position. • Example: ATGCT vs. ACCGCT Match: +2, mismatch: 0, gap: -1 Matching matrix, NOT the dynamical programming matrix!

  29. A T _ A _C A T A C A T A C A _T A C_ Dynamic programming matrix 0 1 2 3 4 5 6 0 1 2 3 4 5 6 7 A A A _ A T A T A _ A T A C

  30. Optimal alignment by traceback • We “traceback” a path that gets us the highest score. If we don't have “end gap” penalties, then takeany path from thelast row or columnto the first. • Otherwise we needto include the top and bottom corners 0 1 2 3 4 5 6 0 1 2 3 4 5 6 7 AT - GCT ACCGCT A - TGCT ACCGCT

  31. Dynamic programming • Global alignment – an alignment of two or more sequences that matches as many characters as possible in all of the sequences. Needleman-Wunch algorithm • Local alignment – an alignment that includes only the best matching, highest-scoring regions in two or more sequences. Smith-Waterman algorithm • Difference – all the scores are kept in the dynamical programming matrix for global alignment; only the positive scores are kept in the dynamical programming matrix for local alignment, the negative ones are converted to zero.

  32. Major issues in genomics • Homology • Annotation • Database • Search • Alignment as an optimization problem • Dynamical programming • BLAST • Tools and examples

  33. Sequence alignment (BLAST) The Basic Local Alignment Search Tool (BLAST) finds regions of local similarity between sequences. The program compares nucleotide or protein sequences to sequence databases and calculates the statistical significance of matches. BLAST can be used to infer functional and evolutionary relationships between sequences as well as help identify members of gene families. http://www.ncbi.nlm.nih.gov/blast/

  34. BLAST – Algorithm Intuition The BLAST algorithm.The BLAST algorithm is a heuristic search method that seeks words of length W (default = 3 in blastp) that score at least T when aligned with the query and scored with a substitution matrix. Words in the database that score T or greater are extended in both directions in an attempt to find a locally optimal ungapped alignment or HSP (high scoring pair) with a score of at least S or an E value lower than the specified threshold. HSPs that meet these criteria will be reported by BLAST, provided they do not exceed the cutoff value specified for number of descriptions and/or alignments to report.

  35. BLAST – Algorithm Intuition Databases are pre-indexed by the words. Without gaps: Altschul, S. F., Gish, W., Miller, W., Myers, E. W., Lipman, D. J., J. Mol. Biol. (1990) 215:403-410 With gaps: Altschul, S. F., Madden, T. L., Schaffer, A. A., Zhang, J., Zhang, Z., Miller, W., Lipman, D. J., Nucleic Acids Research (1997) 25(17):3389-3402 http://www.compbio.dundee.ac.uk/ftp/preprints/review93/Figure10.pdf

  36. BLAST – Scoring Matrices DNA scoring matrix (substitution matrix) A T G C A 5 -4 -4 -4 T -4 5 -4 -4 G -4 -4 5 -4 C -4 -4 -4 5 ATTTAGCCG ACTTGGCCT 5 55 555 Score = 5X6+(-4)X3=18 http://www.ncbi.nlm.nih.gov/BLAST/matrix_info.html#matrix

  37. BLAST – Scoring Matrices • DNA is relatively easy to choose and protein is harder. • PAM (Percent Accepted Mutation) matrices: predicted matrices, most sensitive for alignments of sequences with evolutionary related homologs. The greater the number in the matrix name, the greater the expected evolutionary (mutational) distance, i.e. PAM30 would be used for alignments expected to be more closely related in evolution than an alignment performed using the PAM250 matrix. • BLOSUM (Blocks Substitution Matrix): calculated matrices, most sensitive for local alignment of related sequences, ideal when trying to identify an unknown nucleotide sequence. BLOSUM62 is the default matrix set be the BLAST search tool.

  38. BLAST – Parameters • Word size – for MegaBlast, it works best if w=16, can work to w=8. • Expected – statistical based notion, compare the matched sequence with random sequence (the likelihood). The larger the score, the smaller the expected value, the more significant the result. • Percent Identity, match/mismatch scores.

  39. BLAST – Program Selection • Nucleotide • Quickly search for highly similar sequences (megablast) • Quickly search for divergent sequences (discontiguous megablast) • Nucleotide-nucleotide BLAST (blastn) • Search for short, nearly exact matches • Search trace archives with megablast or discontiguous megablast • Protein • Protein-protein BLAST (blastp) • Position-specific iterated and pattern-hit initiated BLAST (PSI- and PHI-BLAST) • Search for short, nearly exact matches • Search the conserved domain database (rpsblast) • Protein homology by domain architecture (cdart)

  40. Exercises • Question 1 – Gene BLAST • Using MegaBLAST to BLAST the human Rb1 gene using the appropriate BLAST tool. Try the following three different settings: • WordSize 16, Percent Identity, match/mismatch score: None, 1, -2 • WordSize 64, Percent Identity, match/mismatch score: None, 1, -2 • WordSize 16, Percent Identity, match/mismatch score: 95, 1, -3 • Report the following: the total number of hits for each query; the first 20 distinctive matches (not including the query gene) for each query in the sequence they appeared in the hit table (GI number is enough). • Question 2 – Protein BLAST • BLAST the human Rb1 protein using BLASTP. Try the following three different scoring matrices with Word Size 3 and Gap Costs 11/1: PAM30, PAM70, BLOSUM62. Report the first 20 distinctive matches for each query. • Reading – BLAST algorithm • Altschul, S. F., Madden, T. L., Schaffer, A. A., Zhang, J., Zhang, Z., Miller, W., Lipman, D. J., Nucleic Acids Research (1997) 25(17):3389-3402

More Related