Diversidad alfa: Shannon, Simpson, Chao1

r
bioconductor
metagenomica
Qué mide cada índice de diversidad alfa, cuándo aplicar cada uno, por qué rarefying es controvertido y cómo reportar diferencias entre grupos con tests apropiados.

Qué significa “diversidad alfa”

Diversidad alfa = diversidad dentro de una muestra. Es un escalar por muestra. Responde:

¿Cuán rica y variada es esta comunidad?

Distinto de diversidad beta, que mide cuán distintas son dos comunidades entre sí (lo veremos en el siguiente tutorial).

Tres conceptos que conviene tener separados:

Concepto Qué mide Índice típico
Riqueza Cuántas especies hay Observed, Chao1
Equidad Cómo de repartidas están Pielou’s J
Diversidad efectiva Riqueza + equidad en un solo número Shannon, Simpson

Un ejemplo intuitivo: comunidad A con 100 ASVs donde uno es el 95 % y los otros 99 reparten el 5 %. Comunidad B con 100 ASVs todos al 1 %. Misma riqueza, distinta equidad, distinta diversidad efectiva (B es más diversa).

Los índices que importan

Observed (riqueza)

El más simple: número de ASVs presentes (con conteo > 0) en la muestra.

  • Sensible al esfuerzo de muestreo: una muestra con 50k reads detecta más ASVs raros que una con 5k. La comparación entre muestras de distinta profundidad es sesgada.
  • Sirve como sanity check. no es robusto para comparar entre estudios.

Chao1 (riqueza estimada)

Estima la riqueza verdadera (incluidos los ASVs raros no detectados) basándose en cuántos ASVs aparecen una sola vez (singletons) vs dos veces (doubletons).

  • Más robusto que Observed.
  • Pero DADA2 elimina singletons por diseño (modelo de errores los considera ruido). Esto rompe Chao1 técnicamente, la fórmula clásica espera singletons reales. Es decir: Chao1 sobre datos DADA2 da números bajos. No es comparable a Chao1 sobre OTUs clásicos.

Shannon (diversidad efectiva, peso log)

Combina riqueza y equidad. Fórmula: H = −Σ pᵢ · ln(pᵢ).

  • Rango: 0 (una sola especie domina) a ln(S) (todas las especies igualmente abundantes).
  • Pesa más los ASVs raros que Simpson.
  • El más reportado en literatura.

Simpson (diversidad efectiva, peso por dominancia)

D = 1 − Σ pᵢ². Da más peso a las especies dominantes.

  • Rango: 0 a 1.
  • Menos sensible a ASVs raros que Shannon.
  • Algunos papers reportan inverso de Simpson (1/D), que da una “diversidad efectiva” interpretable como “número equivalente de especies igualmente abundantes”.

Pielou’s J (equidad pura)

J = H / ln(S). Diversidad Shannon dividida por el máximo teórico.

  • Rango: 0 a 1.
  • Mide solo equidad (qué tan repartidas están las abundancias).
  • Útil cuando quieres separar el efecto “más especies” del efecto “más equilibrado”.

Calcular alfa diversidad con phyloseq

Continuando desde el ps del tutorial 9:

library(phyloseq)
library(tidyverse)

ps <- readRDS("output/phyloseq.rds")

# Todas las métricas a la vez
alpha_div <- estimate_richness(
    ps,
    measures = c("Observed", "Chao1", "Shannon", "Simpson", "InvSimpson")
)

head(alpha_div)
#          Observed Chao1 se.chao1 Shannon Simpson InvSimpson
# SampleA       243 248.4    3.21    4.12   0.967      30.04
# SampleB       198 201.2    2.78    3.89   0.961      25.71
# ...

estimate_richness devuelve un data.frame con una columna por métrica. Útil para combinar con metadata:

alpha_df <- alpha_div |>
    rownames_to_column("sample") |>
    left_join(
        as.data.frame(sample_data(ps)) |> rownames_to_column("sample"),
        by = "sample"
    )

Ahora puedes comparar entre grupos.

Visualizar diferencias entre grupos

El plot canónico, boxplot + jitter + p-valor:

library(ggpubr)

ggplot(alpha_df, aes(x = condition, y = Shannon, fill = condition)) +
    geom_boxplot(alpha = 0.6, outlier.shape = NA) +
    geom_jitter(width = 0.15, alpha = 0.6, size = 1.4) +
    stat_compare_means(method = "wilcox.test", label = "p.format") +
    scale_fill_manual(values = c("#5F8575", "#A24B4B")) +
    labs(x = NULL, y = "Diversidad Shannon",
         title = "Diversidad alfa por condición") +
    theme_minimal(base_size = 12) +
    theme(legend.position = "none")

stat_compare_means añade automáticamente el p-valor del test Wilcoxon. Para más de dos grupos, usa Kruskal-Wallis (method = "kruskal.test") con post-hoc Dunn.

Por qué Wilcoxon y no t-test: la diversidad alfa rara vez sigue normalidad y los tamaños muestrales son modestos. Wilcoxon es la opción honesta por defecto.

Múltiples métricas en un panel

Para informe, varios índices a la vez:

library(tidyr)

alpha_long <- alpha_df |>
    select(sample, condition, Observed, Chao1, Shannon, InvSimpson) |>
    pivot_longer(cols = Observed:InvSimpson,
                 names_to = "metric", values_to = "value")

ggplot(alpha_long, aes(x = condition, y = value, fill = condition)) +
    geom_boxplot(alpha = 0.6, outlier.shape = NA) +
    geom_jitter(width = 0.15, alpha = 0.5, size = 1) +
    facet_wrap(~ metric, scales = "free_y") +
    scale_fill_manual(values = c("#5F8575", "#A24B4B")) +
    labs(x = NULL, y = NULL) +
    theme_minimal() +
    theme(legend.position = "none")

Cuatro paneles, una vista por métrica. Si las conclusiones son consistentes entre métricas, el efecto es robusto. Si solo aparece en una (típicamente Shannon), ten cuidado al reportar.

El problema del esfuerzo de muestreo

Aquí está el debate clásico del campo. Las muestras tienen profundidades distintas (5k-100k reads). ASVs raros aparecen solo en muestras con mucha profundidad. Esto sesga Observed y Chao1 hacia arriba en las muestras más secuenciadas.

Tres aproximaciones para corregirlo:

Opción A: rarefying (rarefacción)

Submuestrear cada muestra al mínimo común de reads:

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

Lo polémico: McMurdie & Holmes (2014, PLoS Comput Biol) demostraron que rarefying descarta información útil y reduce poder estadístico. Mucha gente lo abandonó.

Lo defendible: para diversidad alfa específicamente, sigue siendo lo más simple y comparable entre estudios. Willis (2019, Microbiome) y otros lo defienden con matices.

Opción B: coverage-based rarefaction

Submuestrear hasta una cobertura estimada de Good (proporción de la comunidad capturada), no a un número fijo de reads. Más justo cuando las muestras tienen diversidad muy distinta. Paquete iNEXT o breakaway.

Opción C: estimadores que corrigen sesgo

breakaway (Willis & Bunge 2015) estima riqueza con intervalos de confianza honestos sin necesidad de rarificar:

library(breakaway)
estimates <- breakaway(otu_table(ps))

Es el approach moderno conservador. Más complejo de explicar a un revisor.

Mi recomendación

  • Para Observed y Chao1: rarifica al mínimo común, o declara explícitamente que no rarificas y muestra que la profundidad no correlaciona con tu variable de interés.
  • Para Shannon y Simpson: el efecto del esfuerzo de muestreo es mucho menor (estos índices son robustos). Suelen poderse comparar sin rarificar siempre que las muestras tengan al menos ~ 5 000 reads.
  • Reporta la decisión en Métodos. Sin disculpas.

Reportar como un grown-up

En la sección Métodos del informe:

La diversidad alfa se calculó con phyloseq v1.46 usando los índices Observed, Chao1, Shannon y InvSimpson. Para evitar el sesgo por profundidad de secuenciación, las muestras se rarificaron al mínimo común (8 542 reads) antes del cálculo. Las diferencias entre grupos se evaluaron con el test no paramétrico de Wilcoxon-Mann-Whitney.

Y en Resultados:

La diversidad alfa de los individuos con obesidad fue significativamente menor que la del grupo control (Shannon mediana: 3.41 vs 4.08, p = 0.003. Wilcoxon).

Esto es lo que se publica.

¿“Más diverso” es “más sano”?

Es la falsa intuición más extendida. No:

  • Microbioma vaginal sano: dominado por Lactobacillus. Diversidad alfa muy baja. Es lo bueno.
  • Microbioma intestinal sano: diversidad alta. Asociada a salud.
  • Bocas adultas: diversidad alta variable, contexto-dependiente.

Lo correcto: la diversidad alfa óptima depende del nicho. Para intestino, suele asociarse a salud. Para vagina o piel, no necesariamente.

Tampoco es necesariamente direccional al envejecimiento: el microbioma del lactante es poco diverso (esperado). El del adulto sano es más diverso. El del centenario vuelve a perder diversidad. Diversidad es una métrica, no un valor.

El patrón idiomático

library(phyloseq); library(tidyverse); library(ggpubr)

ps <- readRDS("output/phyloseq.rds")

# Rarefacción opcional (decisión documentada)
set.seed(42)
ps_rare <- rarefy_even_depth(ps, sample.size = min(sample_sums(ps)),
                              replace = FALSE, rngseed = 42)

# Cálculo
alpha <- estimate_richness(ps_rare,
                            measures = c("Observed", "Shannon", "InvSimpson"))

alpha_df <- alpha |>
    rownames_to_column("sample") |>
    left_join(as.data.frame(sample_data(ps)) |> rownames_to_column("sample"),
              by = "sample")

# Estadística
wilcox <- wilcox.test(Shannon ~ condition, data = alpha_df)
broom::tidy(wilcox)

# Plot
ggplot(alpha_df, aes(condition, Shannon, fill = condition)) +
    geom_boxplot(alpha = 0.6) +
    geom_jitter(width = 0.15) +
    stat_compare_means(method = "wilcox.test") +
    theme_minimal()

Veinte líneas para la sección de diversidad alfa de un paper.

Trampas habituales

  • Comparar Observed entre muestras de profundidad muy distinta sin rarificar ni declarar. Sesgo seguro hacia las muestras más secuenciadas.
  • Reportar solo Shannon “porque es el más bonito”. Reporta 2-3 métricas. Si las conclusiones discrepan entre ellas, hay matiz biológico que investigar.
  • t-test sobre diversidad. Casi nunca cumple normalidad. Wilcoxon es la opción honesta.
  • “Más diverso = mejor” sin justificarlo. Depende del nicho biológico. Asegúrate de defender por qué la dirección que ves es interesante.
  • Cálculo sobre conteos transformados (CLR, relativa). Los índices de diversidad están definidos sobre conteos crudos. Calcula primero, transforma después si vas a otra cosa.

En la siguiente entrega

Tienes diversidad dentro de cada muestra. La siguiente entrega es diversidad entre muestras, distancias (Bray-Curtis, UniFrac), ordenaciones (PCoA, NMDS) y el test estadístico que va con ellas (PERMANOVA). El plot canónico de cualquier paper de microbioma. Lo siguiente.