Merge pairs + remoción de quimeras

r
bioconductor
metagenomica
Cómo combinar R1 y R2 en un amplicón completo con mergePairs(), construir la tabla de secuencias y filtrar quimeras de PCR con removeBimeraDenovo. Decisiones que afectan a la sensibilidad y la precisión finales.

Lo que queda por hacer

Tras dada() tienes ASVs de R1 y R2 por separado. Eso no es lo que quieres. La biología está en el amplicón completo. Toca:

  1. Mergear R1 + R2 en el amplicón completo.
  2. Construir la matriz ASV × muestra final.
  3. Filtrar por longitud esperada (descarta fragmentos sospechosos).
  4. Quitar quimeras (artefactos de PCR).

Estos cuatro pasos son rápidos (segundos a minutos) comparados con dada(), pero las decisiones tienen impacto.

Merge pairs

mergePairs() toma los ASVs de R1, los de R2 y los reads filtrados originales. Para cada pareja R1-R2, busca el overlap y los combina.

mergers <- mergePairs(
    dadaFs, filtFs,
    dadaRs, filtRs,
    verbose = TRUE
)

Output por consola:

36527 paired-reads (in 256 unique pairings) successfully merged out of 39014
   (in 1245 pairings) input.

Lectura: ~ 93 % de los reads pareados pudieron mergear correctamente. Cualquier cosa < 80 % es señal de problema:

  • Overlap insuficiente: trunaste demasiado en el tutorial 5. Revisa truncLen.
  • Reads de baja calidad en el overlap: los R2 al final tienen errores que impiden detectar el solapamiento exacto.
  • Amplicón de longitud variable: si tu amplicón tiene mucha variación de tamaño (poco común en 16S, normal en ITS), puede que tu setup truncLen no acomode los más largos.

Parámetros útiles de mergePairs

  • minOverlap = 12 (default): mínimo de overlap requerido. Aumentar a 20 si quieres ser estricto.
  • maxMismatch = 0 (default): cero mismatches permitidos en el overlap. Subir a 1 si tienes mucha caída de calidad en R2, recupera reads a costa de algo de falsos positivos.

Si tu tasa de merge es baja, prueba:

mergers <- mergePairs(
    dadaFs, filtFs, dadaRs, filtRs,
    minOverlap   = 12,
    maxMismatch  = 1,
    verbose      = TRUE
)

Construir la tabla de secuencias

seqtab <- makeSequenceTable(mergers)
dim(seqtab)
# [1]   66  2543

# 66 muestras × 2543 ASVs (antes de quimeras)

La estructura: filas son muestras, columnas son secuencias ASV. Los valores son conteos de cuántas reads de esa muestra son ese ASV.

seqtab[1:5, 1:3]
#         TACGTAGGGT...  TACGTAGGGA...  TACGGGGGGG...
# SampleA          254             0             88
# SampleB          312            42             71
# ...

Las columnas se llaman por la secuencia exacta. Es eficiente pero pesado de mirar. En el siguiente tutorial las reemplazaremos por IDs cortos al construir phyloseq.

Filtrar por longitud esperada

Inspecciona la distribución de longitudes:

table(nchar(getSequences(seqtab)))
# 396  398  400  402  ... 460  462
#   2    5   17   34  ... 287   12

Para V3-V4, el rango esperado es ~ 430-460 bp. Los ASVs muy cortos (< 380 bp) o muy largos (> 480 bp) suelen ser artefactos:

  • Fragmentos de PCR no específicos.
  • Quimeras parciales que aún no eliminaste.
  • Contaminación de otras regiones amplificadas inespecíficamente.

Filtrado típico:

seqtab2 <- seqtab[, nchar(colnames(seqtab)) %in% 400:460]
dim(seqtab2)
# [1]   66  2401

Para V4 (250 bp esperados), filtrar a 240-260. Para ITS (longitud variable), no filtres por longitud, la variabilidad es biológica.

Remoción de quimeras

Las quimeras de PCR son artefactos donde dos secuencias distintas se recombinaron durante la amplificación. Aparecen como ASVs nuevos que en realidad no existen. Pueden ser el 5-15 % de los ASVs.

DADA2 las detecta con removeBimeraDenovo:

seqtab_nochim <- removeBimeraDenovo(
    seqtab2,
    method      = "consensus",
    multithread = TRUE,
    verbose     = TRUE
)
# Identified 287 bimeras out of 2401 input sequences.

dim(seqtab_nochim)
# [1]   66  2114

Métodos disponibles:

  • "consensus" (recomendado): un ASV se marca como quimera si lo es en al menos la mitad de las muestras donde aparece. Robusto.
  • "pooled": trata todas las muestras como una. Más sensible pero también más agresivo.
  • "per-sample": detecta por muestra. Menos consistente, un ASV puede ser quimera en una muestra y no en otra.

Casi siempre "consensus". "pooled" solo con muestras técnicamente muy homogéneas.

¿Cuántas quimeras esperar?

Lo importante no es el % de ASVs marcados como quimeras (puede ser 10-40 %), sino el % de reads que representan. Quimeras de baja abundancia son normales. Quimeras dominantes son raras.

Calcular el % de reads perdidos:

total_reads_pre  <- sum(seqtab2)
total_reads_post <- sum(seqtab_nochim)
pct_loss <- round((1 - total_reads_post / total_reads_pre) * 100, 1)
pct_loss
# [1] 2.3

Esperable: 1-5 % de reads perdidos en chimera removal. Si pierdes > 15 %, tu PCR fue ruidosa (muchos ciclos, polimerasa de baja fidelidad), revisa el wet-lab.

Tracking final del pipeline

Buena disciplina: una tabla con cuántos reads sobrevivieron en cada paso:

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

track <- cbind(
    out,                                       # input + filtered
    denoised_F = sapply(dadaFs, getN),
    denoised_R = sapply(dadaRs, getN),
    merged     = sapply(mergers, getN),
    nonchim    = rowSums(seqtab_nochim)
)

colnames(track) <- c("input", "filtered",
                     "denoised_F", "denoised_R",
                     "merged", "nonchim")
rownames(track) <- sample_names
head(track)

#         input  filtered  denoised_F  denoised_R  merged  nonchim
# SampleA 48512     45102       43821       42198   38542    37621
# SampleB 51844     48791       47102       46213   42103    41245
# ...

Esto es lo que se reporta en la sección Métodos de cualquier paper de microbioma serio. Lo veremos en el tutorial 12 (caso completo) integrado en el informe Quarto.

Reportar el tracking

library(tidyverse)

track_df <- as.data.frame(track) |>
    rownames_to_column("sample") |>
    mutate(
        pct_retained = round(nonchim / input * 100, 1)
    )

# Plot
track_df |>
    pivot_longer(cols = input:nonchim, names_to = "step", values_to = "reads") |>
    mutate(step = factor(step, levels = c("input", "filtered",
                                          "denoised_F", "denoised_R",
                                          "merged", "nonchim"))) |>
    ggplot(aes(x = step, y = reads, group = sample)) +
    geom_line(alpha = 0.4, color = "#5F8575") +
    geom_point(alpha = 0.6, color = "#3F6354") +
    scale_y_continuous(labels = scales::comma) +
    labs(x = NULL, y = "Reads por muestra",
         title = "Pérdida de reads a través del pipeline DADA2") +
    theme_minimal()

Lo que esperas ver: caída suave en cada paso. No picos en una muestra concreta. Si una muestra cae bruscamente en merged, probablemente sus R2 eran demasiado malos.

Serializar antes de seguir

saveRDS(seqtab_nochim, "output/seqtab_nochim.rds")
saveRDS(track,         "output/track.rds")

seqtab_nochim es la pieza central del análisis. Todo lo siguiente (taxonomía, phyloseq, diversidad, análisis diferencial) parte de este objeto.

Trampas habituales

  • Solapamiento insuficiente y % merge bajo. Si retornas a tu yo del tutorial 5: trunaste demasiado en R2 y ahora no hay overlap. Solución: re-filtrar con menos truncado en R2, aunque pierdas algo de calidad.
  • Filtrar por longitud sin mirar la distribución. El rango “400-460” es típico de V3-V4. Para otra región es otro. Mira table(nchar(...)) primero.
  • Usar method = "pooled" por defecto. Más sensible pero a veces es demasiado agresivo y descarta secuencias reales como quimeras. Quédate en "consensus" salvo razón concreta.
  • No reportar el tracking. Sin saber cuántos reads sobrevivieron, no puedes contestar “¿por qué la muestra X tiene tan pocos ASVs?”. El tracking es la primera tabla que defiende tu análisis.
  • Recalcular cada vez. El pipeline completo desde filterAndTrim hasta removeBimeraDenovo puede ser una hora. Guarda los .rds intermedios y usa freeze: auto en Quarto.

En la siguiente entrega

Tienes ASVs limpios sin quimeras. Pero no sabes qué bacterias son, todavía son secuencias anónimas. La siguiente entrega: asignación taxonómica con SILVA, comparación con GTDB, y cómo extender la asignación hasta especie cuando es posible. Lo siguiente.