Biofisica Computazionale Appunti Personali

42
Biofisica Computazionale 1. L’allineamento di una coppia di sequenze: Per verificare se due sequenze (proteiche o di DNA) siano correlate bisogna allinearle tr correlate quando sono più simili tra loro rispetto a quello che mi aspetto per caso. In qu una correlazione evolutiva tra le due proteine, definendole omologhe, ovvero provenienti d Capire, nonostante l’allineamento possa sembrare in entrambe i casi molto simile, se due d proteine omologhe o meno rappresenta la grande sfida della bioinformatica. La sequenza di della proteina in questione proviene da una serie di mutazioni del suo ancestrale. Prender di mutamenti: le sostituzioni tra i vari amminoacidi (ovvero uno è stato sostituito nella altro cambiandone magari in una certa misura la struttura terziaria) . Le inserzioni e le mancanti o aggiuntivi nella catena di amminoacidi) che insieme chiameremo gaps. La logica che c’è dietro ogni allineamento tra sequenze è quella di assegnare per ogni all (score). Per l’assegnazione dello score bisogna quindi scegliere un modello che ci ottimiz analizza il risultato dello score ottenuto e si constata se è effettivamente maggiore o m potrebbe aspettare dall’allineamento di due sequenze casuali (evidentemente non correlate) Lo score sarà la somma di termini per ogni coppia di amminoacidi delle due catene più i te Nell’interpretazione probabilistica questo corrisponderà al logaritmo del rapporto tra la sequenze sono correlate e la probabilità che non lo siano. Possiamo inoltre avere due sequ m oppure di due lunghezze differenti m,n. Pertanto è possibile richiedere sia un allineame sequenze sia un allineamento locale o a blocchi. Il primo è evidentemente applicato alle lunghezza (ovvero con poche inserzioni e delezioni) mentre il secondo è utilizzato solitam lunghezza differente. (N.B. per sequenze di proteine dovremmo lavorare con un alfabeto di corrispondono ai 20 amminoacidi possibili per ogni posizione, invece per sequenze di DNA u 4 caratteri, A C G T, che corrispondono ai quattro nucleotidi). Come possiamo sapere i termini dovuti alle possibili sostituzioni o conservazioni per ogni due catene allineate? La risposta si trova in quelle che chiamiamo matrici di sostituzione L’idea della matrice di sostituzione sta nel formalizzare matematicamente quanto abbiamo detto sul punteggio compl all’allineamento di due sequenze: Il logaritmo del rapporto delle probabilità del modello modello R (ovvero Mach, probabilità che le sequenze siano correlate e Random, probabilità

Transcript of Biofisica Computazionale Appunti Personali

Biofisica Computazionale

1. Lallineamento di una coppia di sequenze:Per verificare se due sequenze (proteiche o di DNA) siano correlate bisogna allinearle tra loro. Due sequenze sono correlate quando sono pi simili tra loro rispetto a quello che mi aspetto per caso. In questo senso possiamo stabilire una correlazione evolutiva tra le due proteine, definendole omologhe, ovvero provenienti da un progenitore comune. Capire, nonostante lallineamento possa sembrare in entrambe i casi molto simile, se due diversi allineamenti sono di proteine omologhe o meno rappresenta la grande sfida della bioinformatica. La sequenza di amminoacidi analizzati della proteina in questione proviene da una serie di mutazioni del suo ancestrale. Prenderemo in considerazione 3 tipi di mutamenti: le sostituzioni tra i vari amminoacidi (ovvero uno stato sostituito nella medesima posizione da un altro cambiandone magari in una certa misura la struttura terziaria) . Le inserzioni e le delezioni (ovvero pezzi mancanti o aggiuntivi nella catena di amminoacidi) che insieme chiameremo gaps. La logica che c dietro ogni allineamento tra sequenze quella di assegnare per ogni allineamento un punteggio (score). Per lassegnazione dello score bisogna quindi scegliere un modello che ci ottimizza lallineamento. In ultimo si analizza il risultato dello score ottenuto e si constata se effettivamente maggiore o minore di quello che ci si potrebbe aspettare dallallineamento di due sequenze casuali (evidentemente non correlate). Lo score sar la somma di termini per ogni coppia di amminoacidi delle due catene pi i termini di ogni gaps. Nellinterpretazione probabilistica questo corrisponder al logaritmo del rapporto tra la probabilit che le due sequenze sono correlate e la probabilit che non lo siano. Possiamo inoltre avere due sequenze della stessa lunghezza m oppure di due lunghezze differenti m,n. Pertanto possibile richiedere sia un allineamento globale delle due sequenze sia un allineamento locale o a blocchi. Il primo evidentemente applicato alle sequenze della stessa lunghezza (ovvero con poche inserzioni e delezioni) mentre il secondo utilizzato solitamente per sequenze di lunghezza differente. (N.B. per sequenze di proteine dovremmo lavorare con un alfabeto di 20 caratteri, che corrispondono ai 20 amminoacidi possibili per ogni posizione, invece per sequenze di DNA utilizzeremo un alfabeto di 4 caratteri, A C G T, che corrispondono ai quattro nucleotidi). Come possiamo sapere i termini dovuti alle possibili sostituzioni o conservazioni per ogni coppia di amminoacidi delle due catene allineate? La risposta si trova in quelle che chiamiamo matrici di sostituzione. Lidea della matrice di sostituzione sta nel formalizzare matematicamente quanto abbiamo detto sul punteggio complessivo che diamo allallineamento di due sequenze: Il logaritmo del rapporto delle probabilit del modello M con la probabilit del modello R (ovvero Mach, probabilit che le sequenze siano correlate e Random, probabilit che non lo siano)

Ad ogni a e b fissati (coppia di amminoacidi) otteniamo un elemento di matrice, della matrice S di sostituzione che per ogni possibile scambio tra due amminoacidi ci dice quanto questo possa essere frequente (ovvero ci fornisce la probabilit che questo avvenga). Le matrici di sostituzione pi utilizzate sono BLOSTUM e PAM. Entrambe le matrici vengono calcolate sullassunzione che analizzando sequenze correlate filogeneticamente si pu calcolare la probabilit con cui ogni amminoacido subisce un evento di sostituzione. Matrici BLOSTUM: queste matrici sono derivate da blocchi di sequenze appartenenti ad una stessa famiglia con la condizione che debbano essere considerate equivalenti due sequenze, quando la percentuale dei loro residui identici (ovvero conservati) sia superiore di un certo valore percentuale stabilito. Matrici PAM: Per creare le matrici PAM furono calcolate le probabilit di sostituzione per ogni amminoacido con ciascun altro, prevedendo in media una sostituzione ogni 100 amminoacidi. Questa viene chiamata distanza evolutiva 1 PAM. Proprio dallapplicazione delle probabilit di sostituzione degli amminoacidi tra le sequenze correlate definite in PAM1 possibile calcolare matrici PAM di ordine superiore. Tipicamente si utilizzano matrici comprese tra PAM10 (per sequenze filogeneticamente vicine) e PAM240 (per sequenze filogeneticamente lontane). La penalizzazione del punteggio dovuta alla gap piuttosto intuitivo che ogni inserzione o delezione dovranno contribuire nel calcolo dello score totale con termini negativi, dato che allaumentare del numero di gaps durante un allineamento aumenta la probabilit che le due sequenze siano scorrelate (e quindi un punteggio totale pi basso). allo stesso modo intuitivo che il peso della penalizzazione tanto pi alto quanto pi lunga la stringa di amminoacidi appartenenti alla gaps.

Come abbiamo gi accennato in precedenza il secondo step fondamentale per far si che sia attendibile la valutazione finale dellallineamento, chiedere che lo stesso allineamento tra la coppia di sequenze presa in esame sia ottimale. Prenderemo ora in considerazione solamente allineamenti globali (quindi sequenze della stessa lunghezza n). Gli allineamenti globali possibili per queste due sequenze sono:

un numero talmente elevato (si pensi che una proteina ha solitamente un numero maggiore di 100 residui) che non sarebbe computazionalmente possibile (in tempi ragionevoli) esplorare tutte le possibili combinazioni di allineamenti per determinare quella con lo score maggiore, ovvero l allineamento ottimale. Utilizzeremo pertanto degli algoritmi di computazione dinamica (due in particolare) che ci permetteranno di ottimizzare lallineamento. Distingueremo i due algoritmi a seconda se si tratta di un allineamento globale o locale. Algoritmo Needleman Wunsch : allineamento globale 1. Come prima cosa si costruisce una matrice F con dimensione il prodotto delle dimensioni lineari delle due sequenze comparate. Se le due sequenze saranno della stessa lunghezza n allora la matrice F sar nxn. La matrice sar costruita disponendo su un lato la prima sequenza e su laltro la seconda, e sar riempita attraverso i valori delle frequenze di scambi tra amminoacidi controllate su una matrice BLOSTUM o PAM. Riempita la matrice F, si parte dallultimo elemento di matrice in basso a destra per risalire verso lalto attraverso un percorso che segue le seguenti condizioni:

2.

3.

Per far si che lalgoritmo sia giusto bisogna inserire le condizioni iniziali sugli elementi di matrice e considerare i termini di penalizzazione delle possibili gaps.

Secondo questo algoritmo lelemento di matrice finale del percorso individuato lo score migliore che potremmo ottenere dallallineamento delle suddette sequenze. Algoritmo Smith Waterman: allineamento locale Lallineamento locale si applica generalmente per proteine che non hanno la stessa lunghezza, ma di cui si suppone possibile che abbiamo dei domini in comune (quindi delle sottosequenze simili). Un altro caso invece riferito al confronto tra il genoma del DNA. Lalgoritmo che ci permette di ottimizzare lallineamento locale semplicemente un caso particolare di quello globale, ha pertanto la stessa idea di fondo e la stessa struttura:

La differenza sostanziale tra i due algoritmi che questultimo pu terminare in ogni elemento di matrice, al contrario di quello gi visto per lallineamento globale che termina sempre nella prima casella in alto a sinistra (lultima del percorso ricorsivo tracciato secondo le condizioni dellalgoritmo stesso). Metodi di ottimizzazione euristici Gli algoritmi di programmazione dinamica appena descritti sono metodi corretti ma molto lenti, e in molti casi (per sequenze lunghe) la lentezza e la memoria da utilizzare diventano problemi computazionali non trascurabili. Per questo motivo si ricorre alluso di metodi euristici, i quali non godono della stessa precisione ma hanno dei tempi di risposta molto pi brevi. Lo scopo di questi metodi quello di cercare la pi piccola frazione possibile delle celle nella matrice di programmazione dinamica, mentre nello stesso tempo ricercare tutti gli allineamenti che danno un punteggio alto. 1 metodo: BLAST consente di ricercare le similarit tra una sequenza data, chiamata query ed un database di sequenze, chiamate seeds. Il metodo consiste in unapprossimazione dellalgoritmo di Smith Waterman rendendolo in questo modo 50 volte pi veloce. Possiamo dividere lalgoritmo di BLAST in tre fasi: i) Cerca nelle sequenze della banca dati pezzi in comune, esatti e senza gaps, con la sequenza query sopra un certo valore di soglia dello score. La lunghezza della stringa da ricercare cambia a seconda se si stratta di sequenze proteiche (3) o se si stratta di sequenze di DNA (11). Il processo di allineamento estende i seeds iniziali di lunghezza W in entrambe le direzioni nel tentativo di incrementare lo score in modo che scenda al di sotto del valore di soglia definito. Solo ora BLAST opera la costruzione di un allineamento con gap tra la query e quelle della banca dati. Infine vengono mostrati gli allineamenti statisticamente significativi.

ii) iii)

La significativit di un allineamento generato da BLAST stabilito dal valore del parametro E-value che rappresenta il numero di allineamenti differenti con score maggiore o uguale a S (ovvero quello ottenuto per il mio allineamento, che mi aspetto di trovare per caso in una ricerca di sequenza in banca dati)

2 metodo: FASTA che consente di cercare una sequenza query in un database di sequenze, le subject. Anche questo prevede tre fasi: i) ii) Inizializzazione: la query viene divisa in parole di lunghezza pi piccola (2,3,4 caratteri) e si memorizzano tutte le posizioni di inizio parola. Inizia cos lanalisi di tutte le sequenze subject. Ogni volta che il programma trova una parola comune sia alla query che alla subjet viene memorizzata sulla tabella relativa la corrispondente posizione dellindice. Una volta terminata la lettura della sequenza subject vengono estratte le pi lunghe sottosequenze ottenute. Tramite una matrice di sostituzione viene calcolato lallineamento locale. Alla fine di questa fase viene compilata una graduatoria delle migliori similarit ottenute, che saranno quelle con cui si lavorer nelle fasi successive. Raffinamento: il programma tenta di migliorare lallineamento congiungendo le best initial region con le gap. Sulle sequenze che hanno ottenuto il migliore punteggio viene applicata una variante dellalgoritmo di smith-waterman.

iii)

Fase finale dellallineamento: valutazione dello score La valutazione dellefficacia di un metodo un discorso del tutto generale che va affrontato ogni volta che lavoriamo con modelli che seguono una certa statistica. La valutazione di quanto siamo sicuri che due sequenze siano o non siano omologhe sicuramente una tappa fondamentale per far si che le nostre conclusioni siano attendibili. Il primo valore importante nel giudicare un allineamento pertanto il punteggio associato che non ha una validit generale (basti pensare che la sola scelta della matrice di sostituzione utilizzata, dettata dallesperienza, influisce sul risultato finale). Questo perch il problema principale che riguarda lallineamento riguarda la possibilit di ottenere dei falsi positivi, ovvero dalla possibilit di considerare due sequenze omologhe anche se realmente non lo sono (perch lallineamento fornisce un valore dello score sopra il livello di soglia). Attualmente il metodo pi utilizzato per la valutazione dellallineamento quello di considerare gli allineamenti casuali. A tale scopo viene utilizzato un parametro chiamato Z-score ottenuto statisticamente generando un numero sufficientemente grande di allineamenti casuali, sui quali si calcola la media dei punteggi degli allineamenti e la deviazione standard. Se S il punteggio ottenuto dallallineamento che voglio valutare allora lo Z-score associato sar:

tale valore ci dice a che distanza il punteggio da valutare si trova rispetto al punteggio medio degli allineamenti casuali generati . Prima di passare allallineamento multiplo di sequenze soffermiamoci sullimportanza di stabilire se due sequenze siano omologhe o meno. Notiamo che chiedere lomologia significa chiedere la correlazione evolutiva tra due sequenze, questo utile per varie ragioni: 1. 2. 3. Identificazioni di quali sono le regioni conservate. Inferenza funzionale Modellizzazione funzionale

2. Allineamento multiplo di sequenze (MSA)Lo scopo quello di comporre una sequenza con una famiglia di sequenze per individuare una possibile similarit, che indicherebbe (come ormai noto) una possibile omologia. Gli allineamenti multipli sono spesso utilizzati in bioinformatica perch riescono a generare allineamenti pi accurati rispetto agli allineamenti tra coppie di sequenze. Analogamente a quanto discusso in precedenza, per lallineamento multiplo si utilizzer una programmazione dinamica a N dimensioni. Questo corrisponde pertanto alla generalizzazione degli algoritmi visti per lallineamento tra coppie, utilizzando matrici di sostituzione a N dimensioni invece che a 2 dimensioni. Per trovare lallineamento multiplo ottimale tra sequenze come prima cosa generalizziamo la funzione che utilizzavamo per comparare 2 sequenze S(i1,i2,.,in), la inizializziamo S(0,0,0.0) = 0 e infine definiamo:

Dove il secondo termine rappresenta lo score di ogni colonna e dove x un vettore binario diverso da zero:

Quando scegliamo la funzione di score per MSA lapproccio generale quello di misurare la mutua similarit (o conservazione) dei singoli caratteri della colonna. La possibile funzione di score include le seguenti caratteristiche:

Ovvero una colonna ha il punteggio (lo score) alto se ci sono molti residui identici o simili. La seconda propriet il C-G score che misura il punteggio di una colonna rispetto al centro di massa di un set di residui. Questo lamminoacido che minimizza le similarit con tutti gli altri residui della colonna:

Alte conservazioni delle colonne sono una forte indicazione della presenza di un ancestrale comune tra le sequenze. In ultimo, la misura dellentropia riflette lordine osservato allinterno di un set di misurazione ed definita:

Dove P(a) rappresenta la probabilit empirica del carattere a nella colonna. Questo implica che lentropia di una colonna completamente conservata nulla.

Anche la programmazione dinamica a N dimensioni incontra gli stessi problemi incontrati per lallineamento a coppie riguardanti il tempo e la memoria necessaria per lapplicazione degli algoritmi globali e locali. Bisogna anche in questo caso ricorrere a metodi euristici, come lallineamento a stella e lallineamento ad albero. Lapprendimento automatico Lenorme quantit di dati che contengono informazioni biologiche (come nel caso di un allineamento multiplo tra tante sequenze proteiche mediamente lunghe) sono state per un lungo periodo di tempo trattate esclusivamente dagli algoritmi studiati. Tuttavia, tale approccio non pi sempre in grado di studiare con efficienza tutti i dati a disposizione. Si pone cos un problema reale sullanalisi di sequenze nel caso in cui abbiamo a diposizione una quantit di dati superiore alla capacit di analisi degli algoritmi di base. Lapproccio dellapprendimento automatico, come le reti neurali, i modelli nascosti di Markov e le reti bayesiane, adatto a studiare tutti i set di dati di grandi quantit, talvolta anche rumorosi, e dietro i quali non disponiamo di teorie generali per la possibile interpretazione. Lidea fondamentale di questo approccio quella di imparare la teoria direttamente dai dati, e pu portare a risultati soddisfacenti se utilizzato in maniera complementare ai metodi convenzionali. Prima di vedere in dettaglio i metodi dellapprendimento automatico, sar meglio chiarire il concetto di programmazione dinamica di cui abbiamo gi parlato per gli algoritmi utilizzati nellallineamento ottimale di sequenze. Pi generalmente il termine Programmazione Dinamica si riferisce ad una tecnica generale di ottimizzazione che pu essere applicata ad un problema scomponibile in problemi simili tra loro ma di calibro minore. In questo modo la soluzione del problema originario pu essere vista come somma di sottoproblemi pi piccoli. Un tipico problema che riguarda la programmazione dinamica quello di trovare il percorso pi breve tra due nodi di un grafo. abbastanza intuitivo che il percorso pi corto che unisce A con B e che passa per C, sar lunione del percorso pi corto che unisce A con C e il percorso pi corto che unisce C con B. Questo semplice principio prende il nome di principio di Bellman. In questo senso, la soluzione del problema generale originario costruita ricorsivamente unendo le soluzioni di sottoproblemi pi semplici. Allo stesso modo allineare due sequenze di lunghezza N richiede di trovare il percorso pi breve in un grafo di NxN vertici. Poich la programmazione dinamica richiede essenzialmente di visitare TUTTI i vertici del grafo, ne segue che la complessit in tempo dellordine di O(N^2). In generale gli algoritmi di apprendimento possono essere suddivisi in tre categorie a seconda dellobiettivo: 1. 2. 3. La stima completa della distribuzione e della propriet a posteriori P(w|D) La stima dellinsieme ottimo dei parametri w massimizzando P(w|D) La stima dei valori medi rispetto alla probabilit a posteriori.

I modelli nascosti di Markov Spesso, soprattutto quando si ha a che fare con frammenti di sequenze biologiche come il genoma umano, risulta abbastanza complicato il confronto tra la sequenza da analizzare con quelle presenti nellarchivio. I HMM formano una classe di modelli grafici probabilistici molto utilizzata a tale scopo. Infatti questi modelli possono servire a confrontare non coppie di sequenze, ma insiemi di sequenze, per scoprirne le caratteristiche comuni. Definizione: Un HMM discreto del 1 ordine (perch dipende solo dallo stato corrente) un modello probabilistico per le serie temporali definito da: un insieme finito di stati S un alfabeto discreto di simboli A una matrice delle probabilit di transizione T=P(i|j) una matrice delle probabilit di emissione E=P(X|i)

dove P(i|j) indica la probabilit di transizione dal nodo j al nodo i e P(X|i) indica la probabilit che il nodo i emetta il simbolo X. Il sistema evolver casualmente da uno stato allaltro, ogni volta emettendo un simbolo dellalfabeto (dim4 per il DNA e RNA e dim20 per le proteine). I HMM possono essere utilizzati diversamente a seconda di come formulato il nostro problema, ci consentono di: 1. 2. 3. trovare la probabilit di occorrenza di una sequenza a partire da un modello (stima) identificare la sequenza di stati che pi probabilmente ha generato una certa sequenza (decodifica) generare un HMM basandosi su una sequenza di osservazioni (apprendimento)

I modelli di Markov rappresentano semplicemente unestensione delle catene di Markov, ovvero un modello probabilistico in cui la probabilit di un evento dipende solamente dallevento precedente. Una catena discreta si dice di Markov se:

Dato che gli stati sono in numero finito, come anticipato nella definizione di cui sopra, allora una catena di Markov pu essere rappresentata nella forma di un automa stocastico, ovvero di un grafo in cui i nodi corrispondono agli stati (o eventi) e le frecce che collegano i nodi corrispondono alla probabilit di passare da uno stato allaltro.

Il grafo mostrato rappresenta un automa stocastico che genera sequenze di DNA (automa generativo). Calcoliamoci ora analiticamente la probabilit di occorrenza di una sequenza X = X1,..,Xn. Assumiamo che s e t appartengono allinsieme degli indici che costituiscono lalfabeto del caso in questione e poniamo

Allora otteniamo che:

N.B. In una catena di Markov c una corrispondenza biunivoca tra i simboli emessi dallautoma e gli stati corrispondenti. Cosa succede per un HMM? Questa propriet non vale pi perch gli stati interni sono appunto

nascosti, pertanto losservatore pu accedere soltanto ad una sequenza di simboli osservabili. Solo attraverso questa potr ricavare la probabilit degli stati corrispondenti. Esempio esplicativo: un naufrago su un isola si accorge per esperienza dellesistenza di una correlazione tra le previsioni del tempo e lanalisi delle alghe del mare. Se volesse conoscere con una certa probabilit le previsioni metereologiche dovr fare uninferenza sullo stato delle alghe osservate. Vediamo come i gruppi di stati sono due: quelli OSSERVABILI (lo stato delle alghe) e quelli NASCOSTI (le previsioni meteorologiche). I HMM permettono di risolvere questo problema senza necessariamente implementare un algoritmo specifico per ogni caso distinto. Risulta pi chiaro ora che un modello nascosto di Markov costituito da una catena di Markov nascosta (ovvero con stati interni) e da un processo osservabile (ovvero la sequenza di simboli), descritto mediante funzioni probabilistiche. La scelta dellarchitettura per lHMM dipende molto dal tipo di problema. Nel caso di sequenze biologiche, laspetto lineare delle sequenze ben rappresentato dallarchitettura left-right (architettura lineare standard):

Questo tipo di architettura non permette il ritorno in uno stato dopo una transizione da questo. Oltre ai due stati START ed END ci sono altri tre tipi di stati intermedi chiamati MAIN, INSERT e DELETE. I nodi main ed insert emettono un simbolo, mentre i delete sono muti. Questi vengono inseriti, nel caso della modellizzazione di una famiglia di sequenze correlate tra loro, proprio perch la lunghezza delle sequenze non costante tra loro ma sono presenti cancellazioni e inserimenti di amminoacidi.

Applicazione biologica dei HMM: Data una certa sequenza, qual la probabilit che sia stata generata dallHMM dato? (e quindi che appartenga alla famiglia di sequenze che lha generata?) Qual la sequenza pi probabile di match, insert e delete che ha generato la sequenza data? Un esempio di HMM molto utilizzato quello del casin disonesto:

In un casin utilizziamo per la maggior parte del tempo un dado, per poi passare per una frazione di tempo pi piccola ad un dado truccato. Tale dado ha probabilit 0.5 che esca 6 e 0.1 per i numeri che vanno da 1 a 5. Il passaggio dal dato truccato a quello giusto un processo di Markov :

Cosa c di nascosto nel modello dei dadi del casin appena descritto? Se tu puoi vedere solo una sequenza di dati che esce dal lancio dei dadi, non potrai distinguere se si sta utilizzando il dado truccato o quello corretto, perch questa un informazione tenuta segreta dal casin. La domanda che ci poniamo : la mia sequenza di amminoacidi pu essere stata creata da questo schema? Qual la sequenza di eventi pi probabile che lha generata? In generale possiamo generare da questo schema dado giusto dado truccato (ovvero da questo modello di Markov nascosto) una sequenza nel modo seguente: viene scelto uno stato K1 secondo la probabilit a(0,i). In questo stato unosservazione emessa secondo la distribuzione e relativa a questo stato. Un nuovo stato k2 viene scelto secondo la transizione di probabilit a(k1,i) e cos via.

In questo modo viene formata una sequenza random artificiale. P(X) la probabilit che la sequenza x venga generata dal modello. La probabilit congiunta di osservare una successione di simboli x generata da una sequenza di stati k :

Tuttavia nella maggior parte dei casi la sequenza degli stati attraversati non nota a priori, allora ricavare il cammino nel grafo dellHMM che ha portato alle generazione di una certa stringa si pu ridurre semplicemente alla determinazione di quello pi probabile (Most Probable Path, MPP). Ovvero alla sequenza di stati

Che massimizza la probabilit delle varie transizioni. La propriet delle catene di Markov permette di calcolare questo mediante un algoritmo di programmazione dinamica noto come ALGORITMO DI VITERBI. Algoritmo di Viterbi: lidea generale dellalgoritmo quella di sostituire i calcoli relativi a tutti i percorsi con altri che invece riguardano solo un numero piccolo di percorsi (tipicamente 1). Pertanto vengano ignorati tutti i percorsi possibili tranne il pi (o i pi) probabile, sfruttando in questo modo solo parte dellinformazione disponibile che pu offrire lHMM.

Nonostante lapprossimazione fatta per lalgoritmo Viterbi fornisce buoni risultati, non sempre necessario ricorrervi. Sfruttando le propriet di Markov possibile progettare un altro algoritmo di programmazione dinamica simile a questo.

Algoritmo di Forward:

Pertanto i HMM sono ottimi metodi utilizzati che rispondo alle seguenti domande poste allinizio del capitolo:

1. 2.

3. Gli alberi (1 metodo di classificazione)Sono molti gli algoritmi di allineamento multiplo che utilizzano gli alberi. Vediamo nello specifico cos e come si costruisce un albero filogenetico: Tutti gli organismi della terra hanno un ancestrale comune e la relazione tra le specie viventi viene chiamata filogenetica. Utilizzando nello specifico delle sequenze di DNA possiamo dedurre quando c un collegamento genetico e rappresentare tale relazione da un albero. Limportanza di una tale ricostruzione filogenetica stata chiarita anche nelle sezioni precedenti: le relazioni filogenetiche possono aiutare nella predizione di funzione, e questo rappresenta uno dei problemi ancora aperti della bioinformatica. Le procedure per le analisi filogenetiche sono strettamente collegate con quelle per gli allineamenti di sequenze (a coppie o multipli). Riportiamo nella seguente figura le possibilit di evoluzione che permettono le diversificazioni (collegate temporalmente) tra sequenze:

I geni che si sono diversificati a causa di speciazioni si chiamano ORTOLOGHI I geni che si sono diversificati a seguito di una duplicazione genica si chiamano PARALOGHI

Considereremo sempre gli alberi binari (nel senso che da ogni nodo possono partire solamente due rami perch questo quanto ci si aspetta dallevoluzione!) dato che, anche alberi di dimensioni maggiori, possono essere ricondotti ad alberi binari con laggiunta di piccole braccia. Un vero albero genetico ha una radice, cio esiste un progenitore iniziale di tutte le altre sequenze. Questo viene comunemente raffigurato attraverso un albero di tipo root. Esistono algoritmi che riescono a fare una congettura sulla sua posizione. Tuttavia altri non possono stimarla e quindi occorre utilizzare altri criteri. Gli alberi che non hanno un progenitore iniziale definito vengono detti unrooted. Questi sono pertanto alberi i cui nodi non sono collegati temporalmente.

Un albero root deriva da un albero unrooted dopo la scelta della posizione del root principale, ovvero dellinizio dellalbero. La root pu essere collocata in qualsiasi posizione quindi a parit di nodi per ogni albero unrooted ci sono 2n-3 alberi root

Se un albero ha n foglie allora ha (n-1) nodi e (2n-2) rami. Per una corretta definizione di albero filogenetico dovremmo anche definire la lunghezza delle sue braccia. Definiamo pertanto lunghezza di ogni ramo come il numero di cambiamenti di sequenza che avvengono fino al nodo successivo. Alcuni tipi di analisi filogenetica assumono che il tasso di evoluzione nei rami sia costante, ovvero stessa lunghezza dei rami nellalbero). Questa viene chiamata ipotesi di molecular clock, utilizzata anche nellUPGMA (metodo per la costruzione di un albero). La somma della lunghezza di tutti i rami detta lunghezza dellalbero. Come costruire un albero: Lidea di fondo quella di calcolare le distanze tra gli elementi dellinsieme che dovranno essere disposti nellalbero. La distanza tra le sequenze pu essere definita in vari modi: 1. 2. La frazione dei siti in cui gli amminoacidi sono diversi Distanza di Jukes Cantor definita come: -3/4log(1-4f/3) derivata da un modello di Markov la quale tende ad infinito per f = 75%

Metodo UPGMA : Assume la stessa velocit di evoluzione in tutti i rami dellalbero. Inoltre la distanza tra ciascuna coppia di foglie la somma delle lunghezze delle linee (edge) che li connettono (assunzione di additivit). Limportanza nel definire le distanza tra coppie di sequenze risiede nella possibilit di costruire una matrice delle distanze. chiaro che sequenze meno distanti tra loro secondo lindice riportato nella matrice delle distanza si troveranno ad una distanza minore anche nellalbero stesso. Per cluster intendiamo un raggruppamento di pi sequenze in modo tale che la somma fra le distanze di tutte le sequenze presenti nello stesso gruppo sia minima. In pratica: unire in gruppi le sequenze che sono pi vicine tra loro. Definiamo la distanza tra due clusters Ci e Cj come la distanza media tra coppie di sequenze di ogni cluster

Il metodo UPGMA inizia scegliendo allinterno della matrice delle distanze la coppia di sequenze con la distanza minima. Le due sequenze formano in questo modo un primo cluster (se ci sono distanze uguali la scelta sar random). Questo cluster prender il posto delle due sequenze nella nuova matrice delle distanze. Definire un nuovo cluster k da Ck = C1 U Cj e definire d(kl) per tutti gli l tali che:

Definire un nodo k con nodi figli i e j e porre questo ad altezza d(ij)/2 Aggiungere k al cluster corrente.

Riportiamo di seguito lo schema grafico seguito dal metodo UPGMA per la costruzione di un albero.

Prima fase

Seconda fase

Terza fase

Metodo di Neighbour Joining: se si assume per lUPGMA ladditivit ma non la costanza della velocit di evoluzione, allora si pu usare il metodo di NJ. Vediamo in breve quale sia lidea di base di tale algoritmo. Dato un albero T con lunghezze additive d: Trova una coppia di nodi adiacenti i e j Li elimina dalla lista dei nodi e li sostituisce con il nodo k la cui distanza da un nodo n definita:

Il processo continua fino a rimanere con un nodo solo.

Metodo di massima parsimonia: Questo metodo utilizza tuttaltra strategia rispetto ai metodi sulle distanze che abbiamo visto prima. Questo metodo infatti cerca lalbero che pu spiegare le sequenze osservate con il numero minimo di sostituzioni. Pertanto dobbiamo: A) Calcolare il costo di un albero B) Valutare tutti gli alberi per cercare il minimo Esempio: Abbiamo le sequenze riportate e costruiamo gli alberi possibili:

Possiamo associare un costo S(a,b) ad ogni sostituzione. Metodi possibili per la ricerca dellalbero: 1) Invertire a caso coppie di rami e tenere il nuovo albero se il suo costo minore 2) Scegliamo 3 sequenze a caso. Unaltra sequenza viene scelta e aggiunta al ramo che produce lalbero a costo minore e cos via. 3) Branch and Bound: cominciamo a costruire alberi sistematicamente aumentando man mano il numero di foglie ma ci fermiamo appena lalbero incompleto ha un costo maggiore di quello fino ad ora ottenuto per un albero completo. Qualit di un albero: bootstrap Metodo utilizzato per determinare la qualit dellalbero: Dato un set di sequenze allineate, creiamo un insieme di nuovi allineamenti tra sequenze artificiali (tipicamente 100) della stessa dimensione del nostro set iniziale. Selezioniamo colonne a caso dellallineamento anche con rimpiazzo Ripetiamo la costruzione dellalbero con i nuovi data set La frequenza con cui una certa caratteristica dellalbero viene ritrovata presa come misura di confidenza.

Cosa possiamo fare con gli alberi? 3. Ricostruire la filogenia Aiutare nella costruzione di allineamenti multipli Usarli come tecniche di clustering Usarli per identificare coppie di famiglie che co-evolvono E

4. Altri metodi di classificazionePartiamo dal presupposto che le tecniche di classificazione sono il cuore dellanalisi dei dati biologici. Tratteremo in questa sezione altri metodi di classificazione oltre agli alberi che forniscono solitamente una classificazione complessa, con la speranza che sia in qualche modo basata sulla distanza evolutiva. Allo scopo di riuscire a classificare dei dati a noi disposizione, dove per varie ragioni non possibile lanalisi attraverso gli algoritmi tradizionali, utilizzeremo ancora una volta i metodi di apprendimento automatico (Machine Learning), i quali si occupano di realizzare sistemi basati su osservazioni, o esempi noti, per la sintesi di nuova conoscenza (in questo caso la classificazione dei dati stessi). Ogni metodo di clustering pu appartenere ad una delle due categorie: i) Apprendimento supervisionato: cerca la suddivisione che meglio spiega in nostri dati ii) Apprendimento non supervisionato: divide i dati secondo qualche propriet dei dati stessi. A noi interesseranno i metodi di apprendimento automatico supervisionato. Lo scopo sar quello di istruire un sistema informatico in modo da consentirgli di risolvere dei compiti in automatico. Bisogna definire i dati in ingresso, di solito in forma vettoriale, come un insieme di esempi I; Poi si definisce l'insieme dei dati in uscita come insieme di output O (gli output possono essere unetichetta numerica, ad esempio 1); Infine bisogna definire una funzione f che associa ad ogni dato in ingresso (I) la sua risposta corretta (O).

Tutti gli algoritmi di apprendimento supervisionato partono dal presupposto che se forniamo all'algoritmo un numero adeguato di esempi l'algoritmo sar in grado di creare una funzione f1 che approssimer la funzione f. Se l'approssimazione di f risulter adeguata, quando proporremo ad f1 dei dati in ingresso, mai analizzati in precedenza, la funzione dovrebbe essere in grado di fornire delle risposte in uscita simili a quelle fornite da f (quindi corrette entro una certa approssimazione). Si pu facilmente intuire che il buon funzionamento di questi algoritmi dipende in modo significativo dai dati in ingresso; se si forniscono pochi ingressi l'algoritmo potrebbe non aver abbastanza elementi per apprendere, mentre molti dati in ingresso potrebbero renderlo eccessivamente lento, dato che la funzione f1 generata dagli ingressi potrebbe essere molto complicata. Questi algoritmi sono molto sensibili al rumore, anche pochi dati errati potrebbero rendere l'intero sistema non affidabile e condurlo a decisioni errate. Appartengono alla categoria dellapprendimento automatico supervisionato i metodi delle reti neurali , del random forest e del SVM. Per tutti questi metodi useremo sempre il seguente processo: 1. Ci serve un set di dati per i quali conosciamo la risposta e li utilizziamo per addestrare il sistema. Ovvero modificheremo dei parametri in maniera tale che la risposta fornita dal sistema sia quanto pi vicino possibile a quella vera che conosciamo. (training-set) Poi ci serve un set di dati differenti dalla prima (anchessi con risposta nota) e vediamo quanto bene funziona il sistema. In modo tale che siamo sicuri che il sistema funzioni anche quando non conosciamo la risposta sui dati di cui disponiamo (ci che realmente ci interessa). (test-set) Infine serve un terzo set di dati per rivalutare il tutto, aggiustando le differenze riscontrate nel secondo punto, ed avere definitivamente un sistema che con buona probabilit mi potr fornire una risposta sensata al mio set sconosciuto di dati.

2.

3.

N:B: importante che ogni set di dati sia usato una e una sola volta! Altrimenti non possiamo essere sicuri della veridicit del sistema.

Reti Neurali Artificiali Qual la logica di una rete neurale? Nel cervello un neurone funziona cos: Ogni neurone connesso con una serie di altri neuroni. Se la somma dei segnali che un neurone riceve superiore ad una certa soglia che chiameremo B allora il neurone emetter a sua volta un segnale. Altrimenti non succeder nulla.

Ogni neurone pu avere pi segnali di ingresso (input) ed un solo segnale di uscita (output). Questo significa che per una rete a pi layer ogni nodo connesso a tutti i nodi del layer successivo. Ognuno di questi ultimi emetter solo se la somma dei valori in uscita dei nodi del layer precedente superer il valore di soglia. Ogni connessione ha un peso associato e la funzione che ora stiamo considerando come una somma pu anche essere generalmente una funzione diversa che verr integrata per verificare se rientra o supera il valore di soglia stabilito. Gli esempi in cui richiedo una classificazione binaria possono essere tanti: un esempio pu essere proteine extracellulari o intracellulari. Assoceremo questi a casi positivi (1) e a casi negativi (0).

Come prima cosa bisogna addestrare la rete. Quindi gli dar in pasto una certo numeri di casi positivi (magari extracellulari) e mi aspetto che la rete mi restituisca una risposta di dati positivi. Se cos non , allora dovr modificare i singoli pesi relativi ad ogni freccia che collega i nodi della rete in maniera da ottenere la risposta desiderata (ovvero casi positivi). Questo viene chiamato training-set, ovvero se risponde correttamente il maggio numero di volte. Ora devo vedere se la rete ha imparato a generalizzare (test set). Ovvero vediamo quante volte la rete mi da la risposta giusta e quante volte mi da la risposta sbagliata.

Dobbiamo stare attenti alle possibili propriet che classificano i due tipi di proteine perch potrebbero creare confusione alla rete. Per esempio se le proteine extracellulari sono mediamente pi lunghe di quelle intracellulari, la rete pu confondere la classificazione da fare. Questo significa che dobbiamo stare molto attenti al tipo di set di dati che vogliamo dare impasto alla rete. Dobbiamo fare in modo, infatti, che i dati siano tra loro quanto pi simili possibili salvo per la caratteristica per i quali li vogliamo classificare (evitando in questo modo confusione nella rete). Le funzioni che si utilizzando per la divisione tra le due classi di dati possono essere delle pi svariate.

Importante la fase del valutation-set perch se la rete non fornisce in buona percentuale il risultato corretto, ma i parametri di tutti i collegamenti tra i nodi sono gi stati ottimizzati, allora possibile che bisogna cercare la migliore funzione. A questo punto, dopo il cambio della funzione, rifaccio il training-set e se il risultato del successivo test-set migliore allora me lo tengo, altrimenti ricomincio con altre funzioni. (N.B. la nostra valutazione in questo modo non oggettiva, perch nel test set vogliamo che la rete agisca su parametri che non ha mai visto. Ogni set deve essere usato una e una sola volta!). Questo significa che ad ogni passaggio va cambiato il set. Ovvero che nella fase del test la rete veda cose nuove. Esempio: predire la struttura secondaria delle proteine Sappiamo che le proteine sono formate da strutture secondarie che possono essere o alfa eliche o foglietti beta. Voglio che la rete deduca quali sequenze hanno maggiori probabilit di essere alfa elica e quali maggiore probabilit di essere foglietti beta. metodo 1 Prendo strutture con forma nota. Prendo le sequenze di tutte le alfa eliche e di tutti i foglietti beta.

Come facciamo a riportare la sequenza di amminoacidi in termini numerici per poter essere analizzati dalla rete? evidente che non possiamo numerarli 1,2,3,4,20 perch altrimenti la rete vedrebbe distanze diverse tra il primo e lultimo amminoacido rispetto che tra il primo ed il secondo. Allora possiamo pensare ogni amminoacido come un vettore di dimensione 20 e in modo tale che insieme formano una base completa: (1,0,0,.,0), (0,1,0,,0), .., (0,,1) Li do in pasto alla rete e chiedo che me li classifichi. Io gi conosco la risposta e procedo con la massimizzazione dei parametri. Daremo alla rete gruppi di tre amminoacidi e la risposta apparterr solamente allamminoacido centrale (chiedo se appartiene ad una sottosequenza di foglietto beta o alfa elica). Questo un buon metodo per la predizione della struttura delle proteine, ma potrei risolvere il problema partendo da unaltra costatazione.

Metodo 2 Sappiamo che proteine omologhe hanno struttura simile. Allora prendo la mia proteina, faccio una ricerca in banca dati e ottengo un allineamento multiplo di proteine omologhe. Se di una proteina che appartenente allallineamento ne conosco la sua struttura allora posso affermare con una buonissima approssimazione che: se una data regione di questa proteina un alfa elica allora la stessa regione sar un alfa elica in tutte le altre (lo stesso vale per i foglietti beta). Per ogni allineamento multiplo si pu definire un profilo ovvero la tabella della frequenza di ciascun amminoacido per ciascuna posizione. Questo significa che i vettori che costituiscono la tabella del profilo sono 20 e di dim20, in cui in ogni coordinata riporter quanto frequente quellamminoacido in quella data posizione nella famiglia di proteine omologhe identificate in banca dati. Dando in pasto tali vettori nella rete possiamo avere in risposta delle predizioni sugli amminoacidi. Ad ogni sequenza nota quindi possibile ricavare la struttura della proteina seguendo il procedimento sopra descritto. Tenendo conto per dellerrore che ci portiamo dietro dallallineamento multiplo di sequenze, proprio perch la qualit del risultato dipender dalla qualit dellallineamento. SVM Le macchine a supporto vettoriale (Support Vector Machines, SVM) sono un'insieme d i algoritmi per la regressione e la classificazione di pattern. Le SVM possono essere pensate come una tecnica alternativa alle reti neurali che permettono lapprendimento mediante classificatori polinomiali. In parole semplici, una SVM un classificatore binario che apprende il confine fra esempi appartenenti a due classi diverse. Funziona proiettando gli esempi in ingresso, aventi una certa dimensione nello spazio multidimensionale, mediante una funzione che opera una mappatura dei dati nel nuovo spazio, e cercando un iperpiano di separazione ottimale i n questo spazio.

In alto a sinistra vediamo come sono mischiate le due classi di esempi (classe x e classe o)nello spazio di ingresso; mediante la funzione di mapping ,, i due insiemi vengono proiettati in uno spazio di dimensione maggiore dove possibile trovare un iperpiano ottimale di separazione.

L'iperpiano di separazione massimizza la sua distanza (il margine) dagli esempi di training pi vicini. Le propriet pi importante delle SVM quella di saper gestire dati con molte caratteristiche descrittive.

Random Forest tra i tre il metodo pi recente. Lalgoritmo segue la seguente procedura: Supponiamo di avere N casi del mio data set (e li conosco) e M variabili con M