Least Squares para árvores filogenéticas.

Uma coisa que sempre achei complicado de entender, era como árvores filogenéticas eram feitas com verosimilhança.
Pensando em um teste t ou uma regressão linear, sabemos que temos um conjunto de amostras, e ajustamos um modelo a eles, ajustando os parâmetros do modelo de forma a minimizar o erro em todas as amostras.

Talvez a parte complexa é imaginar o que diabos é equivalente a amostra, para usar na construção de uma árvore. Então pensando no método dos mínimos quadrados, temos o seguinte.

Primeiro, imagine que temos quatro sequências qualquer, de quatro linhagens de alguma espécie. Mas vamos criá-las no formato fasta, que é um formato padrão para guardar sequências.

>A AAAAA >B CAAAA >C CGGGG >D GGGGG

Então temos as sequências A,B,C e D, e aqui, mesmo sem análise, da para chutar fácil quem é mais parecido com quem.

Mas tudo bem, o modelo para gerar uma árvore vai ser assim.
Veja que podemos calcular a uma distância entre sequências:

> ex.dna <- read.dna("exemplo.fasta", format = "fasta") > dist.dna(ex.dna,model="raw") A B C B 0.2 C 1.0 0.8 D 1.0 1.0 0.2

O pacote ape tem a função read.dna, que nos permite ler o DNA de várias formas, então como temos um arquivos fasta, lemos o dna no formato fasta, e calculamos a distância entre as sequências, usando um modelo “raw”, que e só o número de bases diferentes dividido pelo tamanho da sequência, esse nem é um modelo usado mesmo na biologia molecular, mas depois vemos isso.
A para B é 0.2 diferente, uma base em cinco, A e C é 1 diferente, porque todas as bases são diferentes, e assim vai.

Agora pense em uma topologia qualquer, uma organização numa árvore, com tamanhos nas arestas, com uma quantidade para a diferença. Vamos criar uma.

arvore1 <- "((A:1,B:1):1,(C:1,D:1):1);" arvore1.tre <- read.tree(text=arvore1) ###Como é árvore plot(arvore1.tre)

Veja que esse formato é relativamente simples, basicamente, uma virgula divide dois indivíduos, até um ancestral.

A:1,B:1

O dois pontos, é o tamanho desse cara até o ancestral comum, ai como temos que fazer mais uniões, colocamos ele em parenteses

(A:1,B:1)

Agora esse é como se fosse um indivíduo, é um ancestral de a com b, mas é um indivíduos, veja que unimos o C e o D

(C:1,D:1)

e depois unimos os ancestrais desses dois caras

((A:1,B:1):1,(C:1,D:1):1)

O ponto e vírgula é para indicar que a arvore terminou.
E essa árvore é a seguinte:

01

Beleza, então temos uma árvore e uma topologia, agora é o seguinte.

Essa topologia é um grafo, onde as folhas são as sequências e a arestas internas são ancestrais, só para esclarecer, na teoria dos grafos, uma árvore é um grafo não direcionado e em que quaisquer dois vértices são conectados por exatamente um caminho. Em outras palavras, qualquer gráfico acíclico conectado é uma árvore. Não existem ciclos.

Captura de tela de 2016-03-20 00-28-56

Agora essas sequências tem uma distância entre elas. Veja que cada aresta tem um valor, e a distância que aparece na matriz de distâncias da sequência tem que ser a soma dessas distâncias. Assim…

Captura de tela de 2016-03-20 00-29-15

Nesse exemplo, a distância filogenética entre a sequência 1 e 2 é o percurso de a->b->c, por isso a distância é a+b+c. Mas a não ser que essa seja uma filogenia perfeita, essa conta nunca fecha. Sempre falta ou sobra um pouquinho.

Captura de tela de 2016-03-20 00-29-39

E que pouquinho que sobra é exatamente o que permite calcular o mínimo quadrado, ou seja, quais valores a gente atribui para a,b,c,d,e que fazem com que tenhamos o mínimo possível de sobra? E veja que essa conta vai ser feita para todos as distâncias, todos os pares de distâncias, logo é minimizar o valor para todos os casos, como se fossem “amostras”.

E assim o mínimo quadrado é aplicado a filogenia.

No pacote phangorn, a gente tem esse método implementado já, e por enquanto vamos usar a função de la, mas pegando aquelas nossas sequências, podemos usar a função pml, que calcula o “Likelihood of a tree”. Então com a árvore e nossas sequências, assumindo um modelo de evolução calculamos um likelihhod, logo, assim como na parcimônia podemos decidir qual a melhor árvore, mas temos um valor de verosimilhança, que diz o quão melhor ela é.

Então pensando numa segunda árvore.

02

Veja que essa árvore é muito ruim, ela junta as sequências A e D, que são completamente diferentes, e usando a função pml vemos o seguinte:

03

Agora se a gente quer fazer uma topologia é so calcular todas as topologia possíveis, com todas as combinações de tamanho de arestas possíveis, e encontramos a melhor. Veja que o valor da melhor topologia é maior (menos negativo). Lembrando que estamos trabalhando com log aqui, mas podemos voltar para os número de probabilidade com exponencial

> format(exp(ajuste1$logLik),scientific = F)
[1] “0.000000000003288564”
> format(exp(ajuste2$logLik),scientific = F)
[1] “0.000000000000593519”

E é isso.

Bem é isso ai, tentando apenas fazer um link entre verossimilhança e árvores, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:
O curso Computational Molecular Evolution

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
rm(list=ls())
##install.packages("phangorn")
library(phangorn)
 
 
###Gerando um arquivo com sequencias
cat(">A",
    "AAAAA",
    ">B",
    "CAAAA",
    ">C",
    "CGGGG",
    ">D",
    "GGGGG",
file = "exemplo.fasta", sep = "\n")
 
###Lendo sequencias e transformando na classe phyDat
ex.dna <- read.dna("exemplo.fasta", format = "fasta")
dist.dna(ex.dna,model="raw")
 
ex.phyDat<-as.phyDat(ex.dna)
 
 
#####################################
## Topologia 1
#####################################
###Criando uma árvore
arvore1 <- "((A:1,B:1):1,(C:1,D:1):1);"
arvore1.tre <- read.tree(text=arvore1)
 
###Como é árvore
plot(arvore1.tre)
###Calculando o ajuste
ajuste1<-pml(arvore1.tre, data=ex.phyDat,model="JC")
ajuste1
 
#####################################
## Topologia 2
#####################################
###Agora uma topologia ruim
arvore2 <- "((C:1,B:1):1,(A:1,D:1):1);"
arvore2.tre <- read.tree(text=arvore2)
 
###Como é árvore
plot(arvore2.tre)
ajuste2<-pml(arvore2.tre, data=ex.phyDat,model="JC")
 
ex.dna <- read.dna("exemplo.fasta", format = "fasta",as.character=T,as.matrix=F)
 
arvore1.tre$tip.label<-sapply(ex.dna,paste,collapse = "")[match(arvore1.tre$tip.label,names(ex.dna))]
arvore2.tre$tip.label<-sapply(ex.dna,paste,collapse = "")[match(arvore2.tre$tip.label,names(ex.dna))]
 
par(mfrow=c(1,2),mar=c(0,0,2,0))
plot(arvore1.tre,main=paste("Verossimilhança:",round(ajuste1$logLik,2)),type="unrooted")
plot(arvore2.tre,main=paste("Verossimilhança:",round(ajuste2$logLik,2)),type="unrooted")
 
format(exp(ajuste2$logLik),scientific = F)