Bioinformatica
Catene di Markov, HMM – GenScan – EasyBack
Dr. Giuseppe Pigola – [email protected]
Catene di Markov
Una catena di Markov del primo ordine è una tripla ( S ,  , A) , dove:
 S={s1,s2,…,sN} è un insieme finito di stati (eventi);
  i  p(q1  si ) 1  i  N è la probabilità iniziale degli stati (è
rappresentato da un vettore);
 A è un insieme di probabilità di transizione tale che:

aij  p(qt  s j | qt 1  si )
1  i, j  N
Dove con qt indichiamo lo stato al tempo t (discreto);
A è rappresentata da una matrice.Valgono le seguenti:
aij  0
N
a
j 1
2
ij
1
Bioinformatica
Catene di Markov

La catene di Markov del primo ordine sono anche dette catene di Markov
memoryless dato che la probabilità che avvenga un evento dipende solo
dall’evento che si è verificato all’istante precedente;

In una catena di Markov di ordine k la probabilità che avvenga un evento
dipende dai k precedenti stati;
3
Bioinformatica
Catene di Markov

In una catena di Markov del primo ordine, se vogliamo calcolare la
probabilità di una serie di stati

Dato che il valore al’istante i dipende solo dall’istante precedente
p(qt | q1 ,..., qt 1 )  p(qt | qt 1 ) 
Ne segue che

4
Bioinformatica
Catene di Markov

In una catena di Markov del primo ordine, se vogliamo calcolare la
probabilità di una serie di stati

Dato che il valore al’istante i dipende solo dall’istante precedente

Ne segue che
5
Bioinformatica
Catene di Markov
Esempio:


P(Sun, Rain, Rain, Rain, Snow, Snow) =
P(Sun) P(Rain | Sun) P(Rain | Rain)P(Rain | Rain) P(Snow | Rain) P(Snow | Snow)
6
Bioinformatica
Catene di Markov
Catena di Markov per stringhe di DNA:

7
Bioinformatica
Catene di Markov
DEFINIZIONE

Se nella diagonale della matrice di transizione troviamo un 1, allora lo stato
contrassegnato da questo valore viene chiamato stato assorbente (Se si arriva a quello
stato non si uscira più);
DEFINIZIONE

Una catena di Markov si dice Finita se è formata da un insieme finito di stati;
DEFINIZIONE

Una catena di Markov si dice Aperiodica se gli stati non vengono osservati mai in
modo periodico;
DEFINIZIONE

Una catena di Markov si dice Irriducibile se tutti gli stati prima o poi vengono raggiunti
(non ci sono zeri nella matrice di transizione);
8
Bioinformatica
Catene di Markov
ESEMPIO: CATENA DI MARKOV PERIODICA

9
Bioinformatica
Catene di Markov

ESEMPIO
10
Bioinformatica
Catene di Markov

ESEMPIO
11
Bioinformatica
Catene di Markov

ESEMPIO: Considerazioni sulla matrice di transizioni.

Moltiplicando la matrice per se stessa ci riferiamo all’istante t+2.
12
Bioinformatica
Catene di Markov

ESEMPIO: Se una persona compra attualmente Coca, quale sarà la
probabilità che tra tre volte comprerà Pepsi?
13
Bioinformatica
Catene di Markov

ESEMPIO: Se una persona compra attualmente Coca, quale sarà la
probabilità che tra tre volte comprerà Pepsi?
14
Bioinformatica
Catene di Markov

Possiamo rappresentare la matrice di transizione attraverso autovalori e
autovettori e sfruttare questi per evitare di fare i prodotti di matrici.

Suponiamo di avere una matrice P = s x s con s autovalori distinti (ottenuti
calcolando det(P-I)=0):

Per le proprietà di autovalori e autovettori, ogni autovalore ha associato un
autovettore destro R e un autovettore sinistro L tali che
15
Bioinformatica
Catene di Markov

AUTOVALORI

AUTOVETTORI
16
Bioinformatica
Catene di Markov

Se vale

La matrice si può scrivere come (operazione Spettrale):

E quindi
17
Bioinformatica
Catene di Markov

Una persona compra una bibita alla settimana;

Inizialmente abbiamo una distribuzione in cui il 60% delle persone compra
Coca e il rimanente 40% compra Pepsi;

Quale è la frazione di persone che comprerà Coca fra tre settimane?

Vogliamo conoscere la distribuzione degli acquisti in un certo istante (X=3)
a partire da una data distribuzione iniziale;
18
Bioinformatica
Catene di Markov

Ci serve

Dobbiamo considerare la colonna relativa alla Coca e moltiplicare per la
distribuzione iniziale
19
Bioinformatica
Catene di Markov


Generalizzando:

D0 : Distribuzione iniziale degli stati;

Di : Distribuzione nella settimana i (Di = D0*Pi);
Proprietà: All’aumentare di i, la distribuzione di probabilità varia fino a
raggiungere valori stazionari (Distribuzione Stazionaria);
20
Bioinformatica
Catene di Markov –
Distribuzione Stazionaria

D0 * P = D 1

D 0 * P2 = D 2

D 0 * Pn = D n

Una distribuzione

Riferendoci al generico elemento della distribuzione si avrà

Mentre riferendoci alla catena di Markov

Si dimostra che la distribuzione stazionaria corrisponde all’autovettore sinistro associato al
primo autovalore (quello dominante), cioè
21
Distribuzione di probabilità degli stati dopo n istanti;
si dice stazionaria se
Bioinformatica
Catene di Markov –

Distribuzione Stazionaria
Consideriamo ora la probabilità di transizione dallo stato i allo stato j. Per
Bayes
Si ha
dove

rappresenta la probabilità di transizione dallo stato i allo stato j in una
catena di Markov identica in cui però la successione di istanti decresce (di
tipo backward) e quindi equivale a
22
Bioinformatica
Catene di Markov –
Distribuzione Stazionaria

Per valori di t sufficientemente elevati vale

quindi

Quest’ultima condizione esprime la condizione necessaria e sufficiente
affinché sia la distribuzione stazionaria di una catena di Markov avente P
come matrice di transizione.
23
Bioinformatica
Catene di Markov –
Distribuzione Stazionaria

Se una catena di Markov è finita, aperiodica e irriducibile, allora essa ha una
distribuzione stazionaria.

Possiamo calcolare la distribuzione stazionaria risolvendo semplicemente il
sistema:
24
Bioinformatica
Catene di Markov –
Markov Chain Montecarlo

Data una catena di Markov finita, aperiodica e irriducibile possiamo
determinare la distribuzione stazionaria associata;

Viceversa, data una distribuzione stazionaria possiamo costruire una catena
di Markov che converga alla distribuzione stazionaria (Markov Chain
Montecarlo):

Metodo Hasting Metropolis;

Gibbs Sampling;
25
Bioinformatica
Catene di Markov –
Hasting Metropolis

In questo caso abbiamo a disposizione una distribuzione stazionaria

Consideriamo delle costanti (definite in modo random) come segue

Definiamo inoltre

e
P

Si dimostra che la matrice di probabilità P è la matrice di transizione di una
catena di Markov che ha per distribuzione stazionaria quella data.
26
Bioinformatica
Catene di Markov –
Gibbs Sampling

Sia x=(x1,x2,…,xn) una variabile casuale e sia y un vettore di k valori scelti a
caso tra x1,x2,…,xn con possibili ripetizioni. Sia Px(x) la distribuzione di
probabilità della variabile x.

Definiamo la catena di Markov i cui stati sono tutti i possibili valori di y.

Siano i vettori i e j due stati della catena. Definiamo la probabilità come:


27
se i vettori i e j differiscono su più componenti;
Se i e j differiscono su al più una componente h:
Bioinformatica
Catene di Markov –

Gibbs Sampling
Si dimostra che Px(x) è la distribuzione stazionaria della catena di Markov
che ha per matrice di transizione la matrice P.
28
Bioinformatica
Catene di Markov –
Allineamento

Supponiamo di avere N sequenze L1,L2,…,Ln in un alfabeto  (ad es.
Sequenze aminoacidiche);

Vogliamo trovare N segmenti di lunghezza w tali da massimizzare la
similarità tra le sequenze
29
Bioinformatica
Catene di Markov –
Allineamento

In ogni sequenza i possibili segmenti di lunghezza w sono

Un algoritmo brute-force dovrebbe considerare tutte le combinazioni di
segmenti di lunghezza w nelle sequenze

Consideriamo invece una catena di Markov con S stati in cui ogni stato
rappresenta una scelta di segmenti nelle N sequenze (e quindi a un
possibile allineamento);
30
Bioinformatica
Catene di Markov –
Allineamento

Scegliamo uno stato iniziale casuale (scegliamo uno dei possibili allineamenti
a caso);

Consideriamo la matrice A NxW dei segmenti scelti;

Rimuoviamo ad ogni passo da A il segmento relativo ad una sequenza e lo
sostituiamo con un altro della stessa sequenza in base alla probabilità
definita in modo opportuno;

In questo modo passiamo da uno stato ad un altro della catena;
31
Bioinformatica
Catene di Markov –
32
Allineamento
Bioinformatica
Catene di Markov –
Allineamento

Ad ogni carattere j presente nell’allineamento è associata una probabilità pj
(frequenza nelle sequenze di input);

Definiamo la probabilità di transizione come

Dove:


cij il numero di volte che il carattere j compare nella colonna i all’interno dell’array
ridotto;
bj è detta probabilità di background e serve a fare in modo che qij non sia nulla. In
generalmente è inizializzata a

33
Bioinformatica
Catene di Markov –
Allineamento

L’obbiettivo ad ogni step è quello di fare una sostituzione di segmenti che
migliori la qualità dell’allineamento;

Sia

Definiamo:
il segmento che deve rimpiazzare quello cancellato;
P prodotto delle probabilità di background
Q probabilità stimata dei caratteri di x nell’insieme
ridotto di n-1 sequenze
34
Bioinformatica
Catene di Markov –
Allineamento

Definiamo il likelyhood_ratio di x come

Allora la probabilità di transizione dallo stato s (quello contenente la
vecchia sequenza) allo stato u (quello contenente x) come


Lh è la sequenza contenente x;
Il denominatore è la somma di tutti i LR di tutti i possibili segmenti di
lunghezza w in Lh;

In altre parole P rappresenta la probabilità di scegliere il segmento x;
35
Bioinformatica
Catene di Markov –
Allineamento

Come scegliamo il segmento x tra le

La scelta più ragionevole sembra essere quella di scegliere il segmento x
con LR massimo: Questo suggerirebbe che la frequenza di caratteri di x
nella popolazione ridotta delle sequenze è molto più alta della frequenza
dei caratteri x nella popolazione totale (mediamente);

In tal caso si avrebbe appunto

e quindi LR alto;
36
possibilità???
Bioinformatica
Catene di Markov –

Allineamento
In realtà è più conveniente scegliere x in modo caso (seguendo la probabilità di
transizione dei possibili stati successivi)!!!

E’ più veloce (non dobbiamo calcolare probabilità P);

E’ più accurato (il metodo precedente potrebbe portare a massimi locali)!!!!

Dietro questa scelta apparentemente senza senso ci sono in realtà dei fondamenti
matematici basati sul concetto di entropia relativa;

Se osserviamo bene come è definita la catena di Markov ci accorgiamo che
rientriamo nella definizione di catena di Markov del Gibbs Sampling;

Ne segue che la Catena di Markov converge a una distribuzione stazionaria (anche
se scegliamo ad ogni istante uno stato in modo casuale tra i possibili);
37
Bioinformatica
Catene di Markov –
Allineamento

Come visto in precedenza

Dove cij(s) è il numero di volte che il carattere j si presenta nella colonna i nello
stato s;

Si dimostra che Il vettore  definito da

Rappresenta la distribuzione stazionare della nostra catena di Markov.
38
Bioinformatica
Hidden Markov Models

Un HMM è definito da una quintupla ( S , V ,  , A, B ) , dove:
 S={s1,s2,…,sN} è un insieme finito di stati nascosti;
 V={v1,v2,…,vM} è un insieme di stati osservabili;
  i  p(q1  si ) 1  i  N è la probabilità iniziale degli stati (è
rappresentato da un vettore);
 A è un insieme di probabilità di transizione tale che
aij  p(qt  s j | qt 1  si )

Dove con qt indichiamo lo stato al tempo t (discreto);
B rappresenta le probabilità di emissione tale che:
b j (k )  p (vk at t | qt  s j )

39
1  i, j  N
1 j  N
1 k  M
Anche B è rappresentata da una matrice.
Bioinformatica
Hidden Markov Models

Esempio:

Data la sequenze di stati osservabili, quale è la sequenza di stati nascosti più
probabile che l’ha generata?
40
Bioinformatica
Catene di Markov

ESEMPIO: Cerchiamo un modello per discriminare due ipotetiche regioni
di DNA (Sopponiamo per semplicità di avere solo due nucleotidi).
I
II
ATTA
TTAT
AAAT
TAAT
…
TTAA
TATA
ATTA
ATAT
ATTT
…
ATAT
Probabilità di transizione dalla
regione I alla II con la sequenza TT
41
Bioinformatica
Catene di Markov

A quale regione appartiene la sequenza TTAT ?

Calcoliamo la probabilità di tutte le possibili sequenze di nucleotidi
appartenenti alle due regioni.

TITIAITI=1.1x10-1
TITIIAITI=1.8x10-3
TIITIAITI=6.0x10-3
TIITIIAITI=9.0x10-3
TITIAITII=8.8x10-3
TITIIAITII=1.4x10-4
TIITIAITII=4.8x10-4
TIITIIAITII=7.2x10-4
TITIAIITI=5.5x10-4
TITIIAIITI=1.0x10-3
TIITIAIITI=3.0x10-5
TIITIIAIITI=5.2x10-3
TITIAIITII=1.4x10-4
TITIIAIITII=8.4x10-3
TIITIAIITII=2.4x10-4
TIITIIAIITII=4.2x10-2
Risulta più probabile che la sequenza appartiene integralmente alla
regione I
42
Bioinformatica
Catene di Markov

ESEMPIO: ISOLE C-G
Esistono delle zone di DNA che evidenziano una presenza superiore d coppie CG;
Nucleotidi appartenenti a isole C-G hanno una diversa probabilità;
43
Bioinformatica
Hidden Markov Models
Dato un HMM λ e la sequenza di osservazioni la sequenza di osservazioni
O=o1,o2,…,oT, possiamo affrontare tre problemi:

Evaluation: Quale è la probabilità di ottenere O nel modello, p(O| λ)?

Decoding: Quale è la corrispondente sequenza di stati nascosti
Q=q1,q2,…,qT che ha generato O?

Learning: Come possiamo “aggiustare” i parametri del modello λ per
massimizzare P(O| λ)?
44
Bioinformatica
Hidden Markov Models
ESEMPIO: Il casinò disonesto.
In un casinò ci sono due dadi di cui uno truccato:
 Dado non truccato
P(1) = P(2) = P(3) = P(4) = P(5) = P(6) = 1/6
 Dado truccato
P(1) = P(2) = P(3) = P(4) = P(5) = 1/10
P(6) = ½
Il croupier passa dal dado non truccato a quello
truccato e viceversa.
45
Bioinformatica
Hidden Markov Models
ESEMPIO: Il casinò disonesto.
Quanto è probabile la seguente osservazione nel nostro modello?
1245526462146146136136661664661636616366163616515615115146123562344
Si tratta di un problema di Evaluation.
46
Bioinformatica
Hidden Markov Models
ESEMPIO: Il casinò disonesto.
Quale porzione di questa sequenza è stata prodotta dal dado truccato e
quale dal dado non truccato?
1245526462146146136136661664661636616366163616515615115146123562344
Si tratta di un problema di Decoding.
47
Bioinformatica
Hidden Markov Models
ESEMPIO: Il casinò disonesto.
Quanto “truccato” era il dado truccato? Quanto spesso i due dadi venivano
cambiati?
1245526462146146136136661664661636616366163616515615115146123562344
Si tratta di un problema di Learning.
48
Bioinformatica
Hidden Markov Models
ESEMPIO: Il casinò disonesto. Modello di Markov Nascosto
0.05
0.95
FAIR
P(1|F) = 1/6
P(2|F) = 1/6
P(3|F) = 1/6
P(4|F) = 1/6
P(5|F) = 1/6
P(6|F) = 1/6
49
Transizioni
0.95
LOADED
0.05
Emissioni
P(1|L) = 1/10
P(2|L) = 1/10
P(3|L) = 1/10
P(4|L) = 1/10
P(5|L) = 1/10
P(6|L) = 1/2
Bioinformatica
Hidden Markov Models
EVALUATION
Vogliamo calcolare la probabilità dell’osservazione O=o1,o2,…,oT dato il modello λ.
Considerata una sequenza di stati nascosti fissata Q=q1,q2,…,qT la probabilità
dell’osservazione O per la sequenza di stati nascosti Q assumendo che le osservazioni
siano indipendenti, è data da:
T
p(O | Q,  )   p(ot | qt ,  )
Quindi
t 1
p(O | Q,  )  bq1 (o1 )bq2 (o2 )    bqT (oT )
La probabilità della sequenza degli stati nascosti è data da:
p(Q |  )   q1 aq1q2 aq2 q3    aqT 1qT
La probabilità che si verifichino contemporaneamente O e Q sarà il prodotto
p(O, Q |  )  p(O | Q,  ) p(Q |  )   q1 bq1 (o1 )aq1q2 bq2 (o2 )aq2 q3    aqT 1qT bqT (oT )
La probabilità di O nel modello sarà allora ottenuta sommando queste probabilità su
tutti i possibili stati nascosti.
p(O |  ) 
p(O | Q,  ) p(Q |  )

q1 ,q1 ,..., qT
50
Bioinformatica
Hidden Markov Models
EVALUATION
Interpretazione della formula:
p(O |  ) 
 p(O | Q,  ) p(Q |  )  
q1 ,q1 ,..., qT
b (o1 )aq1q2 bq2 (o2 )aq2q3    aqT 1qT bqT (oT )
q1 q1
q1 ,q1 ,..., qT
Inizialmente al tempo t=1 ci troviamo nello stato q1con probabilità  q e
generiamo il simbolo o1con probabilità bq (o1 )
1
1
Al tempo t=2 avremo una transizione allo stato q2 con probabilità a q1q2 e
verrà generato il simbolo osservabile o2 con probabilità bq (o2 )
2
Il processo continua fino ad arrivare al tempo T
Complessità: 2TNT calcoli
51
Bioinformatica
Hidden Markov Models
EVALUATION: FORWARD ALGORITHM
Consideriamo la variabile forward:
 t (i)  p(o1o2 ...ot | qt  si ,  )
Essa rappresenta la probabilità di aver osservato al tempo t o1o2…ot e di
trovarci nello stato nascosto si.
Procedendo induttivamente:
1 (i)   i bi (o1 )
1 i  N
Inizializzazione
N
 t ( j )  [  t 1 (i )aij ]b j (ot ) 2  t  T 1  j  N
i 1
N
P(O |  )    T (i )
Induzione
Terminazione
i 1
52
Bioinformatica
Hidden Markov Models
EVALUATION: FORWARD ALGORITHM
Lo stato sj può essere raggiunto al tempo
t+1 dagli N possibili stati al tempo t.
 t 1 (i)rappresenta la probabilità di
aver osservato la sequenza o1o2…ot-1 e di
trovarsi a tempo t-1 nello stato nascosto si.
Il prodotto t 1 (i) aij rappresenta la
probabilità di osservare o1o2…ot-1 e di
raggiungere lo stato sj al tempo t
proveniendo dallo stato si
Sommando su tutti i possibili stati avremo la probabilità di osservare sj al tempo t. Per
tener conto della probabilità di emissione dello stato osservabile ot moltiplichiamo per
la relativa probabilità di emissione
N
 t ( j )  [  t 1 (i )aij ]b j (ot ) 2  t  T 1  j  N
i 1
53
Bioinformatica
Hidden Markov Models
EVALUATION: FORWARD ALGORITHM
Lo step finale ci darà la probabilità cercata
N
P(O |  )    T (i )
i 1
Complessità:
O(N2T) Tempo
O(NT) Spazio
54
Bioinformatica
Hidden Markov Models
EVALUATION: BACKWARD ALGORITHM
In modo speculare consideriamo la variabile backward:
t (i)  p(ot 1ot 2 ...oT , | qt  si ,  )
Essa rappresenta la probabilità della osservazione parziale dal tempo t+1 fino
al tempo T e di trovarci nello stato nascosto si al tempo t.
Procedendo induttivamente:
T (i)  1 1  i  N
N
t (i)   t 1 ( j )aijb j (ot 1 ) t  T , T  1,...,1
Inizializzazione (arbitraria)
1  i  N Induzione
j 1
N
p (O |  )   1 (i )
Terminazione
i 1
55
Bioinformatica
Hidden Markov Models
EVALUATION: BACKWARD ALGORITHM
E’ analogo al caso forward. L’unica differenza
è che in questo caso andiamo a
ritroso.
Complessità:
O(N2T) Tempo
O(NT) Spazio
N
t (i)   t 1 ( j )aijb j (ot 1 ) t  T , T  1,...,1
1 i  N
j 1
56
Bioinformatica
Hidden Markov Models
DECODING: FORWARD-BACKWARD
Se formuliamo il problema di Decoding come scegliere ad ogni passo lo stato qt
che è individualmente il più probabile, possiamo usare gli algoritmi ForwardBackward per risolvere il problema. Definiamo la variabile:
 t (i)  p(qt  si | O,  )
Cioè la probabilità di trovarsi nello stato si al tempo t data l’osservazione O.
Tale espressione può esprimersi in termini di variabili forward-backward
 t (i ) 
 t (i )  t (i )
 (i )  t (i )
 N t
P (O |  )
  t (i )  t (i )
Normalization factor
i 1
Questo perchè t tiene conto dell’osservazione parziale O1O2…Ot e di
trovarsi nello stato si al tempo t, mentre t tiene conto della rimanente
osservazione Ot+1Ot+2…OT trovandosi nello stato si al tempo t.
57
Bioinformatica
Hidden Markov Models
DECODING: FORWARD-BACKWARD
Possiamo allora trovare lo stato individualmente più probabile al tempo t
con (e quindi tutti gli stati al variare di t):
qt 
a rg m a x  t ( i ) 
1 t T
1 i  N
Tale formula non considera la probabilità di una sequenza di stati ma solo quella
dello stato più probabile ad ogni istante.
Se al tempo t teniamo conto della probabilità della sequenza di stati ai passi
precedenti e dello stato con più alta probabilità al passo corrente avremo un
altro modo di risolvere il problema di decoding (avendo definito un nuovo
criterio di “ottimalità”).
58
Bioinformatica
Hidden Markov Models
DECODING: FORWARD-BACKWARD - Posterior Probability
Se consideriamo le probabilità degli stati nascosti ottenuti da
qt 
a rg m a x  t ( i ) 
1 t T
1 i  N
Questi ci daranno una stima di quanto buona è la predizione nel modello λ.
59
Bioinformatica
Hidden Markov Models
DECODING: VITERBI
La soluzione del problema Evaluation ci permette di avere la somma di tutti i possibili
cammini tra stati nell’HMM.
Vogliamo trovare tra tutti i possibili path di stati nascosti quello Q=q1…qT con più alta
probabilità (in base a quanto visto nella slide precedente);
Il procedimento è simile a quello visto per l’algoritmo Forward. Invece di sommare le
probabilità di transizione, calcoliamo quella massima:
N

 t ( j)   t1 (i) aij b j (ot )
i1



t ( j)  max t1 (i) aij b j (ot )
1iN
60
Forward
Viterbi
Bioinformatica
Hidden Markov Models
DECODING: VITERBI
A differenza dell’algoritmo Forward, nel passaggio dal tempo t-I al tempo t, invece di
sommare le probabilità di transizione, prendiamo quella massima;
Se inoltre teniamo traccia dell’indice dello stato migliore ad ogni passo, alla fine
potremo recuperare la sequenza di stati nascosti più probabili;
Definiamo:
 t (i)  max P(q1 , q2 ,..., qt  i, o1 , o2 ,...ot |  )
q
rappresenta la più alta probabilità lungo un path di stati nascosti che termina in qt=Si.
Utilizzeremo invece:
 t (i )
Per mantenere gli indici degli stati nascosti migliori ad ogni passo.
61
Bioinformatica
Hidden Markov Models
DECODING: VITERBI
Inizializzazione
1 (i)   i bi (o1 ) 1  i  N
 1 (i)  0
Ricorsione


t ( j)  max t1 (i) aij b j (ot )
1iN
 t ( j )  arg max ( t 1 (i )aij )
2  t  T ,1  j  N
1i  N
Terminazione
P   max  T (i)
1i  N
q  arg max T (i)

T
62
P* probabilità finale
Stato finale raggiunto
1i  N
Bioinformatica
Hidden Markov Models
DECODING: VITERBI
Inizializzazione
Ottenere la sequenza di stati:
1 (i)   i bi (o1 ) 1  i  N
 1 (i)  0
qt   t 1 (qt1 )
Ricorsione


t ( j)  max t1 (i) aij b j (ot )
1iN
 t ( j )  arg max ( t 1 (i )aij )
t  T  1, T  2,...,1
2  t  T ,1  j  N
1i  N
Terminazione
P   max  T (i)
1i  N
q  arg max T (i)

T
63
P* probabilità finale
Stato finale raggiunto
1i  N
Bioinformatica
Hidden Markov Models
DECODING: VITERBI vs FORWARD
Inizializzazione
1 (i)   i bi (o1 )
1 i  N
Ricorsione


t ( j)  max t1 (i) aij b j (ot )
1iN
Terminazione
P  max  T (i)

1i  N
64

1 (i)   i bi (o1 )
1 i  N
N

 t ( j)   t1 (i) aij b j (ot )
i1

N
P(O |  )    T (i )
i 1
Bioinformatica
Hidden Markov Models
DECODING: VITERBI
b2(V1)
b1(V2)
V1
bk(V3)
V2
b2(Vk)
V3
Vk
Ad ogni passo scegliamo lo stato nascosto di probabilità massima (tenendo
conto anche della probabilità della sequenza di stati ottenuti al passo
precedente).
O(N2T) Tempo O(NT) Spazio
65
Bioinformatica
Hidden Markov Models
LEARNING: Baum-Welch
“Aggiustare” ( , A, B ) in modo da massimizzare la probabilità di una
osservazione, è il problema più difficoltoso.
Non esiste un metodo analitico.
La procedura iterativa di Baum-Welch (o equivalentemente il metodo
Expectation Maximization EM) permette di massimizzare localmente la
probabilità dell’osservazione.
Definiamo:
 t (i, j )  p(qt  si , qt 1  s j | O,  )
Cioè la probabilità di essere nello stato si al tempo t e nello stato sj al tempo
t+1
66
Bioinformatica
Hidden Markov Models
LEARNING : Baum-Welch
Avevamo in precedenza definito:
 t (i)  p(qt  si | O,  )
 t (i )  t (i )
 t (i )  t (i )
 t (i ) 
 N
P (O |  )
  t (i )  t (i )
i 1
Analogamente avremo:
 t (i)aijb j (ot 1 ) t 1 (i)
t (i, j ) 

P(O |  )
 t (i)aijb j (ot 1 ) t 1 (i)
N
N
 (i)a b (o
i 1 j 1
67
t
ij j
t 1
) t 1 (i)
Bioinformatica
Hidden Markov Models
LEARNING : Baum-Welch
Possiamo mettere in relazione:
 t (i)
 t (i, j )
Dato che vale
N
 t (i)    t (i, j )
j 1
Allora:
T 1
  (i)
t 1
t
Numero atteso di transizioni da si
(numero di volte che viene visitato)
T 1
  (i)
t 1
68
t
Numero atteso di transizioni da si a sj
Bioinformatica
Hidden Markov Models
LEARNING : Baum-Welch
Un insieme ragionevole di formule di rivalutazione del modello è

ˆ   1 (i)
la frequenza attesa dello stato si al tempo t= 1
T 1
 aˆij 
  (i, j )
t 1
T 1
t
  (i)
t 1

numero atteso di transizio ni da s i a s j
numero atteso di transizio ni da s i
t
T
 bˆ (k ) 
j
  t ( j)
t 1 s .t .ot  vk
T
  ( j)
t 1
69
numero atteso di volte in cui ci troviamo

nello stato j e osserviamo il simbolo v k
numero atteso di volte in cui ci troviamo nello stato j
t
Bioinformatica
Hidden Markov Models
LEARNING : Baum-Welch
Iterando
  ( A, B,  )
  ( A, B,  )
È stato dimostrato che ad un certo punto si verifica che
 
Oppure
p(O |  )  p(O |  )
L’algoritmo di Baum-Welch è un caso particolare di EM:
E-step: calcolo di  e 
M-step: calcolo di A, B, 
Complessità #Iterazioni*O(N2T)
70
Bioinformatica
Hidden Markov Models
Problemi
Viterbi: Utilizzare la somma dei Log;
Forward-Backward: Riscalare ad ogni step moltiplicando per una costante;
Baum-Welch: Può convergere a massimi locali;
71
Bioinformatica
GenScan -
72
http://genes.mit.edu/GENSCAN.html
Bioinformatica
GenScan -
http://genes.mit.edu/GENSCAN.html

Il tool di gene prediction più utilizzato;

Presenta il miglior compromesso tra Sensibilità e Specificità (sono due
misure di accuratezza);

Largamente utilizzato dal Consorzio Internazionale durante il Progetto
Genoma Umano;

Utilizza come algoritmo di base l’ Hidden Markov Model (generalizzato);
73
Bioinformatica
GenScan -
http://genes.mit.edu/GENSCAN.html
Le coppie di introni/esoni rappresentano i
differenti modi in cui un introne può
interrompere una coding sequence (dopo la
1° base, dopo la 2° o dopo la 3°)
E0
E1
E2
I0
I1
I2
Einit
Esone iniziale e finale
Esngl
5’ UTR
Filamento sense
Eterm
3’ UTR
polyA
P
N
Filamento antisense
74
…………………..
…………………..
Bioinformatica
GenScan -
http://genes.mit.edu/GENSCAN.html
Scelta una caratteristica (es: identificazione esoni). Possiamo definire i
seguenti valori:

1.
2.
3.
4.
TP (true positive) = Numero di esoni predetti, che sono risultati veri esoni.
FP (false positive) = Numero di esoni predetti che sono in realtà dei falsi.
TN (true negative) = Numero di esoni falsi, identificati come tali.
FN (false negative)= Numero di esoni reali, identificati come falsi.
Avremo le seguenti misure:
TP
numero di esoni correttame nte identifica ti

TP  FN
numero totale degli esoni reali
TN
numero di falsi esoni correttame nte identifica ti
Specificit à 

TN  FP
numero totale di esoni identifica ti come falsi
Sensibilit à 
75
Bioinformatica
GenScan -
http://genes.mit.edu/GENSCAN.html
Possiamo calcolare l’accuratezza come il rapporto tra risultati positivi e
l’intera popolazione;

Accuratezza = (TP+TN)= / (TP+FP+FN+TN)
76
Bioinformatica
GenScan -
77
http://genes.mit.edu/GENSCAN.html
Bioinformatica
GenScan -
http://genes.mit.edu/GENSCAN.html
Possibilità di trovare esoni subottimali:
Indica la soglia di score per cui si trova un
esone. Se il valore scende più esoni (magari
meno probabili) verranno dati in output.
Il valore di default sulla pagina web è 1,00, il
che significa che non vengono stampati esono
subottimali.
Per la maggior parte delle applicazioni, un
valore di cutoff di circa 0,10 è raccomandato.
L'impostazione del valore più basso di 0,10
porterà spesso ad una esplosione del numero
di esoni subottimali, la maggior parte dei quali
probabilmente non sarà utile. D'altra parte, se
il valore è impostato molto superiore a 0,10,
gli esoni subottimali potenzialmente
interessanti potrebbero essere persi.
78
Bioinformatica
GenScan -
http://genes.mit.edu/GENSCAN.html
Probabilità
Inizio, Fine e
Score
del
sito
Score della
che Score
Filamento
suldell’
Frame del
lunghezza
Numerazione
Accettore
dieelemento
codinge sequence
l’elemento
quale
viene
fatta
primo codone
elemento
delTipo
Gene
dei
complessivo
Donatore
di
splicing
calcolata
sia un esone
la predizione
dell’elemento
calcolato
suoi riconosciuto
elementi
dell’esone
Proteina predetta sulla base della CDS calcolata
Type:
Init = Initial exon (ATG to 5' splice site)
Intr = Internal exon (3' splice site to 5' splice site)
Term = Terminal exon (3' splice site to stop codon)
Sngl = Single-exon gene (ATG to stop)
Prom = Promoter (TATA box / initation site) PlyA
79
Bioinformatica
GenScan 
http://genes.mit.edu/GENSCAN.html
ESERCIZIO I

Data la seguente sequenza sconosciuta:
>SCONOSCIUTA
ACATTTGCTTCTGACACAACTGTGTTCACTAGCAACCTCAAACAGACACCATGGTGCATCTGACTCCTGA
GGAGAAGTCTGCCGTTACTGCCCTGTGGGGCAAGGTGAACGTGGATGAAGTTGGTGGTGAGGCCCTGGGC
AGGCTGCTGGTGGTCTACCCTTGGACCCAGAGGTTCTTTGAGTCCTTTGGGGATCTGTCCACTCCTGATG
CTGTTATGGGCAACCCTAAGGTGAAGGCTCATGGCAAGAAAGTGCTCGGTGCCTTTAGTGATGGCCTGGC
TCACCTGGACAACCTCAAGGGCACCTTTGCCACACTGAGTGAGCTGCACTGTGACAAGCTGCACGTGGAT
CCTGAGAACTTCAGGCTCCTGGGCAACGTGCTGGTCTGTGTGCTGGCCCATCACTTTGGCAAAGAATTCA
CCCCACCAGTGCAGGCTGCCTATCAGAAAGTGGTGGCTGGTGTGGCTAATGCCCTGGCCCACAAGTATCA
CTAAGCTCGCTTTCTTGCTGTCCAATTTCTATTAAAGGTTCCTTTGTTCCCTAAGTCCAACTACTAAACT
GGGGGATATTATGAAGGGCCTTGAGCATCTGGATTCTGCCTAATAAAAAACATTTATTTTCATTGC
1.
Facendo variare il parametro “Suboptimal exon cut” verificare se GenScan riesce a
trovare esoni;
2.
Utilizzare la proteina predetta da Genscan per fare un BLAST proteico (BLASTP)
per vedere a cosa corrisponde la predizione fatta da Genscan;
80
Bioinformatica
GenScan 
http://genes.mit.edu/GENSCAN.html
ESERCIZIO II

Prelevare da NCBI la sequenza di mRNA relativa a cox4i1 nel topo.

Usare GenScan per predire la struttura del gene.
Quanti esono vengono prodotto con i parametri standard?
Corrispondono a quelli reali?
Utilizzare la sequenza aminocidica predetta per fare un BLAST. Dai risultati ottenuti si
può dire che la proteina è stata predetta da GenScan in modo corretto?



81
Bioinformatica
EasyBack -
82
http://ferrolab.dmi.unict.it/easyback.html#demo
Bioinformatica
EasyBack 
http://ferrolab.dmi.unict.it/easyback.html#demo
Cosa è a Backtranslation?
Permette di ottenere la sequenza codificante a partire da una proteina;
La traduzione di mRNA in proteina è un processo univoco, la
backtranslation è ambigua;
… A T G
M
G C C
T G G
A C T
A
W
T
S
ACT
TCT
ACC
TCC
ACA
TCA
ACG
……
GCT
ATG
GCC
GCA
GCG
83
TGG
T C A …
Traduzione di mRNA in Proteine
Backtranslation
Bioinformatica
EasyBack 
http://ferrolab.dmi.unict.it/easyback.html#demo
METODI PER LA BACKTRANSLATION
La maggior parte dei tools utilizza la specie specificità come principio base
della backtranslation: Si imita l’uso tipico dei codoni in una
determinata specie.
Ci sono essenzialmente due step:
STEP 1 (training): Costruzione della CODON USAGE TABLE
1.
Prendere una famiglia di proteine la cui sequenza codificante è
nota;
2.
Assegnare per ogni aminoacido il codone più freuqente;
STEP2 (decoding):
1.
Backtranslation usando la codon usage table.

TOOL: BBOCUS, LBT,TIP, BackTranseq (da EMBOSS);
84
Bioinformatica
EasyBack -
http://ferrolab.dmi.unict.it/easyback.html#demo

IL NOSTRO APPROCCIO

EasyBack non è basato sull’imitazione dell’uso dei codoni in un organismo, ma
piuttosto sulla similarità delle sequenze (nelle varie specie).

Data la proteina di input, viene costruito un dataset di proteine (e relative CDS)
eseguendo un BLAST su NCBI;

Il training set sarà il più piccolo possibile affinchè l’HMM possa fare una predizione;

HMM:
 Gli stati nascosti saranno I codoni, le le probabilità di transizione sono le
probabilità tra codoni contigui nel training set (probabilità di passare da un
codone ad un altro):
 I simboli osservabili sono gli aminoacidi e la probabilità di emissione rappresenta
la probabilità che un determinato aminoacido sia stato generato da uno specifico
codone nel trainign set
 E’ possibile stimare la qualità dell’output usando le posterior probabilities.
85
Bioinformatica
EasyBack -
http://ferrolab.dmi.unict.it/easyback.html#demo
Probabilità che il codone i segua il codone j
nelle sequenza del training set
Transition Probabilities
GCT
GCC
…
Emission Probabilities
AAT
A
GCT
GCT
GCC
GCC
…
…
AAT
AAT
M
…
S
Start Probabilities
86
…
…
ATG
…
…
0
0
1
0
0
Bioinformatica
EasyBack -
http://ferrolab.dmi.unict.it/easyback.html#demo
Amino acid
input sequence
Run BLAST
N sequences (cDNA and
corresponding peptides)
Reduce the training
set to N/2 sequences
Buld HMM
Yes
Success?
HMM
Try to make
a prediction
No
repeat
Enlarge the training set
to 3/2N sequences
87
Bioinformatica
EasyBack -
http://ferrolab.dmi.unict.it/easyback.html#demo
EasyBackPanels
Desktop
Log
Area
88
Bioinformatica
EasyBack -
89
http://ferrolab.dmi.unict.it/easyback.html#demo
Bioinformatica
EasyBack -
90
http://ferrolab.dmi.unict.it/easyback.html#demo
Bioinformatica
EasyBack -
91
http://ferrolab.dmi.unict.it/easyback.html#demo
Bioinformatica
EasyBack -
92
http://ferrolab.dmi.unict.it/easyback.html#demo
Bioinformatica
EasyBack 
http://ferrolab.dmi.unict.it/easyback.html#demo
ESERCIZIO I

Data la proteina
>gi|4504347|ref|NP_000549.1| alpha 1 globin [Homo sapiens]
MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFD
LSHGSAQVKGHGKKVADALTNA VAHVDDMPNALSALSDLHAHKLRVD
PVNFKLLSHCLLVTLAAHLPAEFTPAVHASLDKFLASVSTVLTSKYR
Effettuare una backtranslation e confrontare il risultato ottenuto con la sua CDS.
Come recuperare la CDS?
93
Bioinformatica