Hidden Markov Modelle

1. Hidden Markov Modelle Sven Schuierer Universität Freiburg

2. Übersicht • Markov-Ketten • CpG-Inseln • Hidden Markov Modelle • Zustandsfolge dekodieren • A posteriori Dekodieren • Parameter schätzen

3. Übersicht - 2 • Profil HMMs • Zustandsfolge dekodieren • A posteriori Dekodieren • Parameter schätzen

4. Motivation Markov Ketten Wichtige Probleme in der Biologie: • Finden von Alignments für lange Zeichenketten • Erkennung von Mustern in Zeichenketten Beispiel: Signifikante Stellen im Genom (CpG-Inseln) • Erkennung von CpG-Inseln. • Entscheidung: Zeichenkette X enthält eine CpG-Insel

5. Markov-Ketten Definition: Eine Markov-Kette ist ein Tupel (Σ,Q,A) mit: Σ - Alphabet Q - endliche Zustandmenge, jeder Zustand entspricht einem Zeichen ausΣ A - {as,t | s, t  Q}, mit as,t ≡ P(xi=t | xi-1=s)(Übergangswahrscheinlichkeiten)

6. Markov-Ketten • X = (x0,…, xL)ist ein Zufallsvektor • jede Zufallsvariable xi ist nur abhängig von ihrem Vorgänger xi-1 ist, d.h:

7. Markov-Ketten • Zwei neue Zustände: • s (Startzustand) • e (Endzustand), • Wahrscheinlichkeitmenge: A0 - {P(x0 = s) = a0,s} : Wahrscheinlichkeit für s  Q Startzustand

8. Beispiel: Markov Ketten A T s e C G

9. A T s e C G Beispiel: Markov Ketten

10. A T s e C G Beispiel: Markov Ketten P(ATGA$)=as,A·aA,T·aT,G·aG,A ·aA,e

11. A T s e C G L-1 i=1 L-1 i=0 Beispiel: Markov Ketten P(X)=p(x1)·Πaxi,xi+1, mit P(x0=s)=a0,s: P(X)=Πaxi,xi+1

12. CpG Inseln • CpG ist Nukleotidpaar CG • CpG's sind selten im Genom • CpG-Inseln: Kurze Teilfolgen (ca. 100 Nukleotide), in denen CpG's vermehrt vorkommen • CpG-Inseln in wichtigen Teilen des Genoms (z.B. Promotoren vieler Gene)

13. Erkennung von CpG-Inseln Gegeben X = (x0,…, xL)  Σ* wobei Σ={A,C,G,T} Frage Ist Xeine CpG-Insel?

14. P(X|CpG-Insel) P(X|nicht CpG-Insel) _______________ axi-1,xi, axi-1,xi, + ____ - L Σ i=1 Erkennung von CpG-Inseln Verfahren • Zwei Markov-Ketten: • eine Markov-Kette mit Wahrscheinlichkeiten für CpG-Insel (+), • eine Markov-Kette für nicht-CpG-Insel (-) Score(X) = log = log • Je grösser Score(X)ist, desto wahrscheinlicher ist XCpG-Insel (Log-likelihood test)

15. Lokalisieren von CpG-Inseln Gegeben X = (x0,…, xL)  Σ* wobei Σ={A,C,G,T} Problem Finden von CpG-Inseln in langen DNA-Sequenzen Zwei Ansätze

16. L l Lokalisieren von CpG-Inseln Ansatz 1 Sliding window ...ACGATACGATAAGTACGATGACCGT... • Ähnlich zur Erkennung von CpG Inseln • Problem: • Laufzeit • Grösse der Insel l nicht bekannt

17. Lokalisieren von CpG-Inseln Ansatz 2 Zwei Markov Ketten in einem Modell → Hidden Markov Modell • Wahrscheinlichkeit für Übergang zwischen „CpG-Insel“ und „nicht-CpG-Insel“

18. Hidden Markov Modell (HMM) Definition Ein Hidden Markov Model (HMM) ist Tupel M= (Σ,Q,Θ) mit: • Σ - Alphabet • Q - Endliche Zustandsmenge • Zustand q gibt Symbol a aus Σmit Wahrscheinlichkeit eq(a) aus.

19. Hidden Markov Model (HMM) Definition (II) Θ - Menge der Wahrscheinlichkeiten: • A:Übergangswahrscheinlichkeiten A:{akl | k,l  Q}, akl= P(πi=l |πi-1=k) • E:Ausgabewahrscheinlichkeiten E:{ek(b) | k  Q, b  Σ}, ek (b) = P(xi=b |πi=k)

20. A+ C+ G+ T+ A- C- G- T- HMM für CpG Inseln in DNA sequenz

21. HMM Übergangswahrscheinlichkeiten p: P(bleibt in CpG-Insel) q: P(b. in nicht CpG-Insel)

22. L i=1 Hidden Markov Model (HMM) WegΠ= (π1,…,πL)ist eine Abfolge von Zuständen πi im ModellM Wahrscheinlichkeit für die Erzeugung der Sequenz X aus Weg Π: P(X,Π)= aπ0,π1·Π eπi (xi) ·aπi,πi+1 π0 = begin, πL+1 =end Gesucht: Weg Π,der X erzeugt hat. Π ist nicht direkt aus X ablesbar (hidden path)

23. π1 π2 π3 πn Hidden Markov Model (HMM) x1 x2 x3 xn Beobachtbar …. Versteckt

24. Dekodierungsproblem Gegeben HMM M= (Σ,Q,Θ),Sequenz X  Σ* Gesucht Wahrschl. Pfad Π* = (π1,…,πL), der X erzeugt: Π*= arg max{P(X,Π)}  Lösung für CpG-Insel Problem Π

25. Viterbi Algorithmus Definition Πk- beliebiger Weg, der im Zust. kendet vk(i) - Wahrsch. des für die Erzeugung von (x1, x2,…, ,xi)wahrscheinlichsten Weges Πk, der im Zust. kendet: vk(i) = maxP(x1,…,xi,Π) Wahrschl. Weg Π*der X erzeugt P(X,Π*)= max {vk(L)· ak,end } {Π | Πi=k} k  Q

26. Viterbi Algorithmus Verfahren (Dynamische Programmierung) Initialisieren: • vbegin(0) = 1 • vk(0) = 0 für allek ≠ begin Für jedes i=0,...L-1: vl(i+1) = el (xi+1) · max{vk(i)·akl} Damit: P(X,Π*) = max{vk(L)·ak,end} kQ kQ

27. Viterbi Algorithmus Qi+1 Qi Q1 v1(i) a1l v2(i) a2l xi+1 el (xi+1) vl(i+1) a|Q|l v|Q|(i)

28. Viterbi Algorithmus Bei Berechnung dervl(i+1) Zeiger auf das (ein) maximale(s)vk(i) Element speichern. Weg Π*: Folge Rückwärtszeigern

29. Viterbi AlgorithmusWertberechnung mit Logarithmus • Produkte von Wahrscheinlichkeitenkönnen sehr klein werden • Möglichkeit von Underflows • Logarithmische Wertberechnung

30. Viterbi Algorithmus Verfahren (Wertberechnung mit Logarithmus) Initialisierung: • vbegin(0) = 0 • vk(0) = - ∞für allek ≠ begin Für jedes i=0,...L-1: vl(i+1) = log el(xi+1) +max{vk(i) + log(akl)} Damit ist: Score(X,Π*) = max{vk(L) + log(ak,end)} kQ kQ

31. Viterbi Algorithmus Komplexität vl(i+1) = el(xi+1) · max{vk(i)·akl} Zeit O(L · |Q|2) • Berechnung einesvk(i) mit O(|Q|) Operationen • vk(i) berechnen k  Q, i ≤ L L · |Q|Einträge Speicherbedarf O(L · |Q|) • vk(i) speichern  kQ, i ≤ L IQ

32. ei(xi) ____ p(xi) Alignmentprofil Definition Profil Pder Länge L ist: Menge der Wahrscheinlichkeiten ei(b),dass Zeichen b an der i-ten Position vorkommt, für Alle b Σ, 1 ≤ i ≤ L Wahrscheinlichkeit der Zeichenkette X=(x1,...,xL) geg. Profil P: P(X|P)=Πei(xi) Wert eines lückenlosen Alignments von X an Profil P (log-likelihood): Score(X| P) = Σ log p(b) – Hintergrundhäufigkeit des Zeichen b L i=1 L i=l

33. Profile Alignment HMM • “Match States“ M1,...,ML entsprechen Übereinstimmung von Zeichen der Zeichenkette mit Profilpositionen • Zustände sequentiell verbunden • ej(b) Ausgabe Wahrscheinlickeit von b im Zustand Mj • ai,j+1=1, für 1≤ j ≤ L:Übergangswahrscheinlichk. • Alignment trivial, da es immer genau einen Nachfolgezustand gibt

34. Profil Alignment M1 M2 M3

35. Profile AlignmentEinfügezustände • I1,...,IL Einfügezustände • Für alle b  Σ, eIj(b) = p(b) • Affine Kostenfunktion einer Lücke der Länge h log(aMj,Ij) + log(aIj,Mj+1) + (h-1)·log(aIj,Ij) Kreieren einer Lücke Erweitern der Lücke

36. Profil Alignment l0 l1 l2 l3 M1 M2 M3

37. Profil Alignment Löschzustände • D1,...,DL Löschzustände • Können kein Zeichen ausgeben: „silent“ • Miteinander sequentiell verbunden • Auch mit „Match States“ verbunden

38. Lokales Alignment D2 D3 D1 l0 l1 l2 l3 M1 M2 M3

39. Komplettes Profil HMM • L Ebenen mit je drei Zuständen Mj,Ij, Dj • Endzuständen und Startzuständen • Übergänge von I zu D oder umgekehrt sind unwahrscheinlich

40. D2 D3 D1 l0 l1 l2 l3 M1 M2 M3 Start Ende Lokales Alignment • Füge neue Zustände ein

41. Dekodieren Viterbi Algorithmus • vj(i) – Wert des wahrschl. Pfades der (x1,...,xi) Präfix von aligniert und im Zustand Zj endet (Z  {M,I,D}) • Viterbi Algorithmus funktioniert wie zuvor, neu: • Höchstens drei Vorgängerzustände • Dikönnen kein Symbol ausgeben Z

42. vj-1(i-1) + log(aMj-1,Mj) vj-1(i-1) + log(aDj-1,Mj) M D M D vj (i-1) + log(aMj,Ij) vj (i-1) + log(aDj,Ij) eMj(xi) eIj(xi) _____ _____ p(xi) p(xi) Viterbi - Berechnung Verfahren vbegin(0)=0 vj(i)= log + max vj-1(i-1) + log(aIj-1,Mj) vj(i)= log + max vj(i-1) + log(aIj,Ij) M I I I

43. vj-1(i) + log(aMj-1,Dj) vj-1(i) + log(aDj-1,Dj) M D vL(m) + log(aDL,end)` vL(m) + log(aML,end)` D M Viterbi - Berechnung Verfahren vj(i) = max vj(i-1) + log(aIj-1,Dj) Score(X|Π*) = max vL(m) + log(aIL,end) D I I

44. Viterbi - Berechnung Komplexität • Es werden O(L·|X|)Werte berechnet • Jede Berechnung braucht O(1) Operationen (nur drei Vorgänger) • O(L·|X|) Zeit und O(L·|X|) Platz

45. n i=1 Parameter schätzen für HMMs Gegeben • X(1),...,X(n)  Σ*(Trainings-Sequenzen) Zeichenketten der Längen L(1),...,L(n), die vom gleichen HMM M generiert wurden • Wahrscheinlichkeiten für Zeichenketten schreiben: P(X(1),...,X(n) |Θ) = Π P(X(i)|Θ)

46. Parameter schätzen für HMMs Gesucht • Maximieren des logarithmischen Werts: Θ* = arg max {Score(X(1),...,X(n) |Θ)} Wobei Score(X(1),...,X(n) |Θ) = log P(X(1),...,X(n) |Θ) = Σlog(P(X(i)|Θ)) Θ n i=1

47. n i=1 Parameter schätzen für HMMs Gegeben • X(1),...,X(n)  Σ*(Trainings-Sequenzen) Zeichenketten der Längen L(1),...,L(n), die vom gleichen HMM M generiert wurden • Wahrscheinlichkeiten für Zeichenketten schreiben: P(X(1),...,X(n) |Θ) = Π P(X(i)|Θ)

48. Ekl Akl _____ _____ ΣEk(σ) ΣAkq q Q q Q Parameter schätzen HMMsZustandsreihenfolge bekannt Verfahren • Akl : # Übergänge von Zustand kzu l • Ek(b) : # Ausgaben von b in Zustand k Man erhält akl = ,ek(b) = Maximum likelihoodSchätzer

49. Parameter schätzen HMMsZustandsreihenfolge bekannt • Um WSK = 0 zu vermeiden, addiere zu Akl , Ek(b) einen Laplace-Korrektorrklbzw rk(b) (z.B. 1)

50. Parameter schätzen HMMsZustandsreihenfolge unbekannt • Wenn Zustandsreihenfolge unbekannt, ist das Problem, die optimalen Parameter für Θ zu finden NP-vollständig • Benutze Baum-Welch-Algorithmus (Optimierungsverfahren) zum heuristischen Finden einer Lösung Notation • fk (i), bk(i)sind Forward- bzw Backward-WSKn für die Zeichenketten X(j) (j) (j)