Gensuche mit Hidden Markov Modellen Zentrum für Bioinformatik der Universität des Saarlandes

1. Gensuche mit Hidden Markov Modellen Zentrum für Bioinformatik der Universität des Saarlandes WS 2002/2003

2. Startkodon Stopkodon ... CAT ATG TTT CCA AGT ACA TGG TAT GTA TAA GGG CAT ... Kodierender Bereich Worum geht es? • In der enormen Datenmenge eines Genoms sollen die kodierenden Regionen bestimmt werden

3. Suche nach Genen • Aufgabenstellung: gegeben eine DNA-Sequenz, klassifiziere jede einzelne Base als Teil eines • Exons (kodierender Bereich) • Introns (nichtkodierender Bereich innerhalb eines Gens) • Zwischengenetischen Bereichs (nichtkodierender Bereich zwischen zwei Genen) Exon 1 Intron 1 Exon 2 Exon 1 Intron 1 Gen 1 Zwischengen. Bereich Gen 2

4. Wie stellen wir das an? • Idee: Suche nach Startkodon/Stopkodon – Paaren. Alles dazwischen ist kodierend. Nachteil: Funktioniert nicht! u.A. werden überlappende Gene und Introns nicht berücksichtigt • Statt dessen: verwende statistische Informationen um Teilsequenzen zu klassifizieren • Analogie: Automatische Erkennung der Sprache eines Textes. In einem typischen deutschen Text macht der Buchstabe ‚e‘ ca. 16,55% aller Buchstaben aus, in einem schwedischen nur ca. 9.77%.  zähle die e‘s im Text, um zu berechnen mit welcher Wahrscheinlichkeit es sich um einen deutschen Text handelt

5. Kurzwiederholung: Eine Markovketteist ein stochastischer Prozess, der nacheinander eine Reihe von Zuständen mit einer gewissen Wahrscheinlichkeit durchläuft. Dabei hängt die Wahrscheinlichkeit für den jeweils nächsten Zustand nur vom aktuellen Zustand ab: P(ti+1|ti, ti-1,...,tj) = P(ti+1|ti) Hidden Markov Modelle In ähnlichen Fragestellungen (z.B. in der Spracherkennung) haben sich Hidden Markov Modelle als sinnvoll erwiesen

6. Hidden Markov Modelle • Ein Hidden Markov Modell besteht aus einer Markovkette, bei der jedoch einige Zustände versteckt sind, d.h. wir können nicht genau angeben, in welchem Zustand sich das System befindet. • In diesem Fall können wir nur von den Effekten die wir beobachten auf die Wahrscheinlichkeit jedes Zustands zurückschließen • Ein Hidden Markov Modell lässt sich als Graph darstellen, dessen Knoten die Zustände und dessen Kanten die Übergänge darstellen. Die Kanten sind mit den Übergangswahrscheinlichkeiten gewichtet.

7. VEIL – The Viterbi Exon-Intron Locator • Ein einfaches Beispiel für ein solches Hidden Markov Modell zur Gensuche ist VEIL • VEIL wurde vorgestellt in: Henderson, Salzberg and Fasman: „Finding Genes in DNA with a Hidden Markov Model“, J. Comp. Biol. (1996) • Der Aufbau von VEIL besteht aus 3 Schritten: • Definition des Modells • Training des Modells mit dem EM-Algorithmus • Klassifizieren von Sequenzen mittels Viterbi-Algorithmus

8. Die Modelldefinition von VEIL • VEIL ist ein modular aufgebautes Hidden-Markov-Modell • Es besteht aus einzelnen Komponenten (selber wieder vollständige HMM‘s), die zu einem Gesamtmodell verdrahtet werden • Jedes Modul repräsentiert eine bestimmte Klassifikation der DNA

9. Downstream Upstream Stopkodon Startkodon Exon Intron 5´Splice Site 3´Splice Site Das Gesamtmodell von VEIL

10. 16 Backedges a a a Anpassung des Reading Frames Anpassung des Reading Frames c c c g g g t t t 3´Splice Site Eine beliebige Folge von Codons 5´Splice Site Das Intron-Modul von VEIL

11. a a c c g g t t a a g g Das Exon-Modul von VEIL Startkodon 16 Backedges a c g t 5´Splice Site 3´Splice Site Downstream

12. a a c c g g t t a a g g Beispiel 1 – ein 1-Exon Gen Startkodon 16 Backedges a c g t 5´Splice Site 3´Splice Site Downstream (Start)

13. a a c c g g t t a a g g Beispiel 1 – ein 1-Exon Gen Startkodon 16 Backedges a c g t 5´Splice Site 3´Splice Site Downstream (Start) A

14. a a c c g g t t a a g g Beispiel 1 – ein 1-Exon Gen Startkodon 16 Backedges a c g t 5´Splice Site 3´Splice Site Downstream (Start) ACC

15. a a c c g g t t a a g g Beispiel 1 – ein 1-Exon Gen Startkodon 16 Backedges a c g t 5´Splice Site 3´Splice Site Downstream (Start) ACC T

16. a a c c g g t t a a g g Beispiel 1 – ein 1-Exon Gen Startkodon 16 Backedges a c g t 5´Splice Site 3´Splice Site Downstream (Start) ACC TAA (Downstream) Stopkodon

17. a a c c g g t t a a g g Beispiel 2 – ein initiales Exon Startkodon 16 Backedges a c g t 5´Splice Site 3´Splice Site Downstream (Start) CTT

18. a a c c g g t t a a g g Beispiel 2 – ein initiales Exon Startkodon 16 Backedges a c g t 5´Splice Site 3´Splice Site Downstream (Start) CTT ACC

19. a a c c g g t t a a g g Beispiel 2 – ein initiales Exon Startkodon 16 Backedges a c g t 5´Splice Site 3´Splice Site Downstream (Start) CTT ACC AT Beliebig, zur Anpassung des Reading Frames

20. a a c c g g t t a a g g Beispiel 2 – ein initiales Exon Startkodon 16 Backedges a c g t 5´Splice Site 3´Splice Site Downstream (Start) CTT ACC AT (Übergang zum 1. Intron) Beliebig, zur Anpassung des Reading Frames

21. Das Training des Modells • Die Topologie des HMM ist nun bekannt. Es fehlen „nur“ noch die Übergangswahrscheinlichkeiten Pij • Für einige Pij haben wir Schätzungen, z.B. aus biologischen Modellen. Diese verwenden wir dann als Ausgangsbelegung für die jeweiligen Kanten. • Die restlichen Kanten belegen wir mit zufälligen Startwerten zwischen 0 und 1. • Anschließend werden mit diesem initialen Modell bekannte Trainingssequenzen untersucht, und die Kantengewichte rekursiv angepasst, um die Vorhersagequalität zu optimieren

22. Der EM-Algorithmus • Klassifiziere eine bekannte Sequenz s1 vom Typ M mit dem HMM. Dies liefert P(s1|M). Wiederhole dies für alle Trainingssequenzen si • Passe die Pij iterativ an, so dass P(S|M) := Pi P(si|M) maximal wird. Dazu werden die Pij in jedem Iterationsschritt so verändert, dass Pnew(si|M)¸ Pold(si|M) i Laufzeit pro Iteration: O(ne), wobei n:= Gesamtlänge aller Trainingssequenzen und e := Anzahl Kanten im HMM

23. Problemstellung: Für eine Eingabesequenz Si der Länge N, bestimme den Pfad qi durch die Zustände der HMM , der Si am wahrscheinlichsten generiert hat, d.h. den Pfad, der P(q1,..,qN|S1,...SN,) maximiert Anwenden des Modells • Da nun das Modell vollständig ist, können wir es auf unbekannte Sequenzen anwenden. • Naiver Brute-Force-Ansatz: berechne alle Pfade der Länge N, überprüfe, welche davon Si generieren können, und wähle denjenigen mit maximaler Wahrscheinlichkeit  exponentielle Laufzeit!

24. Paradigma: Speichern und Wiederverwerten bereits berechneter Informationen Der Viterbi-Algorithmus • Um die Laufzeit zu reduzieren, verwendet man einen Algorithmus, der an dynamische Programmierung angelehnt ist • Dazu bauen wir rekursiv eine Datenstruktur (Trellis) auf, die alle Pfade der Länge i enthält. • Der Trellis besteht aus i Schichten, die jeweils alle N Zustände des HMM enthalten. • Gibt es im HMM eine Kante von i nach j, dann gibt es in jeder Schicht S(t) eine Kante von i zum Knoten j in der Schicht S(t+1)

25. di+1(qj) := max1kN[di(qk)·P(Si+1|qj)·P(qj|qk)] Yi+1(qj):= argmax1kN[di(qk)·P(Si+1|qj)·P(qj|qk)] Der Viterbi-Algorithmus • Wir verwenden den Trellis, um die folgenden Informationen zu berechnen: • di+1(qj) ist die WK des wahrscheinlichsten Pfades, der beim (i+1)ten Eingabezeichen in qj endet. Yi+1(qj) ist der Vorgängerknoten von qj auf diesem wahrscheinlichsten Pfad. • Die Formeln erklärt man so: um bei i+1 in qj zu enden, mußte man zuvor einen anderen Knoten qk erreichen (di(qk)), von qk nach qj übergehen (P(qj|qk)) und schließlich in qj Si+1 ausgeben (P(Si+q|qj))

26. di+1(qj) := max1kN[di(qk)·P(Si+1|qj)·P(qj|qk)] Yi+1(qj):= argmax1kN[di(qk)·P(Si+1|qj)·P(qj|qk)] Der Viterbi-Algorithmus • Initialisierung: - d1(qj) = Anfangswahrscheinlichkeit von qj j - Y1(qj) = 0  j • Terminierung: - P* = maxi[dT(qi)] - q*T = argmaxi[dT(qi)] • Backtracking: - q*t = Yt+1(q*t+1), t = T-1,...,1

27. Der Viterbi-Algorithmus • Nachteil: der Trellis benötigt für ein HMM mit e Kanten und eine Eingabesequenz der Länge n O(ne) Speicherplatz • Idee: benutze die Markoveigenschaft, um den Trellis zu verkleinern! • Dies liegt daran, daß die Übergangswahrscheinlichkeiten einer Markovkette unabhängig von der Geschichte des Pfades sind! Ein Knoten, der auf keinem optimalen Pfad von der Schicht i in die Schicht i+1 liegt, kann aus der Schicht i+1 und allen darauffolgenden Schichten gestrichen werden!

28. Jetzt fehlt nur noch... Euch allen Frohe Weihnachten und erholsame Ferien zu wünschen!