Caso completo: ANCOM-BC + informe Quarto

r
bioconductor
metagenomica
Un análisis 16S de principio a fin sobre el dataset de obesidad y Bacteroides 2 enterotype. DADA2 → phyloseq → diversidad → análisis diferencial composicional con ANCOM-BC y MaAsLin2 → informe Quarto publicable.

El caso

Cerramos la ruta con un análisis real de extremo a extremo. Pregunta biológica:

¿Difiere la composición del microbioma intestinal entre individuos con obesidad y controles, y qué taxones explican la diferencia?

Dataset: estudio de intervención dietética con foco en el enterotipo Bacteroides 2 (Vandeputte et al. 2019, Gut). 16S V3-V4 sobre heces, 66 muestras (obesidad + controles), metadata con BMI, edad, sexo, batch técnico.

Acceso público: SRA BioProject (búsqueda por keywords “Bacteroides 2 enterotype dietary fiber”). El dataset está reanalizado en múltiples papers, usable como referencia.

Este tutorial es el más largo de la ruta (25-30 minutos de lectura, más al ejecutarlo). Reúne todo lo anterior en un caso aplicado, con la complejidad estadística que un análisis real exige.

Estructura del proyecto

caso-obesidad/
├── caso-obesidad.Rproj
├── data/
│   ├── fastq/             ← reads paired-end
│   └── metadata.csv
├── ref/
│   ├── silva_nr99_v138.1_*.fa.gz
│   └── silva_species_assignment_v138.1.fa.gz
├── R/
│   └── funciones.R
├── output/
│   └── *.rds              ← objetos intermedios
└── informe.qmd            ← el documento Quarto final

1. Setup e importación

library(dada2)
library(phyloseq)
library(microbiome)
library(ANCOMBC)
library(Maaslin2)
library(tidyverse)
library(ggpubr)

# Cargar metadata
metadata <- read.csv("data/metadata.csv", row.names = 1)

# Variables clave
metadata$condition <- factor(metadata$condition, levels = c("control", "obese"))
metadata$batch <- factor(metadata$batch)

table(metadata$condition, metadata$batch)
#         A  B  C
# control 12 11 10
# obese   12 11 10

Buen diseño: condiciones balanceadas entre batches. Si no lo fuera, la covariable batch quedaría confundida con condition y limitaría el análisis.

2. Pipeline DADA2 (sumario)

Asumiendo que ejecutaste los tutoriales 5-8 sobre estos FASTQ:

# Filtrado (tutorial 5)
out <- filterAndTrim(
    fnFs, filtFs, fnRs, filtRs,
    truncLen = c(240, 200),
    maxEE = c(2, 2), truncQ = 2, rm.phix = TRUE,
    compress = TRUE, multithread = TRUE
)

# Errores + denoising (tutorial 6)
errF <- learnErrors(filtFs, multithread = TRUE)
errR <- learnErrors(filtRs, multithread = TRUE)
dadaFs <- dada(filtFs, err = errF, pool = "pseudo", multithread = TRUE)
dadaRs <- dada(filtRs, err = errR, pool = "pseudo", multithread = TRUE)

# Merge + quimeras (tutorial 7)
mergers <- mergePairs(dadaFs, filtFs, dadaRs, filtRs)
seqtab <- makeSequenceTable(mergers)
seqtab <- seqtab[, nchar(colnames(seqtab)) %in% 400:460]
seqtab_nochim <- removeBimeraDenovo(seqtab, method = "consensus",
                                      multithread = TRUE)

# Taxonomía (tutorial 8)
taxa <- assignTaxonomy(seqtab_nochim, "ref/silva_nr99_v138.1_wSpecies_train_set.fa.gz",
                       multithread = TRUE)
taxa <- addSpecies(taxa, "ref/silva_species_assignment_v138.1.fa.gz")

# Filtrar a bacterias
keep <- !is.na(taxa[, "Kingdom"]) & taxa[, "Kingdom"] == "Bacteria"
seqtab_nochim <- seqtab_nochim[, keep]
taxa <- taxa[keep, ]

saveRDS(seqtab_nochim, "output/seqtab_final.rds")
saveRDS(taxa, "output/taxa.rds")

Output del tracking: ~88 % retención global, 2 114 ASVs finales en 66 muestras. Razonable para heces humanas.

3. Construir phyloseq y QC final

ps <- phyloseq(
    otu_table(seqtab_nochim, taxa_are_rows = FALSE),
    sample_data(metadata),
    tax_table(taxa)
)

# Renombrar ASVs a IDs cortos (tutorial 9)
dna <- Biostrings::DNAStringSet(taxa_names(ps))
names(dna) <- paste0("ASV", sprintf("%04d", seq_along(dna)))
taxa_names(ps) <- names(dna)
ps@refseq <- dna

# Filtro mínimo: ASV en ≥ 3 muestras con > 5 reads
ps <- filter_taxa(ps, function(x) sum(x > 5) >= 3, prune = TRUE)
ntaxa(ps)
# [1] 1247

# Library size por muestra
qplot(sample_sums(ps), bins = 30) +
    labs(x = "Reads por muestra (post-filtro)", y = "Cuenta") +
    theme_minimal()

Rango típico de 7 000-80 000 reads por muestra. Sin muestras catastróficamente bajas.

4. Composición taxonómica

Vista resumida a nivel de phylum:

ps_rel <- transform_sample_counts(ps, function(x) x / sum(x))
ps_phylum <- tax_glom(ps_rel, taxrank = "Phylum", NArm = TRUE)

plot_bar(ps_phylum, x = "sample", fill = "Phylum") +
    facet_wrap(~ condition, scales = "free_x") +
    scale_fill_brewer(palette = "Set2") +
    labs(x = NULL, y = "Abundancia relativa",
         title = "Composición a nivel de phylum") +
    theme_minimal(base_size = 11) +
    theme(axis.text.x = element_blank(),
          axis.ticks.x = element_blank())

Lo esperable: Bacillota (Firmicutes) y Bacteroidota dominan. Ratio F/B (Firmicutes/Bacteroidetes) suele reportarse en literatura de obesidad, aunque su valor como biomarcador es debatido.

5. Diversidad alfa: ¿hay menos diversidad en obesidad?

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

alpha <- estimate_richness(ps_rare,
                            measures = c("Observed", "Shannon", "InvSimpson"))
alpha_df <- alpha |>
    tibble::rownames_to_column("sample") |>
    left_join(as.data.frame(sample_data(ps)) |>
                  tibble::rownames_to_column("sample"),
              by = "sample")

ggplot(alpha_df, aes(condition, Shannon, fill = condition)) +
    geom_boxplot(alpha = 0.6, outlier.shape = NA) +
    geom_jitter(width = 0.15, alpha = 0.6) +
    stat_compare_means(method = "wilcox.test") +
    scale_fill_manual(values = c("#5F8575", "#A24B4B")) +
    labs(x = NULL, y = "Diversidad Shannon") +
    theme_minimal() +
    theme(legend.position = "none")

Resultado típico de este dataset: diversidad significativamente menor en obesos (Shannon, p ≈ 0.01). Consistente con literatura.

6. Diversidad beta: ¿se separan los grupos?

library(vegan)

ord_bray <- ordinate(ps_rare, method = "PCoA", distance = "bray")
ord_wuf  <- ordinate(ps_rare, method = "PCoA", distance = "wunifrac")

p1 <- plot_ordination(ps_rare, ord_bray, color = "condition") +
    geom_point(size = 4, alpha = 0.85) +
    stat_ellipse(level = 0.8) +
    scale_color_manual(values = c("#5F8575", "#A24B4B")) +
    labs(title = "PCoA — Bray-Curtis") +
    theme_minimal()

p2 <- plot_ordination(ps_rare, ord_wuf, color = "condition") +
    geom_point(size = 4, alpha = 0.85) +
    stat_ellipse(level = 0.8) +
    scale_color_manual(values = c("#5F8575", "#A24B4B")) +
    labs(title = "PCoA — Weighted UniFrac") +
    theme_minimal()

library(patchwork)
p1 + p2

# PERMANOVA + betadisper
dist_bray <- distance(ps_rare, method = "bray")
metadata_df <- as(sample_data(ps_rare), "data.frame")

ad <- adonis2(dist_bray ~ condition + batch + age + bmi,
              data = metadata_df, permutations = 999, by = "margin")
ad
#         Df SumOfSqs    R2     F  Pr(>F)
# condition 1   1.421 0.094  6.92  0.001 ***
# batch     2   0.412 0.027  1.00  0.434
# age       1   0.198 0.013  0.96  0.460
# bmi       1   0.621 0.041  3.02  0.012 *
# Residual 60  12.31  0.815
# Total    65  15.10  1.000

bd <- betadisper(dist_bray, metadata_df$condition)
anova(bd)
# p = 0.621 — dispersiones homogéneas

Reportable:

La composición del microbioma difiere significativamente entre obesos y controles (PERMANOVA Bray-Curtis, R² = 0.094, p = 0.001), con dispersiones beta homogéneas (betadisper p = 0.62). El BMI también contribuye independientemente (R² = 0.041, p = 0.012).

7. El problema composicional

Hasta aquí, exploratorio. Para preguntar qué taxones difieren entre grupos hay que decidir el método estadístico.

¿Por qué no DESeq2? Aunque comparta principios con el análisis diferencial de RNA-seq, DESeq2 asume independencia entre features. En microbioma eso es falso: los datos son composicionales, el total siempre suma 1. Si una bacteria sube en abundancia relativa, las demás bajan aunque su número absoluto no haya cambiado. Aplicar DESeq2 ingenuamente da falsos positivos sistemáticos.

Métodos modernos diseñados para esto:

Método Idea Cuándo
ANCOM-BC Test analítico con corrección por bias muestral Default robusto, dos grupos
MaAsLin2 Modelo lineal mixto sobre log-transform con covariables Datos longitudinales, varias covariables
ALDEx2 Inferencia bayesiana sobre CLR Conservador, intervalos honestos
LinDA Modelo lineal sobre CLR Rápido, bien calibrado

Vamos a usar ANCOM-BC (más extendido) y MaAsLin2 (más flexible con covariables).

8. Análisis diferencial con ANCOM-BC

ANCOM-BC trabaja directamente sobre el phyloseq:

library(ANCOMBC)

# Aglomerar a género (más interpretable que ASV)
ps_genus <- tax_glom(ps, taxrank = "Genus", NArm = TRUE)

# Ejecutar ANCOM-BC
ancom_out <- ancombc2(
    data         = ps_genus,
    assay_name   = "counts",
    tax_level    = "Genus",
    fix_formula  = "condition + batch + age + bmi",
    p_adj_method = "fdr",
    prv_cut      = 0.10,           # ASV en ≥ 10 % de muestras
    lib_cut      = 1000,            # mínimo de reads por muestra
    group        = "condition",
    struc_zero   = TRUE,
    neg_lb       = TRUE,
    n_cl         = 4
)

Los argumentos clave:

  • fix_formula: la fórmula del modelo. Como en DESeq2, la variable de interés (condition) más covariables a ajustar.
  • prv_cut: prevalencia mínima (porcentaje de muestras con count > 0). Filtro de baja prevalencia.
  • lib_cut: muestras con < 1000 reads se descartan.
  • struc_zero = TRUE: detecta y modela ceros estructurales (taxones que están realmente ausentes en un grupo entero).

Interpretar el output de ANCOM-BC

res <- ancom_out$res
head(res)
#       taxon            lfc_(Intercept)  lfc_conditionobese  se_conditionobese
#  ASV0042  Bacteroides            3.42                  1.21               0.34
#  ...

# Filtrar significantes
sig <- res |>
    filter(diff_conditionobese == TRUE) |>
    arrange(p_conditionobese)

nrow(sig)
# [1] 18

diff_conditionobese == TRUE significa que el género difiere significativamente entre obesos y controles tras corrección FDR.

Plotear los hallazgos

sig_top <- sig |>
    arrange(desc(abs(lfc_conditionobese))) |>
    head(20)

ggplot(sig_top, aes(x = lfc_conditionobese,
                    y = reorder(taxon, lfc_conditionobese),
                    fill = lfc_conditionobese > 0)) +
    geom_col() +
    scale_fill_manual(values = c("TRUE" = "#A24B4B", "FALSE" = "#5F8575"),
                      labels = c("TRUE" = "Up en obesos", "FALSE" = "Down en obesos")) +
    labs(x = "Log fold change (obesos vs controles)",
         y = NULL, fill = NULL) +
    theme_minimal()

Lo que esperas ver: Bacteroides y Faecalibacterium menos abundantes en obesos. Prevotella y géneros pro-inflamatorios más abundantes (patrón consistente con literatura).

9. MaAsLin2 con covariables

ANCOM-BC funciona con dos grupos. Para modelar varias covariables continuas y categóricas, MaAsLin2 es más flexible:

library(Maaslin2)

# Preparar inputs
abundance <- as.data.frame(t(otu_table(ps_genus)))
rownames(abundance) <- as.data.frame(tax_table(ps_genus))$Genus

meta_df <- as.data.frame(sample_data(ps))

# Ejecutar
maaslin_out <- Maaslin2(
    input_data    = abundance,
    input_metadata = meta_df,
    output         = "output/maaslin2_obesidad",
    fixed_effects  = c("condition", "age", "bmi", "batch"),
    reference      = c("condition,control"),
    normalization  = "TSS",
    transform      = "LOG",
    analysis_method = "LM",
    min_prevalence = 0.10
)

MaAsLin2 crea automáticamente:

  • significant_results.tsv: taxones con efectos significativos en cualquier variable.
  • all_results.tsv: todos los resultados con coeficientes y p-valores.
  • Plots de boxplot para los más significativos.

Cargar resultados:

results <- read_tsv("output/maaslin2_obesidad/significant_results.tsv")
results |> filter(metadata == "condition") |> arrange(qval)

10. Reconciliación entre métodos

Aquí está el matiz importante: ANCOM-BC y MaAsLin2 pueden dar resultados ligeramente distintos. Esperado, no problema. La buena práctica:

  • Reporta hallazgos consistentes entre métodos como conclusiones robustas.
  • Hallazgos solo en uno pero biológicamente plausibles: menciona como sugestivos.
  • Hallazgos discrepantes: investiga (¿muestra dominada por un outlier? ¿taxón raro?).
ancom_sig <- sig$taxon

maaslin_sig <- results |>
    filter(metadata == "condition", qval < 0.05) |>
    pull(feature)

intersect(ancom_sig, maaslin_sig)
# Los que ambos métodos coinciden — los hallazgos robustos

11. El informe Quarto

Toda esta historia se empaqueta en un único .qmd para presentación. Estructura típica:

---
title: "Análisis diferencial del microbioma en obesidad"
author: "Tu nombre"
date: today
format:
  html:
    toc: true
    code-fold: true
execute:
  echo: false
  warning: false
bibliography: refs.bib
---

## Resumen ejecutivo

Análisis del microbioma intestinal de 66 individuos (33 obesos, 33 controles)
mediante secuenciación 16S V3-V4. Encontramos diversidad alfa reducida en
obesos y diferencias significativas en 18 géneros tras corrección FDR.

## Métodos

Procesamiento DADA2 v1.30, SILVA 138.1, ANCOM-BC v2 y MaAsLin2 v1.18 con R 4.4.
Código y datos en [repositorio].

## Resultados

[Composición, diversidad alfa, diversidad beta, análisis diferencial]

## Discusión

[Comparación con literatura, limitaciones, próximos pasos]

## sessionInfo

Y al final del documento:

#| label: session-info
#| echo: true

sessionInfo()

sessionInfo() es lo que defiende tu análisis ante un reviewer. Captura versiones exactas de todos los paquetes.

12. Lo que aprendes resolviendo este caso

Hablando claro:

  • DADA2 no es trivial. Las decisiones de truncado y los parámetros importan. Diagnostica con plotErrors siempre.
  • El composicional no es opcional. Si reportas análisis diferencial sin método composicional, los reviewers te lo van a pedir.
  • Diversidad alfa baja en obesidad es robusta entre estudios. Diferencias específicas de taxones, menos.
  • La variabilidad inter-individuo en heces humanas es enorme. R² de 0.05-0.15 para condición es lo esperable, no decepcionante.
  • Una sola cohorte rara vez es definitiva. Los hallazgos reproducibles entre cohortes (HMP + tu dataset + Vandeputte) son la verdadera evidencia.

Has terminado la ruta

Con esto cierras Metagenómica con DADA2 y phyloseq. Has visto:

  • Conceptos: 16S vs shotgun, regiones variables, bases públicas.
  • Pipeline DADA2 completo: QC, errores, denoising, merge, quimeras, taxonomía SILVA.
  • phyloseq como contenedor y composición taxonómica.
  • Diversidad alfa (Shannon, Simpson, Chao1) y beta (Bray-Curtis, UniFrac, PCoA, PERMANOVA).
  • Análisis diferencial composicional con ANCOM-BC y MaAsLin2.
  • Un caso real de extremo a extremo.

A partir de aquí, las direcciones naturales: shotgun (Kraken2, MetaPhlAn4, HUMAnN3), análisis funcional inferido (PICRUSt2, Tax4Fun2), estudios longitudinales, microbioma + multi-ómica.

Cuando publiquemos el libro Metagenómica con DADA2 y phyloseq: del FASTQ al informe, profundizaremos en multi-cohorte, batch effects con RUV, modelos mixtos y un caso completo con repositorio reproducible público.

Para terminar de redondear

Si te ha quedado curiosidad sobre QIIME2, la alternativa Python/CLI al pipeline R que hemos seguido, el Anexo la cubre. No es necesario para haber completado la ruta, pero ayuda a entender por qué Rmori adopta el camino R/Bioconductor.