biologia computazionale n.
Download
Skip this Video
Loading SlideShow in 5 Seconds..
Biologia computazionale PowerPoint Presentation
Download Presentation
Biologia computazionale

Loading in 2 Seconds...

play fullscreen
1 / 22

Biologia computazionale - PowerPoint PPT Presentation


  • 148 Views
  • Uploaded on

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. p 8. Hidden Markov Models. Programmazione dinamica in PERL (2)

loader
I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
capcha
Download Presentation

PowerPoint Slideshow about 'Biologia computazionale' - cachez


An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.


- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript
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

obiettivi

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
linee guida

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
slide4

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.

slide5

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’algo-

ritmo che realizza l’iterazione !

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

realizzazione algoritmo viterbi

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
realizzazione algoritmo viterbi1

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
realizzazione algoritmo viterbi2

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

realizzazione algoritmo viterbi3

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

realizzazione algoritmo viterbi4

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

realizzazione algoritmo viterbi5

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.

realizzazione algoritmo viterbi6

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.

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

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)

viterbi iterazione1

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”

viterbi iterazione2

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

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 colonna

Valore > determina stato finale

VITERBI: determinaz. stato finale

Stampa stato finale “vincente” e relativo

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

}

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

viterbi traceback1

$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
viterbi traceback2

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

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 arrays

da hash di array

VITERBI: Output aggiuntivo
viterbi output aggiuntivo1

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