QC y filtrado de reads amplicón

r
bioconductor
metagenomica
El primer paso real de DADA2. Lectura de perfiles de calidad por posición, decisiones de truncado de R1 y R2, eliminación de primers con cutadapt y filtrado con filterAndTrim(), todo en R.

La regla de oro

DADA2 modela los errores reales del secuenciador para distinguir secuencias verdaderas de ruido. Para que el modelo funcione, los reads de entrada deben cumplir tres cosas:

  1. Sin primers (los primers de PCR no son ruido natural, confunden el modelo de errores).
  2. Truncados al punto donde la calidad cae (mejor 5 bases menos que un read entero ruidoso).
  3. Filtrados por número máximo de errores esperados (mejor descartar 5 % de reads que arrastrar basura).

Las tres decisiones tienen impacto directo en la sensibilidad y precisión de los ASVs.

El flujo de este tutorial

library(dada2)

# 1. Localizar archivos FASTQ paired-end
path <- "data/fastq"
fnFs <- sort(list.files(path, pattern = "_R1.fastq.gz", full.names = TRUE))
fnRs <- sort(list.files(path, pattern = "_R2.fastq.gz", full.names = TRUE))

# 2. Visualizar perfiles de calidad
plotQualityProfile(fnFs[1:4])
plotQualityProfile(fnRs[1:4])

# 3. (Si los primers están en los reads) Recortar con cutadapt

# 4. Filtrar y truncar
filtFs <- file.path(path, "filtered", basename(fnFs))
filtRs <- file.path(path, "filtered", basename(fnRs))

out <- filterAndTrim(
    fnFs, filtFs, fnRs, filtRs,
    truncLen = c(240, 200),
    maxN = 0, maxEE = c(2, 2),
    truncQ = 2, rm.phix = TRUE,
    compress = TRUE, multithread = TRUE
)

# 5. Inspeccionar pérdida de reads
head(out)

Vamos paso por paso.

Localizar archivos paired-end

Asumiendo la estructura del tutorial 4:

library(dada2)
library(tidyverse)

path <- "data/fastq"

# Patrones consistentes Illumina: ..._R1.fastq.gz / ..._R2.fastq.gz
fnFs <- sort(list.files(path, pattern = "_R1.fastq.gz", full.names = TRUE))
fnRs <- sort(list.files(path, pattern = "_R2.fastq.gz", full.names = TRUE))

# Verificar simetría
length(fnFs); length(fnRs)
stopifnot(length(fnFs) == length(fnRs))

# Extraer nombres de muestra (todo lo que va antes de _R1)
sample_names <- gsub("_R1.fastq.gz", "", basename(fnFs))
head(sample_names)

Crítico: comprobar que cada R1 tiene su R2 emparejado. Si length(fnFs) != length(fnRs), algo va mal, algún archivo perdido, mal nombrado, o el centro de secuenciación entregó solo single-end por error.

Perfiles de calidad

plotQualityProfile() dibuja la calidad media por posición del read, agregada sobre miles de reads:

plotQualityProfile(fnFs[1:4])  # primeros 4 R1
plotQualityProfile(fnRs[1:4])  # primeros 4 R2

Lectura del plot:

  • Eje X: posición en el read (de 5’ a 3’).
  • Eje Y: calidad Phred (mayor = mejor).
  • Línea verde: calidad media por posición.
  • Línea naranja: cuantiles 25 % y 75 %.
  • Línea roja: fracción de reads que llega a cada posición. En 2×300 MiSeq baja al 100 % hasta el final, ya que todos los reads tienen longitud completa.

Lo esperable en V3-V4 con MiSeq 2×300:

  • R1 (forward): calidad alta (Q30+) hasta posición ~ 230-260, después cae bruscamente.
  • R2 (reverse): degradación notable desde posición ~ 180-220. Esto es normal en Illumina: el R2 siempre es peor que el R1.

Decisiones de truncado: truncLen

truncLen = c(forward, reverse) corta los reads a longitud fija. La decisión clave de DADA2.

Regla práctica:

  • Corta donde la mediana de calidad baja de Q30.
  • Pero asegúrate de que R1 + R2 ≥ longitud del amplicón + 20 bp de solapamiento. Sin solapamiento, mergePairs falla.

Cálculo para V3-V4 (~ 460 bp con primers, ~ 410 bp sin primers):

Sin primers Con 20 bp solap. Mínimo R1+R2
Amplicón V3-V4 410 bp 430 bp 430 bp

Si tus reads son 2×300 y los R2 degradan a posición 200: truncLen = c(240, 200) da 440 bp = solapa ~ 30 bp. OK.

Si los R2 degradan a posición 150: truncLen = c(280, 150) da 430 bp = solapa ~ 20 bp. Mínimo, no recomendable. Mejor descartar R2 muy degradado.

Trampa: con V4 (~ 250 bp), no hace falta tanto solapamiento. truncLen = c(230, 180) ya da 410 bp = 160 de solapamiento. Cómodo.

Decisiones de calidad: maxEE, truncQ, maxN

Aunque trunques, dentro del read siguen quedando errores. Tres parámetros controlan qué se filtra:

maxEE: expected errors máximo

maxEE = c(2, 2) significa: descarta reads con más de 2 errores esperados (sumando todas las posiciones, según el modelo Phred).

Es la decisión más importante del filtrado. Más estricto (maxEE = c(1, 1)) reduce ASVs ruidosos pero pierde reads. Más laxo (maxEE = c(2, 4)) conserva muestras de baja calidad a costa de más ruido downstream.

Default razonable: c(2, 2) para data limpia, c(2, 4) o c(2, 5) si los R2 son malos pero quieres conservarlos.

truncQ: trunca también en cuanto haya una base de calidad baja

truncQ = 2 (default): si encuentra una base con Q ≤ 2 antes de truncLen, trunca ahí mismo. Esto recorta dinámicamente reads con caídas tempranas.

maxN, rm.phix

  • maxN = 0: ningún read con N (base indeterminada) pasa. DADA2 no acepta Ns en el modelo de errores. Sin negociación.
  • rm.phix = TRUE: elimina reads que mapean al control PhiX (spike-in que Illumina añade para calibrar). En MiSeq con baja diversidad (como amplicón), PhiX puede ser 10-20 % de los reads.

Ejecutar filterAndTrim()

Output a una carpeta data/fastq/filtered/:

# Crear paths de salida
filtFs <- file.path(path, "filtered", paste0(sample_names, "_F_filt.fastq.gz"))
filtRs <- file.path(path, "filtered", paste0(sample_names, "_R_filt.fastq.gz"))
names(filtFs) <- sample_names
names(filtRs) <- sample_names

# Filtrar y truncar
out <- filterAndTrim(
    fnFs, filtFs,
    fnRs, filtRs,
    truncLen   = c(240, 200),
    maxN       = 0,
    maxEE      = c(2, 2),
    truncQ     = 2,
    rm.phix    = TRUE,
    compress   = TRUE,
    multithread = TRUE
)

head(out)
#                          reads.in reads.out
# SampleA_R1.fastq.gz         48512     45102
# SampleB_R1.fastq.gz         51844     48791
# ...

La matriz out te dice cuántos reads entraron y cuántos sobrevivieron por muestra. Lo importante:

  • % retención: ideal > 80 %. Si pierdes > 30 %, revisa los criterios.
  • Variación entre muestras: si una muestra retiene 90 % y otra 40 %, esa segunda muestra tiene problemas de calidad, investigar (¿extracción mala? ¿lote distinto?).

Reportar el filtrado

Para informe Quarto:

out_df <- as.data.frame(out) |>
    rownames_to_column("file") |>
    mutate(
        sample = gsub("_R1.fastq.gz", "", file),
        retention = round(reads.out / reads.in * 100, 1)
    ) |>
    select(sample, reads.in, reads.out, retention) |>
    arrange(retention)

knitr::kable(out_df, caption = "Filtrado DADA2: reads por muestra")

Y un plot de retención:

ggplot(out_df, aes(x = reorder(sample, retention), y = retention)) +
    geom_col(fill = "#5F8575") +
    geom_hline(yintercept = 80, linetype = "dashed", color = "grey50") +
    coord_flip() +
    labs(x = NULL, y = "% reads retenidos",
         title = "Retención por muestra tras filterAndTrim") +
    theme_minimal()

Visualizar esto antes de continuar es disciplina básica. Una muestra con 20 % de retención no va a dar ASVs fiables, mejor descartarla ahora que arrastrarla.

¿Y los primers?

Hay dos escenarios:

Escenario A: el centro de secuenciación ya recortó los primers

Lo más habitual hoy. Cada read empieza directamente con la secuencia variable. Salta el siguiente bloque.

Escenario B: los primers están en los reads

Mira el README del centro de secuenciación o haz head a un FASTQ:

zcat SampleA_R1.fastq.gz | head -4

Si las primeras bases coinciden con la secuencia del primer (ej. CCTACGGGNGGCWGCAG para 341F V3-V4), están dentro. Tienes que recortarlos antes de filterAndTrim.

Recomendación: cutadapt (Python, CLI, mucho más robusto que el trimLeft de DADA2):

# Para cada muestra:
cutadapt \
    -g CCTACGGGNGGCWGCAG \
    -G GACTACHVGGGTATCTAATCC \
    -o SampleA_R1_trimmed.fastq.gz \
    -p SampleA_R2_trimmed.fastq.gz \
    SampleA_R1.fastq.gz SampleA_R2.fastq.gz \
    --discard-untrimmed \
    --minimum-length 50

--discard-untrimmed descarta reads donde el primer no se encontró (probablemente reads de PhiX o contaminación).

Para batch sobre muchas muestras, escribir un script bash o usar Herper (paquete R que llama a cutadapt). El libro asociado a esta ruta cubre el batch en serio.

Trampas habituales

  • Truncar agresivo y romper el solapamiento. Si R1 + R2 < amplicón + 12, mergePairs no encontrará overlap y perderás casi todos los reads en el siguiente paso. Calcula primero la longitud esperada del amplicón.
  • maxEE = c(2, 4) “porque es lo que vi en un tutorial”. Los defaults dependen de tu calidad real. Si tus R2 son muy malos, c(2, 5) puede valer. Si son excelentes, c(2, 2) o más estricto.
  • No comprobar PhiX. En MiSeq amplicón, PhiX puede ser 15-20 % de reads. Sin rm.phix = TRUE, esos reads contaminan los ASVs como una especie fantasma.
  • Procesar muestras con calidades muy distintas en el mismo filterAndTrim. Los parámetros deberían ajustarse por lote. Si tienes batches con calidad muy distinta, procesa cada lote por separado.
  • Saltarse plotQualityProfile(). Decidir truncLen “por intuición” sin ver los plots da resultados malos. Es el paso de 30 segundos que más decisiones informa.

En la siguiente entrega

Tienes reads limpios. La siguiente entrega es el corazón estadístico de DADA2: aprender el modelo de errores del secuenciador con learnErrors(), denoising con dada(), y entender cómo el algoritmo distingue una variante real de un error de secuenciación. Lo siguiente.