Diversidad beta y ordenaciones (PCoA, NMDS, UniFrac)

r
bioconductor
metagenomica
Cómo se compara la composición entre comunidades. Distancias Bray-Curtis y UniFrac (weighted/unweighted), ordenaciones PCoA y NMDS, el test PERMANOVA con adonis2 y por qué dispersion puede confundir tu interpretación.

Qué mide la diversidad beta

Diversidad beta = diversidad entre muestras. Responde:

¿Cuán parecidas son dos comunidades microbianas?

A partir de pares de muestras se construye una matriz de distancia (N × N). Esa matriz es la base de:

  • Ordenaciones (PCoA, NMDS): proyectar la matriz a 2-3 dimensiones para visualizar.
  • Tests estadísticos (PERMANOVA): cuantificar si grupos de muestras tienen composiciones distintas.
  • Dendrogramas, heatmaps de distancia: vistas alternativas.

Las distancias que importan

Hay docenas de distancias para microbioma. Cuatro son las que verás casi siempre:

Distancia Considera abundancia Considera filogenia Cuándo
Jaccard No (binaria) No Presencia/ausencia, lo más simple
Bray-Curtis Sí No Default ecológico. Lo más reportado
UniFrac unweighted No (binaria) Sí Filogenia con énfasis en ASVs raros
UniFrac weighted Sí Sí Filogenia con énfasis en abundantes

Reglas rápidas:

  • Bray-Curtis para casi todo. Es el default robusto.
  • UniFrac weighted cuando tu pregunta involucra abundancias y quieres incorporar parentesco filogenético (ej. dos Bacteroides distintos cuentan como “similares” porque son parientes).
  • UniFrac unweighted cuando quieres detectar diferencias en presencia/ausencia de linajes (más sensible a ASVs raros pero relevantes biológicamente).
  • Jaccard raramente la primera elección, pero útil como sanity check.

Antes de las distancias: transformar

Las distancias se calculan sobre conteos. Pero conteos crudos sufren el sesgo de profundidad. Tres aproximaciones:

Opción A: Rarificar y trabajar con conteos rarificados

set.seed(42)
ps_rare <- rarefy_even_depth(ps, sample.size = min(sample_sums(ps)),
                              replace = FALSE, rngseed = 42)

Lo más común en literatura. Polémico (ver tutorial anterior) pero defendible.

Opción B: Abundancia relativa

ps_rel <- transform_sample_counts(ps, function(x) x / sum(x))

Compatible con Bray-Curtis. No con UniFrac (espera conteos).

Opción C: CLR

library(microbiome)
ps_clr <- microbiome::transform(ps, "clr")

Composicional correcta (lo veremos a fondo en el tutorial 12). Sobre CLR, la distancia adecuada es euclídea (vegan::vegdist(..., method = "euclidean")), llamada distancia de Aitchison.

Para este tutorial: rarificar es lo más sencillo y didáctico. Para análisis serio del tutorial 12 usaremos CLR + Aitchison.

Calcular distancias

phyloseq tiene una función única para todas:

library(phyloseq)

ps_rare <- readRDS("output/ps_rare.rds")

# Bray-Curtis (sobre conteos rarificados)
dist_bray <- distance(ps_rare, method = "bray")

# Jaccard
dist_jaccard <- distance(ps_rare, method = "jaccard", binary = TRUE)

# UniFrac (necesita árbol filogenético en el objeto)
dist_uu <- distance(ps_rare, method = "unifrac", weighted = FALSE)
dist_uw <- distance(ps_rare, method = "wunifrac")

Cada distance() devuelve un objeto dist (matriz triangular). N muestras → N×(N−1)/2 distancias.

Comprobación: las distancias deben estar en [0, 1] para Bray-Curtis y Jaccard, [0, 1] (típicamente menos) para UniFrac.

Ordenaciones: PCoA y NMDS

Para visualizar la matriz de distancia en 2D.

PCoA (Principal Coordinates Analysis)

Análogo a PCA pero sobre una matriz de distancia arbitraria. Lineal, captura varianza explicada en cada eje:

ord_pcoa <- ordinate(ps_rare, method = "PCoA", distance = "bray")

plot_ordination(ps_rare, ord_pcoa, color = "condition") +
    geom_point(size = 4, alpha = 0.85) +
    stat_ellipse(aes(group = condition), level = 0.8) +
    labs(title = "PCoA — distancia Bray-Curtis") +
    theme_minimal()

Cada eje muestra el % de varianza que captura. Si PC1 captura > 30 %, hay una fuente de variación dominante (esperable: composición global). Si los grupos se separan a lo largo de un eje, lo señalan.

NMDS (Non-metric Multidimensional Scaling)

No-lineal. Preserva el orden de las distancias en lugar de los valores exactos. Mejor cuando hay no-linealidad:

ord_nmds <- ordinate(ps_rare, method = "NMDS", distance = "bray")
ord_nmds$stress  # < 0.2 es aceptable, < 0.1 es bueno

plot_ordination(ps_rare, ord_nmds, color = "condition") +
    geom_point(size = 4, alpha = 0.85) +
    stat_ellipse(aes(group = condition), level = 0.8) +
    labs(title = "NMDS — distancia Bray-Curtis",
         subtitle = paste0("Stress = ", round(ord_nmds$stress, 3))) +
    theme_minimal()

Stress es la métrica clave de NMDS:

  • < 0.1: excelente representación.
  • 0.1-0.2: aceptable, interpretable con cuidado.
  • 0.2-0.3: cuestionable.

  • 0.3: la representación no es fiable.

PCoA vs NMDS: ¿cuál?

  • PCoA: lineal, te dice cuánta varianza captura cada eje. Default para reportar.
  • NMDS: no-lineal, mejor para patrones complejos. Pero el “% de varianza” no aplica.

Para informe rápido, PCoA. Si el patrón es claro, mantén PCoA. Si NMDS revela algo que PCoA no, reporta ambos.

PERMANOVA: el test estadístico

Mirar el plot y decir “se separan” es subjetivo. Para cuantificar:

PERMANOVA (Permutational Multivariate ANOVA): test no paramétrico que evalúa si los centroides de varios grupos en una matriz de distancia son distintos.

library(vegan)

# Convertir sample_data a data.frame
metadata <- as(sample_data(ps_rare), "data.frame")

# adonis2 = PERMANOVA con la API moderna
adonis_result <- adonis2(
    dist_bray ~ condition + batch,
    data         = metadata,
    permutations = 999,
    by           = "margin"
)

adonis_result
#         Df SumOfSqs      R2      F   Pr(>F)
# condition 1   1.234   0.082  6.21    0.001 ***
# batch     2   0.421   0.028  1.06    0.380
# Residual 62  12.345   0.890
# Total    65  13.99    1.000

Lectura:

  • R²: fracción de la varianza explicada por la variable. Un R² de 0.08 significa que la condición explica el 8 % de la variación total, modesto pero significativo si p < 0.05.
  • Pr(>F): p-valor del test (basado en permutaciones).
  • by = "margin": ajusta cada variable controlando por las demás (Tipo III). Para diseños con covariables, es lo correcto.

El sanity check: betadisper

PERMANOVA detecta diferencias entre centroides. Pero también puede dar significativo si los grupos tienen dispersión muy distinta dentro de sí (un grupo muy concentrado vs otro muy disperso). Eso confunde la interpretación.

betadisper testa la homogeneidad de dispersiones:

bd <- betadisper(dist_bray, metadata$condition)
anova(bd)
#               Df  Sum Sq Mean Sq F value Pr(>F)
# Groups         1 0.00210 0.00210   0.342  0.561
# Residuals     64 0.39231 0.00613

Si p > 0.05: dispersiones homogéneas → PERMANOVA mide diferencias reales de centroide. Bien.

Si p < 0.05: dispersiones distintas → la “significación” de PERMANOVA puede deberse a dispersión, no a desplazamiento. Interpretar con cuidado.

Lo idiomático: reportar ambos:

Las muestras de obesos difirieron significativamente de los controles (PERMANOVA, R² = 0.082, p = 0.001) sin diferencias significativas de dispersión beta (betadisper, p = 0.561).

Ordenación con vectores de variables

Para entender qué variables impulsan la ordenación, puedes superponer envfit:

library(vegan)

ef <- envfit(ord_pcoa$vectors ~ condition + age + bmi,
             data = metadata,
             permutations = 999)
ef
# ***VECTORS / FACTORS
# condition  R² = 0.082, p = 0.001
# age        R² = 0.031, p = 0.067
# bmi        R² = 0.054, p = 0.015

Y para añadir flechas al plot:

plot(ord_pcoa$vectors[, 1:2])
plot(ef, add = TRUE)

Útil para informe. La API de phyloseq + ggplot es algo más verbose para esto. Vegan-base hace el trabajo.

Heatmap de distancias entre muestras

Para datasets pequeños (< 30 muestras), un heatmap de la matriz de distancia es informativo:

library(pheatmap)

mat <- as.matrix(dist_bray)
rownames(mat) <- paste(metadata$condition, metadata$batch, sep = " · ")
colnames(mat) <- NULL

pheatmap(
    mat,
    clustering_distance_rows = dist_bray,
    clustering_distance_cols = dist_bray,
    color = colorRampPalette(c("#3F6354", "white"))(50)
)

Bloques diagonales = grupos cohesivos. Una muestra que se aleja de las demás (línea distinta) = outlier candidato.

El patrón idiomático

library(phyloseq); library(vegan); library(tidyverse)

ps_rare <- readRDS("output/ps_rare.rds")
metadata <- as(sample_data(ps_rare), "data.frame")

# Distancias
dist_bray <- distance(ps_rare, method = "bray")
dist_wuf  <- distance(ps_rare, method = "wunifrac")

# Ordenación
ord <- ordinate(ps_rare, method = "PCoA", distance = "bray")

# Plot
p <- plot_ordination(ps_rare, ord, color = "condition") +
    geom_point(size = 4, alpha = 0.85) +
    stat_ellipse(aes(group = condition), level = 0.8) +
    scale_color_manual(values = c("#5F8575", "#A24B4B")) +
    theme_minimal()

# Estadística
ad <- adonis2(dist_bray ~ condition + batch,
              data = metadata, permutations = 999, by = "margin")
bd <- betadisper(dist_bray, metadata$condition)
ad
anova(bd)

Lo que se incluye en cualquier informe de microbioma.

Trampas habituales

  • Comparar UniFrac entre estudios que usaron árboles distintos. UniFrac depende del árbol. Con árboles diferentes (regiones distintas, métodos distintos) las distancias no son comparables.
  • Interpretar significación de PERMANOVA sin betadisper. Si las dispersiones son distintas, la diferencia “entre grupos” puede ser dispersión y no desplazamiento. Reporta los dos.
  • NMDS con stress > 0.2 y leerlo como un PCoA. NMDS preserva orden, no distancias absolutas. Si el stress es alto, la separación visual puede engañar.
  • Olvidar permutaciones suficientes. permutations = 99 es demasiado bajo para p-valores < 0.01. Default 999 para tests publicables.
  • No considerar batch como covariable. PERMANOVA solo en ~ condition cuando hay batch técnico oculto invariablemente confunde. Modela el batch.

Cierre del Bloque 3

Has cubierto el análisis exploratorio canónico:

  • Composición: plot_bar, top N géneros.
  • Diversidad alfa: Shannon, Simpson, Chao1.
  • Diversidad beta: Bray-Curtis / UniFrac + PCoA / NMDS + PERMANOVA.

Lo que cualquier informe de microbioma incluye en su sección Resultados.

En la siguiente entrega

Hasta ahora has descrito qué microbioma hay. La pregunta final del análisis: ¿qué taxones son distintos entre grupos? Lo siguiente es el análisis diferencial, pero adaptado al microbioma (que es composicional, lo que rompe asunciones de DESeq2). Aprenderás ANCOM-BC y MaAsLin2, los métodos modernos. Y lo aplicaremos a un dataset real (Bacteroides 2 enterotype / obesidad) en un informe Quarto completo. Lo siguiente.