Biologia_computazion.. - Università degli Studi di Milano

UNIVERSITÀ
Docente: Giorgio Valentini
Istruttore: Matteo Re
C.d.l.
DEGLI
STUDI DI MILANO
Biotecnologie Industriali e Ambientali
Biologia
computazionale
A.A. 2010-2011 semestre II
6
Evoluzione e filogenesi - 3
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.
Bio
Massima verosimiglianza
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 ½ .
Bio
Massima verosimiglianza
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
Bio
Massima verosimiglianza
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
filogenetico mediante la tecnica della
verosimiglianza (maximum likelihood)
bisogno innanzitutto di un modello
adatto alle sequenze biologiche.
un albero
massima
abbiamo
evolutivo
Ma come possiamo costruire un tale modello?
Bio
Massima verosimiglianza
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
sequenza all’altra
evoluzione
da
una
Bio
Massima verosimiglianza
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”.
Bio
Massima verosimiglianza
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
Bio
Massima verosimiglianza
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).
Bio
PROCESSO:
Massima verosimiglianza
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:
Bio
PROCESSO:
Massima verosimiglianza
P
• NB: per convenzione sia le righe che le colonne
della matrice corrispondono ai nucleotidi in ordine
alfabetico (quindi: a,c,g,t) *
P ac
Righe sommano
a1
* Alla mutazione a  c è quindi assegnata una probabilità pari a 0.01
Bio
Massima verosimiglianza
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.
Bio
Massima verosimiglianza
ESEMPIO 1 :
likelihood di una sequenza di 1 nt
• 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.
Bio
Massima verosimiglianza
ESEMPIO 2 :
likelihood di una sequenza di 2 nt
• Esempio semplice: 1 sola sequenza,
albero. La sequenza è: ac
2 nt, nessun
Se come composizione utilizziamo le frequenze
nucleotidiche del modello di Jukes-Cantor ( π = [¼ ,
¼ , ¼ , ¼ ] ) allora la likelihood della sequenza “ac”
è:
πa x πc = ¼ x ¼ = 1/16
Bio
Massima verosimiglianza
ESEMPIO 2 :
likelihood di una sequenza di 2 nt
• Esempio semplice: 1 sola sequenza,
albero. La sequenza è: ac
2 nt, nessun
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 π
Bio
Massima verosimiglianza
ESEMPIO 2 :
likelihood di una sequenza di 2 nt
• Esempio semplice: 1 sola sequenza,
albero. La sequenza è: ac
2 nt, nessun
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 π
Bio
Massima verosimiglianza
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:
ccat
ccgt
Per calcolare likelihood ci servono tutte le parti del
modello … sia π che P (P serve quando abbiamo
più di una sequenza)
Bio
Massima verosimiglianza
ESEMPIO 3: likelihood di un albero con un solo ramo
π = [0.1, 0.4 , 0.2 , 0.3]
ccat
ccgt
likelihood
Bio
Massima verosimiglianza
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
Bio
Massima verosimiglianza
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 verosimiglianza
Lunghezza 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.
Bio
Massima verosimiglianza
Lunghezza dei rami:
• Un ramo di lunghezza 1 cde sembra essere un
ramo abbastanza corto.
Valori fuori dalla diagonale bassi:
Poco probabile che un nt muti
in un altro …
Valori sulla diagonale alti:
Molto probabile che un nt non
cambi
Bio
Massima verosimiglianza
Lunghezza dei rami:
• 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.
Bio
Massima verosimiglianza
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 !
Bio
Massima verosimiglianza
Lunghezza dei rami:
Taxon A
ccat
ccgt
A
x ced
B
Taxon B
• 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?
Bio
Massima verosimiglianza
NO ! Esiste un valore massimo:
Likelihood raggiunge un valore massimo in un
punto compreso tra 10 e 20 cde (ced in EN)
Bio
Massima verosimiglianza
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 π .
Bio
Massima verosimiglianza
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
Bio
Massima verosimiglianza
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.
Bio
Massima verosimiglianza
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.
Bio
Massima verosimiglianza
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.
Bio
Massima verosimiglianza
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 fuori
diagonale sommano a 1 (e
quelli sulla diagonale a -1)
Bio
Massima verosimiglianza
Matrici di velocità :
Se moltiplichiamo Q per πdiag (matrice avente i valori di π
sulla diagonale, a volte indicata con Π ) otteniamo
Una matrice in cui i valori fuori
diagonale sommano a 1
L’esponenziale di questa matrice genera una matrice P
utilizzabile per produrre un albero i cui rami hanno
lunghezza espressa in sostituzioni per sito.
Bio
Massima verosimiglianza
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 è:
Bio
Massima verosimiglianza
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.
Bio
Massima verosimiglianza
Interconversione tra P, Q ed R:
NB: i programmi per analisi filogenetiche basati su maximum likelihood
rendono le conversioni tra queste matrici completamente automatiche.
Bio
Massima verosimiglianza
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 le
lunghezze dei rami massimizzando la likelihood
Bio
Massima verosimiglianza
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:
0.1 A
origine
O
0.2
B
Bio
Massima verosimiglianza
Massima verosimiglianza: albero con 2 rami
Ci sono 3 modi di calcolare la likelihood di
quest’albero …
0.1 A
origine
O
0.2
B
Bio
Massima verosimiglianza
Massima verosimiglianza: albero con 2 rami
Modo 1: in un unico passo
A
origine
O
0.3
likelihood
B
Bio
Massima verosimiglianza
Massima verosimiglianza: albero con 2 rami
Modo 2: in 2 passi … da A a O e poi da O a B
Usiamo π perché partiamo da A !
π = [0.1, 0.4 , 0.2 , 0.3]
0.1 A
PROBLEMA:
non conosciamo la sequenza di O !
ccat
????
origine
O
0.2
B
CONSIDERIAMO 1 SOLO SITO:
Le possibilità sono
ca
cc
SOMMIAMO TUTTE LE
cg
PROBABILITA’
ct
Bio
Massima verosimiglianza
Massima verosimiglianza: albero con 2 rami
Modo 2: in 2 passi … da A a O e poi da O a B
PROBLEMA:
non conosciamo la sequenza di O !
0.1 A
origine
O
0.2
ccat
????
likelihood
B
Bio
Massima verosimiglianza
Massima verosimiglianza: albero con 2 rami
Modo 2: in 2 passi … da A a O e poi da O a B
Quando aggiungiamo nel calcolo
il secondo ramo (da O a B) NON
serve includere π … ma solo le
probabilità di arrivo a C partendo
da qualsiasi nt.
0.1 A
origine
O
0.2
B
likelihood
ccat
????
ccgt
Likelihood per 1 sito … se
moltiplico likelihood dei
4 siti ottengo:
0.000177 (come prima)
Bio
Massima verosimiglianza
Massima verosimiglianza: albero con 2 rami
Modo 3: in 2 passi … da O a A + da O a B
PROBLEMA:
non conosciamo la sequenza di O !
0.1 A
origine
O
0.2
B
ccat
????
ccgt
likelihood
Likelihood tot. allineamento: 0.000177
Bio
Massima verosimiglianza
Massima verosimiglianza: albero con 2 rami
3 Modi diversi: stesso valore di likelihood
0.1 A
O
0.2
NB:
Non importa dove mettiamo la radice
… il valore della likelihood
E’ LO STESSO !!!!!
B
Bio
Massima verosimiglianza
Massima verosimiglianza: albero con 3 rami
Allineamento:
A
ccat
B
ccgt
C
gcat
A
0.1
Albero:
O
B
0.3
C
0.2
Consideriamo come origine il nodo interno ed iniziamo da qui il calcolo della
likelihood ( come in Modo 3 dell’esempio precedente)
Bio
Massima verosimiglianza
Massima verosimiglianza: albero con 3 rami
Allineamento:
Albero:
A
0.1
A
ccat
O
0.3
B
ccgt
0.2
B
C
gcat
likelihood
(primo sito)
C
Bio
Massima verosimiglianza
Massima verosimiglianza: albero con 3 rami
Allineamento:
Albero:
A
0.1
A
ccat
O
0.3
B
ccgt
0.2
B
C
gcat
C
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
Bio
Massima verosimiglianza
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)
Bio
Massima verosimiglianza
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
Bio
Massima verosimiglianza
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))
Bio
Massima verosimiglianza
WARNING!
Questo non è codice PERL ma è codice R.
> setwd("/share/home/wim/bin")
> write.dna(seq,"seq.txt", format
="interleaved")
> out <-phymltest("seq.txt",format =
"interleaved", execname ="phyml_linux")
> print(out)
Bio
Massima verosimiglianza
WARNING!
Questo non è codice PERL ma è codice R.
> 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)
Bio
Massima verosimiglianza
WARNING!
Questo non è codice PERL ma è codice R.
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)