Biologia computazionale

Biologia computazionale

A.A. 2010-2011 semestre II

UNIVERSITÀ DEGLI STUDI DI MILANODocente: Giorgio Valentini

Istruttore: Matteo Re

6 Evoluzione e filogenesi - 3

C.d.l. Biotecnologie Industriali e Ambientali

Bio Metodi per costruire alberi filogenetici

Metodi basati su:

• Distanza• Massima parsimonia• Massima verosimiglianza

Questi li abbiamo visti…

Oggi discutiamo questa classe di metodi …

Bio Massima verosimiglianza

Verosimiglianza (likelihood) :

• Probabilità delle osservazioni dato un modello• Quindi è una probailità … perché usare un nome

diverso?• Per porre l’accento sul fatto che non vogliamo

valutare quanto siamo confidenti nell’occorrenza di un determinato evento ma piuttosto valutare quanto i dati sono compatibili con un modello evolutivo che abbiamo scelto.


ESEMPIO : lanciamo una moneta ed otteniamo croce (questo è il dato).

Se dovessi chiedervi qual’è la probabilità dell’evento “osservo croce” probabilmente mi rispondereste ½ .

Questo implica che avete ipotizzato un modello di “moneta onesta” in cui le probabilità di testa e croce sono entrambe uguali a ½ .


ESEMPIO : lanciamo una moneta ed otteniamo croce (questo è il dato).

Supponiamo di definire un modello di moneta con queste caratteristiche: P(testa=1), P(croce=0) … ossia una moneta truccata. I parametri (tutti) del modello li indichiamo complessivamente come Θ

La likelihod dell’osservazione “croce” dato il modello è zero (il che dovrebbe farci venire il dubbio che il modello non è adatto a descrivere i dati osservati)

Se utilizzassimo un modello di moneta truccata (due croci) … la likelihod

dell’osservazione sarebbe uno


Quindi la likelihood è la verosimiglianza di un insieme di osservazioni rispetto ad un modello che dovrebbe descrivere il processo da cui i dati sono stati generati.

Quindi per valutare la verosimiglianza di un albero filogenetico mediante la tecnica della massima verosimiglianza (maximum likelihood) abbiamo bisogno innanzitutto di un modello evolutivo adatto alle sequenze biologiche.

Ma come possiamo costruire un tale modello?


Nel caso dell’evoluzione molecolare i dati sono rappresentati da un allineamento di sequenze ed il modello, in senso molto ampio, è l’albero filogenetico che:

• correla tra di loro le sequenze

• descrive il meccanismo di evoluzione da una sequenza all’altra


L’albero filogenetico ed il modello che descrive il meccanismo attraverso il quale si verificano gli eventi evolutivi, insieme, costituiscono la nostra “ipotesi” rispetto al modo in cui l’evoluzione ha generato le sequenze che stiamo osservando.

Consideriamo le due parti separate : ci riferiamo alle relazioni tra le sequenze (i dati) con il termine “albero filogenetico” mentre ci riferiamo alla parte che descrive il meccanismo evolutivo come “modello”.


L’obiettivo del modello è quello di descrivere il meccanismo attraverso cui le sequenze cambiano nel tempo.

Per semplificare i calcoli ci occuperemo di modelli di sequenze di DNA. Immaginiamo inoltre il modello come diviso in due parti principali:

1) Composizione

2) Processo

descrive le frequenze con cui le parti della sequenza (nt)

cambiano nel tempo


COMPOSIZIONE: π

• Possiamo immaginare un modello in cui ogni nucleotide è presente nelle stesse proporzioni.

• Oppure se vogliamo modellare sequenze che provengono da una isola CpG possiamo immaginare un modello in cui C e G hanno frequenza doppia rispetto ad A e T.

• In alternativa possiamo lasciare che i dati scelgano per noi (nel senso che utilizzeremo delle frequenze nucleotidiche ottenute dai dati che stiamo esaminando).


PROCESSO: P

• Questa parte del modello descrive le frequenze con cui un nucleotide muta in un altro … quindi è una matrice n x n (n = numero possibili nucleotidi). ad esempio:


PROCESSO: P

• NB: per convenzione sia le righe che le colonne della matrice corrispondono ai nucleotidi in ordine alfabetico (quindi: a,c,g,t) *

* Alla mutazione a c è quindi assegnata una probabilità pari a 0.01

P ac

Righe sommanoa 1


ESEMPIO 1 : likelihood di una sequenza di 1 nt

• Esempio semplice: 1 sola sequenza, 1 solo nt, nessun albero. La sequenza è: a

Osservazioni:

Non c’è cambiamento (abbiamo solo una sequenza, quindi non abbiamo bisogno della parte PROCESSO del modello). Ci serve solo la parte COMPOSIZIONE.



• Esempio semplice: 1 sola sequenza, 1 solo nt, nessun albero. La sequenza è: a

Se come composizione utilizziamo le seguenti frequenze π = [1, 0 , 0 , 0 ] allora la likelihood della sequenza “a” è 1. Anche nel caso del vettore delle frequenze l’ordine delle frequenze è, per convenzione, quello dei nucleotidi in ordine alfabetico. La somma dei valori deve essere 1.



• Esempio semplice: 1 sola sequenza, 2 nt, nessun albero. La sequenza è: ac

Se come composizione utilizziamo le frequenze nucleotidiche del modello di Jukes-Cantor ( π = [¼ , ¼ , ¼ , ¼ ] ) allora la likelihood della sequenza “ac” è:

πa x πc = ¼ x ¼ = 1/16



• Esempio semplice: 1 sola sequenza, 2 nt, nessun albero. La sequenza è: ac

Se come composizione utilizziamo le seguenti frequenze nucleotidiche, π = [0.4, 0.1 , 0.2 , 0.3 ] allora la likelihood della sequenza “ac” è:

πa x πc = 0.4 x 0.1 = 0.04

Se calcoliamo la likelihood di tutti i possibili dinucleotidi la somma deve essere uguale a 1. Indipendentemente dal contenuto di π


ESEMPIO 3: likelihood di un albero con un solo ramo

Vogliamo calcolare la likelihood di un albero formato da 1 solo ramo. Questo implica che abbiamo 2 sequenze:

c c a t

c c g t

Per calcolare likelihood ci servono tutte le parti del modello … sia π che P (P serve quando abbiamo più di una sequenza)


ESEMPIO 3: likelihood di un albero con un solo ramo

c c a t π = [0.1, 0.4 , 0.2 , 0.3]

c c g t

likelihood


ESEMPIO 3: Osservazioni

• Le probabilità associate alle colonne (composizione * processo) vengono moltiplicate … assunzione di indipendenza.

• In questo esempio non teniamo conto delle diverse lunghezze dei rami (se avessimo più rami il modello non sarebbe in grado di gestirli separatamente)

likelihood


ESEMPIO 3: Osservazioni

• Come è possibile modificare il modello in modo da ammettere l’esistenza di rami di lunghezza diversa?

• Quale parte del modello descrive i rami?• In cosa differiscono i rami di lunghezze diverse?

likelihood

Bio Massima verosimiglianzaLunghezza dei rami:

• Dipende dalla parte del modello che descrive il processo.

Questa matrice descrive un ramo con una “certa distanza evolutiva” … che non conosciamo. Immaginiamo che corrisponda ad una distanza pari a 1 cde.


• Un ramo di lunghezza 1 cde sembra essere un ramo abbastanza corto.

Valori sulla diagonale alti:Molto probabile che un nt non cambi

Valori fuori dalla diagonale bassi:Poco probabile che un nt muti in un altro …


• Un ramo di lunghezza 1 cde sembra essere un ramo abbastanza corto.

NB: man mano che la lunghezza del ramo cresce i valori nella matrice P diminuiscono lungo la diagonale ed aumentano al di fuori di essa.


Lunghezza dei rami:

• La likelihood calcolata in esempio 3 era per un ramo avente lunghezza pari a 1 unità cde … e se volessimo calcolare la likelihood per un ramo di 2 cde?

MOLTIPLICHIAMO LA MATRICE PER SE’ STESSA !


Lunghezza dei rami:

c c a t

c c g t

• La likelihood calcolata per questo allineamento (branch length = 1 cde) era 0.0000300, per 2 cde sarebbe 0.0000559 (è aumentata), per 3 cde sarebbe 0.000782.

La likelihood cresce indefinitamente?

Taxon A

Taxon B

A Bx ced


NO ! Esiste un valore massimo:Likelihood raggiunge un valore massimo in un punto compreso tra 10 e 20 cde (ced in EN)


Relazione tra π e P:

Se eleviamo la matrice P ad un esponente molto alto, otteniamo delle probabilità tendenti alle frequenze contenute in π !

Quindi π è già “codificato” nella matrice P che descrive il processo (evolutivo) . E’ come se le frequenze di sostituzione codificate in P, dopo un tempo evolutivo infinito, debbano convergere a π .


Matrici di velocità:Se vogliamo calcolare il valore di 54 possiamo calcolarlo come

e4*log(5). Possiamo operare nello stesso modo sulla matrice che rappresenta la parte del modello dedicata al processo:

P4 = e( 4 * log(P) )

Vantaggi:- Possiamo usare esponenti non interi.- Possiamo separare completamente le parti del modello dedicate

alla composizione ed al processo.- Possiamo esprimere lunghezza rami in sost. per sito

Inoltre possiamo usare come lunghezza dei rami qualsiasi numero da 0 a infinito


Matrici di velocità:Il logaritmo della matrice P dei nostri esempi è:

Le righe sommano a 0, la velocità corrisponde ad 1 cde ed e log P restituisce, di nuovo, la matrice P.


Matrici di velocità:

Questa matrice di velocità esprime una velocità di 1 cde … è già un passo avanti ma vorremmo una matrice M il cui esponenziale eM restituisce una matrice corrispondente ad 1 sostituzione per sito.


Matrici di velocità : scalare la matrice che descrive il processo ad una velocità di 1 sost. per sito

Possiamo ottenere questo risultato scalando log P in modo tale che, se moltiplichiamo le sue righe per πrow la SOMMA dei valori al di fuori della diagonale sia 1. In questo modo otteniamo la matrice il cui esponenziale corrisponde a rami da 1 sostituzione per sito.

In generale eQ(v) = P(v) per un ramo di lunghezza v sost. per sito.


Matrici di velocità :

Se scaliamo la matrice log P per un valore v=50 ( 50 sost. sito) otteniamo

Se moltiplichiamo Q per πdiag (matrice avente i valori di π sulla diagonale) otteniamo

Una matrice in cui i valori fuoridiagonale sommano a 1 (e quelli sulla diagonale a -1)


Matrici di velocità :Se moltiplichiamo Q per πdiag (matrice avente i valori di π

sulla diagonale, a volte indicata con Π ) otteniamo

L’esponenziale di questa matrice genera una matrice P utilizzabile per produrre un albero i cui rami hanno lunghezza espressa in sostituzioni per sito.

Una matrice in cui i valori fuoridiagonale sommano a 1


Separazione completa della composizione dalle velocità:

Se dividiamo le colonne di Q per πcol otteniamo la matrice delle velocità R , e separiamo la composizione dalle velocità. L’effetto è che possiamo utilizzare la stessa matrice R per diversi vettori di composizione. La matrice R per gli esempi visti finora è:


Separazione completa della composizione dalle velocità:

Rispetto alla matrice R (matrice velocità):

1. Gli elementi sulla diagonale non contano (trattasi do velocità di sost. e gli elementi sulla diagonale esprimono delle “non sostituzioni”).

2. Lo scaling di Q non ha effetto

3. Se vogliamo un modello reversibile la matrice R dovrebbe essere simmetrica.


Interconversione tra P, Q ed R:

NB: i programmi per analisi filogenetiche basati su maximum likelihood rendono le conversioni tra queste matrici completamente automatiche.


Massima verosimiglianza, lunghezze dei rami in sostituzioni per sito:

La verosimiglianza dell’allineamento di ccat e ccgt a diverse distanze è

Il valore massimo può essere trovato numericamente

mediante approssimazioni successive.

Si trova ad una lunghezza del ramo pari a 0.330614

(valore likelihood: 0.0001777).

Data una topologia è possibile trovare lelunghezze dei rami massimizzando la likelihood


Massima verosimiglianza: albero con 2 rami

Per la matrice Q delle slide precedenti le matrici P corrispondenti a 0.1, 0.2 e 0.3 sostituzioni per sito sono:

AO

B

0.1

0.2

origine



Ci sono 3 modi di calcolare la likelihood di quest’albero …

AO

B

0.1

0.2

origine



Modo 1: in un unico passo

AO

B0.3

origine

likelihood



Modo 2: in 2 passi … da A a O e poi da O a B

π = [0.1, 0.4 , 0.2 , 0.3]A

O

B

0.1

0.2

origine

Usiamo π perché partiamo da A !

PROBLEMA: non conosciamo la sequenza di O !

c c a t ? ? ? ?

CONSIDERIAMO 1 SOLO SITO:Le possibilità sonoc ac cc gc t

SOMMIAMO TUTTE LE PROBABILITA’




AO

B

0.1

0.2

originePROBLEMA: non conosciamo la sequenza di O !

c c a t ? ? ? ?

likelihood




AO

B

0.1

0.2

origineQuando aggiungiamo nel calcoloil secondo ramo (da O a B) NON serve includere π … ma solo le probabilità di arrivo a C partendo da qualsiasi nt.

c c a t ? ? ? ? c c g t

likelihood

Likelihood per 1 sito … se moltiplico likelihood dei

4 siti ottengo:0.000177 (come prima)



Modo 3: in 2 passi … da O a A + da O a B

AO

B

0.1

0.2

originePROBLEMA: non conosciamo la sequenza di O !

c c a t ? ? ? ? c c g t

likelihood

Likelihood tot. allineamento: 0.000177



3 Modi diversi: stesso valore di likelihood

AO

B

0.1

0.2

NB:

Non importa dove mettiamo la radice … il valore della likelihood

E’ LO STESSO !!!!!


Allineamento:

A c c a t

B c c g t

C g c a t

Albero:


A

B

0.1

0.2

C0.3O

Consideriamo come origine il nodo interno ed iniziamo da qui il calcolo dellalikelihood ( come in Modo 3 dell’esempio precedente)


Allineamento: Albero:

A c c a t

B c c g t

C g c a t


A

B

0.1

0.2

C0.3O

likelihood(primo sito)


Allineamento: Albero:

A c c a t

B c c g t

C g c a t


A

B

0.1

0.2

C0.3O

Dopo aver calcolato la likelihood per ognuno dei 4 siti, dato che consideriamo le colonne dell’allineamento indipendenti possiamo moltiplicare per ottenere la likelihood totale:

0.0204 * 0.245 * 0.00368 * 0.166 = 3.04 * 10-6

Fattori che complicano il problema:

1.La selezione agisce su parti diverse delle sequenze (pressione selettiva condivisa da tutti i taxa potrebbe riguardare solo una parte molto ristretta dell’allineamento multiplo)

2.Alcuni siti evolvono velocemente

3.Alcuni siti evolvono molto lentamente (alcuni siti poi non variano del tutto. Questo dipende dalle distanze evolutive tra i taxa e dal gene scelto)


Strumenti free per analisi ML:

PhyML 3.0 :

http://www.atgc-montpellier.fr/phyml/binaries.php

Possiamo interfacciarci a PhyML da R !

NB: per poter effettuare questo test dovete:

1)Scaricare PhyML

2)Posizionarvi nella directory contenente l’eseguibile di PhyML

3)Caricare le librerie R ape e seqinr

4)Utilizzare i comandi che troverete nelle prossime slides


WARNING!

Questo non è codice PERL ma è codice R.

> library(ape);library(seqinr)

> accnr <- paste("AJ5345",26:35,sep="")

> seq <- read.GenBank(accnr)

> names(seq) <- attr(seq, "species")

> dist <- dist.dna(seq, model = "K80")

> plot(nj(dist))


WARNING!


> setwd("/share/home/wim/bin")

> write.dna(seq,"seq.txt", format ="interleaved")

> out <-phymltest("seq.txt",format = "interleaved", execname ="phyml_linux")

> print(out)


WARNING!


> setwd("/share/home/wim/bin")

> write.dna(seq,"seq.txt", format ="interleaved")

> out <-phymltest("seq.txt",format = "interleaved", execname ="phyml_linux")

> print(out)

Tra tutti i modelli testati il migliore è il 27° (GTR+G)


WARNING!


Per stampare l’albero ottenuto (dal 27° modello)

> tr <- read.tree("seq.txt_phyml_tree.txt")

> plot(tr[[27]])

> add.scale.bar(length=0.01)


ATTENZIONE:Questo test è un po’ pesante … i risultati non arrivano in secondi.

(Altri package dedicati in R : phangorn)

Biologia computazionale

Documents

Transcript of Biologia computazionale