Bioinformatica Catene di Markov, HMM – GenScan – EasyBack Dr. Giuseppe Pigola –...

BioinformaticaCatene di Markov, HMM – GenScan – EasyBack

Dr. Giuseppe Pigola – [email protected]

Catene di Markov Una catena di Markov del primo ordine è una tripla ,

dove: S={s1,s2,…,sN} è un insieme finito di stati (eventi); è la probabilità iniziale degli stati (è

rappresentato da un vettore); A è un insieme di probabilità di transizione tale che:

Dove con qt indichiamo lo stato al tempo t (discreto);

A è rappresentata da una matrice. Valgono le seguenti:

Bioinformatica2

1

0

1

N

jij

ij

a

a

),,( AS

Nisqp ii 1 )( 1

Njisqsqpa itjtij ,1 )|( 1

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;

Bioinformatica3

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

Bioinformatica4

)|(),...,|( 111 tttt qqpqqqp

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

Bioinformatica5

Catene di Markov Catena di Markov per stringhe di DNA:

Bioinformatica7

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ù);

DEFINIZIONEUna catena di Markov si dice Finita se è formata da un insieme finito di stati;

DEFINIZIONEUna catena di Markov si dice Aperiodica se gli stati non vengono osservati mai in modo periodico;

DEFINIZIONEUna catena di Markov si dice Irriducibile se tutti gli stati prima o poi vengono raggiunti (non ci sono zeri nella matrice di transizione);

Bioinformatica8

Catene di Markov ESEMPIO: CATENA DI MARKOV PERIODICA

Bioinformatica9

Catene di Markov ESEMPIO

Bioinformatica10

Catene di Markov ESEMPIO

Bioinformatica11

Catene di Markov ESEMPIO: Considerazioni sulla matrice di transizioni.

Moltiplicando la matrice per se stessa ci riferiamo all’istante t+2.

Bioinformatica12

Catene di Markov ESEMPIO: Se una persona compra attualmente Coca, quale

sarà la probabilità che tra tre volte comprerà Pepsi?

Bioinformatica13

Catene di Markov ESEMPIO: Se una persona compra attualmente Coca, quale

sarà la probabilità che tra tre volte comprerà Pepsi?

Bioinformatica14

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

Bioinformatica15

Catene di Markov AUTOVALORI

AUTOVETTORI

Bioinformatica16

Catene di Markov Se vale

La matrice si può scrivere come (operazione Spettrale):

E quindi

Bioinformatica17

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;

Bioinformatica18

Catene di Markov Ci serve

Dobbiamo considerare la colonna relativa alla Coca e moltiplicare per la distribuzione iniziale

Bioinformatica19

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);

Bioinformatica20

Catene di Markov – Distribuzione Stazionaria D0 * P = D1

D0 * P2 = D2

D0 * Pn = Dn Distribuzione di probabilità degli stati dopo n istanti;

Una distribuzione si dice stazionaria se

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è

Bioinformatica21

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

Bioinformatica22

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.

Bioinformatica23

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:

Bioinformatica24

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;

Bioinformatica25

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

Si dimostra che la matrice di probabilità P è la matrice di transizione di una catena di Markov che ha per distribuzione stazionaria quella data.

Bioinformatica26

P

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:

se i vettori i e j differiscono su più componenti;

Se i e j differiscono su al più una componente h:

Bioinformatica27

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.

Bioinformatica28

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

Bioinformatica29


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);

Bioinformatica30


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;

Bioinformatica31


Bioinformatica32


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

Bioinformatica33


L’obbiettivo ad ogni step è quello di fare una sostituzione di segmenti che migliori la qualità dell’allineamento;

Sia il segmento che deve rimpiazzare quello cancellato;

Definiamo:

Bioinformatica34

P prodotto delle probabilità di background

Q probabilità stimata dei caratteri di x nell’insieme ridotto di n-1 sequenze

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;

Bioinformatica35


Come scegliamo il segmento x tra le possibilità???

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;

Bioinformatica36

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);

Bioinformatica37

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.

Bioinformatica38

Hidden Markov Models Un HMM è definito da una quintupla , dove:

S={s1,s2,…,sN} è un insieme finito di stati nascosti;

V={v1,v2,…,vM} è un insieme di stati osservabili; è la probabilità iniziale degli stati (è

rappresentato da un vettore); A è un insieme di probabilità di transizione tale che

Dove con qt indichiamo lo stato al tempo t (discreto); B rappresenta le probabilità di emissione tale che:

Anche B è rappresentata da una matrice.

Bioinformatica39

MkNjsqvpkb jtkj 1 1 )|at t ()(

),,,,( BAVS

Njisqsqpa itjtij ,1 )|( 1

Nisqp ii 1 )( 1

Hidden Markov Models

Bioinformatica40

Esempio:

Data la sequenze di stati osservabili, quale è la sequenza di stati nascosti più probabile che l’ha generata?

Catene di Markov ESEMPIO: Cerchiamo un modello per discriminare due ipotetiche regioni

di DNA (Sopponiamo per semplicità di avere solo due nucleotidi).

Bioinformatica41

I II

ATTATTATAAATTAAT

…TTAA

TATAATTAATATATTT

…ATAT

Probabilità di transizione dalla regione I alla II con la sequenza TT

Catene di Markov A quale regione appartiene la sequenza TTAT ?

Calcoliamo la probabilità di tutte le possibili sequenze di nucleotidi appartenenti alle due regioni.

Bioinformatica42

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

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à;

Bioinformatica43


Bioinformatica44

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| λ)?


Bioinformatica45

ESEMPIO: Il casinò disonesto.

In un casinò ci sono due dadi di cui uno truccato:

Dado non truccatoP(1) = P(2) = P(3) = P(4) = P(5) = P(6) = 1/6

Dado truccatoP(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.


Bioinformatica46


Quanto è probabile la seguente osservazione nel nostro modello?

1245526462146146136136661664661636616366163616515615115146123562344

Si tratta di un problema di Evaluation.


Bioinformatica47


Quale porzione di questa sequenza è stata prodotta dal dado truccato e quale dal dado non truccato?

1245526462146146136136661664661636616366163616515615115146123562344

Si tratta di un problema di Decoding.


Bioinformatica48


Quanto “truccato” era il dado truccato? Quanto spesso i due dadi venivano cambiati?

1245526462146146136136661664661636616366163616515615115146123562344

Si tratta di un problema di Learning.


Bioinformatica49

ESEMPIO: Il casinò disonesto. Modello di Markov Nascosto

FAIR LOADED

0.05

0.05

0.950.95

P(1|F) = 1/6P(2|F) = 1/6P(3|F) = 1/6P(4|F) = 1/6P(5|F) = 1/6P(6|F) = 1/6

P(1|L) = 1/10P(2|L) = 1/10P(3|L) = 1/10P(4|L) = 1/10P(5|L) = 1/10P(6|L) = 1/2

Transizioni

Emissioni


Bioinformatica50

EVALUATIONVogliamo 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:

Quindi

La probabilità della sequenza degli stati nascosti è data da:

La probabilità che si verifichino contemporaneamente O e Q sarà il prodotto

La probabilità di O nel modello sarà allora ottenuta sommando queste probabilità su tutti i possibili stati nascosti.

T

ttt qopQOp

1

),|(),|(

)()()(),|( 21 21 Tqqq obobobQOpT

TT qqqqqqq aaaQp132211

)|(

)()()()|(),|()|,(13222111 21 Tqqqqqqqqqq obaaobaobQpQOpQOp

TTT

Tqqq

QpQOpOp,...,, 11

)|(),|()|(


Bioinformatica51

EVALUATIONInterpretazione della formula:

Inizialmente al tempo t=1 ci troviamo nello stato q1con probabilità e generiamo il simbolo o1con probabilità

Al tempo t=2 avremo una transizione allo stato q2 con probabilità e verrà generato il simbolo osservabile o2 con probabilità

Il processo continua fino ad arrivare al tempo T

Complessità: 2TNT calcoli

T

TTT

T qqqTqqqqqqqqqq

qqq

obaaobaobQpQOpOp,...,,

21,...,, 11

13222111

11

)()()()|(),|()|(

1q

)( 11obq

21qqa

)( 22obq


Bioinformatica52

EVALUATION: FORWARD ALGORITHMConsideriamo la variabile forward:

Essa rappresenta la probabilità di aver osservato al tempo t o1o2…ot e di trovarci nello stato nascosto si.Procedendo induttivamente:

),|...()( 21 ittt sqooopi

Niobi ii 1 )()( 11

N

iT iOP

1

)()|(

Inizializzazione

Induzione

Terminazione

NjTtobaij tj

N

iijtt

1 2 )(])([)(

11


Bioinformatica53

EVALUATION: FORWARD ALGORITHMLo stato sj può essere raggiunto al tempot+1 dagli N possibili stati al tempo t.

rappresenta la probabilità diaver osservato la sequenza o1o2…ot-1 e di trovarsi a tempo t-1 nello stato nascosto si.

Il prodotto rappresenta laprobabilità di osservare o1o2…ot-1 e di raggiungere lo stato sj al tempo tproveniendo 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

NjTtobaij tj

N

iijtt

1 2 )(])([)(

11

)(1 it

ijt ai)(1


Bioinformatica54

EVALUATION: FORWARD ALGORITHMLo step finale ci darà la probabilità cercata

Complessità:

O(N2T) TempoO(NT) Spazio

N

iT iOP

1

)()|(


Bioinformatica55

EVALUATION: BACKWARD ALGORITHMIn modo speculare consideriamo la variabile backward:

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:

),|,...()( 21 itTttt sqooopi

Inizializzazione (arbitraria)

Induzione

Terminazione

NiiT 1 1)(

NiTTtobajiN

jtjijtt

1 1,...,1, )()()(

111

N

i

iOp1

1 )()|(


Bioinformatica56

EVALUATION: BACKWARD ALGORITHME’ analogo al caso forward. L’unica differenzaè che in questo caso andiamo a ritroso.

Complessità:

O(N2T) TempoO(NT) Spazio

NiTTtobajiN

jtjijtt

1 1,...,1, )()()(

111


Bioinformatica57

DECODING: FORWARD-BACKWARDSe formuliamo il problema di Decoding come scegliere ad ogni passo lo stato qt che è individualmente il più probabile, possiamo usare gli algoritmi Forward-Backward per risolvere il problema. Definiamo la variabile:

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

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.

),|()( Osqpi itt

N

itt

ttttt

ii

ii

OP

iii

1

)()(

)()(

)|(

)()()(

Normalization factor


Bioinformatica58

DECODING: FORWARD-BACKWARDPossiamo allora trovare lo stato individualmente più probabile al tempo t con (e quindi tutti gli stati al variare di t):

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à”).

Tttq iNi

t

1 )(m a xa rg1


Bioinformatica59

DECODING: FORWARD-BACKWARD - Posterior ProbabilitySe consideriamo le probabilità degli stati nascosti ottenuti da

Questi ci daranno una stima di quanto buona è la predizione nel modello λ.

Tttq iNi

t

1 )(m a xa rg1


Bioinformatica60

DECODING: VITERBILa 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:

t ( j) t 1(i)aiji1

N

b j (ot )

t ( j) max1iN

t 1(i)aij b j (ot )

Forward

Viterbi


Bioinformatica61

DECODING: VITERBIA 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:

rappresenta la più alta probabilità lungo un path di stati nascosti che termina in qt=Si.

Utilizzeremo invece:

Per mantenere gli indici degli stati nascosti migliori ad ogni passo.

)|,...,,,...,,(max)( 2121 ttq

t oooiqqqPi

)(it


Bioinformatica62

DECODING: VITERBI

t ( j) max1iN

t 1(i)aij b j (ot )

Inizializzazione

Ricorsione

)()( 11 obi ii Ni 10)(1 i

))((maxarg)( 11

ijtNi

t aij

NjTt 1,2

)(max1

iP TNi

)(maxarg1

iq TNi

T

P* probabilità finale

Stato finale raggiunto

Terminazione


Bioinformatica63

DECODING: VITERBI

t ( j) max1iN

t 1(i)aij b j (ot )

Inizializzazione

Ricorsione

)()( 11 obi ii Ni 10)(1 i

))((maxarg)( 11

ijtNi

t aij

NjTt 1,2

)(max1

iP TNi

)(maxarg1

iq TNi

T

P* probabilità finale

Stato finale raggiunto

Terminazione

Ottenere la sequenza di stati:

)( 11

ttt qq 1,...,2,1 TTt


Bioinformatica64

DECODING: VITERBI vs FORWARD

t ( j) t 1(i)aiji1

N

b j (ot )

t ( j) max1iN

t 1(i)aij b j (ot )

Inizializzazione

Ricorsione

)()( 11 obi ii Ni1 Niobi ii 1 )()( 11

N

iT iOP

1

)()|( )(max1

iP TNi

Terminazione


Bioinformatica65

V1 V2 V3 Vk

b2(V1) b1(V2) bk(V3) b2(Vk)

DECODING: VITERBI

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


Bioinformatica66

LEARNING: Baum-Welch“Aggiustare” 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:

Cioè la probabilità di essere nello stato si al tempo t e nello stato sj al tempo t+1

),,( BA

),|,(),( 1 Osqsqpji jtitt


Bioinformatica67

LEARNING : Baum-WelchAvevamo in precedenza definito:

Analogamente avremo:

),|()( Osqpi itt

N

itt

ttttt

ii

ii

OP

iii

1

)()(

)()(

)|(

)()()(

N

i

N

jttjijt

ttjijtttjijtt

iobai

iobai

OP

iobaiji

1 111

1111

)()()(

)()()(

)|(

)()()(),(


Bioinformatica68

LEARNING : Baum-WelchPossiamo mettere in relazione:

Dato che vale

Allora:

),( )( t jiit

N

jt jii

1t ),()(

1

1

( )T

tt

i

1

1

( )T

tt

i

Numero atteso di transizioni da si (numero di volte che viene visitato)

Numero atteso di transizioni da si a sj


Bioinformatica69

LEARNING : Baum-WelchUn insieme ragionevole di formule di rivalutazione del modello è

la frequenza attesa dello stato si al tempo t= 1

i

ji

1

1

1

1

s da ni transiziodi atteso numero

s a s da ni transiziodi atteso numero

)(

),(ˆ

T

tt

T

tt

ij

i

jia

j stato nello troviamoci cuiin voltedi atteso numero

vsimbolo il osserviamo e j stato nello

troviamoci cuiin voltedi atteso numero

)(

)(

)(ˆ k

1

.. 1

T

tt

T

votstt

j

j

j

kb kt

)(ˆ 1 i


Bioinformatica70

LEARNING : Baum-WelchIterando

È stato dimostrato che ad un certo punto si verifica che

Oppure

L’algoritmo di Baum-Welch è un caso particolare di EM:

E-step: calcolo di e

M-step: calcolo di

Complessità #Iterazioni*O(N2T)

),,( BA ),,( BA

)|()|( OpOp

,,BA


Bioinformatica71

Problemi

Viterbi: Utilizzare la somma dei Log;

Forward-Backward: Riscalare ad ogni step moltiplicando per una costante;

Baum-Welch: Può convergere a massimi locali;

GenScan - http://genes.mit.edu/GENSCAN.html

Bioinformatica72


Bioinformatica73

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);


Bioinformatica74

E0 E1 E2

NP polyA

5’ UTR

I0 I1 I2

Esngl

Einit Eterm

Filamento sense

Filamento antisense

3’ UTR

………………….. …………………..

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°)

Esone iniziale e finale


Bioinformatica75

Scelta una caratteristica (es: identificazione esoni). Possiamo definire i seguenti valori:

1. TP (true positive) = Numero di esoni predetti, che sono risultati veri esoni.

2. FP (false positive) = Numero di esoni predetti che sono in realtà dei falsi.

3. TN (true negative) = Numero di esoni falsi, identificati come tali.

4. FN (false negative)= Numero di esoni reali, identificati come falsi.

Avremo le seguenti misure:

falsi come tiidentifica esoni di totalenumero

tiidentifica ntecorrettame esoni falsi di numero

reali esoni degli totalenumero

tiidentifica ntecorrettame esoni di numero

FPTN

TNàSpecificit

FNTP

TPàSensibilit


Bioinformatica76

Possiamo calcolare l’accuratezza come il rapporto tra risultati positivi e l’intera popolazione;

Accuratezza = (TP+TN)= / (TP+FP+FN+TN)Accuratezza = (TP+TN)= / (TP+FP+FN+TN)


Bioinformatica77


Bioinformatica78

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.


Bioinformatica79

Proteina predetta sulla base della CDS calcolata

Numerazione del Gene e dei suoi elementi

Tipo di elemento riconosciuto

Filamento sul quale viene fatta

la predizione

Inizio, Fine e lunghezza dell’

elemento calcolato

Frame del primo codone dell’elemento

Score del sito Accettore e

Donatore di splicing

Score della coding sequence

calcolata

Probabilità che

l’elemento sia un esone

Score complessivo dell’esone

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


Bioinformatica80

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;


Bioinformatica81

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?

EasyBack - http://ferrolab.dmi.unict.it/easyback.html#demo

Bioinformatica82


Bioinformatica83

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;

ATG

GCT

GCC

GCA

GCG

TGG

ACT

ACC

ACA

ACG

TCT

TCC

TCA

……

… A T G G C C T G G A C T T C A …

M A W T S

Traduzione di mRNA in Proteine

Backtranslation


Bioinformatica84

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);


Bioinformatica85

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.


Bioinformatica86

GCT GCC … AAT

GCT

GCC

…

AAT

Transition Probabilities

A M … S

GCT

GCC

…

AAT

Emission Probabilities

… … ATG … …

0 0 1 0 0

Start Probabilities

Probabilità che il codone i segua il codone jnelle sequenza del training set


Bioinformatica87

Amino acid input sequence

Run BLAST

N sequences (cDNA and corresponding peptides) Buld HMM

Try to make a prediction

Success?YesReduce the training

set to N/2 sequences

No

Enlarge the training setto 3/2N sequences

repeat

HMM


Bioinformatica88

EasyBackPanels

Log Area

Desktop


Bioinformatica89


Bioinformatica90


Bioinformatica91


Bioinformatica92


Bioinformatica93

ESERCIZIO I Data la proteina

Effettuare una backtranslation e confrontare il risultato ottenuto con la sua CDS.

Come recuperare la CDS?

>gi|4504347|ref|NP_000549.1| alpha 1 globin [Homo sapiens]MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSHGSAQVKGHGKKVADALTNA VAHVDDMPNALSALSDLHAHKLRVDPVNFKLLSHCLLVTLAAHLPAEFTPAVHASLDKFLASVSTVLTSKYR

Bioinformatica Catene di Markov, HMM – GenScan – EasyBack Dr. Giuseppe Pigola –...

Documents

Transcript of Bioinformatica Catene di Markov, HMM – GenScan – EasyBack Dr. Giuseppe Pigola –...