Dada2 - Remoción de quimeras y seguimiento de lecturas

Las quimeras son secuencias artificiales producidas durante la PCR: una cadena parcialmente extendida se desprende de su molde original y se hibrida con otra molécula diferente, generando una secuencia híbrida que no existe en la naturaleza. Si no se eliminan, inflan artificialmente el número de ASVs y la riqueza de especies observada.

DADA2 elimina quimeras después del denoising y el merge — no antes. Esto es intencional: el denoising corrige errores de secuenciación y el merge ensamble pares, y solo una vez que tenemos secuencias corregidas y completas tiene sentido buscar híbridos entre ellas.

Paso 1 — Remover quimeras con `removeBimeraDenovo()`

Qué hace: compara cada ASV contra todas las demás secuencias de mayor abundancia en el dataset. Si una ASV puede explicarse como el híbrido de dos secuencias parentales más abundantes, se considera quimérica y se elimina.

seqtab.nochim <- removeBimeraDenovo(
  seqtab,
  method      = "consensus",
  multithread = 2,
  verbose     = TRUE
)

Parámetros clave:

method = "consensus" — cada ASV candidata se evalúa de forma independiente en cada muestra; se elimina solo si es quimérica en la mayoría de muestras donde aparece. Es el método recomendado para la mayoría de los estudios.
method = "pooled" — evalúa todas las muestras en conjunto; más sensible para detectar quimeras raras, pero más lento y puede ser más agresivo.

# Proporción de lecturas que sobreviven
freq.nochim <- sum(seqtab.nochim) / sum(seqtab)
message("Lecturas no quiméricas: ", round(freq.nochim * 100, 2), "%")

# Proporción de ASVs únicos que sobreviven
message("ASVs totales antes: ", ncol(seqtab))
message("ASVs despues de quimeras: ", ncol(seqtab.nochim))

Qué chequear: es normal que se elimine una fracción alta de ASVs únicos (secuencias distintas) pero que la proporción de lecturas retenidas sea alta — las quimeras suelen ser raras (pocas lecturas por ASV).

% lecturas retenidas	Interpretación
> 95%	Normal
80–95%	Aceptable — revisar parámetros de PCR
< 80%	Preocupante — revisar calidad del ADN o diseño de primers
< 50%	Grave — revisar el diseño experimental completo

Paso 2 — Distribución de longitudes y filtrado

Qué hace: después de remover quimeras, verifica que todas las ASVs tengan la longitud esperada para el amplicón. Secuencias muy cortas o muy largas pueden ser artefactos de ensamblaje, no quimeras verdaderas pero sí ruido.

Primero visualiza la distribución antes de filtrar:

seq_lengths <- nchar(colnames(seqtab.nochim))
median_len  <- median(seq_lengths)
sd_len      <- sd(seq_lengths)

message("Longitud mediana: ", median_len, " bp")
message("Desviación estándar: ", round(sd_len, 2), " bp")

# Distribución rápida
table(seq_lengths)

El pipeline guarda una gráfica de la distribución antes del filtro:

png(file.path(output_location, "logs",
              paste0(run_name, "_", primer.data$locus_shorthand[i], "_seq_lengths_distribution.png")),
    width = 1800, height = 1200, res = 300)
print(
  ggplot2::ggplot(data = data.frame(Length = seq_lengths), ggplot2::aes(x = Length)) +
    ggplot2::geom_histogram(binwidth = 1, fill = "#69b3a2", color = "#2d5f4f", alpha = 0.8) +
    ggplot2::geom_vline(aes(xintercept = median_len),
                        color = "#e63946", linetype = "dashed", linewidth = 1) +
    ggplot2::geom_vline(aes(xintercept = median_len - sd_len),
                        color = "#f77f00", linetype = "dotted", linewidth = 0.8) +
    ggplot2::geom_vline(aes(xintercept = median_len + sd_len),
                        color = "#f77f00", linetype = "dotted", linewidth = 0.8) +
    ggplot2::labs(title = "Distribución de Longitudes de Secuencias ASV (Sin Filtrar)",
                  subtitle = paste0("n = ", length(seq_lengths), " secuencias"),
                  x = "Longitud (bp)", y = "Frecuencia") +
    ggplot2::theme_minimal()
)
dev.off()

Parámetros clave del filtro: el pipeline calcula automáticamente el rango aceptable como mediana ± 1 desviación estándar. Para 12S (~62 bp) esto es conservador y elimina solo secuencias claramente fuera del rango del amplicón.

min_len <- median_len - sd_len
max_len <- median_len + sd_len

indexes.to.keep <- which(
  nchar(colnames(seqtab.nochim)) <= max_len &
  nchar(colnames(seqtab.nochim)) >= min_len
)

cleaned.seqtab.nochim    <- seqtab.nochim[,  indexes.to.keep, drop = FALSE]
filteredout.seqtab.nochim <- seqtab.nochim[, -indexes.to.keep, drop = FALSE]

message("ASVs antes del filtro: ",   ncol(seqtab.nochim))
message("ASVs después del filtro: ", ncol(cleaned.seqtab.nochim))
message("ASVs filtradas por longitud: ", ncol(filteredout.seqtab.nochim))

Y guarda una segunda gráfica después del filtro para confirmar que el resultado es limpio:

seq_lengths_clean <- nchar(colnames(cleaned.seqtab.nochim))
median_len_clean  <- median(seq_lengths_clean)
sd_len_clean      <- sd(seq_lengths_clean)

png(file.path(output_location, "logs",
              paste0(run_name, "_", primer.data$locus_shorthand[i], "_cleaned_seq_lengths_distribution.png")),
    width = 1800, height = 1200, res = 300)
print(
  ggplot2::ggplot(data = data.frame(Length = seq_lengths_clean), ggplot2::aes(x = Length)) +
    ggplot2::geom_histogram(binwidth = 1, fill = "#06d6a0", color = "#118ab2", alpha = 0.8) +
    ggplot2::geom_vline(aes(xintercept = median_len_clean),
                        color = "#ef476f", linetype = "dashed", linewidth = 1) +
    ggplot2::labs(title = "Distribución de Longitudes de Secuencias ASV (Filtradas)",
                  subtitle = paste0("n = ", length(seq_lengths_clean), " ASVs retenidas"),
                  x = "Longitud (bp)", y = "Frecuencia") +
    ggplot2::theme_minimal()
)
dev.off()

Qué chequear: la gráfica post-filtro debe mostrar un pico estrecho centrado en la longitud esperada del amplicón. Si hay múltiples picos, puede indicar amplificación inespecífica. Las ASVs eliminadas se guardan en final_data/logs/ para revisión.

Paso 3 — Tabla de seguimiento de lecturas (read tracking)

Qué hace: resume cuántas lecturas sobrevivieron cada paso del pipeline. Es la herramienta diagnóstica más importante — una caída abrupta en cualquier paso señala exactamente dónde hay un problema.

getN <- function(x) sum(getUniques(x))

track <- cbind(
  out[exists, ],
  sapply(dadaFs,  getN),
  sapply(dadaRs,  getN),
  sapply(mergers, getN),
  rowSums(seqtab.nochim),
  rowSums(cleaned.seqtab.nochim)
)
colnames(track) <- c("input", "filtered", "denoisedF", "denoisedR",
                     "merged", "nonchim", "length_filter")
rownames(track) <- sample.names

head(track)

# Retenciónes globales por paso
message("Input -> filtered:     ", round(mean(track[,"filtered"]      / track[,"input"]    * 100), 1), "%")
message("Filtered -> merged:    ", round(mean(track[,"merged"]        / track[,"filtered"] * 100), 1), "%")
message("Merged -> nonchim:     ", round(mean(track[,"nonchim"]       / track[,"merged"]   * 100), 1), "%")
message("Nonchim -> length flt: ", round(mean(track[,"length_filter"] / track[,"nonchim"]  * 100), 1), "%")

# Guardar tabla
write.csv(track,
  file.path(output_location, "logs",
            paste0(run_name, "_", primer.data$locus_shorthand[i], "_read_tracking.csv")),
  row.names = TRUE)

Qué significa cada columna:

Columna	Paso	Pérdida esperada por
`input`	Lecturas crudas	—
`filtered`	Tras filtrado de calidad	Reads de baja calidad o muy cortos
`denoisedF/R`	Tras denoising de F y R	Corrección de errores
`merged`	Tras ensamblaje de pares	Reads sin overlap suficiente
`nonchim`	Tras remoción de quimeras	ASVs quiméricas (raras en lecturas)
`length_filter`	Tras filtro de longitud	ASVs fuera del tamaño del amplicón

Valores de referencia por transición:

Transición	Normal	Revisar si
Input → filtered	> 70%	< 50%
Filtered → merged	> 70%	< 50% (truncLen muy agresivo o overlap insuficiente)
Merged → nonchim	> 85%	< 75% (muchas quimeras — revisar PCR)
Total (input → final)	> 50%	< 30%

Paso 4 — Generar tabla ASV en formato largo

Qué hace: convierte la matriz de abundancias (muestras × ASVs) a formato largo (una fila por muestra-ASV con lecturas > 0) para facilitar el análisis downstream con phyloseq y otras herramientas.

ASV_file <- file.path(output_location, "csv_output",
                      paste0(run_name, "_", primer.data$locus_shorthand[i], "_ASV_table.csv"))

current_asv <- as.data.frame(cleaned.seqtab.nochim) %>%
  tibble::rownames_to_column("Sample_name") %>%
  tidyr::pivot_longer(cols = -Sample_name, names_to = "Sequence", values_to = "nReads") %>%
  dplyr::filter(nReads > 0) %>%
  dplyr::mutate(Label = primer.data$locus_shorthand[i]) %>%
  dplyr::relocate(Label, .after = Sample_name)

readr::write_csv(current_asv, ASV_file)
message("Tabla ASV guardada en: ", ASV_file)

El archivo resultante tiene una fila por muestra × ASV con lecturas > 0, listo para análisis downstream con phyloseq.

Para el curso

¿Qué porcentaje de lecturas se retiene después de remover quimeras? ¿Y de ASVs únicos?
¿Cuántos ASVs quedan antes y después del filtro de longitud? ¿La gráfica post-filtro muestra un pico limpio?
Revisa la tabla de tracking — ¿en qué paso se pierde más lecturas?
¿Hay muestras con retención < 30% en algún paso? ¿Cuáles y por qué podrían serlo?

Checkpoint — Antes de asignación taxonómica

Guarda todos los objetos necesarios antes de continuar con la asignación taxonómica. Si el script se interrumpe podrás retomar desde aquí sin repetir ninguno de los pasos anteriores.

save(cleaned.seqtab.nochim, seqtab.nochim, seqtab,
     dadaFs, dadaRs, mergers, track,
     sample.names, primer.data, run_name, out,
     file = file.path(output_location, "rdata_output",
                      paste0(run_name, "_", primer.data$locus_shorthand[i],
                             "_checkpoint_antes_taxonomy.RData")))
message("Checkpoint guardado — listo para asignación taxonómica")

Para reanudar desde este punto en otra sesión:

load(file.path(output_location, "rdata_output",
               paste0(run_name, "_", primer.data$locus_shorthand[i],
                      "_checkpoint_antes_taxonomy.RData")))