Precisaremos de dois pacotes específicos para trabalhar com filogenias. A instalação destes pacotes é bem simples. (Não esqueça de retirar o “#” para realizar o comando)
#install.packages("ape")
#install.packages("geiger")
Você só precisa instalá-los uma vez. Depois basta usar os seguintes comandos:
library(ape)
library(geiger)
A primeira coisa que faremos é simular uma filogenia.
tr<-sim.bdtree(0.1, 0, stop="taxa", n=50)
plot(tr)
O primeiro valor “0.1” corresponde à taxa de especiação, o segundo valor “0” corresponde à taxa de extinção, e taxa.stop diz com quantas espécies a simulação deve terminar. Podemos também simular uma filogenia em que aconteça extinção a uma taxa constante também.
tr<-sim.bdtree(0.1, 0.02, stop="taxa", n=50)
plot(tr)
Vc notará que várias linhagens extintas ainda estão lá. Podemos “podá-las” assim:
tr.podada<-drop.extinct(tr)
plot(tr.podada)
Agora é um momento útil pra aprendermos mais um truque do R. Se quisermos plotar ambas as figuras, uma ao lado da outra, podemos usar os seguintes comandos
layout(matrix(c(1:2), nrow=1, ncol=2))
plot(tr)
plot(tr.podada)
Uma maneira simples de investigarmos a diversificação de linhagens ao longo do tempo é o que chamamos de diagramas de linhagens no tempo (“lineage-through-time plots”). Nestes gráficos, o número cumulativo de linhagens é acompanhado desde a raiz da árvore até o presente, mostrando isso na forma de uma curva.
ltt.plot(tr.podada)
Como sempre estaremos trabalhando com filogenias de espécies atuais, este gráfico sempre vai crescer e é difícil ter uma idéia mais precisa de como as linhagens se diversificaram. Uma forma simples de fazer isso é alterar o eixo Y
ltt.plot(tr.podada, log="y")
Note que, se a taxa de especiação é constante, o log número de espécies aumenta linearmente com o tempo. Isso é uma boa notícia, porque representa uma expectativa simples se essa é a maneira que as espécies de um clado evoluem.
Podemos usar o poder do R e simular várias filogenias ao mesmo tempo:
trees<-list()
for (i in 1:100) {
trees[[i]]<-sim.bdtree(0.1, 0, stop="taxa", n=50)
}
Esse comando cria uma lista e coloca uma árvore em cada lugar da lista. Pra plotarmos uma dessas filogenias, use:
plot(trees[[1]])
Pra termos uma idéia na variabilidade intrínseca neste tipo de simulações, vamos plotar várias filogenias, uma ao lado da outra. Pra ficar mais limpo, vamos tirar os nomes de cada espécie.
layout(matrix(c(1:6), nrow=2, ncol=3))
plot(trees[[1]],show.tip.label=F)
plot(trees[[2]],show.tip.label=F)
plot(trees[[3]],show.tip.label=F)
plot(trees[[4]],show.tip.label=F)
plot(trees[[5]],show.tip.label=F)
plot(trees[[6]],show.tip.label=F)
Isso plotará as primeiras 6 filogenias. Lembre da nossa aula sobre o acaso e note como, mesmo que a taxa de especiação é constante e a mesma pra todas as filogenias, elas frequentemente são bem diferentes umas das outras. Podemos também fazer o LTT de todas as filogenias ao mesmo tempo.
mltt.plot(trees,legend =F)
Como antes, fica mais fácil de visualizar os resultados com o eixo y logaritmizado
mltt.plot(trees,log="y", legend =F)
Note como há muita variabilidade nos dados, mesmo que o modelo de evolução é exatamente o mesmo. Finalmente, vamos explorar a possibilidade de não termos na nossa filogenia todas as espécies. Como o nosso código já está começando a ficar um pouco mais complicado, vamos aprender um pouco mais sobre como organizá-lo. Uma coisa útil é ter um “cabeçalho” onde colocamos nossas opções básicas. Assim, fica muito mais fácil de alterar as opções pra gerar outros resultados.
S <- 200 # número real de espécies
S_missing <- 150 # número de espécies que faltam
nreps <- 100 # número de replicatas
b <- 0.1 # taxa de especiação
# gerando as filogenias completas
v <- list();
for (i in 1:nreps)
v[[i]] <- sim.bdtree(0.1, 0, stop="taxa", n=S)
Agora podamos cada uma das filogenias de acordo com o que especificamos no cabeçalho. Pra mais informações, use “?drop.tip” no console do R pra entender como funciona esse comando.
nomes<-v[[1]]$tip.label
podadas <- lapply(v, drop.tip, sample(nomes, S_missing))
Finalmente, vamos podar os dois LTTs.
layout(matrix(c(1:2), nrow=2))
mltt.plot(v, legend =F)
mltt.plot(podadas, legend =F)