Dada2 - Asignación taxonómica y guardado de resultados

La asignación taxonómica convierte secuencias anónimas (ASVs) en información biológica: géneros, familias, órdenes y, cuando es posible, especies. Es el puente entre el análisis bioinformático y la ecología de comunidades.

¿Cómo funciona `assignTaxonomy()`? — El clasificador Naive Bayesian

assignTaxonomy() implementa el clasificador bayesiano ingenuo (Naive Bayesian Classifier) del RDP (Ribosomal Database Project), adaptado por el equipo de DADA2. El proceso en cuatro pasos:

Fragmentación en k-mers — cada ASV se divide en todos los fragmentos de 8 nucleótidos posibles (8-mers). Una secuencia de 62 bp genera ~55 8-mers distintos.
Comparación con la base de datos — la base de referencia está indexada en 8-mers por grupo taxonómico. El clasificador calcula la probabilidad de que los 8-mers de la ASV pertenezcan a cada grupo en la base.
Asignación jerárquica con bootstrap — la clasificación se hace nivel por nivel: Reino → Filo → Clase → Orden → Familia → Género → Especie. En cada nivel se reporta un valor de bootstrap (0–100) que indica qué tan consistente fue la asignación al remuestrear los 8-mers. El umbral por defecto es 80: si el bootstrap cae por debajo de ese valor en algún nivel, la asignación se trunca ahí y los niveles inferiores quedan como NA.
Por qué bajan las asignaciones en niveles finos — a nivel de género y especie, la base de datos necesita tener representantes muy cercanos a la secuencia de la muestra. Si el organismo no está bien representado, la asignación se queda en familia u orden. Esto no es un error — es una limitación real de la cobertura de la base de datos.

Nota

¿Por qué “naive” (ingenuo)? Porque asume que cada 8-mer es independiente de los demás, lo cual no es biológicamente cierto — las posiciones en una secuencia están correlacionadas. Pero esta simplificación hace el cálculo muy eficiente y, en la práctica, funciona bien para clasificación taxonómica.

Paso 1 — Asignar taxonomía con `assignTaxonomy()`

Qué hace: toma las secuencias de cleaned.seqtab.nochim y las compara contra la base de datos de referencia especificada en primer_data.csv. Devuelve una tabla de taxonomía y otra de valores de bootstrap por nivel.

tax_location <- file.path(metadata_location, primer.data$db_name[i])

if (!file.exists(tax_location)) {
  stop("No se encuentra la base de taxonomía: ", tax_location)
}

at <- assignTaxonomy(
  colnames(cleaned.seqtab.nochim),
  tax_location,
  tryRC            = TRUE,
  multithread      = 2,
  outputBootstraps = TRUE,
  verbose          = TRUE
)

Parámetros clave:

tryRC = TRUE — prueba también el complemento reverso de cada ASV. Importante cuando el amplicón puede amplificarse en cualquier orientación según el par de primers.
outputBootstraps = TRUE — devuelve los valores de confianza (0–100) por nivel taxonómico. Esencial para evaluar la calidad de cada asignación.
multithread = 2 — usa 2 núcleos (límite del servidor con ~30 usuarios concurrentes).

Este paso puede tardar varios minutos dependiendo del número de ASVs y el tamaño de la base de datos.

Paso 2 — Verificar el resultado

Qué hace: separa la tabla de taxonomía de los valores de bootstrap y calcula el porcentaje de ASVs asignadas por nivel taxonómico.

# Separar taxonomía y bootstraps
tax_mat  <- at$tax
boot_mat <- at$boot

tax_df <- as.data.frame(tax_mat) %>%
  tibble::rownames_to_column("Sequence")

head(tax_df)

# Porcentaje de ASVs asignados por nivel taxonómico
ranks <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")
pct_assigned <- colSums(!is.na(tax_mat)) / nrow(tax_mat) * 100
round(pct_assigned, 1)

Qué chequear: es normal que el porcentaje de asignación baje progresivamente de Reino a Especie. Una caída muy brusca en niveles altos (Filo, Clase) puede indicar un problema con la base de datos o contaminación. Ver tabla de referencia al final de la página.

El pipeline guarda automáticamente el porcentaje de ASVs sin asignación por nivel en final_data/logs/:

pct_unassigned <- vapply(ranks, function(rk) {
  col <- tax_df[[rk]]
  mean(is.na(col) | col == "")
}, numeric(1))

readr::write_csv(
  tibble::tibble(level = ranks, pct_unassigned = pct_unassigned),
  file.path(output_location, "logs",
            paste0(run_name, "_", primer.data$locus_shorthand[i], "_pct_unassigned.csv"))
)

Paso 3 — Guardar archivos de taxonomía

Qué hace: genera dos archivos CSV en final_data/csv_output/:

_taxonomy_output.csv — tabla de taxonomía (Sequence + niveles taxonómicos)
_tax_bootstrap.csv — tabla de bootstraps unida a la taxonomía, para revisar la confianza de cada asignación

taxonomy_file  <- file.path(output_location, "csv_output",
                            paste0(run_name, "_", primer.data$locus_shorthand[i], "_taxonomy_output.csv"))
bootstrap_file <- file.path(output_location, "csv_output",
                            paste0(run_name, "_", primer.data$locus_shorthand[i], "_tax_bootstrap.csv"))

# Tabla de taxonomía
joined_old_new_taxa <- tax_df %>%
  dplyr::select(dplyr::any_of(c("Sequence", ranks)))
readr::write_csv(joined_old_new_taxa, taxonomy_file)

# Tabla de bootstraps unida a taxonomía
boot_df <- as.data.frame(boot_mat) %>%
  tibble::rownames_to_column("Sequence")
names(boot_df)[-1] <- paste0("boot.", names(boot_df)[-1])
ready_new_taxa <- dplyr::left_join(boot_df, joined_old_new_taxa, by = "Sequence")
readr::write_csv(ready_new_taxa, bootstrap_file)

message("Taxonomía guardada en: ", taxonomy_file)
message("Bootstraps guardados en: ", bootstrap_file)

Qué chequear en _tax_bootstrap.csv: busca asignaciones con bootstrap < 80 en niveles altos (Family, Order) — esas son las más inciertas. Un bootstrap alto en Genus pero baja cobertura de esa especie en la base de datos puede resultar en asignaciones incorrectas de especie.

Paso 4 — Guardar insumos para phyloseq (.RData)

Qué hace: empaqueta en un único archivo .RData todos los objetos necesarios para el análisis ecológico: la matriz de abundancias, la taxonomía, los metadatos y los parámetros del run. Esto permite continuar el análisis en otra sesión sin repetir ningún paso del pipeline.

# Preparar matrices para phyloseq
asv_mat <- cleaned.seqtab.nochim
storage.mode(asv_mat) <- "integer"

tax_mat <- joined_old_new_taxa %>%
  tibble::column_to_rownames("Sequence") %>%
  as.matrix()

# Cargar metadatos — detecta automáticamente columna sample_id o Sample_name
meta_path <- file.path(metadata_location, "metadata.csv")
if (!file.exists(meta_path)) stop("No se encontró metadata.csv en: ", meta_path)

sample_meta <- readr::read_csv(meta_path, show_col_types = FALSE)

id_col <- dplyr::case_when(
  "sample_id"   %in% names(sample_meta) ~ "sample_id",
  "Sample_name" %in% names(sample_meta) ~ "Sample_name",
  TRUE ~ NA_character_
)
if (is.na(id_col)) stop("metadata.csv debe tener columna 'sample_id' o 'Sample_name'.")

sample_meta <- sample_meta %>%
  dplyr::mutate(!!id_col := trimws(as.character(.data[[id_col]]))) %>%
  dplyr::rename(SampleID = !!id_col) %>%
  tibble::column_to_rownames("SampleID")

# Alinear con asv_mat y advertir si hay muestras sin metadata
faltan    <- setdiff(rownames(asv_mat), rownames(sample_meta))
sobrantes <- setdiff(rownames(sample_meta), rownames(asv_mat))
if (length(faltan))    warning("Faltan en metadata: ", paste(faltan, collapse = ", "))
if (length(sobrantes)) message("En metadata pero no en datos: ", paste(sobrantes, collapse = ", "))
sample_meta <- sample_meta[rownames(asv_mat), , drop = FALSE]

# Parámetros del run (para trazabilidad)
params <- list(
  run_name   = run_name,
  locus      = primer.data$locus_shorthand[i],
  truncLen_F = where_trim_all_Fs,
  truncLen_R = where_trim_all_Rs,
  min_len    = min_len,
  max_len    = max_len,
  date       = as.character(Sys.time())
)

# Guardar todo en un .RData
rdata_path <- file.path(output_location, "rdata_output",
                        paste0(run_name, "_", primer.data$locus_shorthand[i], "_phyloseq_inputs.RData"))

save(asv_mat, tax_mat, sample_meta, track,
     seqtab_raw = seqtab.nochim, seqtab_clean = cleaned.seqtab.nochim,
     taxonomy_df = joined_old_new_taxa, bootstrap_df = ready_new_taxa,
     params,
     file = rdata_path, compress = "xz")

message("RData guardado en: ", rdata_path)

Qué contiene el .RData:

Objeto	Contenido
`asv_mat`	Matriz de abundancias (muestras × ASVs)
`tax_mat`	Taxonomía asignada por DADA2 (sin curación geográfica)
`sample_meta`	Metadatos de muestras alineados con `asv_mat`
`track`	Tabla de seguimiento de lecturas por paso
`params`	Parámetros del run (truncLen, longitudes, fecha)
`taxonomy_df`, `bootstrap_df`	Tablas de taxonomía y bootstraps completas

Exportar secuencias ASV a FASTA

Qué hace: guarda las secuencias de todas las ASVs en formato FASTA con identificadores legibles (ASV_0001, ASV_0002, …). Este archivo se usa en pasos posteriores como el agrupamiento en OTUs con VSEARCH.

# Crear directorio de salida
dir.create(file.path(output_location, "otus_vsearch"), recursive = TRUE, showWarnings = FALSE)

# Exportar secuencias ASV a FASTA
asv_seqs <- colnames(cleaned.seqtab.nochim)
asv_headers <- paste0("ASV_", sprintf("%04d", seq_along(asv_seqs)))
asv_fasta <- c(rbind(paste0(">", asv_headers), asv_seqs))

fasta_path <- file.path(output_location, "otus_vsearch", "fish_ASV.fasta")
writeLines(asv_fasta, fasta_path)

message("FASTA exportado: ", fasta_path)
message("Total ASVs: ", length(asv_seqs))

Causas de baja asignación taxonómica

Causa	Síntoma	Solución
Base de datos desactualizada	Muchos `NA` en Filo	Usar versión más reciente
Marcador poco representado	% bajos en todos los niveles	Considerar base de datos más específica
`minBoot` muy alto (> 80)	Muchos `NA` en niveles medios	Bajar a 60 con precaución
Quimeras no removidas	Asignaciones inconsistentes entre muestras	Verificar paso anterior