Introducción a modelos State Speciaction-Extinción (SSE)

Presentación: Introducción a modelos State Speciaction-Extinción (SSE)

Haz clic en la imagen para ver el PDF de la presentación

Analizaremos si las tasas de diversificación están asociadas al tipo de hospedero en los escarabajos meloidos (Meloidae). Las subfamilias Nemognathinae y Meloinae presentan un ciclo de vida complejo, caracterizado por múltiples metamorfosis y estrategias de parasitoidismo. La mayoría de los géneros y tribus dentro de estas subfamilias son parasitoides de abejas, incluyendo especies foréticas y no foréticas. Sin embargo, dos tribus han adoptado una estrategia distinta, alimentándose de huevos de saltamontes. Estas diferencias en las estrategias de vida están asociadas con marcadas variaciones en la riqueza de especies entre clados (López-Estrada et al., 2019).

Modelo BiSSE (Binary state speciation and extinction, Madisson et al. 2007)

BiSSE el modelo de especiación y extinción dependendiente de un caracter binario y es un proceso estocástico que se forma de la composición de dos procesos de nacimiento y muerte. Los dos procesos se conectan a través de dos tasas de cambio conocidas como las tasas de transición. El supuesto más importante del modelo BiSSE es que cada estado tiene su propia tasa de especiación \(\lambda\) y su propia tasa de extinción \(\mu\) que representan el momento en el cuál un linaje se divide en dos o se extingue. Este es un supuesto muy importante porque la implicación es que la acumulación (o la falta de acumulación) de linajes es el resultado del valor del estado. Pueden encontrar mas información en el material del taller Filo-Bayes - BISSE.

Escarabajos meloidos

Las transiciones entre estados ocurren a tasas \(q_{01}\) (de abeja a saltamontes) y \(q_{10}\) (de saltamontes a abeja). Cada estado tiene su propia tasa de especiación (\(\lambda_0, \lambda_1\)) y extinción (\(\mu_0, \mu_1\)).

Matriz de tasas de transición (\(Q\))

La matriz de transición \(Q\) define la probabilidad de cambio entre estados en un tiempo infinitesimal:

\[ Q = \begin{bmatrix} -q_{01} & q_{01} \\ q_{10} & -q_{10} \end{bmatrix} \]

Donde:

  • \(q_{01}\): es la tasa de transición de hospedero abeja a hospedero saltamontes.

  • \(q_{10}\): es la tasa de transición de hospedero saltamontes a hospedero abeja.

  • La diagonal contiene los valores negativos de las tasas de salida de cada estado (\(-q_{01}\) y \(-q_{10}\)), de modo que cada fila suma cero.

Parámetros en el modelo BiSSE

El modelo incluye seis parámetros fundamentales:

  • \(\lambda_0\): tasa de especiación en el estado 0 (hospedero abeja)

  • \(\lambda_1\): tasa de especiación en el estado 1 (hospedero saltamontes)

  • \(\mu_0\): tasa de extinción en el estado 0

  • \(\mu_1\): tasa de extinción en el estado 1

  • \(q_{01}\): tasa de cambio de hospedero abeja a hospedero saltamontes

  • \(q_{10}\): tasa de cambio de hospedero saltamontes a hospedero abeja

Descarga del Árbol Filogenético y la matriz codificada

📥 Puedes descargar el archivo .tre con el árbol en el siguiente enlace:

📥 Descargar Árbol Filogenético

📥 Puedes descargar el archivo .nex con la matriz codificada en el siguiente enlace:

📥 Descargar Matriz Codificada

Guarda los archivos en la carpeta correspondiente y verifica su ubicación antes de continuar.

Carga de la filogenia y los datos

Crea un script .Rev y copia el codigo de la siguiente sección:

# Leer el arbol
T <- readTrees("../data/ConLyttini.tre")[1]

# Leer la matriz de caracteres
datos <- readCharacterData("../data/Traits_conLyttini_01.nex")

Definir el número de estados en los modelos BISSE

NUM_STATES = 2

Definir vectores para almacenar los movimientos de MCMC y los monitores de resultados.

moves    = VectorMoves()
monitors = VectorMonitors()

Valores de la log-normal

De acuerdo con Nee et al. (1994) el número esperado de linajes en la corona de un clado con n taxones bajo un proceso de nacimiento y muerte en el tiempo t es:

\((\lambda - \mu) = \frac{In(n)/2}{t}\)

En donde: n = es el número de especies y t = es la edad de la filogenia

# Definimos una a priori log normal
rate_mean <- ln(ln(3000.0/2.0) / T.rootAge())

# Con una varianza que es amplia
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))
}

Definir la matriz de tasas de transición entre los estados del carácter

# Calcula la suma de todas las longitudes de las ramas del árbol filogenético
# Dividido por 10: Esto asume que en promedio ocurren 10 cambios de estado a lo largo del árbol
rate_pr := T.treeLength() / 10

Definición de las 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))

¿Por qué usar una exponencial?

  • La distribución exponencial es adecuada para modelar tasas de transición porque impone la condición de que siempre sean positivas.

  • También asigna mayor probabilidad a valores pequeños, lo que significa que los cambios de estado no ocurren demasiado rápido.

Creación de la matriz de tasas de transición

rate_matrix := fnFreeBinary([rate_12, rate_21 ], rescaled=false)

Frecuencias del estado en la raíz

rate_category_prior ~ dnDirichlet(rep(1,NUM_STATES))
moves.append(mvDirichletSimplex(rate_category_prior, tune=true, weight=2))

Fijar la edad de la raíz

root <- T.rootAge()

Definir \(\rho\) : Proporción de especies muestreadas

rho <- T.ntips()/3000

Construcción del árbol evolutivo generado a partir de un proceso de nacimiento y muerte dependiente del carácter

timetree ~ dnCDBDP(rootAge = root, speciationRates = speciation, extinctionRates = extinction, Q = rate_matrix, pi = rate_category_prior, delta = 1.0, rho = rho, condition = "time")

🔹 Parámetros del proceso de nacimiento y muerte (dnCDBDP):

  • rootAge = root → Se fija la edad de la raíz usando root = T.rootAge().

  • speciationRates = speciation → Se usa la tasa de especiación (λ).

  • extinctionRates = extinction → Se usa la tasa de extinción (μ).

  • Q = rate_matrix → Usa la matriz de tasas de transición de caracteres (rate_12, rate_21).

  • pi = rate_category_prior → Usa la distribución a priori de los estados en la raíz.

  • delta = 1.0 → Parámetro de ajuste para la distribución de tiempos de especiación/extinción.

  • rho = rho → Probabilidad de muestreo de especies (rho = T.ntips()/3000).

  • condition = "time" → Condiciona el proceso a un tiempo de raíz fijo.

Fijar (clamp) los datos observados

timetree.clamp(T)
timetree.clampCharData(datos)

Creación del modelo

### workspace model wrapper ###
mymodel = model(rate_matrix)

Configuración de los monitores

Los monitores guardan los resultados en archivos o los imprimen en pantalla.

monitors.append(mnModel(filename="../out/Bisse/mitos_BiSSE_clyttini.log", printgen=1))
monitors.append(mnJointConditionalAncestralState(tree=timetree, cdbdp=timetree, type="Standard", printgen=1, withTips=true, withStartStates=false, filename="../out/Bisse/anc_states_MITOS_BiSSE_clyttini.log"))
monitors.append(mnStochasticCharacterMap(cdbdp=timetree, printgen=10, filename="../out/Bisse/SCHM_BiSSE_clyttini.log"))
monitors.append(mnScreen(printgen=10, rate_12, rate_21, speciation, extinction))

Inicializar y Ejecutar el MCMC

Con el modelo y los monitores listos, creamos el objeto MCMC e iniciamos la simulación.

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

  • mymodel: hace referencia al modelo que estás analizando.

  • monitors: son las variables que se registrarán en el archivo de salida.

  • moves: operadores MCMC que proponen cambios en los parámetros del modelo.

  • nruns=1: solo una cadena MCMC (puedes usar más si deseas una inferencia más robusta).

  • moveschedule="random": ejecuta los operadores en orden aleatorio para evitar sesgos en la exploración del espacio de parámetros.

Ejecutar el MCMC

mymcmc.run(generations=400)

Generar un árbol anotado con los estados ancestrales (MAP)

anc_states = readAncestralStateTrace("../out/Bisse/anc_states_MITOS_BiSSE_clyttini.log")

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

Generar un mapeo estocástico de caracteres (Stochastic Character Mapping)

anc_state_trace = readAncestralStateTrace("../out/Bisse/SCHM_BiSSE_clyttini.log")

characterMapTree(tree=T, ancestral_state_trace_vector=anc_state_trace, character_file="../out/Bisse/mitos_bisse_stoch_map_character.tree", posterior_file="../out/Bisse/mitos_bisse_stoch_map_posterior.tree", burnin=0.2, reconstruction="marginal")

Visualización en R con RevGadgets

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

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

# Definir paleta de colores personalizada
traitcols <- c("#F29494", "#5FD9D9")  

# Generar el gráfico con estados 
plot <- plotAncStatesMAP(p_anc,
        tree_layout = "rect",
        tip_labels_size = 1) +
        # Aplicar la paleta de colores
        scale_color_manual(values = traitcols) +
        scale_fill_manual(values = traitcols) +
        # 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

# Generar el gráfico con ramas coloreadas por estado ancestral
plot_2 <- ggtree(p_anc, layout = "rectangular", aes(color = anc_state_1)) +
  scale_color_manual(name = "Estado", values = traitcols) +  
  theme_minimal() +
  ggtitle("BISSE") +
  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_2

Graficar las tasas de transición

library(tidyverse)

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

# Definir colores para los estados 0 y 1
traitcols <- c("#F29494", "#5FD9D9")  

# 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) +
  theme_classic()

# Mostrar la gráfica
violin_transitions

Graficar las tasas de diversificación neta

# 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) +
  theme_classic()

# Mostrar la gráfica
violin_diversification

Graficar tasas de especiación

# 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) +
  theme_classic()

# Mostrar la gráfica
violin_speciation

Graficar tasas de extinción

# 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) +
  theme_classic()

# Mostrar la gráfica
violin_extinction

Modelo de diversificación con estados escondidos (HISSE)

Uno de los descubrimientos más importantes en el campo de los modelos de diversificación dependientes de estado (SSE) fue hecho por Rabosky y Goldberg (2015). Estos autores encontraron, que BiSSE generaba un elevado error tipo I. En estadística este error se refiere a rechazar la hipótesis nula cuando en general la hipótesis es verdadera.

Modelo HiSSE con dos estados y asociando a dos estados escondidos

  • Si tenemos dos estados 0 y 1, estos se convierten ahora en cuatro estados 0A, 0B, 1A, 1B.

  • Para cada uno de los estados definimos una nueva tasa de especiación \(\lambda_{0A}\), \(\lambda_{0B}\), \(\lambda_{1A}\) y \(\lambda_{1B}\).

  • Para cada uno de los estados definimos una nueva tasa de extinción \(\mu_{0A}\), \(\mu_{0B}\), \(\mu_{1A}\) y \(\mu_{1B}\).

Matriz de tasas de transición Q

\[ Q = \begin{bmatrix} - & q_{0A \to 0B} & q_{0A \to 1A} & 0 \\ q_{0B \to 0A} & - & 0 & q_{0B \to 1B} \\ q_{1A \to 0A} & 0 & - & q_{1A \to 1B} \\ 0 & q_{1B \to 0B} & q_{1B \to 1A} & - \end{bmatrix} \]

La matriz Q muestra cómo los estados observables y ocultos interactúan, ayudando a inferir cómo los cambios de carácter afectan la especiación y extinción.

Carga de la filogenia y los datos

Crea un script .Rev y copia el codigo de la siguiente sección:

# Leer el arbol
T <- readTrees("../data/ConLyttini.tre")[1]

# Leer la matriz de caracteres
datos <- readCharacterData("../data/Traits_conLyttini_01.nex")

Definir el número de estados en los modelos HISSE y paremetros globales

NUM_TOTAL_SPECIES     = 3000.0
NUM_STATES            = 2
NUM_HIDDEN            = 2
NUM_RATES             = NUM_STATES * NUM_HIDDEN
H                     = 0.587405
# Expansión de los Datos para Incluir Estados Ocultos 
data_exp <- datos.expandCharacters(NUM_HIDDEN)

Definir vectores para almacenar los movimientos de MCMC y los monitores de resultados.

moves    = VectorMoves()
monitors = VectorMonitors()

Extraer información del árbol filogenético

taxa <- T.taxa()
tree_length <- T.treeLength()

Definir la Media de las Tasas de Especiación y Extinción

rate_mean <- (NUM_TOTAL_SPECIES - 2) / tree_length

✅ Calcula la media esperada de las tasas de especiación y extinción basada en el número total de especies y la longitud del árbol.

Especiación y Extinción Oculta (Hidden Speciation Rates)

ln_speciation_hidden_mean <- ln(1.0)
ln_extinction_hidden_mean <- ln(1.0)

✅ Define la media logarítmica de las tasas ocultas de especiación como 1.0 en una escala relativa. ✔️ Esto significa que los valores generados serán proporcionales, y no absolutos.

Variación en las Tasas Ocultas de Especiación y Extinción

speciation_hidden_sd ~ dnExponential(1.0 / H)
extinction_hidden_sd ~ dnExponential(1.0 / H)

Movimiento MCMC para Optimizar la Desviación Estándar

moves.append(mvScale(speciation_hidden_sd, lambda=1, tune=true, weight=2.0))
moves.append(mvScale(extinction_hidden_sd, lambda=1, tune=true, weight=2.0))

Definir la Distribución de las Tasas Ocultas

speciation_hidden_unormalized := fnDiscretizeDistribution(dnLognormal(ln_speciation_hidden_mean, speciation_hidden_sd), NUM_HIDDEN)
extinction_hidden_unormalized := fnDiscretizeDistribution(dnLognormal(ln_extinction_hidden_mean, extinction_hidden_sd), NUM_HIDDEN)

✅ Se genera un conjunto de tasas ocultas de especiación a partir de una distribución lognormal.

Normalización de las Tasas Ocultas

speciation_hidden := speciation_hidden_unormalized / mean(speciation_hidden_unormalized)
extinction_hidden := extinction_hidden_unormalized / mean(extinction_hidden_unormalized)

✅ Se normalizan las tasas ocultas dividiéndolas por su media, para que el promedio sea 1.0. ✔️ Esto garantiza que los valores sean relativos y comparables entre estados.

Definir Tasas de Especiación para Estados Observables

for (i in 1:NUM_STATES) {
    speciation_observed[i] ~ dnLoguniform( 1E-6, 1E2)
    speciation_observed[i].setValue( (NUM_TOTAL_SPECIES-2) / tree_length )
    moves.append( mvScale(speciation_observed[i], lambda=1.0, tune=true, weight=3.0) )
        ### Create a loguniform distributed variable for the extinction rate
    extinction_observed[i] ~ dnLoguniform( 1E-6, 1E2)
    extinction_observed[i].setValue( speciation_observed[i] / 10.0 )
    moves.append( mvScale(extinction_observed[i], lambda=1.0, tune=true, weight=3.0) )
}

Cálculo de las Tasas Combinadas

for (j in 1:NUM_HIDDEN) {
    for (i in 1:NUM_STATES) {
        index = i+(j*NUM_STATES)-NUM_STATES
        speciation[index] := speciation_observed[i] * speciation_hidden[j]
        extinction[index] := extinction_observed[i] * extinction_hidden[j]
    }
}

Definir la Matriz de Transición entre Estados Observables

rate_pr := T.treeLength() / 10
for ( i in 1:(NUM_STATES*(NUM_STATES-1)) ) {
    transition_rates[i] ~ dnExp(rate_pr)
    moves.append( mvScale(transition_rates[i], lambda=0.50, tune=true, weight=3.0) )
}

Definir la Matriz de Transición entre Estados Ocultos

hidden_rate ~ dnExponential(rate_pr)
moves.append( mvScale(hidden_rate, lambda=0.5, tune=true, weight=5) )
for (i in 1:(NUM_HIDDEN * (NUM_HIDDEN - 1))) {
    R[i] := hidden_rate
}

Crear la Matriz de Transición Q Combinada

rate_matrix := fnHiddenStateRateMatrix(transition_rates, R, rescaled=false)

Definir la Distribución de los Estados en la Raíz

rate_category_prior ~ dnDirichlet( rep(1,NUM_RATES) )

Optimización con Movimientos MCMC

moves.append( mvBetaSimplex(rate_category_prior, tune=true, weight=2) )
moves.append( mvDirichletSimplex(rate_category_prior, tune=true, weight=2) )

Definir la Edad de la Raíz del Árbol

root <- T.rootAge()

Definir la Probabilidad de Muestreo de Especies (\(\rho\))

rho <- T.ntips()/3000.0

Definir el Árbol Bajo un Proceso de Diversificación Dependiente del Estado

timetree ~ dnCDBDP( rootAge = root, speciationRates = speciation, extinctionRates = extinction, Q = rate_matrix, delta = 1.0, pi = rate_category_prior, rho = rho, condition = "survival")

✔️ Parámetros clave:

  • rootAge = root → La edad de la raíz se fija a la edad del árbol observado.

  • speciationRates = speciation → Tasas de especiación combinadas (λ), que incluyen efectos de los estados ocultos.

  • extinctionRates = extinction → Tasas de extinción combinadas (μ), modeladas de la misma manera.

  • Q = rate_matrix → Matriz de tasas de transición entre estados (0A ↔︎ 1A, 0B ↔︎ 1B, etc.).

  • delta = 1.0 → Modelo estándar en el que las tasas de diversificación dependen completamente del estado.

  • pi = rate_category_prior → Prior de las probabilidades de estado en la raíz.

  • rho = rho → Ajuste por sesgo de muestreo de especies.

  • condition = "survival" → Se modela solo en árboles que han sobrevivido hasta el presente.

Ajustar el Modelo a la Filogenia Observada

timetree.clamp(T)

Ajustar el Modelo a los Datos de Caracteres

timetree.clampCharData(data_exp)

Definir el Modelo en RevBayes

mymodel = model(rate_matrix)

Configurar Monitores para Guardar Resultados

# Monitor Principal del Modelo
monitors.append( mnModel(filename="../out/Hisse/mitos_HiSSE_clyttini_TUTscript.log", printgen=1) )

# Monitor de Estados Ancestrales
monitors.append(mnJointConditionalAncestralState(tree=timetree, cdbdp=timetree, type="NaturalNumbers", printgen=1, withTips=true, withStartStates=false, filename="../out/Hisse/anc_states_mitos_HiSSE_clyttini_TUTscript.log"))

# Monitor de Mapeo Estocástico de Caracteres
monitors.append(mnStochasticCharacterMap(cdbdp=timetree, printgen=10, filename="../out/Hisse/stoch_char_map_mitos_HiSSE_clyttini_TUTscript.log", include_simmap=true))

# Monitor en Pantalla
monitors.append(mnScreen(printgen=10, speciation_observed, extinction_observed))

Configurar el Algoritmo MCMC

### workspace mcmc
mymcmc = mcmc(mymodel, monitors, moves, nruns=1, moveschedule="random", combine="mixed")

Ejecutar el MCMC

mymcmc.run(generations=400)

Cargar la Trazabilidad de los Estados Ancestrales

anc_states2 = readAncestralStateTrace("../out/Hisse/anc_states_mitos_HiSSE_clyttini_TUTscript.log")

Construcción del Árbol con Estados Ancestrales

anc_tree = ancestralStateTree(tree=T, ancestral_state_trace_vector = anc_states2, include_start_states=false, file = "../out/Hisse/anc_states_mitos_HiSSE_clyttini_TUTscript.tree", burnin = 2, summary_statistic = "MAP", site=1)

Definir Parámetros para la Resumir Estados Ancestrales

burnin=10
n_time_slices = 500

✅ burnin = 25 → Descarta las primeras 25 iteraciones del MCMC antes de hacer el resumen.
✅ n_time_slices = 500 → Divide el árbol en 500 segmentos de tiempo para interpolar los cambios de estado.

Leer la Historia del Carácter Inferida por Mapeo Estocástico

anc_states_SCHM = readAncestralStateTrace("../out/Hisse/stoch_char_map_mitos_HiSSE_clyttini_TUTscript.log")

Crear el Árbol de Mapeo Estocástico de Caracteres

char_map_tree = characterMapTree(tree=T, 
                 ancestral_state_trace_vector=anc_states_SCHM, 
                 character_file="../out/Hisse/stoch_mitos_HiSSE_clyttini_TUTscriptmarginal_char.tree", 
                 posterior_file="../out/Hisse/stoch_mitos_HiSSE_clyttini_TUTscriptmarginal_post.tree", 
                 burnin=burnin, 
                 num_time_slices=n_time_slices)

Visualización en R con RevGadgets

# Cargar el archivo de estados ancestrales
hisse_file <- paste0("../docs/u1_PatDiv/output/hisse/anc_states_mitos_HiSSE_clyttini_TUTscript.tree")

p_anc_hisse <- processAncStates(hisse_file)

# Definir paleta de colores personalizada
traitcols <- c("#F29494", "#5FD9D9", "#9DD962", "#D48DF2")  # Rojo para estado 0, azul para estado 1

# Generar el gráfico con estados ancestrales
plot_hisse <- plotAncStatesMAP(p_anc_hisse,
        tree_layout = "rect",
        tip_labels_size = 1) +
        # Aplicar la paleta de colores
        scale_color_manual(values = traitcols) +
        scale_fill_manual(values = traitcols) +
        # 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_hisse

# Generar el gráfico con ramas coloreadas por estado ancestral
plot_hisse_2 <- ggtree(p_anc_hisse, layout = "rectangular", aes(color = anc_state_1)) +
  scale_color_manual(name = "Estado", values = traitcols) +  
  theme_minimal() +
  ggtitle("HISSE") +
  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_hisse_2

Graficar las tasas de diversificación neta

# Leer el archivo de salida de RevBayes
hisse <- read.table("../docs/u1_PatDiv/output/hisse/mitos_HiSSE_clyttini_TUTscript.log", header = TRUE)

# Definir paleta de colores personalizada
traitcols <- c("#F29494", "#5FD9D9", "#9DD962", "#D48DF2")  # Rojo para estado 0, azul para estado 1

# Calcular la tasa de diversificación neta: especiación - extinción
netdiversification_rates <- data.frame(
  dens = c(hisse$speciation_observed.1. - hisse$extinction_observed.1.,
           hisse$speciation_observed.2. - hisse$extinction_observed.2.,
           hisse$speciation_hidden.1. - hisse$extinction_hidden.1.,
           hisse$speciation_hidden.2. - hisse$extinction_hidden.2.),
  rate = rep(c("0A", "1A", "0B", "1B"), each = length(hisse$speciation_observed.1.))
  )

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

# 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) +
  theme_classic()

# Mostrar la gráfica
violin_diversification

Graficar tasas de especiación

# Crear un dataframe con tasas de especiación para cada estado
speciation_rates <- data.frame(
  dens = c(hisse$speciation_observed.1.,  hisse$speciation_observed.2.,
           hisse$speciation_hidden.1., hisse$speciation_hidden.2.), 
  rate = rep(c("0A", "1A", "0B", "1B"), each = length(hisse$speciation_observed.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) +
  theme_classic()

# Mostrar la gráfica
violin_speciation

Graficar tasas de extinción

# Crear un dataframe con tasas de extinción para cada estado
extinction_rates <- data.frame(
  dens = c(hisse$extinction_observed.1., hisse$extinction_observed.2.,
           hisse$extinction_hidden.1., hisse$extinction_hidden.2.), 
  rate = rep(c("0A", "1A", "0B", "1B"), each = length(hisse$extinction_observed.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) +
  theme_classic()

# Mostrar la gráfica
violin_extinction