Hidden Markov Modelle
Sven Schuierer
Universität Freiburg
Übersicht
• Markov-Ketten
• CpG-Inseln
• Hidden Markov Modelle– Zustandsfolge dekodieren
• A posteriori Dekodieren
• Parameter schätzen
Übersicht - 2
• Profil HMMs– Zustandsfolge dekodieren– A posteriori Dekodieren– Parameter schätzen
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
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)
Markov-Ketten
• X = (x0,…, xL) ist ein Zufallsvektor
• jede Zufallsvariable xi ist nur abhängig von ihrem Vorgänger xi-1 ist, d.h:
iss ,...,1
ii ssiiii asxsxP ,11 1)(
),...,( 1111 iiii sxsxsxP
Markov-Ketten
• Zwei neue Zustände: – s (Startzustand) – e (Endzustand),
• Wahrscheinlichkeitmenge:
A0 - {P(x0 = s) = a0,s} :
Wahrscheinlichkeit für s Q Startzustand
Beispiel: Markov Ketten
A T
C G
s e
Beispiel: Markov Ketten
A T
C G
s e
Beispiel: Markov Ketten
P(ATGA$)=as,A·aA,T ·aT,G·aG,A ·aA,e
A T
C G
s e
Beispiel: Markov Ketten
P(X)=p(x1)·Π axi,xi+1
,
mit P(x0=s)=a0,s:
P(X)=Π axi,xi+1
A T
C G
s e
L-1
i=1
L-1
i=0
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)
Gegeben
X = (x0,…, xL) Σ* wobei Σ={A,C,G,T}
Frage
Ist X eine CpG-Insel?
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 X CpG-Insel (Log-likelihood test)
Erkennung von CpG-Inseln
P(X|CpG-Insel)P(X|nicht CpG-Insel)_______________ L
i=1Σ ____axi-1,xi,
axi-1,xi,
+
-
Gegeben
X = (x0,…, xL) Σ* wobei Σ={A,C,G,T}
Problem
Finden von CpG-Inseln in langen DNA-Sequenzen
Zwei Ansätze
Lokalisieren von CpG-Inseln
l
L
Ansatz 1
Sliding window
...ACGATACGATAAGTACGATGACCGT...
– Ähnlich zur Erkennung von CpG Inseln– Problem:
• Laufzeit• Grösse der Insel l nicht bekannt
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“
Lokalisieren von CpG-Inseln
Hidden Markov Modell (HMM)
DefinitionEin Hidden Markov Model (HMM) ist Tupel M = (Σ,Q,Θ) mit:
• Σ - Alphabet• Q - Endliche Zustandsmenge • Zustand q gibt Symbol a aus Σ mit Wahrscheinlichkeit eq(a) aus.
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)
HMM für CpG Inseln in DNA sequenz
A+
A- C- G- T-
C+ G+ T+
HMM Übergangswahrscheinlichkeiten
πi/πi+1 A+ C+ G+ T+ A- C- G- T-
A+
C+
G+
T+
0.180p 0.274p 0.426p 0.120p
0.171p 0.368p 0.274p 0.188p
0.161p 0.339p 0.375p 0.125p
0.079p 0.355p 0.384p 0.182p
A-
C-
G-
T-
0.300q 0.205q 0.285q 0.210q
0.322q 0.298q 0.078q 0.302q
0.248q 0.246q 0.298q 0.208q
0.177q 0.239q 0.292q 0.292q4
1
4
1
4
1
4
14
1
4
1
4
1
4
14
1
4
1
4
1
4
14
1
4
1
4
1
4
1
qqqq
qqqq
qqqq
qqqq
p:P(bleibt in CpG-Insel)
q:P(b. in nicht CpG-Insel)
4
1
4
1
4
1
4
14
1
4
1
4
1
4
14
1
4
1
4
1
4
14
1
4
1
4
1
4
1
pppp
pppp
pppp
pppp
Hidden Markov Model (HMM)Weg Π = (π1,…,πL) ist eine Abfolge von Zuständen πi im Modell M
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)
i=1
L
Hidden Markov Model (HMM)
Beobachtbar
π1 π2 π3πn…. Versteckt
x1 x2 x3 xn
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
Π
Viterbi AlgorithmusDefinitionΠk - beliebiger Weg, der im Zust. k endet
vk(i) - Wahrsch. des für die Erzeugung von (x1, x2,…, ,xi) wahrscheinlichsten Weges Πk, der im Zust. k endet:
vk(i) = max P(x1,…,xi,Π)
Wahrschl. Weg Π* der X erzeugt
P(X,Π*) = max {vk(L)· ak,end }
{Π | Π i=k}
k Q
Viterbi AlgorithmusVerfahren (Dynamische Programmierung)Initialisieren:• vbegin(0) = 1• vk(0) = 0 für alle k ≠ 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}
kQ
kQ
Viterbi Algorithmus
Q1 Qi Qi+1
v1(i)
vl(i+1)
v2(i)
v|Q|(i)
a1l
a2l
a|Q|l
el (xi+1)
xi+1
Viterbi Algorithmus
Bei Berechnung der vl(i+1) Zeiger auf das (ein) maximale(s) vk(i) Element speichern.
Weg Π*: Folge Rückwärtszeigern
Viterbi Algorithmus Wertberechnung mit Logarithmus
• Produkte von Wahrscheinlichkeiten können sehr klein werden– Möglichkeit von Underflows Logarithmische Wertberechnung
Viterbi AlgorithmusVerfahren (Wertberechnung mit Logarithmus)
Initialisierung:
• vbegin(0) = 0
• vk(0) = - ∞ für alle k ≠ 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)}
kQ
kQ
Viterbi Algorithmus
Komplexität
vl(i+1) = el(xi+1) · max {vk(i)· akl}
Zeit O(L · |Q|2)– Berechnung eines vk(i) mit O(|Q|) Operationen
– vk(i) berechnen k Q, i ≤ L L · |Q| Einträge
Speicherbedarf O(L · |Q|)– vk(i) speichern kQ, i ≤ L
IQ
AlignmentprofilDefinition
Profil P der 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
i=l
L
L
i=1
____ei(xi)p(xi)
Profile AlignmentHMM
• “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
Profil Alignment
M1 M2 M3
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)
Erweitern der LückeKreieren einer Lücke
Profil Alignment
M1 M2 M3
l0 l1 l2 l3
Profil Alignment Löschzustände• D1,...,DL Löschzustände
• Können kein Zeichen ausgeben: „silent“• Miteinander sequentiell verbunden• Auch mit „Match States“ verbunden
Lokales Alignment
M1 M2 M3
l0 l1 l2 l3
D1 D2 D3
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
Lokales Alignment
• Füge neue Zustände ein
Ende
M1 M2 M3
l0 l1 l2 l3
D1 D2 D3
Start
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
– Di können kein Symbol ausgeben
Z
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)
vj-1(i-1) + log(aMj-1,Mj)M
vj-1(i-1) + log(aDj-1,Mj)D
I
I
vj (i-1) + log(aMj,Ij)
M
vj (i-1) + log(aDj,Ij)
D
_____eMj(xi)
p(xi)
_____eIj(xi)
p(xi)
M
I
Viterbi - Berechnung
Verfahren
vj(i) = max vj(i-1) + log(aIj-1,Dj)
Score(X|Π*) = max vL(m) + log(aIL,end)
vj-1(i) + log(aMj-1,Dj)M
ID
vj-1(i) + log(aDj-1,Dj)D
vL(m) + log(aDL,end)`D
vL(m) + log(aML,end)`M
I
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
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)|Θ)n
i=1
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
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)|Θ)n
i=1
Parameter schätzen HMMsZustandsreihenfolge bekanntVerfahren
• Akl : # Übergänge von Zustand k zu l
• Ek(b) : # Ausgaben von b in Zustand k
Man erhält
akl = ,ek(b) =
Maximum likelihood Schätzer
Σ Akqq Q
_____Akl
q QΣ Ek(σ)_____Ekl
Parameter schätzen HMMsZustandsreihenfolge bekannt• Um WSK = 0 zu vermeiden, addiere zu Akl ,
Ek(b) einen Laplace-Korrektor rkl bzw rk(b) (z.B. 1)
Parameter schätzen HMMsZustandsreihenfolge unbekannt• Wenn Zustandsreihenfolge unbekannt, ist dasProblem, die optimalen Parameter für Θ zufinden 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)
Parameter schätzen HMMs Baum-Welch Training
Verfahren
Initialisierung:• Wähle beliebige Werte für Θ und
Erwartungswert:• Wahrscheinlichkeit für Übergang von k nach
l:
P(πi=k, πi+1=l | X,Θ) = . fk(i)·akl·el(xi+1)·bl(i+1)
P(X)___________________
Parameter schätzen HMMs Baum-Welch Training
Erwartungswerte für Übergänge:
1. Akl = Σ · Σ fk(i) ·akl ·el(xi+1) ·bl(i+1)
2. Ek(b) = Σ · Σ fk(i) ·bk(i)
3. Maximieren:• Berechne Θ wie zuvor aus Akl und Ek(b) und
ermittle Score(X(i),…, X(n)| Θ)
• Wiederhole 2. und 3. bis Score(X(i),…, X(n)| Θ) sich um weniger als ε verbessert
i=1 P(X(j))_____1 L(j)
P(X(j))_____1
j=1
n (j) (j) (j)
(j) (j)
(j){i|xi=b}j=1
n
Parameter Schätzen bei Profil HMMsGegebenMultiples Alignment von X(1),...,X(n)
GesuchtWerte akl, ek(b) für Profil HMMIdeeWie zuvor Werte aus Akl, Ek(b) berechnen
akl= , ek(b) =
Wie zuvor um bei kleinen TrainingsmengenWSKs von 0 zu vermeiden Laplace-Korrektor hinzuaddieren
____Akl
ΣAkqqQ
____Ek(b)Σpk(σ)
σΣ
Profil HMM
Beispiel
ATTAAAAGTTCAGTTACATCTCGCGCCACACCTATC
A TTA AAAGTT CA GTTA CATCTCG C GCCA CACCT ATC1 2 34 5 67
D1 D2 D3 D4 D5 D7D6
l1 l2 l3 l4 l5 l7l6
A:B:G:T:
A:B:G:T:
A:B:G:T:
A:B:G:T:
A:B:G:T:
A:B:G:T:
A:B:G:T:
Start End
Multiple Alignments
GegebenMenge S von Zeichenketten X(1),...,X(n)
GesuchtMultiples AlignmentIdeeBenutze Profil HMMZwei Fälle
– Profil HMM bekannt– Profil HMM unbekannt
Multiple Alignments
Profil HMM bekannt
• Aligniere jede Zeichenkette X(i) einzeln mit Profil HMM
• Bilde aus erhaltenen Alignments Multiples Alignment
Multiple Alignments
Profil HMM unbekannt• Wähle Länge L des Profil HMM• Initialisiere Ausgabegangs und Übergangs
Wahrscheinlichkeiten• Trainiere das HMM mit Baum-Welch
Algorithmus auf gesamter Trainingsmenge• Man erhält nun ein Multiple Alignment aus
dem entstandenen Profil HMM, wenn man wie auf Folie zuvor verfährt
A posteriori DekodierenproblemGegebenHMM M = (Σ,Q,Θ), Zeichenkette X
Gesucht• Wahrscheinlichkeit des Zustandes k als i-ter
Zustand, für alle i {1,…,L}, für alle k QP(πi=k | X)
• Zur Berechnung werden zwei weitere Wahrscheinlichkeiten genutzt:– Forward- und – Backward-Wahrscheinlichkeiten
A posteriori DekodierungForward Algorithmus
Definition
• fk (i) = P(x1 ,... , xi , πi = k), Wahrscheinlichkeit, dass man die Zeichenkette X=(x1,...,xi) ausgibt und im
Zustand πi = k endet
A posteriori DekodierungForward algorithmVerfahrenInitialisierung:
• fbegin(0) = 1• fk(0) = 0 für alle k ≠ beginFür alle i {1,...,L-1} Für alle l Q
fl (i+1) = el(xi+1) · Σ fk(i) · akl
Damit ist:• P (X) = Σ fk(L) · ak,end
k Q
k Q
A posteriori DekodierungBackward algorithm
Definition
• bk(i) = P(xi+1,... ,xL | πi = k), Wahrscheinlichkeit, dass man in Schritt i in Zustand πi = k ist und dann die
Zeichenkette X=(xi+1,...,xL) ausgibt
A posteriori DekodierungBackward algorithm
Verfahren
Initialisierung:
• bk(L) = ak,end k Q
Für alle i {1,...,L-1} für alle k Q
bk(i) = Σ akl · el(xi+1) · bl(i+1)
Damit ist:
• P(X) = Σ abegin,l · el(x1) · bl(1)
l Q
l Q
A posteriori Dekodierung
Komplexität
• Berechnung der fk (i) ´ s und bk(i) ´ s in
• O(L · |Q|2) Zeit
• Platz O(L · |Q|)
A posteriori DekodierungBackward algorithm
Logarithmische Wertberechnung• Wegen Summe nicht trivial
→ verwende Exponentialfunktion
– fbegin(0) = 1
– fk(0) = -∞ für alle k ≠ begin
– fl (i+1) = log[el(xi+1)] + log[Σ akl · exp(fk(i))]
– P(X) = log Σ ak,end · exp( fk(L))
kQ
k Q
A posteriori Dekodierung
• Mit fk (i) und bk(i) kann man die Werte
von P(πi=k|X) berechnen
• Abhängigheit nur vom letzten Zustand:
• P(X,πi=k)
=P(x1,…, xi,πi=k)·P(xi+1,…, xL|x1,…, xi,πi=k)
=P(x1,…, xi,πi=k)·P(xi+1,…, xL|πi=k)
=fk(i)·bk(i)
A posteriori Dekodierung
• Mit der Definition der bedingten Wahrscheinlichkeit:
• P(πi=k|X) = =
• P(X) kann mit der forward oder backward Wahrscheinlichkeit berechnet werden, z.B.
• P(X) = Σ abegin,l· el(x1)·bl(1)
P(X,πi=k)P(X)
________ _______fk(i)·bk(i),
P(X),
lQ
A posteriori Dekodierung Anwendung• Nützlich, wenn verschiedene Wege fast die
gleichen Wahrscheinlichkeiten haben
Definition
• Πi**=arg max {P(Πi=k|X)}
• Aber möglicherweise sind einige Übergänge nicht gültig
k
A posteriori DekodierungAnwendung• Andere Eigenschaften: definiere
Funktion g(k) und betrachte:• G(i|X) = Σ {P(Π i=k|X) · g(k)}• Wenn g(k)=1 für Teilmenge S und g(k)=0
in anderen Fall, G(i|X) ist die Wahrscheinlichkeit dass x aus Zustand S erzeugt wird → A posteriori Wahrscheinlichkeit für CpG
Insel
k
Profil HMMS Forward Algorithmus Verfahren
fbegin(0) = 1
fj (i) = eMj(xi)·[ fj-1(i-1) · aMj-1,Mj
+
fj-1(i-1) · aIj-1,Mj +
fj-1(i-1) · aDj-1,Mj]
fj (i) = eIj(xi)·[ fj (i-1) · aMj-1,Ij
+
fj (i-1) · aIj-1, Ij +
fj (i-1) · aDj-1, Ij]
fj (i) = fj-1(i) · aMj-1, Dj + fj-1(i) · aIj-1, Dj
+ fj-1(i) · aDj-1, Dj
M M
I
D
M
I
D
M I D
I
D
Profil HMMs A posteriori DekodierungFinden echt Werte für Übergangwahrsch. Und
Ausgangwahrsch. Für einen Profil HMM
Gegeben X = (x1,...,xm), Zj {Mj, Ij , Dj}Forward Algorithmus
• fj (i) = P(x1,... ,xi, πi = Zj), WSK, dass man die Zeichen (x1,...,xi) ausgibt und im Zustand Zj endet
Backward Algorithmus
• bj (i) = P(xi+1,... ,xm | Zj), WSK, dass man die Zeichen (xi+1,...,xm) ausgeben wird, wenn man in Zustand Zj ist
Z
Z
Backward Algorithmus
Verfahren
bL(m) = aML,end
bj (i) = bj+1(i+1) · aZj,Mj+1 · eMj+1
(xi+1) +
bj (i+1) · aZj,Ij · eIj
(xi+1) +
bj+1(i) · aZj,Dj+1
für alle Z {M, I, D} Q
Z
Z
I
D
M
Baum-Welch
Schätzen der WSK nach Baum-Welch
• Erwartete Ausgabe WSK für Z {M, I}
EZk(a) = Σ fk (i) bk (i)
• Erwartete Übergangs WSK für Z {M, I, D}
AZk,Mk+1 = Σ fk (i) aZk,Mk+1
eMk+1(xi+1) bk+1(i+1)
AZk,Ik = Σ fk (i) aZk,Ik
eIk(xi+1) bk(i+1)
AZk,Dk+1 = Σ fk (i) aZk,Dk+1
bk+1(i)
____1p(X) i|xi=a
____1p(X)
____1p(X)
____1p(X)
i
i
i
Z Z
Z
Z
Z
M
I
D
Casino Problem
Beispiel
1: 1/10
2: 1/10
3: 1/10
4: 1/10
5: 1/10
6: 1/2
0.05
0.1
0.95 0.9
1: 1/6
2: 1/6
3: 1/6
4: 1/6
5: 1/6
6: 1/6
Fair Loaded
aFF = 0.95aFL = 0.05aLF = 0.1aLL = 0.9
1≤i≤6 eF(i) = 1/6
1≤i≤5 eL(i) = 0.1 eL(6) = 0.5
Viterbi Algorithmus Diagramm
Top Related