Análises de diversificação utilizando BAMM e RPANDA
library("BAMMtools")
library("RPANDA")

BAMM

Baixando e instalando o BAMM

BAMM (Bayesian Analysis of Macroevolutionary Mixtures) é um software que utiliza rjMCMC para modelar dinâmicas complexas de especiação, extinção e evolução fenotípica em filogenias. No entanto, o tutorial a seguir irá focar somente nas estimativas de taxas de especiação e extinção.

O primeiro passo para usar o BAMM é baixar o arquivo executável disponível em website do programa. Para instalar o programa você deve seguir as instruções específicas ao seu sistema operacional no site.

O jeito ideal de trabalhar com uma análise do BAMM é criar um folder para cada run/projeto/análise. Dessa maneira, você terá todos os arquivos necessários no local certo e manterá seus resultados organizados.

Existem três arquivos principais para rodar o BAMM:

  • BAMM executable file

  • Control file

  • Tree file (no formato newick)

Filogenia de baleias

Usaremos uma filogenia de baleias como exemplo. Você pode obter a árvore nesse link, e o control file (“divcontrol.txt”) pode ser obtido neste link. Na pasta que você criou para a sua análise, salve os arquivos principais mencionados anteriormente: arquivo executável do BAMM, control file e tree file.

O próximo passo é editar o control file.

Control file

O control file contém todas as informações que você precisará para rodar o BAMM, incluindo o nome dos arquivos de entrada que você quer analisar. O control file possui muitas diferentes opções, algumas previamente definidas, e todas as opções tem comentários no arquivo descrevendo sua função. Você pode obter mais detalhes sobre cada um no site do BAMM.

Rodando o BAMM

Em seu terminal, navegue para a pasta contendo o arquivo executável do BAMM, o control file e a filogenia.

./bamm -c divcontrol.txt ## copie e cole no terminal

Pacote BAMMtools

No R:

install.packages("BAMMtools")
library(BAMMtools)

tree <- read.tree("whaletree.tre")
plot(tree, cex = 0.35);axisPhylo()

Leia o arquivo “eventdata.txt”, que contém a informação que o BAMMtools precisa para analisar os resultados do BAMM.

events <- read.csv("event_data.txt")

No entanto, para explorar algumas das funções do BAMMtools, vamos usar um eventfile resultante de uma análise do BAMM que rodou por mais tempo (com um “chain length mais longo”).

data(whales, events.whales)

Aqui iremos utilizar a função “getEventData” para gerar um objeto da classe “bammdata”. O objeto “bammdata” é uma estrutura de dados complexa que inclui uma filogenia e o mapeamento de todos os parâmetros macroevolutivos amostrados pelo BAMM. Muitos dos métodos no BAMMtools operam diretamente em objetos da classe “bammdata”.

ed <- getEventData(whales, events.whales, burnin=0.25)
head(ed$eventData)

bamm.whales <- plot.bammdata(ed, lwd=2, labels = T, cex = 0.5)
addBAMMshifts(ed, cex = 2)
addBAMMlegend(bamm.whales, location = 'topleft')

Rate through time plots

Também podemos plotar as taxas ao longo do tempo estimadas pelo BAMM

plot.new()
plotRateThroughTime(ed, intervalCol = "red", avgCol = "red", ylim = c(0, 1), cex.axis = 2)
text(x = 30, y = 0.8, label = "All whales", font = 4, cex = 2.0, pos = 4)

Ou plotar as taxas ao longo do tempo estimadas para um dado clado:

plotRateThroughTime(ed, intervalCol = "blue", avgCol = "blue", node = 141, ylim = c(0, 1), cex.axis = 1.5)
text(x = 30, y = 0.8, label = "Dolphins only", font = 4, cex = 2.0, pos = 4)

Ou excluir esse dado clado e, então, plotar as taxas ao longo do tempo para o clado de fundo (“background clade”).

plotRateThroughTime(ed, intervalCol = "darkgreen", avgCol = "darkgreen", node = 141, nodetype = "exclude", ylim = c(0, 1), cex.axis = 1.5)
text(x = 30, y = 0.8, label = "Non-dolphins", font = 4, cex = 2.0, pos = 4)

Extraindo as taxas dos tips

O Bamm estima as taxas de especiação e extinção para cada ramo da filogenia. Aqui iremos extrair as taxas estimadas para cada tip e plotar de diferentes formas.

tip.rates <- getTipRates(ed)
str(tip.rates)

Você pode explorar as taxas dos tips através de um histograma:

hist(tip.rates$lambda.avg, xlab = "Lambda Médio")
hist(tip.rates$mu.avg, xlab = "Mu Médio")

Ou você pode separar diferentes grupos para comparar as taxas dos tips:

boxplot(tip.rates$lambda.avg[53:87], tip.rates$lambda.avg[1:52], col = c("red", "blue"), names = c("dolphins", "other cetaceans"))

RPANDA

RPANDA é outro pacote do R recentemente desenvolvido que nos permite ajustar diferentes modelos de variação das taxas ao longo do tempo e selecionar o melhor modelo através de análise de maximum likelihood.

O RPANDA possui duas diferenças básicas em relação ao BAMM:

  1. No RPANDA, o usuário deve informar quais modelos serão testados, ao passo que no BAMM, o programa irá calcular uma média ponderada das taxas para cada regime de diversificação;

  2. No RPANDA, o usuário deve saber a priori qual são os clados que devem ter regimes de diversificação específicos, ao passo que no BAMM, mudanças no regime de diversificação (“rate shifts”) são posicionados e os diferentes regimes estimados através de um algoritmo rjMCMC.

O pacote também contém algumas funções de simulação bem como conjuntos de dados. Você pode encontrar mais detalhes em CRAN.

Seleção de modelos

Para este exercício, iremos testar quatro diferentes cenários de diversificação, com todas as combinações de taxas de especiação e extinção constantes e/ou variáveis (baseado em Morlon et al. 2011, PNAS):

lambda.cst <- function(x,y){y}
lambda.var <- function(x,y){y[1] * exp(y[2] * x)}
mu.cst <- function(x,y){y}
mu.var <- function(x,y){y[1] * exp(y[2] * x)}

As variações das taxas para as quatro possíveis combinações se parecem com isso:

par(mfrow = c(2,2))
plot(rep(1, length(seq(0, 25, 0.1))) ~ seq(0, 25, 0.1), type = "l", ylim = c(0, 1.2), xlab = "Tempo", ylab = "Taxa", main = "BOTHcst", lwd = 2)
legend("topleft", legend = c("Especiação", "Extinção"), col = c(1, 2), lty = 1, cex = 0.4, lwd = 2)
abline(h = 0.2, col = 2, lwd = 2)
plot(lambda.var(seq(0, 25, 0.1), c(1, -0.1)) ~ seq(0, 25, 0.1), type = "l", xlab = "Tempo", ylab = "Taxa", main = "SPvar", lwd = 2)
legend("topright", legend = c("Especiação","Extinção"), col = c(1, 2), lty = 1, cex = 0.4, lwd = 2)
abline(h = 0.2, col = 2, lwd = 2)
plot(mu.var(seq(0, 25, 0.1), c(0.2, 0.075)) ~ seq(0, 25, 0.1), type = "l", col = 2, xlab = "Tempo", ylab = "Taxa", main = "EXvar", lwd = 2)
legend("topleft", legend = c("Especiação", "Extinção"), col = c(1, 2), lty = 1, cex = 0.4, lwd = 2)
abline(h = 1, lwd = 2)
plot(lambda.var(seq(0, 25, 0.1),c(1, -0.1)) ~ seq(0, 25, 0.1), type = "l", xlab = "Tempo", ylab = "Taxa", main = "BOTHvar", lwd = 2)
legend("topright", legend = c("Especiação", "Extinção"), col = c(1,2 ), lty = 1, cex = 0.4, lwd = 2)
lines(mu.var(seq(0, 25, 0.1), c(0.1, 0.075)) ~ seq(0, 25, 0.1), type = "l", col = 2, lwd = 2)

Tendo decidido quais modelos gostariamos de testar, o próximo passo é extrair clados cujas espécies devem compartilhar um mesmo regime de diversificação. Como dito anteriormente, o RPANDA precisa que o usuário indique a priori quais são esses clados. Neste exemplo, usaremos as quatro principais famílias como clados e o restante dos cetaceos como um quinto clado.

balaenidae <- c(whales$tip.label[grep("Balaena", whales$tip.label)], whales$tip.label[grep("Eubalaena", whales$tip.label)])
balaenopteridae <- c(whales$tip.label[grep("Balaenoptera", whales$tip.label)], whales$tip.label[grep("Megaptera", whales$tip.label)])
delphinidae <- c(whales$tip.label[grep("Delphinus", whales$tip.label)], whales$tip.label[grep("Cephalorhynchus", whales$tip.label)], whales$tip.label[grep("Feresa", whales$tip.label)], whales$tip.label[grep("Globicephala", whales$tip.label)], whales$tip.label[grep("Lagenodelphis", whales$tip.label)], whales$tip.label[grep("Lagenorhynchus", whales$tip.label)], whales$tip.label[grep("Lissodelphis", whales$tip.label)], whales$tip.label[grep("Orcaella", whales$tip.label)], whales$tip.label[grep("Orcinus", whales$tip.label)], whales$tip.label[grep("Peponocephala", whales$tip.label)], whales$tip.label[grep("Pseudorca", whales$tip.label)], whales$tip.label[grep("Sotalia", whales$tip.label)], whales$tip.label[grep("Sousa", whales$tip.label)], whales$tip.label[grep("Stenella", whales$tip.label)], whales$tip.label[grep("Steno", whales$tip.label)], whales$tip.label[grep("Tursiops", whales$tip.label)], whales$tip.label[grep("Grampus", whales$tip.label)])
phocoenidae <- c(whales$tip.label[grep("Neophocaena", whales$tip.label)], whales$tip.label[grep("Phocoena", whales$tip.label)], whales$tip.label[grep("Phocoenoides", whales$tip.label)])
ziphidae <- c(whales$tip.label[grep("Berardius", whales$tip.label)], whales$tip.label[grep("Hyperoodon", whales$tip.label)], whales$tip.label[grep("Indopacetus", whales$tip.label)], whales$tip.label[grep("Mesoplodon", whales$tip.label)], whales$tip.label[grep("Tasmacetus", whales$tip.label)], whales$tip.label[grep("Ziphius", whales$tip.label)])
balaenidae.tree <- drop.tip(whales, whales$tip.label[-match(balaenidae, whales$tip.label)])
balaenopteridae.tree <- drop.tip(whales, whales$tip.label[-match(balaenopteridae, whales$tip.label)])
delphinidae.tree <- drop.tip(whales, whales$tip.label[-match(delphinidae, whales$tip.label)])
phocoenidae.tree <- drop.tip(whales, whales$tip.label[-match(phocoenidae, whales$tip.label)])
ziphidae.tree <- drop.tip(whales, whales$tip.label[-match(ziphidae, whales$tip.label)])
othercetaceans.tree <- drop.tip(whales, c(balaenopteridae, delphinidae,  phocoenidae,  ziphidae))
fit.multi.rpanda <- function(tree, par)
    {
        bcstdcst <- fit_bd(tree,  max(branching.times(tree)),  f.lamb = lambda.cst,  f.mu = mu.cst,  lamb_par = par[[1]][1], mu_par = par[[1]][2], cst.lamb = TRUE, cst.mu = TRUE, cond = "crown", f = 87/89, dt = 1e-3)
        bvardcst <- fit_bd(tree,  max(branching.times(tree)),  f.lamb = lambda.var,  f.mu = mu.cst,  lamb_par = par[[2]][c(1, 2)], mu_par = par[[2]][3], expo.lamb = TRUE, cst.mu = TRUE, cond = "crown", f = 87/89, dt = 1e-3)
        bcstdvar <- fit_bd(tree,  max(branching.times(tree)),  f.lamb = lambda.cst,  f.mu = mu.var,  lamb_par = par[[3]][1], mu_par = par[[3]][c(2, 3)], cst.lamb = TRUE, expo.mu = TRUE, cond = "crown", f = 87/89, dt = 1e-3)
        bvardvar <- fit_bd(tree,  max(branching.times(tree)),  f.lamb = lambda.var,  f.mu = mu.var,  lamb_par = par[[4]][c(1, 2)], mu_par = par[[4]][c(3, 4)], expo.lamb = TRUE, expo.mu = TRUE, cond = "crown", f = 87/89, dt = 1e-3)
        return(list("bcstdcst" = bcstdcst, "bvardcst" = bvardcst, "bcstdvar" = bcstdvar, "bvardvar" = bvardvar))
    }
whales.par <- list(c(0.4, 0), c(0.4, -0.05, 0), c(0.4, 0.1, 0.05), c(0.4, -0.05, 0.1, 0.05))

Estimativa dos parâmetros dos modelos

Tendo todos os modelos e clados, podemos finalmente estimar os parâmetros de todos os quatro modelos para cada um dos cinco clados e criar uma tabela de AICc para selecionarmos qual deles é o melhor modelo que descreve as mudanças nas taxas de diversificação (especiação e extinção). (esta parte do código leva um tempo para rodar).

results <- list("balaenopteridae.res" = fit.multi.rpanda(balaenopteridae.tree, whales.par), 
                "delphinidae.res" = fit.multi.rpanda(delphinidae.tree, whales.par), 
                "phocoenidae.res" = fit.multi.rpanda(phocoenidae.tree, whales.par), 
                "ziphidae.res" = fit.multi.rpanda(ziphidae.tree, whales.par), 
                "othercetaceans.res" = fit.multi.rpanda(othercetaceans.tree, whales.par))
aic.table <- matrix(nrow = 4, ncol = 5, NA)
for(i in 1:5)
    {
        for(j in 1:4)
            {
                aic.table[j, i] <- results[[i]][[j]]$aicc
            }
    }
colnames(aic.table) <- c("Balaenopteridae", "Delphinidae", "Phocoenidae", "Ziphidae", "Other Cetaceans")
rownames(aic.table) <- c("bcstdcst", "bvardcst", "bcstdvar", "bvardvar")
aic.table
par.table <- data.frame("Balaenopteridae" = c(results[[1]]$bcstdcst$lamb_par[1:2], results[[1]]$bcstdcst$mu_par[1:2]), "Delphinidae" = c(results[[2]]$bvardcst$lamb_par[1:2], results[[2]]$bvardcst$mu_par[1:2]), "Phocoenidae" = c(results[[3]]$bcstdcst$lamb_par[1:2], results[[3]]$bcstdcst$mu_par[1:2]), "Ziphidae" = c(results[[4]]$bcstdcst$lamb_par[1:2], results[[4]]$bcstdcst$mu_par[1:2]), "Other Cetaceans" = c(results[[5]]$bcstdvar$lamb_par[1:2], results[[5]]$bcstdvar$mu_par[1:2]))
par.table

Plotando a diversidade no tempo (diversity through time)

Depois de selecionar qual modelo melhor ajusta às filogenias, podemos estimar a trajetória de diversidade para cada um dos cinco clados.

# A funçãoo calcula a riqueza de espécies para um dado ponto no tempo
div.times <- c(max(branching.times(balaenopteridae.tree)), max(branching.times(delphinidae.tree)), max(branching.times(phocoenidae.tree)), max(branching.times(ziphidae.tree)), max(branching.times(othercetaceans.tree)))

# Funçãoo modificada a partir de plot_dtt do pacote RPANDA
plotdtt <- function (fit.bd, tot_time, N0, col = 1, add = FALSE, div.time, xlim, ylim)
{
  if (!inherits(fit.bd, "fit.bd"))
    stop("object \"fit.bd\" is not of class \"fit.bd\"")
  t <- seq(tot_time-div.time, tot_time, 0.01)
  if ("f.mu" %in% attributes(fit.bd)$names) {
    r <- function(t) {
      -fit.bd$f.lamb(t) + fit.bd$f.mu(t)
    }
    R <- function(s) {
      RPANDA:::.Integrate(Vectorize(r), 0, s)
    }
    N <- N0 * exp(Vectorize(R)(t))
                    #dev.new()
    if(add == FALSE)
      {
    plot(-t, N, type = "l", xlab = "time", ylab = "Number of species", 
       main = "Diversity Through Time", col = col, xlim = xlim, ylim = ylim)
  }
    else
      {
        lines(-t, N, type = "l", xlab = "time", ylab = "Number of species", 
           main = "Diversity Through Time", col = col, xlim = xlim, ylim = ylim)
      }
  }
  else {
    r <- function(t) {
      -fit.bd$f.lamb(t)
    }
    R <- function(s) {
      RPANDA:::.Integrate(Vectorize(r), 0, s)
    }
    N <- N0 * exp(Vectorize(R)(t))
                    #dev.new()
    if(add == FALSE)
      {
    plot(-t, N, type = "l", xlab = "time", ylab = "Number of species", 
       main = "Diversity Through Time", col = col, xlim = xlim, ylim = ylim)
  }
    else
      {
        lines(-t, N, type = "l", xlab = "time", ylab = "Number of species", 
           main = "Diversity Through Time", col = col, xlim = xlim, ylim = ylim)
      }
  }
}

par(mfrow = c(1, 1))
plotdtt(results$balaenopteridae$bcstdcst, div.times[1], N0 = Ntip(balaenopteridae.tree), xlim = c(-max(div.times), 0), ylim = c(0, 150), div.time = div.times[1])
plotdtt(results$delphinidae$bvardcst, div.times[2], N0 = Ntip(delphinidae.tree), col = 6, add = TRUE, xlim = c(-max(div.times), 0), ylim = c(0, 150), div.time = div.times[2])
plotdtt(results$phocoenidae$bcstdcst, div.times[3], N0 = Ntip(phocoenidae.tree), col = "goldenrod", add = TRUE, xlim = c(-max(div.times), 0), ylim = c(0, 150), div.time = div.times[3])
plotdtt(results$ziphidae$bcstdcst, div.times[4], N0 = Ntip(ziphidae.tree), col = 4, add = TRUE, xlim = c(-max(div.times), 0), ylim = c(0, 150), div.time = div.times[4])
plotdtt(results$othercetaceans$bcstdvar, div.times[5], N0 = Ntip(othercetaceans.tree), col = "darkred", add = TRUE, xlim = c(-max(div.times), 0), ylim = c(0, 150), div.time = div.times[5])
legend("topleft", legend = c("Balaenopteridae", "Delphinidae", "Phocoenidae", "Ziphidae", "Other Cetaceans"), text.col = c(1, 6, "goldenrod", 4, "darkred"))

Desafio

  1. Ajuste os modelos de diversificação usando o RPANDA para a filogenia inteira de baleias, e plote a diversidade no tempo usando partes dos códigos anteriores.

  2. Rode o BAMM na filogenia de Anolis (disponível no arquivo de exemplo) e plote a filogenia com as taxas de diversificação para cada ramo.