Este tutorial te guía a través de los pasos para realizar un análisis filogenético espacial en R, utilizando el paquete de R phylospatial, el cual fue desarrollado recientemente. Este tutorial se basa en gran medida en las viñetas de este paquete, escritas por Matthew Kling. Puedes encontrar más detalles sobre las funciones y el paquete visitando el repositorio https://matthewkling.github.io/phylospatial/. Una de las novedades más importantes de este paquete es su capacidad para gestionar diferentes tipos de datos de ocurrencia, desde datos de presencia/ausencia, abundancias, hasta probabilidades de ocurrencia derivadas de modelos de distribución de especies (SDM).

Introducción

La diversidad filogenética (PD) mide el parentesco del conjunto de taxa presentes en una región. En términos prácticos, medir diversidad filogenética implica sumar las longitudes de las ramas que conectan un conjunto de taxones en una filogenia. Al igual que con otras medidas de diversidad, cuantificar la PD a lo largo del espacio geográfico nos permite comprender los patrones de diversidad. Por ejemplo, podemos abordar pregunatas como: ¿cuánto varía la PD entre los píxeles (subconjuntos de igual tamaño) de una región dada? ¿Existen píxeles con una PD significativamente mayor o menor? La respuesta a estas preguntas puede contribuir tanto a estudios sobre los procesos que generan dichos patrones como a aplicaciones directas en conservación.

La filogenética espacial

Si bien el estudio de la PD se ha utilizado en múltiples contextos, en este tutorial nos centraremos en un enfoque específico denominado “filogenética espacial”. Este enfoque se centra en el estudio de los patrones de PD y otras medidas derivadas de PD en el espacio geográfico, tratándolo dicho espacio como un conjunto de píxeles.

las briofitas de California

En este tutorial, utilizaremos un conjunto de datos de los musgos de California que es un subconjunto de los datos utilizados en la publicación reciente titulada “Spatial Phylogenetics with Continuous Data: an Application to California Bryophytes” donde estudiamos los patrones de diversidad filogenética de las briofitas en California, y desarrollamos algunos enfoques novedosos para analizar este tipo de datos.

Objetivos de aprendizaje

  • Comprender el tipo de preguntas que podemos responder con el enfoque de la filogenética espacial.
  • Aplicar análisis estándar del enfoque de la filogenética espacial en un ejemplo.
  • Familiarizarse con los datos, la estructura de los objetos, y el código para realizar análisis filogenéticos espaciales en R, utilizando el paquete ‘phylospatial’.

Manos a la obra!

Primero instalemos los paquetes que necesitamos

# Let's make sure we have the libraries we need
# if you haven't installed them, uncomment the lines below
# install.packages('phylospatial')
library(phylospatial)
library(terra)
library(ape)
library(sf)
library(tmap)

1) Los datos

Utilizaremos el conjunto de datos distribuido con la biblioteca phylospatial. Primero, leeremos una raster que contiene los datos de la comunidad.

moss_comm <- rast(system.file("extdata", "moss_comm.tif", package = "phylospatial"))

# checa la estructura de los datos descomentando y corriendo la linea de abajo
# moss_comm

# plotealo!
plot(moss_comm)

Ahora leemos la filogenia y la visualizamos.

moss_tree <- read.tree(system.file("extdata", "moss_tree.nex", package = "phylospatial"))

# checa la estructura de los datos descomentando y corriendo la linea de abajo
# moss_tree

# plotealo!
plot(moss_tree, type = "fan", cex = 0.15)

Ahora, usemos la función phylospatial() para crear un objeto filoespacial. Esta función requiere dos argumentos: el conjunto de datos de la comunidad y una filogenia de los terminales (usualmente especies) de dicho conjunto de datos.

ps <- phylospatial(moss_comm, moss_tree)
plot(ps, "comm")

Preguntas A
  1. ¿Qué tipo de objetos son comm, moss_tree y ps?
  2. ¿Cuál es la estructura del objeto filoespacial ps? Usa la consola para explorar qué información se almacena en este objeto.

Notas importantes sobre la estructura de los datos

  • Los objetos filoespaciales permiten el uso de datos de comunidad expresados como probabilidades (por ejemplo, el resultado de modelos de distribución de especies), abundancias, o datos binarios (presencia-ausencia). La función phylospatial() identificará el tipo de datos proporcionados, pero puedes especificarlo explícitamente declarando el argumento data_type= al construir el objeto filoespacial.
  • Ten en cuenta que los terminales del árbol filogenético y los taxones del conjunto de datos de la comunidad deben ser idénticos.

2) Patrones de diversidad filogenética alfa

Ahora que tenemos nuestro objeto del tipo spatiaphylo, el cual contiene información sobre la ocurrencia de taxones en cada píxel y una filogenia con información sobre las relaciones filogenéticas de dichos taxones, es natural que queramos saber cuál es la diversidad filogenética (PD) de cada píxel. El paquete phylospatial calcula todas las medidas de diversidad filogenética alfa con la función ps_diversity(). En el argumento metric= podemos especificar qué medidas queremos calcular. Hoy, calcularemos la diversidad de taxones (TD), la diversidad filogenética (PD) y el endemismo filogenético (PE) como se describe a continuación (pero ten en cuenta que existen muchas otras medidas que se pueden calcular; puedes explorarlas consultando la documentación de la función ?ps_diversity()):

  • TD: Diversidad Terminal, es decir, riqueza de taxones terminales.
  • TE: Endemismo Terminal, es decir, diversidad total ponderada por endemismo de taxones terminales, también conocido como “endemismo ponderado”.
  • PD: Diversidad Filogenética, es decir, longitud total de las ramas presentes en un sitio.
  • PE: Endemismo Filogenético, es decir, PD ponderada por endemismo.
# calculamos las medidad de diversidad filogenética alfa
div <- ps_diversity(ps, metric = c("TD","TE", "PD", "PE"))

Usamos la consola para explorar la estructura del objeto div. ¡Y vamos a plotearlo!

# plot diversity metrics

#PD
TD_plot <- tm_shape(div$TD) + 
      tm_raster(palette = "inferno", style = "cont") +
      tm_layout(legend.outside = TRUE) + tm_title("Taxon diversity, TD")

TE_plot <- tm_shape(div$TE) + 
      tm_raster(palette = "inferno", style = "cont") +
      tm_layout(legend.outside = TRUE) + tm_title("Taxon endemism, TE")

PD_plot <- tm_shape(div$PD) + 
      tm_raster(palette = "inferno", style = "cont") +
      tm_layout(legend.outside = TRUE) + tm_title("Phylogenetic diversity, PD")


PE_plot <- tm_shape(div$PE) + 
      tm_raster(palette = "inferno", style = "cont") +
      tm_layout(legend.outside = TRUE) + tm_title("Phylogenetic endemism, PE")

# plots
tmap_arrange(TD_plot, TE_plot, PD_plot, PE_plot, nrow = 2, ncol = 2)

Preguntas B
  1. ¿Cómo se interpreta cada medida?
  2. ¿Qué representan las escalas? Consulta la documentación de la función para encontrar las fórmulas para calcular cada medida.
  3. ¿Están relacionadas las diferentes medidas?

En la figura anterior, observamos que la diversidad filogenética está altamente correlacionada con la diversidad de taxones. Cuantos más terminales haya en un píxel, mayor será la longitud de las ramas que añadimos al calcular la PD. Si nos interesa saber si un píxel tiene ramas relativamente más largas, pero utilizando una correción del número de taxones en ese píxel, podemos utilizar “medidas relativas”. Visualicemos cómo se ven la PD y la diversidad filogenética relativa en nuestro conjunto de datos.

# calculamos pd y rpd
pd_rpd <- ps_diversity(ps, metric = c("PD","RPD"))

# las visualizamos

pd_plot  <- tm_shape(pd_rpd$PD) + 
            tm_raster(palette = "inferno", style = "cont") +
            tm_layout(legend.outside = TRUE) + tm_title("Phylogenetic diversity, PD")

rpd_plot <- tm_shape(pd_rpd$RPD) + 
            tm_raster(palette = "inferno", style = "cont") +
            tm_layout(legend.outside = TRUE) + tm_title("Relative phylogenetic diversity, RPD")

tmap_arrange(pd_plot, rpd_plot, ncol = 2)

Preguntas C
  1. ¿Qué diferencias observas entre los patrones espaciales de estas dos medidas?
  2. ¿Cómo interpretas estas diferencias?

3) Pruebas de significancia y análisis categórico del neoendemismo y paleoendemismo

Hasta ahora hemos explorado los patrones espaciales de las medidas de diversidad. Nuestras observaciones nos permiten comparar la diversidad entre píxeles. Sin embargo, probablemente queremos identificar regiones con valores significativamente diferentes a los esperados; es decir, queremos asignar un valor de significancia a nuestras observaciones. Para ello, utilizaremos una aleatorización. Una prueba de aleatorización implica comparar los datos observados con una expectativa. Por lo tanto, necesitaremos generar una expectativa de cuáles serían los valores de las métricas de diversidad para cada píxel. La función ps_rand() aplica internamente el algoritmo quantize() para generar una hipótesis nula basada en la ejecución de múltiples aleatorizaciones estratificadas.

Para ejecutar esta función, necesitamos proporcionarle como argumento un objeto filoespacial, indicar el número de aleatorizaciones que queremos usar para crear nuestra expectativa nula, y las medidas de diversidad que queremos usar. Ten en cuenta que este paso puede tardar un tiempo en correr.

# Corremos la randomizacion para nuestros datos, puedes cambiar el numero de randomizaciones modificando el argymento n_rand
rand <- ps_rand(ps, n_rand = 50, progress = F,
                metric = c("PD", "PE", "CE", "RPE"))

Si indagas el resultado de esta función, observarás que consta de un ráster por métrica. Grafiquemos uno de estos rásters.

tm_shape(rand$qPE) + 
      tm_raster(palette = "inferno", style = "cont") +
      tm_layout(legend.outside = TRUE)

Los valores graficados indican, para cada celda, la proporción de aleatorizaciones en las que el valor observado de la medida es mayor que el valor de la aleatorización. Por ejemplo, en la figura anterior, si una celda tiene un valor de 0.8, significa que en el 80 % de las aleatorizaciones el valor observado de PE fue mayor que el valor simulado. En ciencias, solemos utilizar valores de significancia de alfa = 0.5. Si aplicasemos este criterio en este escenario, solo consideraríamos a las celdas con valores superiores a 0.95 como celdas con PE significativamente mayor al esperado.

CANAPE

Los investigadores pueden estar interesados ​​en evaluar la importancia de cierta medidas específicas de diversidad filogenética, y pueden usar variantes del código anterior para ello. Sin embargo, un análisis que se ha convertido en estándar en filogenética espacial es CANAPE (análisis categórico de neoendemismo y paleoendemismo, descrito en Mishler et al. 2014). CANAPE consta de dos pasos. Primero, se identifican los píxeles con un endemismo filogenético (PE) significativamente alto. Y segundo, cada uno de estos píxeles con un PE significativamente alto se asigna a categorías de endemismo basadas en el RPE, que, como se mencionó anteriormente, refleja la longitud relativa de las ramas de un píxel en comparación con otros píxeles.

# run CANAPE
# nota que podemos modificar el valor de alpha
cp <- ps_canape(rand, alpha = .05)

# ploteamos los results
terra::plot(cp)

Las categorías de endemismo significativo a las que se pueden asignar los píxeles son:

  1. Neoendemismo: píxeles con un endemismo significativamente alto debido a ramas significativemente cortas, indicadas por RPE.

  2. Paleoendemismo: píxeles con un endemismo significativamente alto debido a ramas significativamente largas, indicadas por RPE.

  3. Endemismo mixto: píxeles con un endemismo significativamente alto, pero sin valores significativos en RPE, lo que se interpreta como una mezcla de ramas cortas y largas.

  4. Superendemismo: es una subcategoría del endemismo mixto con un endemismo significativo particularmente alto.

Preguntas D

  1. ¿Cuál es la importancia de usar pruebas de aleatorización?

  2. ¿Cómo interpreta cada tipo de categoría de endemismo? ¿Se te ocurre algún proceso evolutivo que pueda explicar el hallazgo de tales valores de significancia?

4) Patrones de diversidad filogenética beta

Otro aspecto interesante de los estudios de diversidad es comprender la distribución espacial de los ensamblajes de especies, respondiendo a preguntas como ¿qué tan similares son las comunidades a lo largo del espacio geográfico? ¿Existe un recambio evidente en la composición de las comunidades? También abordaremos esta cuestión desde la perspectiva de la diversidad filogenética. Para ello, comenzaremos calculando un índice de disimilitud filogenética para cada par de pixeles, utilizando la distancia de Bray-Curtis, también conocida como disimilitud cuantitativa de Sorensen. La función ps_add_dissim() añadirá una matriz de disimilitud a nuestro objeto filoespacial.

ps <- ps_add_dissim(ps, method = "sorensen", endemism = TRUE, normalize = TRUE)
ps
## `phylospatial` object
##   - 884 lineages across 1116 sites
##   - community data type: probability 
##   - spatial data class: SpatRaster 
##   - dissimilarity data: sorensen

En el código anterior, al indicar endemism = TRUE, los valores de ocurrencia de cada linaje se rescalan para que sumen 1 en todos los sitios, otorgando mayor peso a los taxones con distribución estrecha. Al establecer normalize = TRUE, el valor de ocurrencia total de cada sitio se escala para que sume 1 en todos los taxones, lo que genera una matriz de distancia que enfatiza las diferencias proporcionales en la composición en lugar de los gradientes de diversidad alfa.

Una vez calculada la matriz de disimilitud filogenética, podemos visualizar la similitud filogenética entre píxeles. Para ello, utilizaremos la función ps_rbg(). Internamente, esta función realizará un análisis de ordenación (en este caso, utilizaremos un CMDS: escalamiento multidimensional clásico) que reduce la variación en las comunidades a tres componentes, que luego se utilizan para definir los colores (utilizando las bandas RBG del espacio de color). Como resultado, cuanto más similares sean los colores, más similares filogenéticamente serán las comunidades.

ps %>% # toma el objeto filoespacial
      ps_rgb(method = "cmds") %>% # ddeclaramos el algoritmo de ordenacion
      tm_shape() + # hacemos un mapa
      tm_rgb(max.value = 1, interpolate = FALSE) +
      tm_layout(title = "Phylogenetic community ordination")

Si bien esta visualización ayuda a observar la similitud, puede resultar difícil decidir dónde dividir el mapa si intentamos realizar una regionalización geográfica. Para ello, podemos utilizar un análisis de agrupamiento (clustering) con la función ps_regions(). Esta función requiere que decidamos el número de regiones en las que queremos agrupar el espacio. Si bien no existe una respuesta única para definir el número de regiones y la decisión tiene aspectos algo arbitrarios, un enfoque común en los métodos de agrupamiento consiste en seleccionar la menor cantidad de categorías k que expliquen la mayor parte de la variación de los datos. Para determinar el valor de k, utilizaremos la función ps_regions_eval().

ps_regions_eval(ps, k = 1:15, plot = TRUE, method = "average")

El gráfico anterior sugiere que K=3 es un valor adecuado para nuestro análisis de agrupamiento. Así pues, ahora podemos ejecutar una regionalización con k=3 y, posteriormente, representarla gráficamente.

ps %>%
      ps_regions(k = 3, method = "average") %>% # regionalization
      tm_shape() + # plot
      tm_raster(style = "cat", palette = "Dark2",
                title = "phylogenetic region") +
      tm_layout(legend.outside = TRUE)

Es importante mencionar que existen múltiples algoritmos que se pueden utilizar para generar los clusters, y la regionalización puede variar dependiendo de cuál algoritmo se seleccione, recomiendo leer más sobre los métodos de clusterización antes de tomar una decisión sobre qué algoritmo utilizar.

Preguntas E
  1. Experimenta con los diferentes valores de K. Dado que diferentes algoritmos de regionalización pueden generar resultados distintos, ¿cómo abordarías la regionalización de tu sistema de estudio?

NOTA FINAL

Espero que este tutorial se ade ayuda para comenzar a utiliar este paquete de R. Este tutorial no es exhaustivo, ya que hay muchas otras funcionalidades por explorar. Les animo a checar las viñetas del paquete para aprender más. En particular, un aspecto que podría ser de interés para muchos es la aplicación en conservación con el algoritmo de priorización incluido en el paquete.