Biologia computazionale

Università degli studi di milano Docente: Giorgio Valentini Istruttore: Matteo Re C.d.l. Biotecnologie Industriali e Ambientali Biologia computazionale A.A. 2010-2011 semestre II p8 HiddenMarkovModels

Programmazionedinamica in PERL (2) • Implementazione di un algoritmocheutilizzitecniche di programmazionedinamicaed un modellogenerativo • Hash di hash • Hash di array • Biologiacomputazionale • Decoding diunasequenzadiosservazionidato un HMM • ImplementazionealgoritmoVITERBI • Addestramentodi un HMM a partiredaunacollezionedisequenzebiologiche a statinoti Obiettivi

Il livello di complessità di questaesercitazione è alto • Cercate di risolvereilproblemadopoaverlosuddiviso in sottoproblemi • Indipendentemente dal fattoche lo script Perl funzioni o menol’esercizioNON verràvalutatose, insiemeallo script, non verràinviatoanche lo pseudocodice. • Modalità di svolgimentodell’esercitazione: • Scaricare dal sito web del corsoil file HMM_exampleVIT_ls.pl • Questo script è incompleto • La posizionedelleparti da aggiungere è evidenziata da questocommento: ########### description # FILL IN ################## description: descrizionedell’operazione da svolgere • Unadifferenzaimportantetra le esercitazioniprecedenti e questa è che, dopo aver resoutilizzabile lo script dell’algoritmodiViterbi, dovretescriveresenzanessun template, degli script Perl che vi permettanodistimareiparametridell’HMMdainserire in HMM_example_VIT_ls.pl per effettuareil decoding diunasequenzaproteica. I file contenenti le collezionidisequenzesonoreperibilisulsito del corso. Linee guida

3 1 5 1 1 4 6 6 6 6 6 2 4 (osservazioni) + Input: HMM DECODING (VITERBI): 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.

VITERBI: algoritmo Input: 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* = Ptri (i) E’ di fondamentale importanza progettare bene la parte dell’algoritmo che realizza l’iterazione ! Lo stato finale è determinato dallo stato a valore massimo dell’ultimo step di iterazione.

Acquisizione sequenza di osservazioni • Definizione parametri del modello HMM • Inizializzazione Viterbi • Iterazione Viterbi • Determinazione ultimo stato • Traceback (costruzione sequenza di stati) • Stampa output: sequenza stati Realizzazionealgoritmo VITERBI

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 ; Realizzazionealgoritmo VITERBI

Inizializzazione Viterbi # Initialization: my %v = ( 'F' => [0], 'L' => [-0.51] ); Obs. F 0 St. 3 1 5 1 1 6 6 6 FF 0.64 LF -1.61 FL -2.30 LL 0.59 L -0.51 Realizzazionealgoritmo VITERBI V_F 0 Assumiamo uguali le probabilità Iniziali di transire in stato F o L. F = E_f(0) + 0 + 0 = 0 L = F_l(0) + 0 + 0 = -0.51 V_L -0.51 Questo passaggio corrisponde alla prima colonna nella Matrice di progr. dinamica. p_F p_L

Definizione della procedura di iterazione prob. diessere in stato kalla posizionei Per ogni statok, e per una posizione fissata i nella seq., Vk(i) =max{1… i-1} P[x1…xi-1, 1, …, i-1, xi, i = k] Pos. i Emissione osservazione in pos. i Massimizzazione probabilità seq. di stati da pos. 1 a pos. i-1 Osservazioni da posizione 1 a posizione i-1 Seq. Stati da posizione 1 a posizione i-1 Realizzazionealgoritmo VITERBI ITERAZIONE Assumendo di essere in stato k Pos. i+1 Vl(i+1)= el(xi+1)maxkaklVk(i) emissione Viterbi stato k posizione i transizione

Emissione_F(5) = 0 • Definizione della procedura di iterazione transizione 1 3 5 FF 0.64 F 0.64 F F -2.12 L L L max Realizzazionealgoritmo VITERBI Vl(i+1)= el(xi+1)maxkaklVk(i) emissione Viterbi stato k posizione i transizione VF(i+1) = 0 + 0.64 + 0.64 = 1.28

Emissione_F(5) = 0 • Definizione della procedura di iterazione transizione 1 3 5 FF 0.64 F 0.64 F F -2.12 L L L HMM 1° ordine max Realizzazionealgoritmo VITERBI Vl(i+1)= el(xi+1)maxkaklVk(i) emissione Viterbi stato k posizione i transizione Questo calcolo va eseguito, per ogni coppia di simboli (i, i-1), per ogni possibile transizione.

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'); foreach my $state1 qw(F L) { 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 Vl(i+1)= el(xi+1)maxk aklVk(i)

foreach my $state2 qw(F L) { my ($max_v, $max_prev_state, $max_path) = (-1,'X','XX'); foreach my $state1 qw(F L) { my $tmp_trans = $state1.$state2; ####### FILL IN # # 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; } symbol (seq. cycle) idx state1 state2 VITERBI: iterazione F 0.64 F F NB: Salvataggio variabili VITERBI per lo step successivo ? (if) -2.12 L Scelta miglior percorso “locale”

Lunghezza seq. : L 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'); foreach my $state1 qw(F L) { 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; } N° stati : K N° stati : K VITERBI: iterazione Complessità temporale = K * K * L = K2L

## 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 colonna Valore > determina stato finale VITERBI: determinaz. stato finale Stampa stato finale “vincente” e relativo valore VITERBI

## 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; } TRANSIZIONE: Qui $prev_state di L diventa F !!! 6 0.64 1.92 2.56 3.20 4.48 5.12 1.28 3.84 VITERBI: Traceback F 0 -2.12 -1.96 -1.88 -1.80 0.39 2.08 -2.12 -1.72 St. 3 1 5 1 1 6 6 6 -2.81 -1.53 -0.89 -0.25 2.64 3.28 -2.81 2.00 L -0.51 -0.43 -0.27 -0.19 -0.11 3.69 5.38 -0.43 1.58

$last_state (ultima col.) ## 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; } Si parte da $last_state dell’ultima colonna Ma questo viene ripetutamente sostituito dal valore del puntatore della colonna precedente. Array stati prodotto da Traceback (ordine inverso) necessario reverse … 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; } ## Print most probable path ## print "\n[Output]\n"; print join("",@obs_list),"\n"; print join("",reverse @state_seq),"\n"; VITERBI: Traceback OUTPUT VITERBI: sequenza di stati che ha prodotto l’intera serie di osservazioni con probabilità maggiore.

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 arrays da 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) VITERBI: Output aggiuntivo

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

Biologia computazionale