Biologia computazionale

1. 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 7 HiddenMarkovModels

2. IDENTIFICAZIONE DI ELEMENTI FUNZIONALI IN UNA SEQUENZA (seq. nucleotidiche o proteiche) Bio PROBLEMA: Abbiamoricevuto un file FASTA contenente le sequenzecromosomicheassemblateprodottenell’ambito di un progetto di sequenziamentogenomico. • Questi file non contengonoaltrainformazione se non la sequenzanucleotidica • Perchè la sequenzasia utile dobbiamoannotarlaidentificandoglielementifunzionali in essacontenutie la loroposizione. Come possiamoannotareilgenomaappenaricevuto?

3. IDENTIFICAZIONE DI ELEMENTI FUNZIONALI IN UNA SEQUENZA Bio IPOTESI DI LAVORO: • E’ disponibileilgenomaannotato di una specie evolutivamentevicina.  allineamento al nuovogenoma di tutte le proteinedella specie giàannotata • Sequenziamentodituttiitrascrittiseguitodaallineamentodellelorosequenze al genoma OSSERVAZIONI: • Non è dettochesiadisponibileilgenomaannotatodiuna specie evolutivamentevicina • Non tuttiglielementifunzionalisonogeni(quindiil “trucco” deitrascritti non è sempreapplicabile) I geni sono comunque tra i primi elementi funzionali da cercare in una seq. genomica non annotata

4. IDENTIFICAZIONE DI ELEMENTI FUNZIONALI IN UNA SEQUENZA Bio IPOTESI DI LAVORO (II): • Supponiamodiaverealcuneannotazioni per ilgenoma (ad esempiodisponiamodelleannotazioni del 20% (stima) deigenidell’organismo) e divoleridentificare le posizionideglialtrigeni. POSSIBILE SOLUZIONE: • Utilizzodeglielementinoti per costruiremodelli • Utilizzodeimodelliprodotti per effettuareunascansionedellasequenzagenomicachepermettadirilevare la presenzadiregionisimili al modelloconsiderato

5. IDENTIFICAZIONE DI ELEMENTI FUNZIONALI IN UNA SEQUENZA Bio DIFFICOLTA’ CHE CARATTERIZZANO QUESTO APPROCCIO: Supponiamo di voler creare dei modelli di elementi funzionali che permettano di predire la presenza (e la struttura) dei geni. L’obiettivo è quello di creare UN MODELLO che permetta di predire tutti I geni presenti in un genoma … MA I GENI DI UN ORGANISMO POSSONO AVERE CARATTERISTICHE MOLTO DIVERSE … E’ DIFFICILE PER UN UNICO MODELLO ADATTARSI A TUTTI I GENI!

6. IDENTIFICAZIONE DI ELEMENTI FUNZIONALI IN UNA SEQUENZA Bio

7. IDENTIFICAZIONE DI ELEMENTI FUNZIONALI IN UNA SEQUENZA Bio

8. DETERMINAZIONE STRUTTURA DI UN GENE Bio • Approccio sperimentale: • Estrazione degli RNA e sequenziamento. • Confronto con dati pubblici: sequenze di cDNA assemblate e sequenze EST. • Limitazioni: Trascritti rari,specificità cellulare e tissutale, espressione condizionale. NON E’ possibile estrarre RNA da tutti i tipi cellulari, a tutti gli stadi di sviluppo e in tutte le possibili condizioni ambientali da un organismo complesso. • Approccio Computazionale: • Predizione ab initio basata su caratteristiche della sequenza in esame. • Approccio ibrido: • Predizione ab initio integrata con dati ottenuti mediante sequenziamento di RNA parziali e confronto con strutture di geni omologhi noti in altri organismi. (similarità a livello di sequenza proteica o nucleotidica)

9. IDENTIFICAZIONE DI ELEMENTI FUNZIONALI IN UNA SEQUENZA Bio Problema pratico: La complessità dei modelli probabilistici è proporzionale alla complessità strutturale degli elementi funzionali che vogliamo modellare. Questo implica che un modello adatto a descrivere le caratteristiche di un gene eucariotico ha già un livello di complessità abbastanza elevato. Per introdurre I modelli basati su catene di Markov e gli Hidden Markov Models (HMM) utilizzeremo esempi più semplici (e quindi associati a elementi genomici aventi una struttura meno complessa di quella dei geni).

10. IDENTIFICAZIONE DI ELEMENTI FUNZIONALI IN UNA SEQUENZA Bio CS CARATTERISTICHE DEL PROBLEMA: • I “modelli” di cui abbiamobisognodovrannoesserebasatiunicamentesullasequenzagenomica. • Hanno alcuniaspetti in comune con imodellicheabbiamovistonella parte del corsodedicataallafilogenesi • Sonomodelliprobabilistici • Sono in gradodiGENERARE SEQUENZE “CONFORMI AL MODELLO” (per questovengonodettimodelligenerativi)

11. ISOLE CpG Bio Nel genoma umano i dinucleotidi CpG sono più rari di quanto saremmo portati ad aspettarci sulla base delle probabilità indipendenti di osservare una C o una G a causa di un motivo biologico. Le coppie CpG vanno incontro ad un processo di metilazione che modifica la citosina: C  T Le regioni promotore sono ricche in CpG in quanto i promotori non vengono metilati e quindi, in essi, la frequenza di mutazione C  T è molto minore che nel resto del genoma. Se siamo in grado, data una sequenza genomica in formato FASTA, di identificare al suo interno una o più isole CpG queste possono essere viste come l’indicazione che, nelle “immediate” vicinanze, è presente un gene.

12. IDENTIFICAZIONE DI ELEMENTI FUNZIONALI IN UNA SEQUENZA Bio Stat CATENE DI MARKOV per la ricerca di isole CpG

13. A T C G IDENTIFICAZIONE DI ELEMENTI FUNZIONALI IN UNA SEQUENZA Bio Stat • Una catena di Markov è una • collezionedistati e transizioni. • Seguendo le transizioninel • diagrammasipossono • generaretutte le sequenze • osservabili. • Questo è ilclassicodiagrammadi catena di Markov per la generazionedisequenzedi DNA.

14. IDENTIFICAZIONE DI ELEMENTI FUNZIONALI IN UNA SEQUENZA Bio Stat • Ad ognifreccianellafigura è associato un valorecherappresenta la probabilitàche ad un certont ne segua un altro. • Aat = P(xL=T | xL-1=A) Probabilità che A sia seguita da T in una sequenza di DNA. A T C G

15. CATENE DI MARKOV : PROBABILITA’ Bio Stat La probabilitàdiunasequenzaxpuòesseredefinita come: P(x) = P(xL,xL-1,…,x1) Applicando P(X,Y) = P(X|Y)P(Y) più volte: = P(xL|xL-1,…,x1)P(xL-1|xL-2,…,x1)…P(x1) REGOLA DI MARKOV: La probabilitàdiosservare un datoeventonondipendadatutta la sequenzadeglieventiprecedenti ma solo daNeventiprecedenti (nelcasopiùsemplice N=1)

16. CATENE DI MARKOV : PROBABILITA’ Bio Stat ES. CATENA DI MARKOV DEL PRIMO ORDINE : L’osservazione in posizione i-esima dipende unicamente da quello che osserviamo in posizione i-1.

17. CATENE DI MARKOV : PROBABILITA’ Bio Stat ES. CATENA DI MARKOV DEL SECONDO ORDINE : L’osservazione in posizione i-esima dipende unicamente da quello che osserviamo in posizione i-1 ed in posizione i-2.

18. CATENE DI MARKOV : PROBABILITA’ Bio Stat P(X) = P(xL|xL-1,…,x1)P(xL-1|xL-2,…,x1)…P(x1) Regoladimarkov: Datoche la probabilitàdiogniosservazionedipende solo dall’osservazioneprecedente, possiamosemplificarein: P(X) = P(xL|xL-1)P(xL-1|xL-2)…P(x2|x1)P(x1) = X = seq. di lunghezza L simbolo osservato in pos. i simbolo osservato in pos. i-1

19. CATENE DI MARKOV : isole CpG Bio Stat Proviamo ad usare2 catenedi Markov per trovareisoleCpG in unasequenzagenomica. Ogni catena di Markov costituisceun modellochedeveesserebasatosullaprobabilitàdiosservareildinucleotideCpG in sequenzeNOTE. Dopo aver preparato due set disequenze (CpG islands e NON_CpG islands) sicreano due tabelledicontingenzariportanti le probabilità (osservate) che un qualsiasint y segua un nucleotide x. Le tabelleottenutesarannoalla base dellecatenedi Markov dautilizzare per ilnostroesperimentodidiscriminazione.

20. CATENE DI MARKOV : isole CpG Bio Stat Il modello + (CpG): si basa sulla matrice di transizioni += (a+st), in cui: a+st = (probabilità che t segua s in un’isola CpG) Xi Xi-1

21. CATENE DI MARKOV : isole CpG Bio Stat Il modello - (noCpG): si basa sulla matrice di transizioni -= (a-st), in cui: a-st = (probabilità che t segua s al di fuori di un’isola CpG) Xi Xi-1

22. CATENE DI MARKOV : isoleCpG (confrontotabelleditransizione) Bio Stat Xi + Xi-1 Xi - Xi-1

23. CATENE DI MARKOV : P cheuna seq. derividaun’isolaCpG (I) Bio Stat Data la sequenza GCG, calcolare la probabilità che essa appartenga ad un’isola CpG o meno. P(G+C+G+) P(G-C-G-) Calcoliamo la tabella riportante il logaritmo dei rapporti di frequenza: +aa = 0.180 log2(0.180/0.300) = -0.740 -aa = 0.300

24. CATENE DI MARKOV : P cheuna seq. derividaun’isolaCpG (II) Bio Stat Le catene di Markov che stiamo utilizzando come modello sono catene del primo ordine e quindi la probabilità di osservare un nucleotide in posizione i dipende UNICAMENTE da ciò che si osserva in posizione i-1. Calcoliamo P(Ci|Gi-1)P(Gi|Ci-1) nelle due catene (Dato che siamo in spazio logaritmico i prodotti diventano somme) 0.461+1.812 = LogRatio(Ci|Gi-1)+LogRatio(Gi|Gi-1) = 2.273 Ora torniamo alle probabilità relative: P+(GCG)/P-(GCG) = 22.273 = 4.461 In definitivaP+(GCG) =4.461P-(GCG) Dato che gli stati sono solo 2 (+/-) e la somma delle probabilità è 1… P(+) + P(-) = 1 quindi 4.461 P(-) + (P-) = 1 P(-)= 17% , P(+)= 83%

25. CATENE DI MARKOV : PROBLEMI PRATICI Bio Stat Il nostro test è basato sul calcolo di un rapporto di probabilità di una sequenza osservata IN STATI NOTI … quando stiamo analizzando una sequenza genomicaNON ANNOTATA … P(G+C+G+) P(G-C-G-) Dove inizia e dove finisce la regione da testare ? NON SAPPIAMOIN CHE STATO CI TROVIAMO !!! ? start ? end genome

26. CATENE DI MARKOV : PROBLEMI PRATICI Bio Stat Questo equivale a dire che quando testiamo una regione genomica come potenziale isola CpG siamo NOI a scegliere la regione da valutare (le sue estremità e la sua estensione). Quello che vogliamo ottenere, invece, è un modello che valuti TUTTA la sequenza del genoma e decida al posto nostro che dalla posizione a alla posizione b è presente (potenzialmente) un’isola CpG. In definitiva, per quanto statisticamente corretta la soluzione appena presentata NON E’ applicabile per la predizione delle isole CpG in una sequenza genomica non annotata!

27. CATENE DI MARKOV : PROBLEMI PRATICI Bio Stat Possibile soluzione: scegliere una dimensione w FISSA (es. 1000 nt) ed utilizzare una finestra scorrevole sul genoma valutando il risultato del test ad ogni passo. Questo approccio è affetto da diversi problemi: Le isole CpG hanno lunghezze diverse e quindi una finestra di estensione fissa potrebbe essere inadatta per esprimere un giudizio. - Se scegliamo una dimensione della finestra w troppo piccola, tenderemo a giudicare ogni occorrenza del dinucleotide CG come proveniente da un’isola CpG. - Se scegliamo una dimensione di w troppo grande non avremo abbastanza potere di discriminazione.

28. CATENE DI MARKOV : PROBLEMI PRATICI Bio Stat Soluzione alternativa : Incorporare i due modelli (le due catene di Markov) in un UNICO modello.

29. CATENE DI MARKOV : PROBLEMI PRATICI Bio Stat Soluzione alternativa : Vantaggi: Questo modello riflette meglio la realtà. Non c’è più dipendenza dalla dimensione della finestra dato che il modello prevede ad ogni step la possibilità di transire da un qualsiasi nucleotide in stato + (isola CpG) ad un nucleotide in stato – (non isola CpG). ATTENZIONE: La maggior flessibilità del modello ha un costo!

30. HIDDEN MARKOV MODELS ( HMM ) Bio Stat Soluzione alternativa : Svantaggi: Il modello contiene non uno ma DUE possibili stati ( + / - ) per ogni possibile nucleotide! Osservando la sequenza non sappiamo più quale stato ha emesso un certo nt. In definitiva gli stati sono NASCOSTI ! Da qui il nome : HIDDEN Markov Models

31. HIDDEN MARKOV MODELS ( HMM ) Bio Stat CARATTERISTICHE degli HMM : EVALUATION: Siamo liberi di definire il modello in modo da adattarlo il più possibile al tipo di elementi funzionali che vogliamo modellare ma dovremo comunque rispondere a questa domanda: Data una sequenza qual è la PROBABILITA’che essa sia stata generata dal modello? DECODING: Data una sequenza decidere, per ogni posizioneQUAL’E’ LO STATO (NASCOSTO) PIU’ PROBABILE. TRAINING: Data una collezione di sequenze annotate come possiamo utilizzarle per addestrare un HMM ?

32. HIDDEN MARKOV MODELS ( HMM ) Bio Stat CARATTERISTICHE degli HMM : Il modello su cui stiamo ragionando ha una caratteristica particolare … come le catene di Markov è composto UNICAMENTE DA PROBABILITA’ DI TRANSIZIONE. In questo modo siamo OBBLIGATI a creare un NUMERO DI STATI PARI a N * dim(alfabeto), ossia Numero_Stati * Numero_Nucleotidi (nel nostro caso 2 * 4 = 8). E’ molto scomodo. Gli Hidden Markov Models, per risolvere questo problema, prevedono l’esistenza di DUE TIPI DI PROBABILITA’ : probabilità di transizione e probabilità di emissione. P_transz = probabilità di transizione tra STATI P_em = probabilità di emettere un dato simbolo quando ci troviamo in un certo stato

33. 1: 1/6 2: 1/6 3: 1/6 4: 1/6 5: 1/6 6: 1/6 1: 1/10 2: 1/10 3: 1/10 4: 1/10 5: 1/10 6: 1/2 HIDDEN MARKOV MODELS ( HMM ) Bio Stat • Consideriamo un casinò in cui viene utilizzato, per la maggior parte del tempo, un dado onesto ma, occasionalmente, si passa ad un dado truccato. • Il dado truccato ha una probabilità pari a 0.5 di fare 6 e 0.1 per gli altri numeri da uno a cinque. Il casinò passa dal dado onesto al dado truccato con una probabilità di 0.05 e torna al dado onesto con una probabilità doppia a quella della transizione inversa (0.1). • Possiamo esprimere questa situazione attraverso il seguente modello: 0.9 0.95 In un dato momento non sappiamo quale dado sta usando il casinò. 0.05 0.1 Fair Loaded

34. HIDDEN MARKOV MODELS ( HMM ) Bio Stat Il casinò disonesto: DECODING • Domanda: • Data unaseriediosservazioni (e un HMM) siamo in gradodideciderequaliderivanodall’utilizzo del dado onesto e qualidal dado truccato? • Seriediosservazioni: • 3 1 5 1 1 6 6 6 6 3 6 2 4

35. HIDDEN MARKOV MODELS ( HMM ) Bio CS Il casinò disonesto: DECODING • DECODING: • Decodificareunasequenzasignificatrovare, tratuttequellepossibili, la seriedistati in gradodigenerare, con probabilitàmassima, le osservazioni . • Un HMM assume che la sequenzaosservatasiastatageneratada un processomarkoviano. A causadellaregoladi Markov questiprocessisonodetti memory-less (si “ricordano” solo di N passiprecedenti a quelloconsiderato). • Il problemachevogliamorisoverepuòesseresuddiviso in un numerofinitodisottoproblemiripetitivi • Puòessererisoltomedianteprogrammazionedinamica

36. HIDDEN MARKOV MODELS ( HMM ) Bio CS Diagramma TRELLIS: Il diagramma trellis permettedirappresentareunaserietemporaleditransizionitra un set finitodistati. Supponiamodiavere2 possibilistatidi un HMM edunaseriediosservazioni. Essesipossonorappresentarecosì: osservazioni Ad ogni step vengono rappresentate tutte le possibili transizioni (4).

37. HIDDEN MARKOV MODELS ( HMM ) Bio CS DECODING definizione del problema: Dato un HMM edunasequenzadiosservazionitrovare la seriedistati (s*) del modello in gradodigenerarecon probabilitàmassima la seriediosservazioni. Data una sequenza osservabilex = (x1,…,xL), La più probabile sequenza si stati s* =(s*1,…,s*L) è quella che massimizzap(s|x). NB: molto spesso la seq. di stati si indica con la lettera greca  ( ad indicare “path” … percorso)

38. HIDDEN MARKOV MODELS ( HMM ) Bio CS DECODING : AlgoritmodiViterbi Un problemapraticodiquestoalgoritmo, comune a tuttiglialgoritmichemoltiplicanoprobabilità, è ilrischiodi underflow. Per evitarequestoproblemadobbiamoportare le probabilitàdell’HMM in spaziologaritmico e poi sostituire le moltiplicazioni con somme. Come possiamoconvertire le probabilità in log(probabilità)?

39. 1: 1/6 2: 1/6 3: 1/6 4: 1/6 5: 1/6 6: 1/6 1: 1/10 2: 1/10 3: 1/10 4: 1/10 5: 1/10 6: 1/2 HIDDEN MARKOV MODELS ( HMM ) Bio CS DECODING : P in spaziologaritmico P transizione: La sommaditutte le probabilitàditransizione in uscitadaunostatodeveesserepari a 1. Nelcaso in cui tutte le possibilitransizionefosseroequiprobabili le loroprobabilitàsarebberotutteuguali a 1/K (k = numerodeglistati) 0.9 0.95 FL FF LL 0.05 0.1 LF Fair Loaded

40. HIDDEN MARKOV MODELS ( HMM ) Bio CS DECODING : P in spaziologaritmico P transizione: La sommaditutte le probabilitàditransizione in uscitadaunostatodeveesserepari a 1. Nelcaso in cui tutte le possibilitransizionefosseroequiprobabili le loroprobabilitàsarebberotutteuguali a 1/K (k = numerodeglistati) Se fossero equiprob. varrebbero ½ … Se fossero equiprob. varrebbero ½ …

41. HIDDEN MARKOV MODELS ( HMM ) Bio CS DECODING : P in spaziologaritmico P transizione (matrice A): A(FF) = ln(0.95 / 0.5) = 0.6418539 A(FL) = ln(0.05 / 0.5) = - 2.302585 A(LF) = ln(0.10 / 0.5) = - 1.609438 A(LL) = ln(0.90 / 0.5) = 0.5877867 somma = 1 somma = 1 A

42. HIDDEN MARKOV MODELS ( HMM ) Bio CS DECODING : P in spaziologaritmico P emissione (FAIR) : matrice EF EF(1) = ln( 1/6 / 1/6 ) = 0 EF(2) = ln( 1/6 / 1/6 ) = 0 EF(3) = ln( 1/6 / 1/6 ) = 0 EF(4) = ln( 1/6 / 1/6 ) = 0 EF(5) = ln( 1/6 / 1/6 ) = 0 EF(6) = ln( 1/6 / 1/6 ) = 0 somma = 1

43. HIDDEN MARKOV MODELS ( HMM ) Bio CS DECODING : P in spaziologaritmico P emissione (LOADED) : matrice EL EL(1) = ln( 1/10 / 1/6 ) = -0.51 EL(2) = ln( 1/10 / 1/6 ) = -0.51 EL(3) = ln( 1/10 / 1/6 ) = -0.51 EL(4) = ln( 1/10 / 1/6 ) = -0.51 EL(5) = ln( 1/10 / 1/6 ) = -0.51 EL(6) = ln( 1/2 / 1/6 ) = 1.10 somma = 1

44. HIDDEN MARKOV MODELS ( HMM ) Bio CS DECODING : P in spaziologaritmico 0.59 0.64 FF FL LL 1: 0 2: 0 3: 0 4: 0 5: 0 6: 0 1: -0.51 2: -0.51 3: -0.51 4: -0.51 5: -0.51 6: 1.10 -2.30 -1.61 LF Fair Loaded

45. HIDDEN MARKOV MODELS ( HMM ) Bio CS • DECODING : AlgoritmodiViterbi • E’ basatosull’ideadirisolvere un problemaintrattabilesuddividendolo in tantisottoproblemitrattabili (progr. dinamica). • Risolveilproblema del decoding: data unasequenzadiosservazionicalcola la seriedistatiche ha generato la sequenzadiosservazioni con probabilitàmassima. • E’ difondamentaleimportanzadefinireisottoproblemi in mododapermetteredirisolverli in manieraiterativa. Mano a manochel’algoritmoprocede I nuovisottoproblemiverrannorisoltisfruttando le soluzionideisottoproblemiprecedentementeincontrati.

46. HMM : DECODING (VITERBI, dyn.prog.) Bio CS Dato che 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] prob. di essere in stato k alla posizione i Emissione osservazione in pos. i Osservazioni da posizione 1 a posizione i-1 Seq. Stati da posizione 1 a posizione i-1 Massimizzazione probabilità seq. di stati da posizione 1 a posizione i-1 Assumendo di essere in stato k Come calcolare Vl(i+1) ?

47. HMM : DECODING (VITERBI, dyn.prog.) Bio CS Dato che 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] Per definizione, Vl(i+1) = max{1… i}P[ x1…xi, 1, …, i, xi+1, i+1 = l ] = max{1… i}P(xi+1, i+1 = l | x1…xi, 1,…, i) P[x1…xi, 1,…, i] = max{1… i}P(xi+1, i+1 = l | i ) P[x1…xi-1, 1, …, i-1, xi, i] = maxk [P(xi+1, i+1 = l | i=k) max{1… i-1}P[x1…xi-1,1,…,i-1, xi,i=k]] = maxk [ P(xi+1 | i+1 = l )P(i+1 = l | i=k)Vk(i)] = el(xi+1)maxkaklVk(i) Come calcolare Vl(i+1) ? emissione Viterbi stato k posizione i transizione

48. HMM : DECODING (VITERBI, dyn.prog.) Bio CS Vk(i) = max{1… i-1} P[x1…xi-1, 1, …, i-1, xi, i = k] Vl(i+1) = max{1… i}P[ x1…xi, 1, …, i, xi+1, i+1 = l ] = el(xi+1)maxkaklVk(i) RICORSIONE emissione Viterbi stato k posizione i transizione

49. HMM : DECODING (VITERBI, dyn.prog.) Bio CS Input: x = x1……xL 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)

50. HMM : DECODING (VITERBI, dyn.prog.) Bio CS DECODING : calcolo Prima didescrivereipassaggidell’algoritmodiViterbidobbiamorisolvere un problema: qual’è la strutturadatichepermettedieffettuareicalcoli in manierapiùefficientepossibile? Abbiamogiàconsiderato un problemacherichiedeval’utilizzodellaprogrammazionedinamica (allineamentosequenze) e, in quelcaso, la strutturadatifondamentale era unamatrice. Anche in questocaso è possibileutilizzareunamatrice?