Ejemplo BiSSE con ranas de la familia Dendrobatidae

En este tutorial vamos a buscar entender si hay una relación entre la coloración de ranas dardo y las tasas de diversificación de estas ranas.

Para esto, vamos a partir de dos datos de entrada: La filogenia de este grupo (que para este caso ASUMIMOS como fija) y una matriz donde codifiquemos la coloración, en este caso como “aposematico” y criptico.

Code

# Leer el arbol
T <- readTrees("./trees/dendobatidae_tree.nex")[1]

# Leer la matriz de caracteres
data <- readCharacterData("./data/dendrobatidae_coloration.nex")
## 1 para especies con colores conspicuos y 0 para especies con colores no conspicuos
##Numero de estados
NUM_STATES = 2

##Edad de la raíz
root <- T.rootAge()

##Proporción de esécies muestreadas
rho <- T.ntips()/220

Despues, vamos a definir los vectores de movimientos y monitores, esta parte es importante, porque a traves de estos nuestra MCMC va a poder avanzar (movimientos) y va a poder evaluar y registrar la probabilidad posterior de los nuevos cambios en la cadena (monitores)

Code

##Definir vectores para movimientos y monitores

moves    = VectorMoves()
monitors = VectorMonitors()

Ahora, vamos a definir nuestros priors, definiendo distribuciones con base en los conocimientos que tenemos de nuestro grupo de estudio.

Code

#  Tasa promedio de diversificación

rate_mean <- ln(ln(220.0/2.0) / T.rootAge())

# varianza
rate_sd <- 0.587405

for (i in 1:NUM_STATES) {
# lognormales de las especiaciones
  log_speciation[i] ~ dnNormal(mean=rate_mean, sd=rate_sd) 
  speciation[i] := exp(log_speciation[i])
  moves.append(mvSlide(log_speciation[i], delta=0.20, tune=true, weight=3.0))

# lognormales de las extinciones
  log_extinction[i] ~ dnNormal(mean=rate_mean, sd=rate_sd)  
  extinction[i] := exp(log_extinction[i])
  moves.append(mvSlide(log_extinction[i], delta=0.20, tune=true, weight=3.0))
}

##Matriz de transición entre estados del carácter

# Cantidad de transiciónes en el árbol
rate_pr := T.treeLength() / 10

##Tasas de transición

# Tasa de transición de estado 1 → estado 2 
rate_12 ~ dnExp(rate_pr)

# Tasa de transición de estado 2 → estado 1 
rate_21 ~ dnExp(rate_pr)

# Definición de movimientos 
moves.append(mvScale(rate_12, weight=2))
moves.append(mvScale(rate_21, weight=2))

##Crear matriz de tasas de transición
rate_matrix := fnFreeBinary([rate_12, rate_21 ], rescaled=false)

##Frecuencias de estado en la raíz
rate_category_prior ~ dnDirichlet(rep(1,NUM_STATES))
moves.append(mvDirichletSimplex(rate_category_prior, tune=true, weight=2))

El caso de la filogenia es especial, ya que definimos el tipo de proceso que la generó, pero que al no estimarla la “clampeamos”, es decir que la tomamos como el resultado observado de ese proceso. Igualmente, tomamos la coloración de las especies actuales como el resultado “observado” del proceso de evolución de esta característica.

Code

##Construcción del árbol a partir de un proceso BD dependiente de caracter
timetree ~ dnCDBDP(rootAge = root, speciationRates = speciation, extinctionRates = extinction, Q = rate_matrix, pi = rate_category_prior, delta = 1.0, rho = rho, condition = "time")

##Fijar árbol y data "observados"

timetree.clamp(T)
timetree.clampCharData(data)

Ahora, vamos a configurar la cadena MCMC para iniciar, primero vamos a “envolver nuestro modelo” (Especificamos una variable y se van a tomar todos los nodos del modelo conectadas directa o indirectamente a esta), para despues definir los archivos de salida de nuestro análisis.

Code

## "Envolver" modelo

mymodel = model(rate_matrix)

##Configuración de los monitores

monitors.append(mnModel(filename="./output/BiSSE_Dendrobatidae_coloration.log", printgen=1))
monitors.append(mnJointConditionalAncestralState(tree=timetree, cdbdp=timetree, type="Standard", printgen=1, withTips=true, withStartStates=false, filename="./output/anc_states_BiSSE_Dendrobatidae_coloration.log"))
monitors.append(mnStochasticCharacterMap(cdbdp=timetree, printgen=10, filename="./output/SCHM_BiSSE_Dendrobatidae_coloration.log"))
monitors.append(mnScreen(printgen=10, rate_12, rate_21, speciation, extinction))

Ahora, vamos a ejecutar el MCMC

Code

## Iniciar y ejecutar el MCMC

mymcmc = mcmc(mymodel, monitors, moves, nruns=2, moveschedule="random")

#mymcmc.burnin(generations = 4000, tuningInterval=100)
mymcmc.run(generations=50000)

Despues de los análisis vamos a generar archivos con las historias del caracter reconstruidas sobre el arbol filogenético.

Code


##Arbol anotado con estados ancestrales (MAP)
anc_states = readAncestralStateTrace("./output/anc_states_BiSSE_Dendrobatidae_coloration_run_1.log")

anc_tree = ancestralStateTree(tree=T, ancestral_state_trace_vector=anc_states, include_start_states=false, file="./output/BiSSE_anc_states_Dendrobatidae_results.tree", burnin=0.2, summary_statistic="MAP", site=1)

#Mapeo estocastico de caracteres

anc_state_trace = readAncestralStateTrace("./output/SCHM_BiSSE_Dendrobatidae_coloration_run_1.log")

characterMapTree(tree=T, ancestral_state_trace_vector=anc_state_trace, character_file="./ouput/BiSSE_Dendrobatidae_stoch_map_character.tree", posterior_file="./output/BiSSE_Dendrobatidae_stoch_map_posterior.tree", burnin=0.2, reconstruction="marginal")

Finalmente, vamos a gráficar los resultados del análisis BiSSE utilizando los paquetes ggplot2 y ggtree de R.

Code

library(tidyverse)
library(RevGadgets)
library(ggtree)

##Cargar el archivo de estados ancestrales
bisse_file <- paste0("./output/BiSSE_anc_states_Dendrobatidae_results.tree")
p_anc <- processAncStates(bisse_file)


  |                                              
  |                                        |   0%
  |                                              
  |========================================| 100%

Code

# Definir paleta de colores personalizada
traitcols <- c("#634138","#FF620B")  

# Generar el gráfico con estados 
plot <- plotAncStatesMAP(p_anc,
        tree_layout = "circular",
        tip_labels_size = 1) +
        # Aplicar la paleta de colores
        scale_color_manual(name = "Estado", values = traitcols, 
                           breaks = c("A", "B"), labels = c("Críptico", "Aposemático")) +
        scale_fill_manual(name = "Estado", values = traitcols, 
                           breaks = c("A", "B"), labels = c("Críptico", "Aposemático")) +
        # Modificar leyenda y otros elementos de ggplot2
        theme(legend.position = c(0.1, 0.85),
              legend.key.size = unit(0.3, 'cm'), # Tamaño de la clave de la leyenda
              legend.title = element_text(size = 6), # Tamaño del título de la leyenda
              legend.text = element_text(size = 4))  # Tamaño del texto de la leyenda

# Mostrar el gráfico
plot

Code

# Mostrar el gráfico con probabilidades posteriores en forma de Pie

plot_pie<- plotAncStatesPie(p_anc,layout = "circular",tip_labels_size = 1,pie_color = c("A" = "#634138",
                                            "B" = "#FF620B")) +
  scale_color_manual(name = "Estado", values = traitcols, 
                           breaks = c("A", "B"), labels = c("Criptico","Aposematico")) +
  theme(legend.position = c(0.1, 0.85),
        legend.key.size = unit(0.3, 'cm'), # Tamaño de la clave de la leyenda
        legend.title = element_text("Estado",size = 6), # Tamaño del título de la leyenda
        legend.text = element_text(size = 4))  # Tamaño del texto de la leyenda
                                     
plot_pie

Code

# Generar el gráfico con ramas coloreadas por estado ancestral
plot_2 <- ggtree(p_anc, layout = "circular", aes(color = anc_state_1)) +
  scale_color_manual(name = "Estado", values = traitcols, 
                           breaks = c("A", "B"), labels = c("Criptico","Aposematico" )) +  
  ggtitle("Dendrobatidae aposematism") +
  theme(plot.title = element_text(hjust = 0,vjust = 1.0),
        legend.position = c(0.0, 0.9),
        legend.key.size = unit(0.3, 'cm'), # Tamaño de la clave de la leyenda
        legend.title = element_text(size = 6), # Tamaño del título de la leyenda
        legend.text = element_text(size = 4))  # Tamaño del texto de la leyenda


plot_2

Code

# Leer el archivo de salida de RevBayes
bisse <- read.table("./output/BiSSE_Dendrobatidae_coloration_run_1.log", header = TRUE)

# Definir colores para los estados 0 y 1
 traitcols <- c("#634138","#FF620B")    

# Crear un dataframe con solo las tasas de transición
transition_rates <- data.frame(
  dens = c(bisse$rate_12, bisse$rate_21), 
  rate = rep(c("0", "1"), each = length(bisse$rate_12))
)

# Crear un gráfico de violín para las tasas de transición
violin_transitions <- ggplot(transition_rates, aes(x = rate, y = dens, fill = rate)) +
  geom_violin(trim = FALSE) +
  labs(title = "Tasas de Transición", x = "Estado", y = "Tasa") +
  scale_fill_manual(name = "Estado", values = traitcols,labels = c("Criptico","Aposematico" )) +
  theme_classic()

# Mostrar la gráfica
violin_transitions

Code

# Calcular la tasa de diversificación neta: especiación - extinción
netdiversification_rates <- data.frame(
  dens = c(bisse$speciation.1. - bisse$extinction.1.,
           bisse$speciation.2. - bisse$extinction.2.),
  rate = rep(c("0", "1"), each = length(bisse$speciation.1.))
)

# Convertir la variable de estado en factor
netdiversification_rates$rate <- factor(netdiversification_rates$rate, levels = c("0", "1"))

# Crear el gráfico de violín para tasas de diversificación neta
violin_diversification <- ggplot(netdiversification_rates, aes(x = rate, y = dens, fill = rate)) +
  geom_violin(trim = FALSE) +
  labs(title = "Tasas de Diversificación Neta", x = "Estado", y = "Tasa") +
  scale_fill_manual(name = "Estado", values = traitcols,labels = c("Criptico","Aposematico" )) +
  theme_classic()

# Mostrar la gráfica
violin_diversification

Code

# Crear un dataframe con tasas de especiación para cada estado
speciation_rates <- data.frame(
  dens = c(bisse$speciation.1.,  bisse$speciation.2.), 
  rate = rep(c("0", "1"), each = length(bisse$speciation.1.))
)

# Crear gráfico de violín para las tasas de especiación
violin_speciation <- ggplot(speciation_rates, aes(x = rate, y = dens, fill = rate)) +
  geom_violin(trim = FALSE) +
  labs(title = "Tasas de Especiación", x = "Estado", y = "Tasa de Especiación") +
  scale_fill_manual(name = "Estado", values = traitcols,labels = c("Criptico","Aposematico" )) +
  theme_classic()

# Mostrar la gráfica
violin_speciation

Code

# Crear un dataframe con tasas de extinción para cada estado
extinction_rates <- data.frame(
  dens = c(bisse$extinction.1., bisse$extinction.2.), 
  rate = rep(c("0", "1"), each = length(bisse$extinction.1.))
)

# Crear gráfico de violín para tasas de extinción
violin_extinction <- ggplot(extinction_rates, aes(x = rate, y = dens, fill = rate)) +
  geom_violin(trim = FALSE) +
  labs(title = "Tasas de Extinción", x = "Estado", y = "Tasa de Extinción") +
  scale_fill_manual(name = "Estado", values = traitcols,labels = c("Criptico","Aposematico" )) +
  theme_classic()

# Mostrar la gráfica
violin_extinction