„„Fortgeschrittene Fortgeschrittene algorithmische Bioinformatik“algorithmische Bioinformatik“

Thema: Profile HMMsThema: Profile HMMs

ein Vortrag von Gunar Maiwaldein Vortrag von Gunar Maiwald

GliederungGliederung

Einführung Biologischer Hintergrund Multiples Sequenzalignment

Profile HMMs in der Theorie Grundidee Parameterabschätzung Suche mit Profile HMMs

Profile HMMs in der Praxis PFAM

Fazit

Biologischer HintergrundBiologischer Hintergrund

verschiedene Organismen haben Proteine mit ähnlichen Funktionen

Funktionen in Merkmalsregionen (Domänen) konserviert

Domänen: Alpha-Helices und Beta-Faltblätter

50 Aminosäuren und mehr

in Domänen-Familien gruppierbar

Frage: Gegeben ein Protein - lassen sich Domänen entdecken (... und damit Funktionen ableiten) ?

Biologischer HintergrundBiologischer Hintergrund

verschiedene Organismen haben Proteine mit ähnlichen Funktionen

Proteine mit ähnlichen Funktionen in Protein-Familien gruppierbar

Frage: Gegeben ein Protein - kann man es einer Protein-Familie zuordnen ?

Profile HMMs - BegriffsklärungProfile HMMs - Begriffsklärung

Definition: probabilistisches Modell zur Charakterisierung von Merkmalsregionen

Bestandteile: 3 verschiedene Zustände (M, I, D) an jeder Position

sowie Start- und Endzustand Emissionswahrscheinlichkeiten Übergangswahrscheinlichkeiten

Voraussetzung: korrektes MSA dient der Generierung des Profile HMMs

Multiples SequenzalignmentMultiples Sequenzalignment

MSA ist (optimale) Anordnung mehrerer (Protein)sequenzen

MSA zeigt Merkmalsregion(en) einer Domäne

Merkmalsregionen sind stark konserviert mit wenigen Gaps (Helices, Faltblätter)

dazwischen Regionen mit vielen Gaps (Loops)

111111111111111111HBA_HUMAN ---------------VLSPADKTNVKAAWGKVG--HBB_HUMAN --------------VHLTPEEKSAVTALWGKV---GLB5_PETMA -----PIVDTGSV-APLSAAEKTKIRSAWAPVY--MYG_PHYCA ---MACRCEPHALUSVLSEGEWQLVLHVWAKVE--GLB1_GLYDI ---------BLDWRMGLSAAQRQVIAATWKDIAGAGLB3_CHITP THUMMIPIGERMIDGELSADQISTVQASFDKVK--LGB2_LUPLU ---------LUPINGALTESQAALVKSSWEEFN-- *: : . : .:

222222222222222233333333333HBA_HUMAN AHAGEYGAEALERMFLSFPTTKTYFPHF-DLSH-- HBB_HUMAN -NVDEVGGEALGRLLVVYPWTQRFFESFGDLSTPD GLB5_PETMA STYETSGVDILVKFFTSTPAAQEFFPKFKGLTTAD MYG_PHYCA ADVAGHGQDILIRLFKSHPETLEKFDRFKHLKTEAGLB1_GLYDI DNGAGVGKDCLIKFLSAHPQMAAVFGFSGASDP-- GLB3_CHITP ----GDPVGILYAVFKADPSIMAKFTQFAGKDLES LGB2_LUPLU ANIPKHTHRFFILVLEIAPAAKDLFSFLKGTSEVP . : . * * . .

Profile HMMs – GrundideeProfile HMMs – Grundidee

Ideal: MSA ohne Gaps Sequenz x „matcht“ mit Konsensussequenz an

allen Positionen x1...xn

Länge des HMM = Länge d. Konsensussequenz

Emissonswahrscheinlichkeiten für alle AS a: eMj (a)

Transitionswahrscheinlichkeiten: tMjMj+1 = 1

BeginBegin MM11 EndEndMMnnMMjj MMj+1j+1

Profile HMMs – GrundideeProfile HMMs – Grundidee

Insert: in Sequenz x wird Buchstabe xj an Position j eingefügt

Emissonswahrscheinlichkeiten für alle AS a: eIj(a)

Transitionswahrscheinlichkeiten: tMjIj tIjIj

tIjMj+1

BeginBegin EndEndMMjj MMj+1j+1

IIjj

Profile HMMs – GrundideeProfile HMMs – Grundidee

Delete: Teile aus Sequenz x werden gelöscht Realisierung durch Sprung von Mj nach Mj+k

Problem: zu viele Übergänge nötig

BeginBegin EndEndMMjj

Profile HMMs – GrundideeProfile HMMs – Grundidee

„silent states“: Zustände ohne Emission ermöglichen Gaps variabler Länge in Sequenz x Transitionswahrscheinlichkeiten: tMj-1Dj

tDjDj+1 tDjMj+1

BeginBegin EndEndMMjj MMj+1j+1

DDjj

Profile HMMs – GrundideeProfile HMMs – Grundidee

MM00

BeginBeginMML+1L+1

EndEndMMjj

IIjj

DDjj

Insgesamt:

Profile HMMs – GenerierungProfile HMMs – Generierung

bat V E V - - L

rat - E - E V L

cat L E V – E -

gnat L - - L E L

goat - E V - - L

. 1 2 . . 3

Profile HMMs – GenerierungProfile HMMs – Generierungmatch emissions:

0 1 2 3

E 4 - -

L - - 4

V - 3 -

... - - -

insert emissions:

0 1 2 3

E - - 3 -

L 2 - 1 -

V 1 - 1 -

... - - -

state transitons:

0 1 2 3

M - M 2 3 2 4

M – D - 1 -

M – I 3 - 1 -

0 1 2 3

I – M 2 - 2 -

I – D 1 - 1

I – I - - 2 -

0 1 2 3

D– M - - 1

D– D 1 -

D – I - 2 -

Profile HMMs – GenerierungProfile HMMs – Generierungmatch probabilities: eMj

(xj)

0 1 2 3

E 1 - -

L - - 1

V - 1 -

... - - -

insert probabilities: eIj (xj)

0 1 2 3

E - - 0.6 -

L 0.67 - 0.2 -

V 0.33 - 0.2 -

... - - -

state transitons: tZj-1Zj

0 1 2 3

M - M 0.4 0.75 0.67 1

M – D - 0.25 -

M – I 0.6 - 0.33 -

0 1 2 3

I – M 0.67 - 0.4 -

I – D 0.33 - 0.2

I – I - - 0.4 -

0 1 2 3

D– M - - 1

D– D 1 -

D – I - 1 -

Profile HMMs – GenerierungProfile HMMs – Generierung

MM00

BeginBeginMM44

EndEndMM22

II22

DD22DD11 DD33

MM11 MM33

II00L 0.67V 0.33...

E 0.6L 0.2V 0.2...

E 1.0...

V 1.0...

L 1.0...

0.4

0.6

0.67

0.33

ParameterabschätzungParameterabschätzung

Frage: Wie werden Emissions- und Transitions-wahrscheinlichkeiten abgeschätzt ?

einfacher Ansatz: „Maximum Likelihood“:

eMj (a) = Vorkommen der AS a an Position j

Anzahl aller AS b an Position j

Problem: Kommt AS a’ an Position j im MSA nicht vor, so ist eMj

(a‘) = 0

Grund: Trainingsdaten überdecken nicht alle in der Realität existierenden Fälle

ParameterabschätzungParameterabschätzung

Pseudocounts: häufig verwendet Laplace (k=1)

eMj (a) = Vorkommen von AS a an Position j + 1*k

Anzahl aller AS b an Position j + 20*k

kleine Trainingsmenge: Wahrscheinlichkeit nicht gesehener Ereignisse überschätzt

grosse Trainingsmenge: Angleichung an Maximum Likelihood Werte

Problem: grosser Aufwand nötig, um k gut abzuschätzen (50 und mehr Beispiele)

ParameterabschätzungParameterabschätzung

Mixmodell: Berechnung der Pseudocounts durch Einbeziehen einer Substitutionsmatrix

Umrechnung von Matrixeintrag s(b,a) nach P(a|b) positionsspez.Pseudocount: αja = beMj

(b)*P(a|b)

eMj (a) = Vorkommen von AS a an Position j + αja

Anzahl aller AS b an Position j + αjb

Problem: heuristisches Modell ohne statistisch fundierte Erklärung der Herangehensweise

Suche mit Profile HMMsSuche mit Profile HMMs

Suche: Hauptanwendung von Profile HMMs 1. Entdecken von Merkmalsregionen in Proteinen 2. Zuordnung von Proteinen zu Familien zwei unterschiedliche Algorithmen:

Viterbi-Algorithmus Forward-Algorithmus

dynamische Programmierung

Suche mit Profile HMMsSuche mit Profile HMMs Gegeben: Sequenz x1...xn, Profile-HMM λ Frage: Wie wahrscheinlich ist es, dass x1...xn

durch λ modelliert wird ? Brute-Force:

Durchlaufe alle potentiellen Pfade π1 ... πm für x1...xn und berechne die Wahrscheinlichkeiten p1 ... pm

Summiere alle Wahrscheinlichkeiten auf Wenn Schwellwert überschritten, dann Treffer

Problem: # potientieller Pfade: m >> 3n

# Rechenschritte pro Pfad: 2n

Suche mit Profile HMMsSuche mit Profile HMMs

Viterbi: ermittelt die wahrscheinlichste Abfolge π* von

versteckten Zuständen gegeben eine Beobachtungsfolge x und ein HMM λ

Beobachtungsfolge ist die Sequenz des Proteins versteckte Zustände sind Mj, Ij und Dj

Falls P( x, π* | λ ) einen Schwellwert übersteigt, gehört x der durch λ beschriebenen Familie an

Hier: Variante des Viterbi-Algorithmus speziell für Profile-HMMs

Suche mit Profile HMMsSuche mit Profile HMMs

dynamische Programmierung: Sei M0 = Anfangszustand

mit einem „Viterbi-Score“ V0M(0) = 0

Sei ML+1 = Endzustand

mit einem „Viterbi-Score“ VL+1M(n),für einen optimalen

Pfad von Zuständen z0,...,zL+1 mit der Ausgabe x0,...,xn

Suche mit Profile HMMsSuche mit Profile HMMs

dynamische Programmierung:

Sei z0,...,zj-1 eine „optimale“ Zustandsfolge für die Ausgabe x1...xi-1

VjM(i) ist der „Viterbi-Score“ für die Zustandsfolge

z0...zj-1,Mj mit der Ausgabe x1...xi-1,xi

VjI(i) ist „Viterbi-Score“ für z0...zj-1,Ij und x1...xi-1,xi

VjD(i-1) ist „Viterbi-Score“ für z0...zj-1,Dj undx1...xi-1

VjM(i ) = log

eMj (xi)

qxi+ max

Vj-1M(i-1) + log tMj-1Mj

Vj-1I(i-1) + log tIj-1Mj

Vj-1D(i-1) + log tDj-1Mj

VjI(i ) = log

eIj (xi)

qxi+ max

VjM(i-1) + log tMjIj

VjI(i-1) + log tIjIj

VjD(i-1) + log tDjIj

VjD(i) = max

Vj-1D(i) + log tMj-1Dj

Vj-1D(i) + log tIj-1Dj

Vj-1D(i) + log tDj-1Dj

Suche mit Profile HMMsSuche mit Profile HMMs

Laufzeit: # möglicher „Viterbi-Scores“: 3i*j # Rechenschritte pro „Viterbi-Score“: 4

Platzkapazität: Backtracking erfordert die Speicherung aller Viterbi-

Scores

Suche mit Profile HMMsSuche mit Profile HMMs

Forward: ermittelt für jeden Buchstaben xj aus Sequenz x den

wahrscheinlichsten Zustand

Zustandsfolge = Aneinanderreihung der wahrscheinlichsten Zustände und eventueller Zwischenzustände

Viterbi:wahrscheinlichste Abfolge von Zuständen

Forward: Abfolge wahrscheinlichster Zustände

PFAMPFAM

DB mit Vielzahl an MSAs und Profile HMMs

analysiert Proteine ermöglicht Domänen-Organisation von

Proteinen zu betrachten 75% alle Proteine mit mind.1 Match in

PFAM

FAZITFAZIT

Profile HMMs aus MSA erzeugbar Wahrscheinlichkeiten für Emission und

Transition werden abgeschätzt Suche findet Proteindomänen und -familien

Viterbi- und Forward-Algorithmus mit dynamischer Programmierung

Realisierung in PFAM