DEC_tutorial
2026-01-27
Creado por Vanessa Parra (Enero 2026) para el taller “Métodos Comparativos para Naturalistas”
Introducción
En el tutorial “Introducción al modelo DEC”, revisamos como configurar la matriz de tasas instantáneas y las probabilidades de transición cladogenética para un modelo DEC sencillo. En este tutorial, vamos a realizar un análisis biogeográfico utilizando un conjunto de datos: lagartijas (Liolaemus) de Sudamérica
Un análisis DEC sencillo
El siguiente tutorial tiene como objetivo estimar las áreas de distribución ancestrales de las lagartijas Liolaemus de Sudamérica, que es uno de los géneros de tetrápodos mas diversos de América, con alrededor de 300 especies descritas. Aunque el género Liolaemus se encuentra ampliamente distribuido en el sur de América del sur, su diversidad se concentra principalmente en los Andes y la Patagonia, regiones en las que se ha propuesto el origen del clado. La orogenia andina ha desempeñado un papel clave en la evolución de este grupo, por lo tanto, se convierten en un sistema ideal para explorar conceptos de biogeografía histórica y filogenia. Para una lectura mas detallada, se recomienda consultar Esquerré et al (2019).
Dentro del genero Liolaemus se reconocen dos clados principales, que han sido clasificados como subgéneros: Liolaemus y Eulaemus. Para este tutorial nos centraremos exclusivamente en Eulaemus. Para comenzar, utilizaremos cuatro áreas —A (Altiplano-Atacama), B (Andes centrales-Chile central), C (Patagonia), D (tierras bajas orientales)—
Para empezar
Para este tutorial se proporcionan los siguientes archivos de datos:
Análisis
Estructura del directorio: Primero, creeamos variables para la gestión de archivos de entrada y salida
range_fn = "data/eulaemus_range.nex"
tree_fn = "data/eulaemus.tre"
out_fn = "output/eulaemus_DEC"Filogenia:
tree <- readTrees(tree_fn)[1]Áreas de distribución:
Leer datos de áreas de distribución binarios (01) de presencia-ausencia
dat_range_01 = readDiscreteCharacterData(range_fn)Convertir rangos binarios en números naturales
dat_range_n = formatDiscreteCharacterData(dat_range_01, "DEC")Parámetros globales
Dimensión de los datos:
Registrar el número de áreas (caracteres) a partir del objeto de datos de caracteres discretos
n_areas = dat_range_01.nchar()
max_areas = 2Aquí podemos visualizar los datos por taxón para observar cómo los caracteres están codificados tanto en un formato legible de presencia-ausencia
dat_range_01[1] Liolaemus_abaucan:
1001como en un formato legible utilizando números naturales
dat_range_n[1]
Liolaemus_abaucan:
8Para la correcta interpretación y visualización de los rangos ancestrales, es necesario registrar la correspondencia entre los estados de rango numéricos y las etiquetas de rango que representan las combinaciones de áreas. Esta información será utilizada posteriormente al generar figuras de estimación de rangos ancestrales.
Para esto almacenaremos el vector que contiene las descripciones de cada estado de rango, el cual asocia cada estado numérico con la combinación de áreas correspondiente
state_desc = dat_range_n.getStateDescriptions()A continuación, esta información se convierte en una cadena de texto con formato tabular y se guarda en un archivo externo. Este archivo funcionará como una tabla de referencia que vincula cada estado numérico con su rango geográfico, facilitando la interpretación de los resultados y la generación de figuras
state_desc_str = "state,range\n"
for (i in 1:state_desc.size())
{
state_desc_str += (i-1) + "," + state_desc[i] + "\n"
}
write(state_desc_str, file=out_fn+".state_labels.txt")Movimientos y monitores:
moves = VectorMoves()
monitors = VectorMonitors()
n_gen = 20000Modelo DEC
Eventos anagenéticos
A continuación, construiremos la matriz de tasas anagenéticas para el modelo DEC. En su forma más simple, esta matriz de tasas requiere un parámetro de dispersión y un parámetro de extirpación.
Para este análisis, se asume que todos los pares de áreas comparten la misma tasa relativa de dispersión y que todas las áreas comparten la misma tasa relativa de extirpación. Con el objetivo de facilitar la exploración de la sensibilidad a los priors, reparametrizaremos la matriz de tasas del modelo DEC para expresar las tasas relativas de los eventos de dispersión frente a los de extirpación. Además, para que las tasas de los eventos anagenéticos puedan medirse en una escala temporal absoluta (por ejemplo, en millones de años), introduciremos un parámetro de tasa biogeográfica.
En primer lugar, se crea un parámetro para la tasa de ocurrencia de eventos anagenéticos de cambio en el rango geográfico a lo largo de las ramas del filogenético.
rate_bg ~ dnLoguniform(1E-4,1E2)
moves.append( mvSlide(rate_bg, weight=2) )Tasa de dispersión:
dispersal_rate ~ dnBeta(1.1, 20)
moves.append( mvSlide(dispersal_rate, weight=2, delta=0.2) )Ahora creamos la matriz para la tasa de dispersión
for (i in 1:n_areas) {
for (j in 1:n_areas) {
dr[i][j] := dispersal_rate
}
}Tasa de extirpación:
A continuación, se define un parámetro para la tasa relativa de extirpación (o tasa de extinción por área). Esta tasa se modela mediante una distribución prior lognormal, parametrizada de modo que tenga una media de 1.
log_sd <- 0.1
log_mean <- ln(1) - 0.5*log_sd^2
extirpation_rate ~ dnLognormal(mean=log_mean, sd=log_sd)
extirpation_rate.setValue(0.01)
moves.append( mvScale(extirpation_rate, weight=1) )Ahora creamos la matriz para la tasa de extirpación
for (i in 1:n_areas) {
for (j in 1:n_areas) {
er[i][j] <- 0.0
}
er[i][i] := extirpation_rate
}¡La matriz Q!
Es importante señalar que la matriz de extirpación (er) es una matriz diagonal, lo que refleja que la extirpación representa la pérdida de una sola área dentro del rango de una especie. Los valores de la diagonal están determinados (:=) por un único parámetro estocástico, extirpation_rate, compartido por todas las áreas. Esto implica que, aunque todas las áreas presentan la misma tasa de extirpación, la pérdida de cada área ocurre de manera independiente a lo largo del tiempo.
Con las matrices de dispersión y extirpación definidas, se construye
la matriz de tasas relativas del modelo DEC (Q_DEC) utilizando la
función fnDECRateMatrix.
Q_DEC := fnDECRateMatrix(dispersalRates=dr,
extirpationRates=er,
maxRangeSize=2)Eventos cladogenéticos
El siguiente paso es construir la matriz de probabilidades cladogenéticas. A diferencia de los procesos anagenéticos, las probabilidades de los eventos cladogenéticos se describen mediante una matriz de probabilidades de transición, y no mediante una matriz de tasas.
Primero vamos a definir un vector que indica que se consideraran eventos de cladogénesis por simpatría de subconjuntos (s) y alopatría (a).
clado_event_types <- [ "s", "a" ]
clado_events ~ dnDirichlet([1,1])
moves.append(mvBetaSimplex(clado_events, alpha=10))
moves.append(mvElementSwapSimplex(clado_events))Probabilidades de simpatría de subconjunto (s) y alopatría (a)
p_sympatry := clado_events[1]
p_allopatry := clado_events[2]Distribucion de las probabilidades
clado_event_probs := simplex(p_sympatry, p_allopatry)Finalmente, construiremos la matriz de probabilidades cladogenéticas
(P_DEC) con fnDECCladoProbs.
P_DEC := fnDECCladoProbs(eventProbs=clado_event_probs,
eventTypes=clado_event_types,
numCharacters=n_areas,
maxRangeSize=2)Integración del modelo DEC
Finalmente juntamos todo los elementos y construimos el modelo DEC.
Una vez definidas la matriz de tasas anagenéticas (Q_DEC) y la matriz
de probabilidades cladogenéticas (P_DEC), ambos componentes se integran
en una cadena de Markov continua sobre el árbol filogenético mediante
dnPhyloCTMCClado.
eulaemus_dec ~ dnPhyloCTMCClado(tree=tree,
Q=Q_DEC,
cladoProbs=P_DEC,
branchRates=rate_bg,
type="NaturalNumbers",
nSites=1)Asociar los datos de áreas de distribución al modelo:
Finalmente, se asocian los rangos observados al modelo. Es importante utilizar los caracteres de rango codificados como números naturales (dat_range_n) y no los caracteres de presencia–ausencia (dat_range_01).
eulaemus_dec.clamp(dat_range_n)Monitores … sí otra vez
A continuación, se definen los monitors, que permiten registrar y visualizar los parámetros y resultados del modelo a lo largo de la inferencia bayesiana.
monitors.append( mnScreen(printgen=100, rate_bg, dispersal_rate, p_sympatry, p_allopatry) )
monitors.append( mnModel(file=out_fn+".model.log", printgen=10) )
monitors.append( mnFile(tree, filename=out_fn+".tre", printgen=10) )
monitors.append( mnJointConditionalAncestralState(tree=tree,
ctmc=eulaemus_dec,
type="NaturalNumbers",
withTips=true,
withStartStates=true,
filename=out_fn+".states.log",
printgen=10) )
monitors.append( mnStochasticCharacterMap(ctmc=eulaemus_dec,
filename=out_fn+".stoch.log",
printgen=100) )¡A correr!
Definimos nuestro modelo
mymodel = model(eulaemus_dec)Creamos el objeto MCMC a partir de las variables del modelo, los movimientos, y los monitores.
mymcmc = mcmc(mymodel, monitors, moves, nruns=2, moveschedule="random")A continuación, se define una regla de parada (stopping rule) basada en el tamaño efectivo de muestra (ESS). Esta regla detiene automáticamente la cadena MCMC cuando el ESS alcanza un valor mínimo predefinido, lo que indica que el muestreo ha sido suficiente y que la cadena ha explorado adecuadamente el espacio de parámetros. En este caso, se establece un ESS mínimo de 250, el cual se evalúa periódicamente a partir del archivo de salida del modelo.
stopping_rules[1] = srMinESS(250, file = out_fn+".model.log", freq = 10000)Finalmente corremos el análisis MCMC
mymcmc.run(n_gen,checkpointInterval=1000,checkpointFile="output/eulaemus_Dec.state", rules = stopping_rules)Resultados
A continuación, obtenemos la reconstrucción de las áreas ancestrales a partir del archivo simple.states.log
Descargar estos archivos:
eulaemus_anc.tre: filogenia con estados biogeográficos ancestrales estimados. Filogenia_anc
eulaemus_DEC_states.log: Muestras posteriores de estados (rangos) geográficos ancestrales del modelo DEC. Datos_anc
tree = readTrees("data/eulaemus.tre")[1]
state_trace = readAncestralStateTrace(file="output/eulaemus_DEC.states_run_1.log")
anc_tree = ancestralStateTree(tree= tree,
ancestral_state_trace_vector=state_trace,
include_start_states=true,
file="output/eulaemus_ancstates_run_1.tre",
burnin=0.25,
site=1)El archivo generado mediante la función ancestralStateTree puede ser revisado en FigTree.
Visualización en R
Para visualizar la reconstrucción de estados ancestrales obtenida bajo el modelo DEC, utilizaremos el paquete de R RevGadgets.
En primer lugar, se especifica el archivo del árbol anotado y se define un diccionario (state_labels) que vincula los identificadores numéricos de estado con las combinaciones de áreas correspondientes. En este ejemplo, las áreas se definen como: A (Altiplano andino + desierto de Atacama), B (Chile central + Andes centrales chilenos), C (Patagonia) y D (tierras bajas orientales). Finalmente, se genera la figura con probabilidades posteriores por nodo y punta.
Para este tutorial necesitaremos los siguientes paquetes
library(RevGadgets)
library(ggplot2)Cargamos el Archivo del árbol anotado con estados ancestrales
dec_eulaemus_file = "output/eulaemus_ancstates_run_1.tre"Etiquetas de estados: números (codificación) -> rangos (combinación de áreas).
labs <-c("1"="A",
"2"="B",
"3"="C",
"4"="D",
"5"="AB",
"6"="AC",
"7"="BC",
"8"="AD",
"9"="BD",
"10"="CD"
)Leer y procesar el árbol con probabilidades posteriores de estados ancestrales
dec_eulaemus <- processAncStates(dec_eulaemus_file,
state_labels = labs)## | | | 0% | |========================================| 100%Podemos seguir modificando el codigo para por ejemplo generar figura con “pies” de probabilidad por nodo y por punta
pie <- plotAncStatesPie(
t = dec_eulaemus,
cladogenetic = TRUE, # incluir eventos cladogenéticos (si están anotados)
tip_labels_states = FALSE, # mostrar etiqueta del estado en las puntas
tip_labels_offset = 0.35, # mover etiquetas de especies para dar espacio al pie
tip_pie_nudge_x = 0.05, # desplazar pies de puntas hacia la derecha
tip_pie_size = 0.70,
node_pie_size = 1,
shoulder_pie_size = 0.75,
state_transparency = 0.85,
timeline = TRUE, # mostrar escala temporal si el árbol está fechado
geo = FALSE,
time_bars = TRUE
)
pie <- pie +
labs(fill = "State") +
theme(
legend.position = c(0.85, 0.5),
legend.title = element_text(size = 10),
legend.text = element_text(size = 8),
legend.key.size = unit(0.4, "cm"),
legend.background = element_rect(
fill = alpha("white", 0.8),
color = NA
)
)
pie
ggsave(
filename = "DEC_ancestral_states.png",
plot = pie,
width = 14, # ancho en pulgadas (clave para árboles)
height = 8, # alto proporcional
dpi = 300, # calidad publicación
bg = "white"
)
Referencias
Esquerré, D., Brennan, I. G., Catullo, R. A., Torres-Pérez, F., & Keogh, J. S. (2019). How mountains shape biodiversity: The role of the Andes in biogeography, diversification, and reproductive biology in South America’s most species-rich lizard radiation (Squamata: Liolaemidae). Evolution, 73(2), 214–230. https://doi.org/10.1111/evo.13657