phyloseq: contenedor y composición taxonómica

r
bioconductor
metagenomica
El contenedor central del microbioma. Cómo se construye un phyloseq con otu_table + sample_data + tax_table + phy_tree desde el output de DADA2, cómo se manipula sin desincronizar y los plots de composición taxonómica que abren cualquier informe.

La idea

Tienes tres piezas que tienen que ir juntas:

  • seqtab_final: tabla ASV × muestra (conteos).
  • taxa: tabla ASV × niveles taxonómicos (Kingdom → Species).
  • metadata: tabla muestra × covariables (condición, batch, edad…).

Y opcionalmente una cuarta:

  • tree: árbol filogenético sobre las ASVs (necesario para UniFrac).

Si trabajas con esas cuatro tablas sueltas, te vas a desincronizar: filtrar muestras y olvidar filtrar metadata, reordenar y desalinear, perder un ASV al filtrar quimeras tarde y olvidar actualizarlo en la taxonomía.

phyloseq resuelve esto: un único objeto que mantiene las cuatro sincronizadas. Filtras, transformas o subseteas, y todo se mueve a la vez. Es el mismo principio que SummarizedExperiment en RNA-seq, adaptado a microbioma.

Anatomía del objeto

phyloseq
├── otu_table   (matriz: filas/cols = ASVs vs muestras)
├── tax_table   (matriz: filas = ASVs, cols = Kingdom...Species)
├── sample_data (DataFrame: filas = muestras, cols = covariables)
└── phy_tree    (objeto ape::phylo, opcional)

Los nombres de fila/columna deben coincidir: las ASVs del otu_table deben estar en tax_table y phy_tree. Las muestras del otu_table deben estar en sample_data. phyloseq lo valida al construir.

Construir el objeto

A partir de las salidas de DADA2 (tutoriales 5-8):

library(phyloseq)
library(tidyverse)

# Cargar piezas
seqtab_final <- readRDS("output/seqtab_final.rds")  # ASVs × muestras
taxa         <- readRDS("output/taxa.rds")          # ASV × taxonomía
metadata     <- read.csv("data/metadata.csv", row.names = 1)

# Verificar alineación muestras
stopifnot(all(rownames(metadata) %in% rownames(seqtab_final)))

# Construir phyloseq
ps <- phyloseq(
    otu_table(seqtab_final, taxa_are_rows = FALSE),
    sample_data(metadata),
    tax_table(taxa)
)

ps
# phyloseq-class experiment-level object
# otu_table()   OTU Table:         [ 2114 taxa and 66 samples ]
# sample_data() Sample Data:       [ 66 samples by 8 sample variables ]
# tax_table()   Taxonomy Table:    [ 2114 taxa by 7 taxonomic ranks ]

Detalles:

  • taxa_are_rows = FALSE porque el output de DADA2 tiene muestras en filas, ASVs en columnas. Si tu matriz está transpuesta, pon TRUE.
  • stopifnot() antes de construir es la línea crítica, si los IDs de muestra no coinciden, phyloseq lo rechaza con un error claro.

Renombrar ASVs a IDs cortos

Las ASVs vienen identificadas por su secuencia completa (TACGTAGGGTGCAAGCG...). Es preciso pero ilegible. Conviene renombrarlas a IDs cortos (ASV0001, ASV0002, …) y guardar las secuencias en un slot separado:

library(Biostrings)

# Guardar las secuencias originales
dna <- DNAStringSet(taxa_names(ps))
names(dna) <- paste0("ASV", sprintf("%04d", seq_along(dna)))

# Renombrar el phyloseq
taxa_names(ps) <- names(dna)

# Adjuntar las secuencias como un slot extra
ps@refseq <- dna

ps
# otu_table()   OTU Table:         [ 2114 taxa and 66 samples ]
# sample_data() Sample Data:       [ 66 samples by 8 sample variables ]
# tax_table()   Taxonomy Table:    [ 2114 taxa by 7 taxonomic ranks ]
# refseq()      DNAStringSet:      [ 2114 reference sequences ]

Ahora las filas de los plots dicen ASV0042 en vez de una secuencia de 450 bp. Las secuencias siguen accesibles via refseq(ps) cuando hace falta (BLAST, árbol, etc.).

Añadir el árbol filogenético

Para UniFrac (lo veremos en el tutorial 11) hace falta un árbol filogenético sobre las ASVs. Se construye desde sus secuencias:

library(phangorn)
library(DECIPHER)

# 1. Alinear las secuencias (DECIPHER es mejor que muscle/mafft para 16S)
alignment <- AlignSeqs(refseq(ps), anchor = NA, verbose = FALSE)

# 2. Construir árbol neighbour-joining inicial
phang_align <- phyDat(as(alignment, "matrix"), type = "DNA")
dm <- dist.ml(phang_align)
treeNJ <- NJ(dm)

# 3. Refinar con maximum likelihood (opcional, mejor pero más lento)
fit <- pml(treeNJ, data = phang_align)
fitGTR <- update(fit, k = 4, inv = 0.2)
fitGTR <- optim.pml(fitGTR, model = "GTR", optInv = TRUE, optGamma = TRUE,
                    rearrangement = "stochastic",
                    control = pml.control(trace = 0))

# Pegar al phyloseq
phy_tree(ps) <- fitGTR$tree

El proceso completo tarda minutos a horas según el número de ASVs. Para datasets exploratorios el NJ puro vale. Para análisis publicable, ML refinado.

Acceso a los componentes

phyloseq tiene un montón de getters/setters:

ntaxa(ps)            # 2114
nsamples(ps)         # 66

sample_names(ps)     # nombres de muestras
taxa_names(ps)       # nombres de ASVs (ASV0001, ASV0002, ...)
rank_names(ps)       # niveles taxonómicos disponibles

otu_table(ps)[1:5, 1:3]
sample_data(ps)
tax_table(ps)[1:5, ]

Para trabajar con tibbles (más cómodo):

library(tibble)

# Tabla ASV plana
otu_df <- as.data.frame(otu_table(ps)) |>
    rownames_to_column("sample")

# Taxonomía como tibble
tax_df <- as.data.frame(tax_table(ps)) |>
    rownames_to_column("asv")

# Metadata
meta_df <- as.data.frame(sample_data(ps)) |>
    rownames_to_column("sample")

Filtros básicos

Los más comunes:

Filtrar muestras

# Sólo muestras del grupo "case"
ps_case <- subset_samples(ps, condition == "case")

# Muestras con > 5000 reads totales
keep <- sample_sums(ps) > 5000
ps <- prune_samples(keep, ps)

Filtrar ASVs

# Solo ASVs presentes en al menos 3 muestras con > 5 reads
filter_taxa(ps, function(x) sum(x > 5) >= 3, prune = TRUE)

# Solo ASVs con > 0.01 % de abundancia relativa global
ps_filt <- filter_taxa(ps, function(x) mean(x) > 1e-4, prune = TRUE)

# Por taxonomía
ps_firmicutes <- subset_taxa(ps, Phylum == "Bacillota")  # antes "Firmicutes"

Lo importante: cualquier subset_* o prune_* mantiene la sincronía entre las cuatro tablas. No vas a quedarte con un tax_table que tenga ASVs que ya no existen en otu_table.

Aglomerar a un nivel taxonómico

Para análisis a nivel de phylum, family, genus…

ps_genus  <- tax_glom(ps, taxrank = "Genus",  NArm = TRUE)
ps_phylum <- tax_glom(ps, taxrank = "Phylum", NArm = TRUE)

ntaxa(ps_genus)
# [1] 287

NArm = TRUE descarta ASVs sin asignación a ese nivel (los NA que vimos en el tutorial 8). Útil cuando quieres barras “limpias” de phylum.

Transformaciones

Los conteos crudos no son lo que quieres pintar. Transforma siempre:

# Abundancia relativa (suma 1 por muestra)
ps_rel <- transform_sample_counts(ps, function(x) x / sum(x))

# Log10 con pseudocount
ps_log <- transform_sample_counts(ps, function(x) log10(x + 1))

# CLR (la opción composicional correcta, lo veremos en el tutorial 12)
library(microbiome)
ps_clr <- microbiome::transform(ps, "clr")

Para los plots de composición de este tutorial, usar transform_sample_counts(ps, function(x) x / sum(x)), abundancia relativa. Es lo legible.

Composición taxonómica: bar plots

El plot canónico del microbioma:

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") +
    labs(x = NULL, y = "Abundancia relativa") +
    theme_minimal(base_size = 11) +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))

Con fill = "Phylum" cada barra (muestra) se rellena con la abundancia de cada filo. En microbioma humano esperas dominancia de Bacillota (antes Firmicutes) + Bacteroidota.

Para agrupar visualmente por condición:

plot_bar(ps_phylum, x = "sample", fill = "Phylum") +
    facet_wrap(~ condition, scales = "free_x") +
    labs(x = NULL, y = "Abundancia relativa") +
    theme_minimal()

Composición a nivel de género (los top N)

A nivel de género hay cientos de taxones. Mostrarlos todos da un arcoíris ilegible. La estrategia: top 15 + Other.

library(dplyr)

ps_genus <- tax_glom(ps_rel, taxrank = "Genus", NArm = TRUE)

# Identifica los top 15 géneros por abundancia media
top15 <- names(sort(taxa_sums(ps_genus), decreasing = TRUE))[1:15]

ps_top <- prune_taxa(top15, ps_genus)

# El resto se agrupa como "Other" sumando lo descartado
other_abund <- 1 - sample_sums(ps_top)

# Plot
plot_bar(ps_top, x = "sample", fill = "Genus") +
    labs(x = NULL, y = "Abundancia relativa (top 15)") +
    scale_fill_brewer(palette = "Set3") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))

Esta es la estética de la mayoría de papers de microbioma. Set3 o paletas custom de ColorBrewer funcionan bien.

Heatmap de los top taxa

Otra vista útil: heatmap de los top N taxones × muestras.

library(pheatmap)

# Matriz de abundancia relativa de los top 30 géneros
top30 <- names(sort(taxa_sums(ps_genus), decreasing = TRUE))[1:30]
ps_h <- prune_taxa(top30, ps_genus)

# Construir matriz: filas = géneros, columnas = muestras
mat <- as.matrix(otu_table(ps_h))
if (!taxa_are_rows(ps_h)) mat <- t(mat)

# Renombrar filas a nombres de género
rownames(mat) <- tax_table(ps_h)[, "Genus"]

# Anotación por condición
annot <- as.data.frame(sample_data(ps_h))[, "condition", drop = FALSE]

pheatmap(log10(mat + 1e-5),
         annotation_col = annot,
         show_colnames = FALSE,
         color = colorRampPalette(c("#F8F5F0", "#5F8575", "#3F6354"))(50))

Es la vista canónica para ver si una condición tiene un perfil de microbioma claramente distinto.

microViz: plots más modernos

microViz (instalado en el tutorial 4) ofrece una API más limpia sobre phyloseq. Una alternativa concisa:

library(microViz)

ps |>
    tax_fix() |>                          # rellena NAs con el nivel padre
    tax_transform("compositional") |>     # abundancia relativa
    comp_barplot(
        tax_level = "Genus",
        n_taxa    = 12,
        bar_outline_colour = NA,
        palette = distinct_palette(n = 12)
    ) +
    facet_wrap(~ condition, scales = "free_x") +
    coord_flip() +
    labs(x = NULL)

Vale la pena explorar microViz cuando llegues a la fase de informe, su API es mucho más concisa que el plot_bar original.

El patrón idiomático

El flujo típico desde DADA2 a un primer plot:

library(phyloseq); library(tidyverse); library(Biostrings)

# Cargar piezas
seqtab_final <- readRDS("output/seqtab_final.rds")
taxa <- readRDS("output/taxa.rds")
metadata <- read.csv("data/metadata.csv", row.names = 1)

# Construir
ps <- phyloseq(
    otu_table(seqtab_final, taxa_are_rows = FALSE),
    sample_data(metadata),
    tax_table(taxa)
)

# Renombrar ASVs y guardar secuencias
dna <- DNAStringSet(taxa_names(ps))
names(dna) <- paste0("ASV", sprintf("%04d", seq_along(dna)))
taxa_names(ps) <- names(dna)
ps@refseq <- dna

# Filtro mínimo: ASVs presentes en al menos 3 muestras
ps <- filter_taxa(ps, function(x) sum(x > 5) >= 3, prune = TRUE)

# Guardar
saveRDS(ps, "output/phyloseq.rds")

# Primer plot
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, fill = "Phylum") +
    facet_wrap(~ condition, scales = "free_x") +
    theme_minimal()

Diez líneas para tener el objeto principal del análisis listo y un primer plot.

Trampas habituales

  • Construir phyloseq con metadata cuyos IDs no coinciden con colnames(seqtab_final). phyloseq tira un error, pero el mensaje confunde. El stopifnot() previo lo cazará antes.
  • Olvidar taxa_are_rows = FALSE. Es la trampa #1 al venir de RNA-seq donde genes (filas) son la unidad. En DADA2, las muestras son filas y las ASVs son columnas.
  • Plotear conteos crudos en lugar de abundancia relativa. Una muestra con 100k reads aparecerá 10× más alta que una con 10k, sin que la biología cambie. Transforma siempre antes de plotear.
  • Mostrar 200 géneros en una sola barra apilada. Ilegible. Top 12-15 + “Other” es lo idiomático.
  • tax_glom con NArm = FALSE que arrastra una categoría enorme “NA” en cada nivel. Casi siempre NArm = TRUE.

En la siguiente entrega

Tienes el contenedor y la composición. La siguiente pieza es diversidad alfa, Shannon, Simpson, Chao1: cuándo aplicar cada índice, qué dice de la biología y por qué “más diverso ≠ mejor microbioma”. Lo siguiente.