Asignación taxonómica con SILVA

r
bioconductor
metagenomica
De ASVs anónimos a nombres bacterianos. Cómo funciona assignTaxonomy() con el clasificador Bayesiano de RDP, comparación SILVA vs GTDB, asignación a nivel de especie con addSpecies() y los límites del 16S.

De secuencia a nombre

Hasta ahora tienes ASVs como TACGTAGGGTGCAAGCGTTAATC..., secuencias verdaderas pero anónimas. La taxonomía les da contexto biológico:

TACGTAGGGTGCAAGCGTT... → Bacteria; Bacteroidetes; Bacteroidia; Bacteroidales;
                          Bacteroidaceae; Bacteroides; Bacteroides fragilis

Sin ese paso, las gráficas y análisis siguientes son ilegibles. Con él, los resultados se conectan con literatura existente.

El clasificador detrás

DADA2 implementa el clasificador Bayesiano de RDP (Ribosomal Database Project), el estándar de facto para asignación taxonómica de 16S. La idea:

  1. Toma k-mers de tu ASV (por defecto kmer = 8).
  2. Compara con todos los k-mers de cada taxón en la base de referencia.
  3. Calcula la probabilidad de que el ASV pertenezca a cada taxón.
  4. Por bootstrap (default 100 iteraciones), estima la confianza de cada asignación a cada nivel.

Solo se reporta una asignación si el bootstrap es ≥ 50 % (default). Lo veremos abajo.

Bases de referencia: SILVA vs GTDB

Dos opciones modernas:

SILVA 138: el estándar histórico

  • Cobertura: la más amplia para 16S/18S. Mejor para muestras “todo”.
  • Curación: alineamiento multi-criterio cuidadoso.
  • Acceso: las bases formateadas para DADA2 están en Zenodo.

Es lo que cubrimos en este tutorial. Default razonable para cualquier muestra 16S.

GTDB r214: la alternativa moderna

  • Cobertura: foco en bacterias y arqueas con genomas completos disponibles.
  • Taxonomía: re-estructurada con criterios filogenómicos (más precisa biológicamente, pero a veces “rompe” géneros familiares).
  • Acceso: data.gtdb.ecogenomic.org/releases/ o paquete DADA2GTDBR (R-universe).

Cuándo usarla: shotgun, MAGs, análisis comparativo con genomas. Para amplicón clínico estándar, SILVA es lo más reportado.

Asignación con SILVA

Asumiendo que descargaste silva_nr99_v138.1_wSpecies_train_set.fa.gz en ref/ en el tutorial 4:

library(dada2)

# Cargar tabla de secuencias sin quimeras (del tutorial 7)
seqtab_nochim <- readRDS("output/seqtab_nochim.rds")

# Asignación hasta género
taxa <- assignTaxonomy(
    seqtab_nochim,
    refFasta    = "ref/silva_nr99_v138.1_wSpecies_train_set.fa.gz",
    multithread = TRUE,
    minBoot     = 50
)

head(taxa)
#                       Kingdom    Phylum         Class             ...  Genus
# TACGTAGGGTGCAAGCG...  "Bacteria" "Bacteroidota" "Bacteroidia"     ...  "Bacteroides"
# TACGTAGGTGGCAAGCG...  "Bacteria" "Firmicutes"   "Clostridia"      ...  "Faecalibacterium"
# ...

assignTaxonomy() tarda algunos minutos a varias horas según el dataset:

  • 2 000 ASVs en 8 cores: ~ 10 min.
  • 10 000 ASVs: ~ 30-60 min.

Salida: matriz con 7 columnas (Kingdom → Species en el caso de la base con especies, hasta Genus si no).

Extender a especie con addSpecies()

assignTaxonomy() solo es seguro hasta género. Para especies, DADA2 usa exact matching contra otra base SILVA específica:

taxa_species <- addSpecies(
    taxa,
    refFasta = "ref/silva_species_assignment_v138.1.fa.gz"
)

head(taxa_species[, c("Genus", "Species")])
#                       Genus              Species
# TACGTAGGGTGCAAGCG...  "Bacteroides"      "fragilis"
# TACGTAGGTGGCAAGCG...  "Faecalibacterium" NA
# ...

Para un ASV se asigna especie solo si su secuencia coincide exactamente con la de SILVA. Muchos ASVs quedarán como NA en Species, eso es honesto: el 16S no resuelve todas las especies.

Porcentaje típico de ASVs con especie asignada: 30-60 % en heces humanas, < 20 % en suelos. Es lo que es.

Inspeccionar los resultados

Conviene mirar varios checks rápidos:

1. ¿Cuántos ASVs por nivel taxonómico?

library(tidyverse)

taxa_df <- as.data.frame(taxa_species) |>
    rownames_to_column("asv_seq")

# Cuántos ASVs identificados a cada nivel
taxa_df |>
    summarise(
        with_kingdom = sum(!is.na(Kingdom)),
        with_phylum  = sum(!is.na(Phylum)),
        with_class   = sum(!is.na(Class)),
        with_order   = sum(!is.na(Order)),
        with_family  = sum(!is.na(Family)),
        with_genus   = sum(!is.na(Genus)),
        with_species = sum(!is.na(Species)),
        total        = n()
    )

Esperable en heces humanas:

Nivel % con asignación
Kingdom > 99 %
Phylum > 95 %
Class > 90 %
Order > 85 %
Family > 80 %
Genus 65-85 %
Species 30-60 %

Si Kingdom es bajo (< 95 %), tu setup tiene problemas: ASVs muy ruidosos, primers mal recortados o base de referencia equivocada.

2. Detección de no-bacterianos

Si assignTaxonomy() clasifica algo como Eukaryota o Archaea cuando esperabas bacterias, hay contaminación. Para 16S V3-V4 con primers bacterianos:

table(taxa_df$Kingdom, useNA = "ifany")
# Bacteria  Archaea  Eukaryota  <NA>
#     2089       12          4    9

12 arqueas es plausible (en heces, microbiota arqueal real). 4 eucariotas suelen ser contaminación de DNA del host (humano) o de comida (planta), filtrarlas antes de análisis:

# Quitar no-bacterianos y NAs en Kingdom
keep <- taxa_df$Kingdom == "Bacteria" & !is.na(taxa_df$Kingdom)
taxa_filtered <- taxa[keep, ]
seqtab_filtered <- seqtab_nochim[, keep]

3. ¿Hay un Phylum mayoritario sospechoso?

taxa_df |>
    count(Phylum, sort = TRUE) |>
    head(10)

En heces humanas esperas Bacillota (antes Firmicutes), Bacteroidota, Actinobacteriota, Proteobacteria. Si dominan filos exóticos, algo va mal: muestras mal procesadas o contaminación.

Bootstrap y minBoot

minBoot = 50 (default): solo reporta asignación si el bootstrap pasa el 50 %.

Aumentar a minBoot = 80 da asignaciones más conservadoras:

taxa_strict <- assignTaxonomy(
    seqtab_nochim,
    refFasta    = "ref/silva_nr99_v138.1_wSpecies_train_set.fa.gz",
    multithread = TRUE,
    minBoot     = 80
)

Resultado: más NA en niveles inferiores. Útil para reports que necesitan certeza alta. Para EDA general, 50 está bien.

Asignación con DECIPHER (alternativa)

DECIPHER::IdTaxa() es una alternativa más moderna a assignTaxonomy(). Más rápida con datasets grandes, soporta GTDB bien. Útil si tu pipeline lo necesita:

library(DECIPHER)

# Cargar base SILVA pre-formateada para DECIPHER
load("ref/SILVA_SSU_r138_2019.RData")  # crea objeto `trainingSet`

# Convertir ASVs a DNAStringSet
asv_seqs <- DNAStringSet(getSequences(seqtab_nochim))

# Clasificar
ids <- IdTaxa(asv_seqs, trainingSet, strand = "top", processors = 8)

DECIPHER da algoritmo distinto pero conceptualmente equivalente. La elección entre assignTaxonomy() y IdTaxa() raramente importa para 16S amplicón típico.

El patrón idiomático

library(dada2)

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

taxa <- assignTaxonomy(
    seqtab_nochim,
    refFasta    = "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"
taxa <- taxa[keep, ]
seqtab_nochim <- seqtab_nochim[, keep]

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

Con esto cerramos el pipeline DADA2 puro. Tienes seqtab_final (ASVs × muestras) y taxa (ASV × niveles). Es todo lo que necesitas para construir un phyloseq.

Limitaciones del 16S: recordatorio honesto

Antes de pasar a análisis, conviene recordar:

  • Cepa no se resuelve. Dos cepas de E. coli clínicamente muy distintas pueden tener V3-V4 idéntico.
  • Especie no siempre. Especialmente géneros como Lactobacillus, Streptococcus: el 16S no distingue rutinariamente.
  • Asignación incompleta es lo normal. Que un 40 % de tus ASVs no tengan especie asignada no es bug, es el límite de la técnica.
  • Una asignación “errónea” puede ser exacta. Si SILVA actualizó la taxonomía de un género, tu asignación con SILVA 138.1 puede no coincidir con un paper de 2018. No es error: es que la taxonomía es viva.

Trampas habituales

  • Bases de SILVA descargadas de silva.de en lugar del Zenodo de DADA2. La estructura es distinta. Solo la del Zenodo funciona con assignTaxonomy(). Error común al empezar.
  • Olvidar filtrar no-bacterianos. Si secuenciaste con primers bacterianos, las arqueas y eucariotas que aparezcan son ruido (excepto en estudios específicos de arqueas). Filtra antes de phyloseq.
  • Reportar “Species: NA” como “no detectada”. NA significa “no asignable con confianza”. El organismo SÍ está. No confundir con ausencia biológica.
  • Recalcular assignTaxonomy() en cada render. Es lento. saveRDS + freeze: auto.
  • Cambiar de SILVA a GTDB sin documentar. Los resultados ASV-a-taxonomía cambian. Si publicas, especifica versión exacta de la base usada (SILVA 138.1, GTDB r214…).

Cierre del Bloque 2

Has completado el pipeline DADA2 puro:

FASTQ → filterAndTrim() → learnErrors() → dada()
     → mergePairs() → makeSequenceTable() → removeBimeraDenovo()
     → assignTaxonomy() → addSpecies()

Tres objetos serializados: seqtab_final.rds (ASVs × muestras), taxa.rds (taxonomía), y track.rds (auditoría del pipeline).

En la siguiente entrega

Tienes todo el material crudo. Lo siguiente es empaquetarlo en phyloseq: el contenedor del microbioma que combina tabla ASV + metadata + taxonomía + árbol filogenético. Una vez en phyloseq, los análisis se simplifican mucho. Lo siguiente.