Biologia computazionale

1
Biologia computazionale
Università degli studi di milano
Docente Giorgio Valentini Istruttore
Matteo Re
C.d.l. Biotecnologie Industriali e Ambientali

• A.A. 2010-2011 semestre II

6
Evoluzione e filogenesi - 3
2
Metodi per costruire alberi filogenetici

• Metodi basati su
• Distanza
• Massima parsimonia
• Massima verosimiglianza

Questi li abbiamo visti
Oggi discutiamo questa classe di metodi
3
Massima verosimiglianza

• Verosimiglianza (likelihood)
• Probabilità delle osservazioni dato un modello
• Quindi è una probailità perché usare un nome
  diverso?
• Per porre laccento sul fatto che non vogliamo
  valutare quanto siamo confidenti nelloccorrenza
  di un determinato evento ma piuttosto valutare
  quanto i dati sono compatibili con un modello
  evolutivo che abbiamo scelto.

4
Massima verosimiglianza

• ESEMPIO lanciamo una moneta ed otteniamo
  croce (questo è il dato).
• Se dovessi chiedervi qualè la probabilità
  dellevento 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 ½ .

5
Massima verosimiglianza

• ESEMPIO lanciamo una moneta ed otteniamo
  croce (questo è il dato).
• Supponiamo di definire un modello di moneta con
  queste caratteristiche P(testa1), P(croce0)
  ossia una moneta truccata. I parametri (tutti)
  del modello li indichiamo complessivamente come
  T
• La likelihod dellosservazione 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 dellosservazione
  sarebbe uno

6
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 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?

7
Massima verosimiglianza

• Nel caso dellevoluzione molecolare i dati sono
  rappresentati da un allineamento di sequenze ed
  il modello, in senso molto ampio, è lalbero
  filogenetico che
• correla tra di loro le sequenze
• descrive il meccanismo di evoluzione da una
  sequenza allaltra

8
Massima verosimiglianza

• Lalbero 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
  levoluzione 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.

9
Massima verosimiglianza

• Lobiettivo 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
10
Massima verosimiglianza

• COMPOSIZIONE p
• 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).

11
Massima verosimiglianza

• 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

12
Massima verosimiglianza

• PROCESSO 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 a 1
Alla mutazione a ? c è quindi assegnata una
probabilità pari a 0.01
13
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.

14
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 p 1, 0 , 0 , 0 allora la
  likelihood della sequenza a è 1. Anche nel
  caso del vettore delle frequenze lordine delle
  frequenze è, per convenzione, quello dei
  nucleotidi in ordine alfabetico. La somma dei
  valori deve essere 1.

15
Massima verosimiglianza

• ESEMPIO 2 likelihood di una sequenza di 2 nt
• Esempio semplice 1 sola sequenza, 2 nt, nessun
  albero. La sequenza è ac
• Se come composizione utilizziamo le frequenze
  nucleotidiche del modello di Jukes-Cantor ( p
  ¼ , ¼ , ¼ , ¼ ) allora la likelihood della
  sequenza ac è
• pa x pc ¼ x ¼ 1/16

16
Massima verosimiglianza

• ESEMPIO 2 likelihood di una sequenza di 2 nt
• Esempio semplice 1 sola sequenza, 2 nt, nessun
  albero. La sequenza è ac
• Se come composizione utilizziamo le seguenti
  frequenze nucleotidiche, p 0.4, 0.1 , 0.2 ,
  0.3 allora la likelihood della sequenza ac
  è
• pa x pc 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 p
17
Massima verosimiglianza

• ESEMPIO 2 likelihood di una sequenza di 2 nt
• Esempio semplice 1 sola sequenza, 2 nt, nessun
  albero. La sequenza è ac
• Se come composizione utilizziamo le seguenti
  frequenze nucleotidiche, p 0.4, 0.1 , 0.2 ,
  0.3 allora la likelihood della sequenza ac
  è
• pa x pc 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 p
18
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
• c c a t
• c c g t
• Per calcolare likelihood ci servono tutte le
  parti del modello sia p che P (P serve quando
  abbiamo più di una sequenza)

19
Massima verosimiglianza

• ESEMPIO 3 likelihood di un albero con un solo
  ramo
• c c a t p
  0.1, 0.4 , 0.2 , 0.3
• c c g t

likelihood
20
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
21
Massima verosimiglianza

• ESEMPIO 3 Osservazioni
• Come è possibile modificare il modello in modo da
  ammettere lesistenza di rami di lunghezza
  diversa?
• Quale parte del modello descrive i rami?
• In cosa differiscono i rami di lunghezze diverse?

likelihood
22
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.
23
Massima verosimiglianza

• Lunghezza dei rami
• 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
24
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.
25
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 !
26
Massima verosimiglianza
Taxon A
A
B
x ced

• 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 B
27
Massima verosimiglianza
NO ! Esiste un valore massimo Likelihood
raggiunge un valore massimo in un punto compreso
tra 10 e 20 cde (ced in EN)
28
Massima verosimiglianza

• Relazione tra p e P
• Se eleviamo la matrice P ad un esponente molto
  alto, otteniamo delle probabilità tendenti alle
  frequenze contenute in p !

Quindi p è 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
p .
29
Massima verosimiglianza

• Matrici di velocità
• Se vogliamo calcolare il valore di 54 possiamo
  calcolarlo come e4log(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
30
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.

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

32
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 prow 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.

33
Massima verosimiglianza

• Matrici di velocità
• Se scaliamo la matrice log P per un valore v50 (
  50 sost. sito) otteniamo
• Se moltiplichiamo Q per pdiag (matrice avente i
  valori di p sulla diagonale) otteniamo

Una matrice in cui i valori fuori diagonale
sommano a 1 (e quelli sulla diagonale a -1)
34
Massima verosimiglianza

• Matrici di velocità
• Se moltiplichiamo Q per pdiag (matrice avente i
  valori di p sulla diagonale, a volte indicata con
  ? ) otteniamo
• Lesponenziale 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 fuori diagonale
sommano a 1
35
Massima verosimiglianza

• Separazione completa della composizione dalle
  velocità
• Se dividiamo le colonne di Q per pcol otteniamo
  la matrice delle velocità R , e separiamo la
  composizione dalle velocità. Leffetto è che
  possiamo utilizzare la stessa matrice R per
  diversi vettori di composizione. La matrice R per
  gli esempi visti finora è

36
Massima verosimiglianza

• Separazione completa della composizione dalle
  velocità
• Rispetto alla matrice R (matrice velocità)
• Gli elementi sulla diagonale non contano
  (trattasi do velocità di sost. e gli elementi
  sulla diagonale esprimono delle non
  sostituzioni).
• Lo scaling di Q non ha effetto
• Se vogliamo un modello reversibile la matrice R
  dovrebbe essere simmetrica.

37
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.
38
Massima verosimiglianza

• Massima verosimiglianza, lunghezze dei rami in
  sostituzioni per sito
• La verosimiglianza dellallineamento 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
39
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

A
0.1
origine
O
B
0.2
40
Massima verosimiglianza

• Massima verosimiglianza albero con 2 rami
• Ci sono 3 modi di calcolare la likelihood di
  questalbero

A
0.1
origine
O
B
0.2
41
Massima verosimiglianza

• Massima verosimiglianza albero con 2 rami
• Modo 1 in un unico passo

A
origine
O
B
0.3
likelihood
42
Massima verosimiglianza

• Massima verosimiglianza albero con 2 rami
• Modo 2 in 2 passi da A a O e poi da O a B
• p
  0.1, 0.4 , 0.2 , 0.3

Usiamo p perché partiamo da A !
A
0.1
origine
O
PROBLEMA non conosciamo la sequenza di O !
c c a t ? ? ? ?
B
0.2
CONSIDERIAMO 1 SOLO SITO Le possibilità sono c ?
a c ? c c ? g c ? t
SOMMIAMO TUTTE LE PROBABILITA
43
Massima verosimiglianza

• Massima verosimiglianza albero con 2 rami
• Modo 2 in 2 passi da A a O e poi da O a B

A
0.1
origine
PROBLEMA non conosciamo la sequenza di O !
c c a t ? ? ? ?
O
B
0.2
likelihood
44
Massima verosimiglianza

• Massima verosimiglianza albero con 2 rami
• Modo 2 in 2 passi da A a O e poi da O a B

A
0.1
origine
Quando aggiungiamo nel calcolo il secondo ramo
(da O a B) NON serve includere p ma solo le
probabilità di arrivo a C partendo da qualsiasi
nt. c c a t ? ? ? ?
c c g t
O
B
0.2
likelihood
Likelihood per 1 sito se moltiplico likelihood
dei 4 siti ottengo 0.000177 (come prima)
45
Massima verosimiglianza

• Massima verosimiglianza albero con 2 rami
• Modo 3 in 2 passi da O a A da O a B

A
0.1
origine
PROBLEMA non conosciamo la sequenza di O !
c c a t ? ? ? ? c c g t
O
B
0.2
likelihood
Likelihood tot. allineamento 0.000177
46
Massima verosimiglianza

• Massima verosimiglianza albero con 2 rami
• 3 Modi diversi stesso valore di likelihood

A
0.1
O
B
0.2
NB Non importa dove mettiamo la radice il
valore della likelihood E LO STESSO !!!!!
47
Massima verosimiglianza

• Massima verosimiglianza albero con 3 rami
• Allineamento
• A c c a
  t
• B c c g
  t
• C g c a
  t
• Albero

A
0.1
C
0.3
O
0.2
B
Consideriamo come origine il nodo interno ed
iniziamo da qui il calcolo della likelihood (
come in Modo 3 dellesempio precedente)
48
Massima verosimiglianza

• Massima verosimiglianza albero con 3 rami
• Allineamento
  Albero
• A c c a t
• B c c g t
• C g c a t

A
0.1
C
0.3
O
0.2
B
likelihood (primo sito)
49
Massima verosimiglianza

• Massima verosimiglianza albero con 3 rami
• Allineamento
  Albero
• A c c a t
• B c c g t
• C g c a t

A
0.1
C
0.3
O
0.2
B
Dopo aver calcolato la likelihood per ognuno dei
4 siti, dato che consideriamo le colonne
dellallineamento indipendenti possiamo
moltiplicare per ottenere la likelihood totale
0.0204 0.245 0.00368 0.166 3.04
10-6
50
Massima verosimiglianza

• Fattori che complicano il problema
• La selezione agisce su parti diverse delle
  sequenze (pressione selettiva condivisa da tutti
  i taxa potrebbe riguardare solo una parte molto
  ristretta dellallineamento multiplo)
• Alcuni siti evolvono velocemente
• Alcuni siti evolvono molto lentamente (alcuni
  siti poi non variano del tutto. Questo dipende
  dalle distanze evolutive tra i taxa e dal gene
  scelto)

51
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
• Scaricare PhyML
• Posizionarvi nella directory contenente
  leseguibile di PhyML
• Caricare le librerie R ape e seqinr
• Utilizzare i comandi che troverete nelle prossime
  slides

52
Massima verosimiglianza

• WARNING!
• Questo non è codice PERL ma è codice R.
• gt library(ape)library(seqinr)
• gt accnr lt- paste("AJ5345",2635,sep"")
• gt seq lt- read.GenBank(accnr)
• gt names(seq) lt- attr(seq, "species")
• gt dist lt- dist.dna(seq, model "K80")
• gt plot(nj(dist))

53
Massima verosimiglianza

• WARNING!
• Questo non è codice PERL ma è codice R.
• gt setwd("/share/home/wim/bin")
• gt write.dna(seq,"seq.txt", format "interleaved")
• gt out lt-phymltest("seq.txt",format
  "interleaved", execname "phyml_linux")
• gt print(out)

54
Massima verosimiglianza

• WARNING!
• Questo non è codice PERL ma è codice R.
• gt setwd("/share/home/wim/bin")
• gt write.dna(seq,"seq.txt", format "interleaved")
• gt out lt-phymltest("seq.txt",format
  "interleaved", execname "phyml_linux")
• gt print(out)
• Tra tutti i modelli testati il migliore è il 27
  (GTRG)

55
Massima verosimiglianza

• WARNING!
• Questo non è codice PERL ma è codice R.
• Per stampare lalbero ottenuto (dal 27 modello)
• gt tr lt- read.tree("seq.txt_phyml_tree.txt")
• gt plot(tr27)
• gt add.scale.bar(length0.01)

ATTENZIONE Questo test è un po pesante i
risultati non arrivano in secondi.
(Altri package dedicati in R phangorn)