Biologia computazionale A.A. 2010-2011 semestre II U NIVERSITÀ DEGLI STUDI DI MILANO Docente:...

Biologia computazionale

A.A. 2010-2011 semestre II

UNIVERSITÀ DEGLI STUDI DI MILANODocente: Giorgio Valentini

Istruttore: Matteo Re

7 Hidden Markov Models

C.d.l. Biotecnologie Industriali e Ambientali

PROBLEMA: Abbiamo ricevuto un file FASTA contenente le sequenze cromosomiche assemblate prodotte nell’ambito di un progetto di sequenziamento genomico.

•Questi file non contengono altra informazione se non la sequenza nucleotidica•Perchè la sequenza sia utile dobbiamo annotarla identificando gli elementi funzionali in essa contenuti e la loro posizione.

Come possiamo annotare il genoma appena ricevuto?

IDENTIFICAZIONE DI ELEMENTI FUNZIONALI IN UNA SEQUENZA (seq. nucleotidiche o proteiche)Bio

IPOTESI DI LAVORO:•E’ disponibile il genoma annotato di una specie evolutivamente vicina. allineamento al nuovo genoma di tutte le proteine della specie già annotata

•Sequenziamento di tutti i trascritti seguito da allineamento delle loro sequenze al genoma

OSSERVAZIONI:

1)Non è detto che sia disponibile il genoma annotato di una specie evolutivamente vicina

2)Non tutti gli elementi funzionali sono geni (quindi il “trucco” dei trascritti non è sempre applicabile)

IDENTIFICAZIONE DI ELEMENTI FUNZIONALI IN UNA SEQUENZA Bio

I geni sono comunque tra i primi elementi funzionali da cercare in una seq. genomica non annotata

IPOTESI DI LAVORO (II):

•Supponiamo di avere alcune annotazioni per il genoma (ad esempio di sponiamo delle annotazioni del 20% (stima) dei geni dell’organismo) e di voler identificare le posizioni degli altri geni.

POSSIBILE SOLUZIONE:

1)Utilizzo degli elementi noti per costruire modelli

2)Utilizzo dei modelli prodotti per effettuare una scansione della sequenza genomica che permetta di rilevare la

presenza di regioni simili al modello considerato


DIFFICOLTA’ CHE CARATTERIZZANO QUESTO APPROCCIO:

Supponiamo di voler creare dei modelli di elementi funzionali che permettano di predire la presenza (e la struttura) dei geni.

L’obiettivo è quello di creare UN MODELLO che permetta di predire tutti I geni presenti in un genoma …

MA I GENI DI UN ORGANISMO POSSONO AVERE CARATTERISTICHE MOLTO DIVERSE … E’ DIFFICILE PER UN UNICO MODELLO ADATTARSI A TUTTI I GENI!



Prodotto del gene Dimensione gene (kb)

Numero di esoni

Dimensione media esoni (bp)

Dimensione media introni (kb)

tRNA-tyr 0.1 2 50 0.02

Insulin 1.4 3 155 0.48

β-globin 1.6 3 150 0.49

Class I HLA 3.5 8 187 0.26

Serum albumin 18 14 137 1.1

Type VII collagen 31 114 77 0.19

Complement C3 41 29 122 0.9

Phenylalanine-hydroxylase

90 26 96 3.5

Factor VIII 186 26 375 7.1

CTFR (cystic fibrosis)

250 27 227 9.1

Dystrophin 24000 79 180 30.0

DETERMINAZIONE STRUTTURA DI UN GENE Bio

Approccio sperimentale:Estrazione degli RNA e sequenziamento.Confronto con dati pubblici: sequenze di cDNA assemblate e sequenze EST.Limitazioni: Trascritti rari,specificità cellulare e tissutale, espressione condizionale. NON E’ possibile estrarre RNA da tutti i tipi cellulari, a tutti gli stadi di sviluppo e in tutte le possibili condizioni ambientali da un organismo complesso.

Approccio Computazionale:Predizione ab initio basata su caratteristiche della sequenza in esame.

Approccio ibrido: Predizione ab initio integrata con dati ottenuti mediante sequenziamento di RNA parziali e confronto con strutture di geni omologhi noti in altri organismi. (similarità a livello di sequenza proteica o nucleotidica)


Problema pratico:

La complessità dei modelli probabilistici è proporzionale alla complessità strutturale degli elementi funzionali che vogliamo modellare.

Questo implica che un modello adatto a descrivere le caratteristiche di un gene eucariotico ha già un livello di complessità abbastanza elevato.

Per introdurre I modelli basati su catene di Markov e gli Hidden Markov Models (HMM) utilizzeremo esempi più semplici (e quindi associati a elementi genomici aventi una struttura meno complessa di quella dei geni).

CARATTERISTICHE DEL PROBLEMA:

•I “modelli” di cui abbiamo bisogno dovranno essere basati unicamente sulla sequenza genomica.•Hanno alcuni aspetti in comune con i modelli che abbiamo visto nella parte del corso dedicata alla filogenesi• Sono modelli probabilistici•Sono in grado di GENERARE SEQUENZE “CONFORMI AL MODELLO” (per questo vengono detti modelli generativi)

IDENTIFICAZIONE DI ELEMENTI FUNZIONALI IN UNA SEQUENZA Bio CS

ISOLE CpGBio

Nel genoma umano i dinucleotidi CpG sono più rari di quanto saremmo portati ad aspettarci sulla base delle probabilità indipendenti di osservare una C o una G a causa di un motivo biologico. Le coppie CpG vanno incontro ad un processo di metilazione che modifica la citosina:

C T

Le regioni promotore sono ricche in CpG in quanto i promotori non vengono metilati e quindi, in essi, la frequenza di mutazione C T è molto minore che nel resto del genoma.

Se siamo in grado, data una sequenza genomica in formato FASTA, di identificare al suo interno una o più isole CpG queste possono essere viste come l’indicazione che, nelle “immediate” vicinanze, è presente un gene.


CATENE DI MARKOVper la ricerca di isole CpG


Una catena di Markov è una

collezione di stati e transizioni.

Seguendo le transizioni nel

diagramma si possono

generare tutte le sequenze

osservabili.

Questo è il classico diagramma di catena di Markov per la generazione di sequenze di DNA.

A

C G

T


• Ad ogni freccia nella figura è associato un valore che rappresenta la probabilità che ad un certo nt ne segua un altro.

• Aat = P(xL=T | xL-1=A)

A

C G

T

Probabilità che A sia seguita da T in una sequenza diDNA.

CATENE DI MARKOV : PROBABILITA’Bio

La probabilità di una sequenza x può essere definita come:

P(x) = P(xL,xL-1,…,x1)

Applicando P(X,Y) = P(X|Y)P(Y) più volte:

= P(xL|xL-1,…,x1)P(xL-1|xL-2,…,x1)…P(x1)

REGOLA DI MARKOV:

La probabilità di osservare un dato evento non dipenda da tutta la sequenza degli eventi precedenti ma solo da N eventi precedenti (nel caso più semplice N=1)


ES. CATENA DI MARKOV DEL PRIMO ORDINE :

L’osservazione in posizione i-esima dipendeunicamente da quello che osserviamo in posizione i-1.


ES. CATENA DI MARKOV DEL SECONDO ORDINE :

L’osservazione in posizione i-esima dipendeunicamente da quello che osserviamo in posizione i-1 ed in posizione i-2.


P(X) = P(xL|xL-1,…,x1)P(xL-1|xL-2,…,x1)…P(x1)

Regola di markov:

Dato che la probabilità di ogni osservazione dipende solo dall’osservazione precedente, possiamo semplificare in:

P(X) = P(xL|xL-1)P(xL-1|xL-2)…P(x2|x1)P(x1)

=

L

iix xax

i2

1 1)(P

X = seq. di lunghezza L

simbolo osservato in pos. i-1

simbolo osservato in pos. i

CATENE DI MARKOV : isole CpGBio

Proviamo ad usare 2 catene di Markov per trovare isole CpG in una sequenza genomica.

Ogni catena di Markov costituisce un modello che deve essere basato sulla probabilità di osservare il dinucleotide CpG in sequenze NOTE.

Dopo aver preparato due set di sequenze (CpG islands e NON_CpG islands) si creano due tabelle di contingenza riportanti le probabilità (osservate) che un qualsiasi nt y segua un nucleotide x.

Le tabelle ottenute saranno alla base delle catene di Markov da utilizzare per il nostro esperimento di discriminazione.


Xi-1

Xi

A C G T

A 0.18 0.27 0.43 0.12

C 0.17 p+(C | C) 0.274 p+(T|C)

G 0.16 p+(C|G) p+(G|G) p+(T|G)

T 0.08 p+(C |T) p+(G|T) p+(T|T)

Il modello + (CpG): si basa sulla matrice di transizioni + = (a+

st), in cui: a+st = (probabilità che t segua s in un’isola CpG)


Il modello - (noCpG): si basa sulla matrice di transizioni - = (a-

st), in cui: a-st = (probabilità che t segua s al di fuori di un’isola

CpG)

Xi-1

Xi

A C G T

A 0.3 0.2 0.29 0.21

C 0.32 p-(C|C) 0.078 p-(T|C)

G 0.25 p-(C|G) p-(G|G) p-(T|G)

T 0.18 p-(C|T) p-(G|T) p-(T|T)

CATENE DI MARKOV : isole CpG (confronto tabelle di transizione)Bio

Xi-1

Xi

A C G T

A 0.180 0.270 0.426 0.120

C 0.171 0.368 0.274 0.188

G 0.161 0.339 0.375 0.125

T 0.079 0.335 0.384 0.182

+

Xi-1

Xi

A C G T

A 0.300 0.205 0.285 0.210

C 0.322 0.298 0.078 0.302

G 0.248 0.246 0.298 0.208

T 0.177 0.239 0.292 0.292

-

CATENE DI MARKOV : P che una seq. derivi da un’isola CpG (I)Bio

Data la sequenza GCG, calcolare la probabilità che essa appartenga ad un’isola CpG o meno.

P(G+C+G+) P(G-C-G-)

Log-odds ratios

A C G T

A -0.740 0.419 0.580 -0.803

C -0.913 0.302 1.812 -0.685

G -0.624 0.461 0.331 -0.730

T -1.169 0.573 0.393 -0.679

Calcoliamo la tabella riportante il logaritmo dei rapporti di frequenza:

+aa = 0.180 log2(0.180/0.300) = -0.740-aa = 0.300

CATENE DI MARKOV : P che una seq. derivi da un’isola CpG (II)Bio

Le catene di Markov che stiamo utilizzando come modello sono catene del primo ordine e quindi la probabilità di osservare un nucleotide in posizione i dipende UNICAMENTE da ciò che si osserva in posizione i-1.

Calcoliamo P(Ci|Gi-1)P(Gi|Ci-1) nelle due catene(Dato che siamo in spazio logaritmico i prodotti diventano somme)

0.461+1.812 = LogRatio(Ci|Gi-1)+LogRatio(Gi|Gi-1) = 2.273

Ora torniamo alle probabilità relative:

P+(GCG)/P-(GCG) = 22.273 = 4.461In definitiva P+(GCG) = 4.461P-(GCG)

Dato che gli stati sono solo 2 (+/-) e la somma delle probabilità è 1…

P(+) + P(-) = 1 quindi 4.461 P(-) + (P-) = 1

P(-)= 17% , P(+)= 83%

CATENE DI MARKOV : PROBLEMI PRATICIBio

Il nostro test è basato sul calcolo di un rapporto di probabilità di unasequenza osservata IN STATI NOTI … quando stiamo analizzandouna sequenza genomica NON ANNOTATA …

P(G+C+G+) P(G-C-G-)

Dove inizia e dove finisce la regione da testare ?

NON SAPPIAMONON SAPPIAMO IN CHE

STATO CI TROVIAMO !!!

genome

? start

? end


Questo equivale a dire che quando testiamo una regione genomica come potenziale isola CpG siamo NOI a scegliere la regione da valutare (le sue estremità e la sua estensione).

Quello che vogliamo ottenere, invece, è un modello che valuti TUTTA la sequenza del genoma e decida al posto nostro che dalla posizione a alla posizione b è presente (potenzialmente) un’isola CpG.

In definitiva, per quanto statisticamente corretta la soluzione appena presentata NON E’ applicabile per la predizione delle isole CpG in una sequenza genomica non annotata!


Possibile soluzione: scegliere una dimensione w FISSA (es. 1000 nt) ed utilizzare una finestra scorrevole sul genoma valutando il risultato del test ad ogni passo.

Questo approccio è affetto da diversi problemi:

Le isole CpG hanno lunghezze diverse e quindi una finestra di estensione fissa potrebbe essere inadatta per esprimere un giudizio.- Se scegliamo una dimensione della finestra w troppo piccola, tenderemo a giudicare ogni occorrenza del dinucleotide CG come proveniente da un’isola CpG.- Se scegliamo una dimensione di w troppo grande non avremo abbastanza potere di discriminazione.


Soluzione alternativa :

Incorporare i due modelli (le due catene di Markov) in un UNICO modello.



Vantaggi:Questo modello riflette meglio la realtà. Non c’è più dipendenza dalla dimensione della finestra dato che il modello prevede ad ogni step la possibilità di transire da un qualsiasi nucleotide in stato + (isola CpG) ad un nucleotide in stato – (non isola CpG).

ATTENZIONE:La maggior flessibilità del modello

ha un costo!

HIDDEN MARKOV MODELS ( HMM )Bio


Svantaggi:Il modello contiene non uno ma DUE possibili stati ( + / - ) per ogni possibile nucleotide! Osservando la sequenza non sappiamo più quale stato ha emesso un certo nt. In definitiva gli stati sono NASCOSTI !

Da qui il nome :

HIDDEN Markov Models


CARATTERISTICHE degli HMM :

EVALUATION:Siamo liberi di definire il modello in modo da adattarlo il più possibile al tipo di elementi funzionali che vogliamo modellare ma dovremo comunque rispondere a questa domanda: Data una sequenza qual è la PROBABILITA’ che essa sia stata generata dal modello?

DECODING:Data una sequenza decidere, per ogni posizione QUAL’E’ LO STATO (NASCOSTO) PIU’ PROBABILE.

TRAINING:Data una collezione di sequenze annotate come possiamo utilizzarle per addestrare un HMM ?

HIDDEN MARKOV MODELS ( HMM )Bio CARATTERISTICHE degli HMM :

Il modello su cui stiamo ragionando ha una caratteristica particolare … come le catene di Markov è composto UNICAMENTE DA PROBABILITA’ DI TRANSIZIONE.In questo modo siamo OBBLIGATI a creare un NUMERO DI STATI PARI a N * dim(alfabeto), ossia Numero_Stati * Numero_Nucleotidi (nel nostro caso 2 * 4 = 8). E’ molto scomodo. Gli Hidden Markov Models, per risolvere questo problema, prevedono l’esistenza di DUE TIPI DI PROBABILITA’ : probabilità di transizione e probabilità di emissione.

P_transz = probabilità di transizione tra STATIP_em = probabilità di emettere un dato simbolo quando ci troviamo in un certo stato


• Consideriamo un casinò in cui viene utilizzato, per la maggior parte del tempo, un dado onesto ma, occasionalmente, si passa ad un dado truccato.

• Il dado truccato ha una probabilità pari a 0.5 di fare 6 e 0.1 per gli altri numeri da uno a cinque. Il casinò passa dal dado onesto al dado truccato con una probabilità di 0.05 e torna al dado onesto con una probabilità doppia a quella della transizione inversa (0.1).

• Possiamo esprimere questa situazione attraverso il seguente modello:

1: 1/62: 1/63: 1/64: 1/65: 1/66: 1/6

1: 1/102: 1/103: 1/104: 1/105: 1/106: 1/2

Fair Loaded

0.05

0.1

0.95 0.9

In un dato momento non sappiamo quale dado sta

usando il casinò.


Domanda:

Data una serie di osservazioni (e un HMM) siamo in grado di decidere quali derivano dall’utilizzo del dado onesto e quali dal dado truccato?

Serie di osservazioni:

3 1 5 1 1 6 6 6 6 3 6 2 4

Il casinò disonesto: DECODING


DECODING:

Decodificare una sequenza significa trovare, tra tutte quelle possibili, la serie di stati in grado di generare, con probabilità massima, le osservazioni .

Un HMM assume che la sequenza osservata sia stata generata da un processo markoviano. A causa della regola di Markov questi processi sono detti memory-less (si “ricordano” solo di N passi precedenti a quello considerato).

Il problema che vogliamo risovere può essere suddiviso in un numero finito di sottoproblemi ripetitivi

Può essere risolto mediante programmazione dinamica

Il casinò disonesto: DECODING


Diagramma TRELLIS:

Il diagramma trellis permette di rappresentare una serie temporale di transizioni tra un set finito di stati. Supponiamo di avere 2 possibili stati di un HMM ed una serie di osservazioni. Esse si possono rappresentare così:

osservazioni

Ad ogni step vengono rappresentatetutte le possibili transizioni (4).


DECODING definizione del problema:

Dato un HMM ed una sequenza di osservazioni trovare la serie di stati (s*) del modello in grado di generare con probabilità massima la serie di osservazioni.

Data una sequenza osservabile x = (x1,…,xL),

La più probabile sequenza si stati s* = (s*1,…,s*

L) è quella che massimizza p(s|x).

1( ,..., )

* *1 1 1* ( ,..., ) ( ,..., | ,..., )maxarg

Ls s

L L Ls s s p s s x x

NB: molto spesso la seq. di stati si indica con la lettera greca ( ad indicare “path” … percorso)


DECODING : Algoritmo di Viterbi

Un problema pratico di questo algoritmo, comune a tutti gli algoritmi che moltiplicano probabilità, è il rischio di underflow. Per evitare questo problema dobbiamo portare le probabilità dell’HMM in spazio logaritmico e poi sostituire le moltiplicazioni con somme.

Come possiamo convertire le probabilità in log(probabilità)?

HIDDEN MARKOV MODELS ( HMM )BioDECODING : P in spazio logaritmico

P transizione:

La somma di tutte le probabilità di transizione in uscita da uno stato deve essere pari a 1. Nel caso in cui tutte le possibili transizione fossero equiprobabili le loro probabilità sarebbero tutte uguali a 1/K (k = numero degli stati)

1: 1/62: 1/63: 1/64: 1/65: 1/66: 1/6

1: 1/102: 1/103: 1/104: 1/105: 1/106: 1/2

Fair Loaded

0.05

0.1

0.95 0.9FL

LF

FF LL

HIDDEN MARKOV MODELS ( HMM )BioDECODING : P in spazio logaritmico

P transizione:

La somma di tutte le probabilità di transizione in uscita da uno stato deve essere pari a 1. Nel caso in cui tutte le possibili transizione fossero equiprobabili le loro probabilità sarebbero tutte uguali a 1/K (k = numero degli stati)

F L

F 0.95 0.05

L 0.10 0.90

Se fossero equiprob. varrebbero ½ …

Se fossero equiprob. varrebbero ½ …


DECODING : P in spazio logaritmico

P transizione (matrice A):

A(FF) = ln(0.95 / 0.5) = 0.6418539

A(FL) = ln(0.05 / 0.5) = - 2.302585

A(LF) = ln(0.10 / 0.5) = - 1.609438

A(LL) = ln(0.90 / 0.5) = 0.5877867

F L

F 0.64 -2.30

L -1.61 0.59

somma = 1

somma = 1A



P emissione (FAIR) : matrice EF

EF(1) = ln( 1/6 / 1/6 ) = 0

EF(2) = ln( 1/6 / 1/6 ) = 0

EF(3) = ln( 1/6 / 1/6 ) = 0

EF(4) = ln( 1/6 / 1/6 ) = 0

EF(5) = ln( 1/6 / 1/6 ) = 0

EF(6) = ln( 1/6 / 1/6 ) = 0

somma = 1



P emissione (LOADED) : matrice EL

EL(1) = ln( 1/10 / 1/6 ) = -0.51

EL(2) = ln( 1/10 / 1/6 ) = -0.51

EL(3) = ln( 1/10 / 1/6 ) = -0.51

EL(4) = ln( 1/10 / 1/6 ) = -0.51

EL(5) = ln( 1/10 / 1/6 ) = -0.51

EL(6) = ln( 1/2 / 1/6 ) = 1.10

somma = 1



1: 02: 03: 04: 05: 06: 0

1: -0.512: -0.513: -0.514: -0.515: -0.516: 1.10

Fair Loaded

-2.30

-1.61

0.64 0.59

FL

LF

FFLL


DECODING : Algoritmo di Viterbi

• E’ basato sull’idea di risolvere un problema intrattabile suddividendolo in tanti sottoproblemi trattabili (progr. dinamica).

• Risolve il problema del decoding: data una sequenza di osservazioni calcola la serie di stati che ha generato la sequenza di osservazioni con probabilità massima.

• E’ di fondamentale importanza definire i sottoproblemi in modo da permettere di risolverli in maniera iterativa. Mano a mano che l’algoritmo procede I nuovi sottoproblemi verranno risolti sfruttando le soluzioni dei sottoproblemi precedentemente incontrati.

HMM : DECODING (VITERBI, dyn.prog.)Bio

Dato che per ogni stato k, e per una posizione fissata i nella

seq.,

Vk(i) = max{1… i-1} P[x1…xi-1, 1, …, i-1, xi, i = k]

prob. di esserein stato k allaposizione i

Massimizzazione probabilità seq. distati da posizione 1a posizione i-1

Osservazioni da posizione 1a posizione i-1

Seq. Stati da posizione 1a posizione i-1

Emissione osservazionein pos. i

Assumendo di essere instato k

Come calcolare Vl(i+1) ?

HMM : DECODING (VITERBI, dyn.prog.)Bio Dato che per ogni stato k,

e per una posizione fissata i nella seq.,


Per definizione,

Vl(i+1) = max{1… i}P[ x1…xi, 1, …, i, xi+1, i+1 = l ] = max{1… i}P(xi+1, i+1 = l | x1…xi, 1,…, i) P[x1…xi, 1,…, i]

= max{1… i}P(xi+1, i+1 = l | i ) P[x1…xi-1, 1, …, i-1, xi, i] = maxk [P(xi+1, i+1 = l | i=k) max{1… i-1}P[x1…xi-1,1,…,i-1, xi,i=k]] = maxk [ P(xi+1 | i+1 = l ) P(i+1 = l | i=k) Vk(i) ] = el(xi+1) maxk akl Vk(i)

Come calcolare Vl(i+1) ?

emissione transizione Viterbi stato k posizione i



Vl(i+1) = max{1… i}P[ x1…xi, 1, …, i, xi+1, i+1 = l ]

= el(xi+1) maxk akl Vk(i)

emissione transizione Viterbi stato k posizione i

RICORSIONE


Input: x = x1……xL

Inizializzazione:V0(0) = 1 (0 è una prima posizione fittizia)Vk(0) = 0, per ogni k > 0

Iterazione:Vj(i) = ej(xi) maxk akj Vk(i – 1)

Ptrj(i) = argmaxk akj Vk(i – 1)

Terminazione:P(x, *) = maxk Vk(L)

Traceback: L* = argmaxk Vk(L) i-1* = Ptri (i)


DECODING : calcolo

Prima di descrivere i passaggi dell’algoritmo di Viterbi dobbiamo risolvere un problema: qual’è la struttura dati che permette di effettuare i calcoli in maniera più efficiente possibile?

Abbiamo già considerato un problema che richiedeva l’utilizzo della programmazione dinamica (allineamento sequenze) e, in quel caso, la struttura dati fondamentale era una matrice.

Anche in questo caso è possibile utilizzare una matrice?


DECODING : calcolo

Le probabilità rappresentate da questo diagramma trellis possono essere associate alle celle di una matrice K x L, dove K = numero di possibili transizioni e L = numero di simboli che compongono la sequenza di osservazioni.

K = 4 (FF,FL,LF,LL)L = 6


Dato un HMM

E una sequenza: 3 1 5 1 1 6 6 6 6 6 6 2 4

1: 02: 03: 04: 05: 06: 0

1: -0.512: -0.513: -0.514: -0.515: -0.516: 1.10

Fair Loaded

-2.30

-1.61

0.64 0.59

FL

LF

FFLL

HMM : DECODING (VITERBI, dyn.prog.)PRIMA COLONNA

BioSt. 3 1 5 1 1 6 6 6

FF

0.64

LF

-1.61

FL

-2.30

LL

0.59

F

L

Obs.

0

-0.51

V_F 0V_L -0.51

p_F Bp_L B

Assumiamo uguali le probabilità Iniziali di transire in stato F o L.

F = E_f(0) + 0 + 0 = 0L = F_l(0) + 0 + 0 = -0.51

BioSt. 3 1 5 1 1 6 6 6

FF

0.64

LF

-1.61

FL

-2.30

LL

0.59

F

L

Obs.

0

-0.51

V_F 0V_L -0.51

p_F Bp_L B

HMM : DECODING (VITERBI, dyn.prog.)el(xi+1) maxk akl Vk(i)

F F(1) = 0 + 0.64 + 0 = 0.64L F(1) = 0 -1.61 -0.51 = -2.12F L(1) = -0.51 -2.30 + 0 = -2.81L L(1) = -0.51 0.59 -0.51 = -0.43

Emiss. 1 in stato F

Max F step preced.

0.64

-2.12

-2.81

-0.43

Max path to F : FF (valore: 0.64)

BioSt. 3 1 5 1 1 6 6 6

FF

0.64

LF

-1.61

FL

-2.30

LL

0.59

F

L

Obs.

0

-0.51

V_F 0 0.64V_L -0.51 -0.43

p_F B FFp_L B LL


0.64

-2.12

-2.81

-0.43

F F F

L L L

el(1) = 0

3B 1

FF0.64

0

BioSt. 3 1 5 1 1 6 6 6

FF

0.64

LF

-1.61

FL

-2.30

LL

0.59

F

L

Obs.

0

-0.51

V_F 0 0.64V_L -0.51 -0.43

p_F B FFp_L B LL


F F(5) = 0 + 0.64 + 0.64 = 1.28L F(5) = 0 -1.61 -0.43 = -2.04F L(5) = -0.51 -2.30 +0.64 = -2.17L L(5) = -0.51 +0.59 -0.43 = -0.35

Emiss. 5 in stato F

Max F step preced.

0.64

-2.12

-2.81

-0.43

1.28

-2.12

-2.81

-0.43

BioSt. 3 1 5 1 1 6 6 6

FF

0.64

LF

-1.61

FL

-2.30

LL

0.59

F

L

Obs.

0

-0.51

V_F 0 0.64 1.28V_L -0.51 -0.43 -0.35

p_F B FF FFp_L B LL LL


F F(5) = 0 + 0.64 + 0.64 = 1.28L F(5) = 0 -1.61 -0.43 = -2.04F L(5) = -0.51 -2.30 +0.64 = -2.17L L(5) = -0.51 +0.59 -0.43 = -0.35

0.64

-2.12

-2.81

-0.43

1.28

-2.04

-2.17

-0.35

F F F

L L L

13 5

FF0.640.64

-2.12

el(5) = 0

Bio

St. 3 1 5 1 1 6 6 6

F

L

Obs.

0

-0.51

p_F B FF FF FF FF FF FF FFp_L B LL LL LL LL LL FL LL

HMM : DECODING (VITERBI, dyn.prog.)Termination: Max prob ultima colonna

0.64

-2.12

-2.81

-0.43

1.28

-2.12

-2.81

-0.43

1.92

-1.96

-1.53

-0.27

2.56

-1.88

-0.89

-0.19

3.20

-1.80

-0.25

-0.11

3.84

-1.72

2.00

1.58

4.48

0.39

2.64

3.69

6

5.12

2.08

3.28

5.38

NBB: STATO FINALE = L !

Bio

St. 3 1 5 1 1 6 6 6

F

L

Obs.

0

-0.51

p_F B FF FF FF FF FF FF FFp_L B LL LL LL LL LL FL LL LL

HMM : DECODING (VITERBI, dyn.prog.)TRACEBACK

0.64

-2.12

-2.81

-0.43

1.28

-2.12

-2.81

-0.43

1.92

-1.96

-1.53

-0.27

2.56

-1.88

-0.89

-0.19

3.20

-1.80

-0.25

-0.11

3.84

-1.72

2.00

1.58

4.48

0.39

2.64

3.69

6

5.12

2.08

3.28

5.38

PARSE : 3 1 5 1 1 6 6 6 6 F F F F F F L L L


Algoritmo di Viterbi: riepilogo (I)

Input: modello , sequenza

Init:

semplice, ma deve contenere uno stato BEGIN che permetta di modellare l’inizio della sequenza

Iteration:

emission(state2_pos_i) +

transition(state1state2) +

maxVITk(pos_i-1)


Algoritmo di Viterbi: riepilogo (II)

Iteration (per ogni transizione):

emission(state2_pos_i) +

transition(state1state2) +

maxVITk(pos_i-1) SALVARE IL VALORE MASSIMO PEROGNI STATO k (serve per avere a dispo-sizione il valore maxVITk(pos_i-1) al passo successivo)

SALVARE LO STATO DA CUI SI PARTE(in posizione i-1) PER ARRIVARE AL VALORE MASSIMO PER OGNI STATO k

Serve per il TRACEBACK,Ossia per produrre la sequenza di stati che haGenerato le osservazioniCon prob. Massima.

HMM : DECODING (VITERBI, dyn.prog.)BioAlgoritmo di Viterbi: riepilogo (III)

Caratteristiche dell’algoritmo:

• Basato su programmazione dinamica• Complessità temporale: O(K2L)• Complessità spaziale: O(KL)

Biologia computazionale A.A. 2010-2011 semestre II U NIVERSITÀ DEGLI STUDI DI MILANO Docente:...

Documents

Transcript of Biologia computazionale A.A. 2010-2011 semestre II U NIVERSITÀ DEGLI STUDI DI MILANO Docente:...