1 / 35

Sequence Alignment Cont’d

Sequence Alignment Cont’d. Linear-space alignment. Iterate this procedure to the left and right!. k *. N-k *. M/2. M/2. The Four-Russian Algorithm. Main structure of the algorithm: Divide N N DP matrix into K K log 2 N-blocks that overlap by 1 column & 1 row For i = 1……K

fedora
Download Presentation

Sequence Alignment Cont’d

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. Sequence AlignmentCont’d

  2. Linear-space alignment • Iterate this procedure to the left and right! k* N-k* M/2 M/2

  3. The Four-Russian Algorithm Main structure of the algorithm: • Divide NN DP matrix into KK log2N-blocks that overlap by 1 column & 1 row • For i = 1……K • For j = 1……K • Compute Di,j as a function of Ai,j, Bi,j, Ci,j, x[li…l’i], y[rj…r’j] Time: O(N2 / log2N) times the cost of step 4 t t t

  4. Heuristic Local AlignersBLAST, WU-BLAST, BlastZ, MegaBLAST, BLAT, PatternHunter, ……

  5. State of biological databases Sequenced Genomes: Human 3109 Yeast 1.2107 Mouse 2.7109  12 different strains Rat 2.6109 Neurospora 4107 14 more fungi within next year Fugu fish 3.3108 Tetraodon 3108 ~250 bacteria/viruses Mosquito 2.8108 Next year: Drosophila 1.2108 Dog, Chimpanzee, Chicken Worm 1.0108 2 sea squirts  1.6108Current rate of sequencing: Rice 1.0109 4 big labs  3 109 bp /year/lab Arabidopsis 1.2108 10s small labs

  6. State of biological databases • Number of genes in these genomes: • Vertebrate: ~30,000 • Insects: ~14,000 • Worm: ~17,000 • Fungi: ~6,000-10,000 • Small organisms: 100s-1,000s • Each known or predicted gene has an associated protein sequence • >1,000,000 known / predicted protein sequences

  7. Some useful applications of alignments • Given a newly discovered gene, • Does it occur in other species? • How fast does it evolve? • Assume we try Smith-Waterman: Our new gene 104 The entire genomic database 1010 - 1011

  8. Some useful applications of alignments • Given a newly sequenced organism, • Which subregions align with other organisms? • Potential genes • Other biological characteristics • Assume we try Smith-Waterman: Our newly sequenced mammal 3109 The entire genomic database 1010 - 1011

  9. BLAST (Basic Local Alignment Search Tool) Main idea: • Construct a dictionary of all the words in the query • Initiate a local alignment for each word match between query and DB Running Time: O(MN) However, orders of magnitude faster than Smith-Waterman query DB

  10. BLAST  Original Version …… Dictionary: All words of length k (~11) Alignment initiated between words of alignment score  T (typically T = k) Alignment: Ungapped extensions until score below statistical threshold Output: All local alignments with score > statistical threshold query …… scan DB query

  11. BLAST  Original Version A C G A A G T A A G G T C C A G T Example: k = 4, T = 4 The matching word GGTC initiates an alignment Extension to the left and right with no gaps until alignment falls < 50% Output: GTAAGGTCC GTTAGGTCC C C C T T C C T G G A T T G C G A

  12. Gapped BLAST A C G A A G T A A G G T C C A G T Added features: • Pairs of words can initiate alignment • Extensions with gaps in a band around anchor Output: GTAAGGTCCAGT GTTAGGTC-AGT C T G A T C C T G G A T T G C G A

  13. Gapped BLAST A C G A A G T A A G G T C C A G T Added features: • Pairs of words can initiate alignment • Nearby alignments are merged • Extensions with gaps until score < T below best score so far Output: GTAAGGTCCAGT GTTAGGTC-AGT C T G A T C C T G G A T T G C G A

  14. Variants of BLAST • MEGABLAST: • Optimized to align very similar sequences • Works best when k = 4i  16 • Linear gap penalty • PSI-BLAST: • BLAST produces many hits • Those are aligned, and a pattern is extracted • Pattern is used for next search; above steps iterated • WU-BLAST: (Wash U BLAST) • Optimized, added features • BlastZ • Combines BLAST/PatternHunter methodology

  15. Example Query:gattacaccccgattacaccccgattaca (29 letters) [2 mins] Database: All GenBank+EMBL+DDBJ+PDB sequences (but no EST, STS, GSS, or phase 0, 1 or 2 HTGS sequences) 1,726,556 sequences; 8,074,398,388 total letters >gi|28570323|gb|AC108906.9|Oryza sativa chromosome 3 BAC OSJNBa0087C10 genomic sequence, complete sequence Length = 144487 Score = 34.2 bits (17), Expect = 4.5 Identities = 20/21 (95%) Strand = Plus / Plus Query: 4 tacaccccgattacaccccga 24 ||||||| ||||||||||||| Sbjct: 125138 tacacccagattacaccccga 125158 Score = 34.2 bits (17), Expect = 4.5 Identities = 20/21 (95%) Strand = Plus / Plus Query: 4 tacaccccgattacaccccga 24 ||||||| ||||||||||||| Sbjct: 125104 tacacccagattacaccccga 125124 >gi|28173089|gb|AC104321.7| Oryza sativa chromosome 3 BAC OSJNBa0052F07 genomic sequence, complete sequence Length = 139823 Score = 34.2 bits (17), Expect = 4.5 Identities = 20/21 (95%) Strand = Plus / Plus Query: 4 tacaccccgattacaccccga 24 ||||||| ||||||||||||| Sbjct: 3891 tacacccagattacaccccga 3911

  16. Example Query: Human atoh enhancer, 179 letters [1.5 min] Result: 57 blast hits • gi|7677270|gb|AF218259.1|AF218259 Homo sapiens ATOH1 enhanc... 355 1e-95 • gi|22779500|gb|AC091158.11| Mus musculus Strain C57BL6/J ch... 264 4e-68 • gi|7677269|gb|AF218258.1|AF218258 Mus musculus Atoh1 enhanc... 256 9e-66 • gi|28875397|gb|AF467292.1| Gallus gallus CATH1 (CATH1) gene... 78 5e-12 • gi|27550980|emb|AL807792.6| Zebrafish DNA sequence from clo... 54 7e-05 • gi|22002129|gb|AC092389.4| Oryza sativa chromosome 10 BAC O... 44 0.068 • gi|22094122|ref|NM_013676.1| Mus musculus suppressor of Ty ... 42 0.27 • gi|13938031|gb|BC007132.1| Mus musculus, Similar to suppres... 42 0.27 gi|7677269|gb|AF218258.1|AF218258 Mus musculus Atoh1 enhancer sequence Length = 1517 Score = 256 bits (129), Expect = 9e-66 Identities = 167/177 (94%), Gaps = 2/177 (1%) Strand = Plus / Plus Query: 3 tgacaatagagggtctggcagaggctcctggccgcggtgcggagcgtctggagcggagca 62 ||||||||||||| ||||||||||||||||||| |||||||||||||||||||||||||| Sbjct: 1144 tgacaatagaggggctggcagaggctcctggccccggtgcggagcgtctggagcggagca 1203 Query: 63 cgcgctgtcagctggtgagcgcactctcctttcaggcagctccccggggagctgtgcggc 122 |||||||||||||||||||||||||| ||||||||| |||||||||||||||| ||||| Sbjct: 1204 cgcgctgtcagctggtgagcgcactc-gctttcaggccgctccccggggagctgagcggc 1262 Query: 123 cacatttaacaccatcatcacccctccccggcctcctcaacctcggcctcctcctcg 179 ||||||||||||| || ||| |||||||||||||||||||| ||||||||||||||| Sbjct: 1263 cacatttaacaccgtcgtca-ccctccccggcctcctcaacatcggcctcctcctcg 1318

  17. BLAT: Blast-Like Alignment Tool Differences with BLAST: • Dictionary of the DB (BLAST: query) • Initiate alignment on any # of matching hits (BLAST: 1 or 2 hits) • More aggressive stitching-together of alignments • Special code to align mature mRNA to DNA Result: very efficient for: • Aligning many seqs to a DB • Aligning two genomes

  18. PatternHunter Main features: • Non-consecutive position words • Highly optimized Consecutive Positions Non-Consecutive Positions 6 hits 7 hits 5 hits 7 hits 3 hits 3 hits On a 70% conserved region: ConsecutiveNon-consecutive Expected # hits: 1.07 0.97 Prob[at least one hit]: 0.30 0.47

  19. Advantage of Non-Consecutive Words 11 positions 11 positions 10 positions

  20. 1 2 2 1 1 1 1 … 2 2 2 2 … K … … … … x1 K K K K x2 x3 xK … Hidden Markov Models

  21. Outline for our next topic • Hidden Markov models – the theory • Probabilistic interpretation of alignments using HMMs Later in the course: • Applications of HMMs to biological sequence modeling and discovery of features such as genes

  22. Example: The Dishonest Casino A casino has two dice: • Fair die P(1) = P(2) = P(3) = P(5) = P(6) = 1/6 • Loaded die P(1) = P(2) = P(3) = P(5) = 1/10 P(6) = 1/2 Casino player switches back-&-forth between fair and loaded die once every 20 turns Game: • You bet $1 • You roll (always with a fair die) • Casino player rolls (maybe with fair die, maybe with loaded die) • Highest number wins $2

  23. Question # 1 – Evaluation GIVEN A sequence of rolls by the casino player 1245526462146146136136661664661636616366163616515615115146123562344 QUESTION How likely is this sequence, given our model of how the casino works? This is the EVALUATION problem in HMMs

  24. Question # 2 – Decoding GIVEN A sequence of rolls by the casino player 1245526462146146136136661664661636616366163616515615115146123562344 QUESTION What portion of the sequence was generated with the fair die, and what portion with the loaded die? This is the DECODING question in HMMs

  25. Question # 3 – Learning GIVEN A sequence of rolls by the casino player 1245526462146146136136661664661636616366163616515615115146123562344 QUESTION How “loaded” is the loaded die? How “fair” is the fair die? How often does the casino player change from fair to loaded, and back? This is the LEARNING question in HMMs

  26. The dishonest casino model 0.05 0.95 0.95 FAIR LOADED P(1|F) = 1/6 P(2|F) = 1/6 P(3|F) = 1/6 P(4|F) = 1/6 P(5|F) = 1/6 P(6|F) = 1/6 P(1|L) = 1/10 P(2|L) = 1/10 P(3|L) = 1/10 P(4|L) = 1/10 P(5|L) = 1/10 P(6|L) = 1/2 0.05

  27. Definition of a hidden Markov model Definition: A hidden Markov model (HMM) • Alphabet = { b1, b2, …, bM } • Set of states Q = { 1, ..., K } • Transition probabilities between any two states aij = transition prob from state i to state j ai1 + … + aiK = 1, for all states i = 1…K • Start probabilities a0i a01 + … + a0K = 1 • Emission probabilities within each state ei(b) = P( xi = b | i = k) ei(b1) + … + ei(bM) = 1, for all states i = 1…K 1 2 K …

  28. A HMM is memory-less At each time step t, the only thing that affects future states is the current state t P(t+1 = k | “whatever happened so far”) = P(t+1 = k | 1, 2, …, t, x1, x2, …, xt) = P(t+1 = k | t) 1 2 K …

  29. 1 1 1 1 … 2 2 2 2 … … … … … K K K K … A parse of a sequence Given a sequence x = x1……xN, A parse of x is a sequence of states  = 1, ……, N 1 2 2 K x1 x2 x3 xK

  30. 1 2 2 1 1 1 1 … 2 2 2 2 … K … … … … x1 K x2 K x3 K K xK … Likelihood of a parse Given a sequence x = x1……xN and a parse  = 1, ……, N, To find how likely is the parse: (given our HMM) P(x, ) = P(x1, …, xN, 1, ……, N) = P(xN, N | N-1) P(xN-1, N-1 | N-2)……P(x2, 2 | 1) P(x1, 1) = P(xN | N)P(N | N-1) ……P(x2 | 2)P(2 | 1)P(x1 | 1)P(1) = a01 a12……aN-1Ne1(x1)……eN(xN)

  31. Example: the dishonest casino Let the sequence of rolls be: x = 1, 2, 1, 5, 6, 2, 1, 6, 2, 4 Then, what is the likelihood of  = Fair, Fair, Fair, Fair, Fair, Fair, Fair, Fair, Fair, Fair? (say initial probs a0Fair = ½, aoLoaded = ½) ½  P(1 | Fair) P(Fair | Fair) P(2 | Fair) P(Fair | Fair) … P(4 | Fair) = ½  (1/6)10  (0.95)9 = .00000000521158647211 = 0.5  10-9

  32. Example: the dishonest casino So, the likelihood the die is fair in all this run is just 0.521  10-9 OK, but what is the likelihood of = Loaded, Loaded, Loaded, Loaded, Loaded, Loaded, Loaded, Loaded, Loaded, Loaded? ½  P(1 | Loaded) P(Loaded, Loaded) … P(4 | Loaded) = ½  (1/10)8  (1/2)2 (0.95)9 = .00000000078781176215 = 7.9  10-10 Therefore, it is after all 6.59 times more likely that the die is fair all the way, than that it is loaded all the way

  33. Example: the dishonest casino Let the sequence of rolls be: x = 1, 6, 6, 5, 6, 2, 6, 6, 3, 6 Now, what is the likelihood  = F, F, …, F? ½  (1/6)10  (0.95)9 = 0.5  10-9, same as before What is the likelihood  = L, L, …, L? ½  (1/10)4  (1/2)6 (0.95)9 = .00000049238235134735 = 0.5  10-7 So, it is 100 times more likely the die is loaded

  34. The three main questions on HMMs • Evaluation GIVEN a HMM M, and a sequence x, FIND Prob[ x | M ] • Decoding GIVEN a HMM M, and a sequence x, FIND the sequence  of states that maximizes P[ x,  | M ] • Learning GIVEN a HMM M, with unspecified transition/emission probs., and a sequence x, FIND parameters  = (ei(.), aij) that maximize P[ x |  ]

  35. Let’s not be confused by notation P[ x | M ]: The probability that sequence x was generated by the model The model is: architecture (#states, etc) + parameters  = aij, ei(.) So, P[ x |  ], and P[ x ] are the same, when the architecture, and the entire model, respectively, are implied Similarly, P[ x,  | M ] and P[ x,  ] are the same In the LEARNING problem we always write P[ x |  ] to emphasize that we are seeking the  that maximizes P[ x |  ]

More Related