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

p8

Hidden Markov Models

C.d.l. Biotecnologie Industriali e Ambientali

Obiettivi

• Programmazione dinamica in PERL (2)• Implementazione di un algoritmo che utilizzi tecniche di

programmazione dinamica ed un modello generativo• Hash di hash• Hash di array

• Biologia computazionale• Decoding di una sequenza di osservazioni dato un HMM• Implementazione algoritmo VITERBI• Addestramento di un HMM a partire da una collezione di

sequenze biologiche a stati noti

Linee guida• Il livello di complessità di questa esercitazione è alto

• Cercate di risolvere il problema dopo averlo suddiviso in sottoproblemi• Indipendentemente dal fatto che lo script Perl funzioni o meno l’esercizio

NON verrà valutato se, insieme allo script, non verrà inviato anche lo pseudocodice.

• Modalità di svolgimento dell’esercitazione:• Scaricare dal sito web del corso il file HMM_exampleVIT_ls.pl• Questo script è incompleto• La posizione delle parti da aggiungere è evidenziata da questo

commento:

########### description # FILL IN ##################

description: descrizione dell’operazione da svolgere

• Una differenza importante tra le esercitazioni precedenti e questa è che, dopo aver reso utilizzabile lo script dell’algoritmo di Viterbi, dovrete scrivere senza nessun template, degli script Perl che vi permettano di stimare i parametri dell’HMM da inserire in HMM_example_VIT_ls.pl per effettuare il decoding di una sequenza proteica. I file contenenti le collezioni di sequenze sono reperibili sul sito del corso.

Permette di identificare: La sequenza di stati che, con probabilità massima, ha prodotto l’intera sequenza di osservazioni (ma non ci dice nulla rispetto alla probabilità di essere in un certo stato ad una determinata posizione nella sequenza).

E’ basato su: tecniche di programmazione dinamica.

DECODING (VITERBI):

Input:

3 1 5 1 1 4 6 6 6 6 6 2 4(osservazioni)

HMM

+

VITERBI: algoritmoInput: x = x1……xL + HMM

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)

E’ di fondamentale importanzaprogettare bene la parte dell’algo-ritmo che realizza l’iterazione !

Lo stato finale è determinato dallo stato a valore massimo dell’ultimo step di iterazione.

Realizzazione algoritmo VITERBI

1) Acquisizione sequenza di osservazioni

2) Definizione parametri del modello HMM

3) Inizializzazione Viterbi

4) Iterazione Viterbi

5) Determinazione ultimo stato

6) Traceback (costruzione sequenza di stati)

7) Stampa output: sequenza stati


• Inizializzazione Viterbi

# Initialization:

my %v = ( 'F' => [0], 'L' => [-0.51] );

NB: hash di arrays .

Una variabile dichiarata in questo modo è un normale hash ma i suoi valori sono arrays (come potete indovinare dalle parentesi quadre) . I valori inseriti si trovano in POSIZIONE 0 (primo elemento) degli array. I valori possono essere aggiunti così:

$v{chiave_hash}->[indice_array] = valore ;


• Inizializzazione Viterbi

# Initialization:

my %v = ( 'F' => [0], 'L' => [-0.51] );

St. 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 p_L

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

Questo passaggio corrispondealla prima colonna nella Matrice di progr. dinamica.


• Definizione della procedura di iterazione

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]

Massimizzazione probabilità seq. di stati da pos. 1a pos. 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

prob. di esserein stato k allaposizione i

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

emissione transizione Viterbi stato k posizione i

Pos. i

Pos. i+1

ITERAZIONE





F F F

L L L

13 5

FF0.640.64

-2.12

Emissione_F(5) = 0

transizione

max

VF(i+1) = 0 + 0.64 + 0.64 = 1.28





F F F

L L L

13 5

FF0.640.64

-2.12

Emissione_F(5) = 0

transizione

max

Questo calcolo va eseguito, per ogni coppia di simboli (i, i-1), per ogni possibile transizione.

HMM 1°ordine

VITERBI: iterazione

my $idx = 1;

foreach my $tmp_obs (@obs_list) {

my $max_state = '';

my $max_state_v = 0;

foreach my $state2 qw(F L) {

my ($max_v, $max_prev_state, $max_path) = (-1,'X','XX');


my $tmp_trans = $state1.$state2;

####### FILL IN #### VITERBI ITERATION STEP ######

# my $tmp_v =

if( $tmp_v > $max_v ) {$max_v = $tmp_v;

$max_prev_state = $state1;

$max_path = $state1.$state2;

}

}

$path{$idx}->{$state2} = $max_prev_state;

$v{$state2}->[$idx] = $max_v;

}

$idx += 1;

}


Per ogni simbolo (colonna)

Per ogni transizione (doppio foreach)

VITERBI: iterazione foreach my $state2 qw(F L) {




####### FILL IN #

# my $tmp_v =




}

}



}

F F F

L

0.64

-2.12

Scelta miglior percorso “locale”

state2state1

idxsymbol(seq. cycle)

NB: Salvataggio variabili VITERBIper lo step successivo

? (if)

VITERBI: iterazione

my $idx = 1;

foreach my $tmp_obs (@obs_list) {

my $max_state = '';

my $max_state_v = 0;





####### FILL IN #### VITERBI ITERATION STEP ######

# my $tmp_v =




}

}



}

$idx += 1;

}

Complessità temporale = K * K * L = K2L

Lunghezza seq. : L

N° stati : K

N° stati : K

VITERBI: determinaz. stato finale

## Determine the last state ##

my $last_state = 'X';

if( $v{'F'}->[$idx-1] > $v{'L'}->[$idx-1] ){

$last_state = 'F';

}else{

$last_state = 'L';

}

print "Final state: $last_state (v= $v{$last_state}->[$idx-1])\n";

Controllo valori ultima colonnaValore > determina stato finale

Stampa stato finale “vincente” e relativovalore VITERBI

VITERBI: Traceback ## Traceback ##my @state_seq = ();

for(my $i=$idx-1;$i>0;$i--) {

my $prev_state = $path{$i}->{$last_state};

push(@state_seq,$prev_state);

$last_state = $prev_state;

}

St. 3 1 5 1 1 6 6 6

F

L

0

-0.51

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

TRANSIZIONE:Qui $prev_statedi L diventa F !!!

VITERBI: Traceback ## Traceback ##my @state_seq = ();

for(my $i=$idx-1;$i>0;$i--) {




}

Si parte da $last_state dell’ultima colonna

Ma questo viene ripetutamente sostituito dal valore del puntatore della colonna precedente.

$last_state(ultima col.)

Array stati prodotto da Traceback (ordine inverso)necessario reverse …

VITERBI: Traceback

## Traceback ##my @state_seq = ();

for(my $i=$idx-1;$i>0;$i--) {




}

## Print most probable path ##

print "\n[Output]\n";

print join("",@obs_list),"\n";

print join("",reverse @state_seq),"\n";

OUTPUT VITERBI: sequenza di stati che ha prodottol’intera serie di osservazioni con probabilità maggiore.

VITERBI: Output aggiuntivo

my $jseq = join(" \t", @obs_list);

print "\tB\t$jseq\n";

my $arrayref = $v{"F"};

my @Vf = @{ $arrayref };

my $Fvit=join("\t", @Vf);

print "F\t$Fvit\n";

$arrayref = $v{"L"};

my @Vl = @{ $arrayref };

my $Lvit=join("\t", @Vl);

print "L\t$Lvit\n\n";

STAMPA MATRICE VITERBI

Estrazione arraysda hash di array

VITERBI: Output aggiuntivo

my $rHoH = \%path;

my $curpos=0;

for my $k1 ( sort {$a <=> $b} (keys %$rHoH) ) {

print "Step: $k1 \tEmission: $obs_list[$curpos]\n";

for my $k2 ( keys %{$rHoH->{ $k1 }} ) { print "\tmax path to $k2 : $rHoH->{ $k1 }{ $k2 }\t";

if($k2 ne $rHoH->{ $k1 }{ $k2 }){

print "*\n";

}else{

print "\n";

}

}

$curpos++;

print "\n";

}

STAMPA MATRICE TRACEBACK (path pointers)

Esercizi

• Rendere lo script funzionante (3 pt) (COMMENTARE IN MANIERA DETTAGLIATA)• Scaricare i seguenti files dal sito del corso: soluble_sequences.txt,

transmembrane_sequences.txt e state_sequences. Essi contengono, rispettivamente, sequenze di porzioni solubili di proteine di lievito, sequenze di porzioni transmembrana di proteine di lievito e sequenze annotate (T= transmembrana, S = solubile) di proteine di lievito. Scrivete 2 script Perl. Il primo servirà a leggere I files con le sequenze amminoacidiche e a generare le probabilità di emissione per gli stati S (solubile) e T (transmembrana). Il secondo script leggerà il file delle annotazioni e stimerà le probabilità di transizione tra I due stati e le probabilità iniziali per gli stati S e T. (3 pt)

• Osservate le pobabilità di emissione dei diversi aa negli stati S e T. Trovate corrispondenza con le proprietà chimico-fisiche degli aa che li rendono più o meno adatti a far parte di una regione transmembrana? (commentate la risposta) (5 pt)

• Modificare lo script HMM_exampleVIT_ls.pl in modo da permettervi di predire se gli aa di una sequenza proteica fanno parte di una regione transmembrana. Provate poi a predire la serie di stati che, con maggior probabilità, ha emesso questa sequenza: KKIIFFFFL. (4 pt)

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