Learn errors + denoising

r
bioconductor
metagenomica
El corazón estadístico de DADA2. Cómo se aprende el modelo de errores del secuenciador con learnErrors(), cómo dada() lo aplica para distinguir variantes reales de ruido, y qué pinta tienen los plots diagnósticos que confirman que el modelo va bien.

Por qué este paso es distinto

La mayoría de pipelines viejos (mothur, USEARCH) agrupaban reads similares y llamaban OTUs a los clusters. La asunción: reads con < 3 % de diferencia son la misma “especie operacional”.

DADA2 hace algo radicalmente distinto: modela los errores reales del secuenciador y, con ese modelo, calcula la probabilidad de que dos reads casi idénticos sean dos secuencias biológicamente distintas o el mismo organismo con un error de secuenciación.

Resultado: resolución de un único nucleótido. Donde mothur veía un OTU, DADA2 puede ver tres ASVs distintos si la estadística los soporta.

Esto pide dos pasos:

  1. learnErrors(): aprende cómo se equivoca tu secuenciador en tu corrida concreta.
  2. dada(): usa ese modelo para inferir secuencias reales.

El modelo de errores

Cada secuenciador comete errores distintos según:

  • La base correcta (A, C, G, T).
  • La base que llamó (A → C, A → G, A → T, …). 12 transiciones posibles.
  • La calidad Phred que asignó.

learnErrors() ajusta un modelo paramétrico: para cada (base correcta, base llamada, Q-score) estima la probabilidad de error. Con esto sabe cuándo creer una calidad reportada y cuándo no.

Necesita un volumen mínimo de reads para que la estimación sea estable. El default es ~ 100 millones de bases (nbases = 1e8). Por eso learnErrors puede tardar varios minutos.

Ejecutar learnErrors()

Para R1 y R2 por separado, sus perfiles de error son distintos:

library(dada2)

# filtFs y filtRs vienen del filterAndTrim del tutorial anterior
errF <- learnErrors(filtFs, multithread = TRUE)
errR <- learnErrors(filtRs, multithread = TRUE)

Output típico en consola:

105712320 total bases in 440468 reads from 12 samples will be used for learning the error rates.

DADA2 ha tomado reads de las primeras 12 muestras hasta acumular ~ 100M bases. Si tu dataset es pequeño (pocas muestras o pocos reads), usará todos.

Plots diagnósticos del modelo

Lo siguiente, y crítico, es verificar visualmente que el modelo es razonable:

plotErrors(errF, nominalQ = TRUE)
plotErrors(errR, nominalQ = TRUE)

Cada plot muestra 12 paneles (una por transición A2A, A2C, A2G, A2T, C2A, …, T2T):

  • Eje X: calidad Phred (Q-score).
  • Eje Y: tasa de error (log).
  • Puntos negros: tasas observadas en los datos.
  • Línea negra: ajuste del modelo paramétrico.
  • Línea roja: tasa esperada según la calidad Phred nominal (la “promesa” del secuenciador).

Lo que esperas ver:

  • La línea negra (modelo) baja a medida que sube Q (más calidad → menos error). Bien.
  • Los puntos siguen la línea, el modelo casa con los datos.
  • La línea roja nominal suele estar por debajo de la observada, el secuenciador es algo más ruidoso de lo que promete (normal).

Lo que sospechar:

  • Línea negra horizontal (el modelo no aprendió relación calidad-error), algo va mal con los reads.
  • Puntos muy por encima del modelo a Qs altas, calidades sobreestimadas. El filtrado de filterAndTrim fue demasiado laxo.
  • Una transición concreta (ej. C2A) mucho más alta que las demás, sesgo de química Illumina específico, normalmente no problema.

Estos plots son lo que separa el análisis con criterio del análisis “a ciegas”. Mirar los dos (R1 y R2) cada vez.

Denoising con dada()

Con el modelo aprendido, infiere ASVs:

dadaFs <- dada(filtFs, err = errF, multithread = TRUE)
dadaRs <- dada(filtRs, err = errR, multithread = TRUE)

Esto puede tardar bastante (segundos a minutos por muestra). Es el cuello de botella del pipeline.

Inspección de una muestra:

dadaFs[[1]]
# dada-class: object describing DADA2 denoising results
# 145 sequence variants were inferred from 12814 input unique sequences.
# Key parameters: OMEGA_A = 1e-40, OMEGA_C = 1e-40, BAND_SIZE = 16

145 ASVs inferidos de 12 814 secuencias únicas. Lo que ves es una matriz de N muestras × K ASVs por muestra.

Pooling: tres modos

dada() tiene un argumento clave: pool. Controla cómo trata las muestras conjuntamente:

Modo Cuándo Qué hace
pool = FALSE (default) Por defecto Cada muestra procesada por separado. Rápido. Pierde ASVs raros entre muestras.
pool = "pseudo" Recomendado para casi todo Procesa una muestra, comparte info con las demás, vuelve a procesar. Buen compromiso.
pool = TRUE Estudios pequeños (< 50 muestras) Todas las muestras juntas. Más sensible para ASVs raros. Mucho más lento.

Por defecto Benjamin Callahan (autor de DADA2) recomienda pool = "pseudo" para análisis comparativos. Lo dejas en:

dadaFs <- dada(filtFs, err = errF, pool = "pseudo", multithread = TRUE)
dadaRs <- dada(filtRs, err = errR, pool = "pseudo", multithread = TRUE)

Si los datasets son grandes (cientos de muestras), pool = FALSE es el único viable por tiempo.

¿Cuántos ASVs deberías esperar?

Depende del tipo de muestra, pero rangos típicos:

  • Heces humanas sanas: 200-800 ASVs por muestra, 2 000-5 000 ASVs en el dataset completo.
  • Suelo: 1 000-10 000 ASVs por muestra (diversidad altísima).
  • Sangre, biopsia (baja biomasa): < 100 ASVs por muestra. Sospechar contaminación si > 200.
  • Probiótico definido (cultivo): 5-20 ASVs (las cepas que hay + algo de ruido).

Si los números se desvían mucho de estos rangos, revisa: ¿filtraste bien? ¿usaste la región correcta? ¿la muestra es lo que crees que es?

Diagnósticos rápidos tras dada()

Tres checks que conviene automatizar:

1. Reads únicos vs ASVs

data.frame(
    sample = sample_names,
    uniques_F = sapply(dadaFs, function(x) length(getUniques(x))),
    asvs_F    = sapply(dadaFs, function(x) sum(getUniques(x) > 0))
) |> head()

La diferencia entre “uniques” (secuencias distintas vistas) y “ASVs” (variantes biológicas inferidas) es lo que DADA2 colapsó como errores.

2. Distribución del tamaño de ASV

Para V3-V4, esperas ASVs de ~ 400-460 bp. Cualquier cosa más corta o larga es sospechosa:

seqs <- getUniques(dadaFs[[1]]) |> names()
hist(nchar(seqs), breaks = 30,
     main = "Longitud de ASVs en SampleA",
     xlab = "Longitud (bp)")

ASVs muy cortos (< 250 bp) suelen ser fragmentos no específicos. Los filtraremos en el siguiente tutorial al construir la tabla de secuencias.

3. Tracking de reads a través del pipeline

Es la métrica que tienes que mantener desde el principio:

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

track <- data.frame(
    sample   = sample_names,
    input    = out[, "reads.in"],
    filtered = out[, "reads.out"],
    denoised_F = sapply(dadaFs, getN),
    denoised_R = sapply(dadaRs, getN)
)
head(track)

Esto te dice cuántos reads quedan tras cada paso. La caída entre filtered y denoised_* típica es de 5-15 %. Si pierdes > 30 %, algo va mal, quizá el modelo de errores no funcionó (revisar plotErrors).

El patrón idiomático

El flujo de los dos pasos juntos:

library(dada2)

# Aprende errores (10-30 min para 50 muestras)
errF <- learnErrors(filtFs, multithread = TRUE)
errR <- learnErrors(filtRs, multithread = TRUE)

# Diagnóstico — siempre
plotErrors(errF, nominalQ = TRUE)
plotErrors(errR, nominalQ = TRUE)

# Denoising (más lento; reserva tiempo)
dadaFs <- dada(filtFs, err = errF, pool = "pseudo", multithread = TRUE)
dadaRs <- dada(filtRs, err = errR, pool = "pseudo", multithread = TRUE)

# Cuántos ASVs
sum(sapply(dadaFs, function(x) length(getUniques(x))))

Cache: guarda el output

dada() puede tardar 30-60 minutos en datasets de cientos de muestras. Guarda el output para no recalcular en cada render:

saveRDS(dadaFs, "output/dadaFs.rds")
saveRDS(dadaRs, "output/dadaRs.rds")
saveRDS(errF, "output/errF.rds")
saveRDS(errR, "output/errR.rds")

Y en el chunk Quarto correspondiente:

#| cache: true
#| dependson: ["learn-errors", "filter-trim"]

O mejor, en proyecto Quarto: freeze: auto en _quarto.yml (ver tutorial de cache/freeze).

Trampas habituales

  • Ignorar plotErrors(). Es el único punto donde puedes detectar que el modelo de errores no se ajustó bien. Saltarlo es procesar a ciegas.
  • Usar pool = TRUE con cientos de muestras. Va a tardar horas o días. Si tienes > 100 muestras, pool = "pseudo" casi siempre.
  • Reusar errF entre datasets distintos. El modelo es específico de la corrida de secuenciación. Una corrida = un learnErrors. Si tienes dos runs juntas, procesa por separado y combina al final.
  • No serializar el output. Si después tu render Quarto recalcula dada() cada vez, te tomará 30 minutos por iteración. Cache con freeze: auto o saveRDS explícito.
  • Asumir que más ASVs es mejor. Más ASVs en una muestra puede indicar más diversidad real, o indicar que el filtrado fue laxo y se colaron ruidos. Cruza con la diversidad alfa (lo veremos en el tutorial 10).

En la siguiente entrega

Tienes ASVs inferidos por R1 y R2 independientemente. Pero todavía no tienes la secuencia completa del amplicón. La siguiente entrega: merge pairs (combinar R1 y R2 en una sola secuencia) y remoción de quimeras (los artefactos de PCR que confunden taxonomía). Lo siguiente.