Objetivos

Paquetes

library(dada2)
library(Biostrings)

  • dada2: Herramienta principal para el análisis de datos de secuenciación de amplicones (ej. 16S, ITS). Permite filtrar lecturas, corregir errores y asignar taxonomía, obteniendo variantes de secuencia exactas (ASVs).
    📖 Documentación
    📑 Artículo (Nature Methods, 2016)

  • Por información más detallada de cada función asociada a los paquetes utilice help(“función”) en RStudio (por ej. help(“learnErrors”))

1. Filtrado de lecturas

Asegúrate de trabajar en la carpeta donde están los archivos. Si no, usa setwd(“ruta/a/los/archivos”).

Ruta y listado de archivos

path <- setwd()

Verificar si la ruta es correcta y si todos los archivos están presentes

list.files(path)

Leer archivos fastq forward y reverse

fnFs <- sort(list.files(path, pattern="_1.fastq", full.names = TRUE))
fnRs <- sort(list.files(path, pattern="_2.fastq", full.names = TRUE))

Función sort()

En el flujo de trabajo, la función sort() se utiliza para ordenar alfabéticamente los nombres de los archivos FASTQ listados con list.files().
Esto asegura que las lecturas forward (_1.fastq) y reverse (_2.fastq) queden organizadas de manera correcta antes de ser procesadas, evitando desajustes entre las muestras.

Extraer nombres de muestra, asumiendo que los archivos tienen el formato: MUESTRANOMBRE_XXX.fastq

sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)

Leer metadatos de las muestras

sample.metadata <-read.csv("sampleData_alimentos.csv")

Verificar calidad de las lecturas por base

Lecturas forward

plotQualityProfile(fnFs[1:2])

Lecturas reverse

plotQualityProfile(fnRs[1:2])

Función plotQualityProfile()

La función plotQualityProfile() del paquete dada2 permite visualizar la calidad de las lecturas en archivos FASTQ.
En este caso, se aplica a las dos primeras muestras (fnFs[1:2]) para inspeccionar de forma preliminar la distribución de calidad a lo largo de las posiciones de las secuencias.

Este gráfico es muy útil porque:
- Muestra la caída típica de la calidad hacia el extremo derecho de las lecturas.
- Permite decidir hasta qué posición recortar (trimming) las secuencias.
- Ayuda a identificar problemas en la calidad de los datos de secuenciación.

Notas:

Las lecturas reverse suelen tener una calidad mucho peor que las forward, especialmente hacia el final.

Truncar las lecturas aproximadamente donde el puntaje de calidad baja por debajo de 30–35.

Crear archivos filtrados en el subdirectorio filtered/

filtFs <- file.path(path, "filtered", paste0(sample.names, "_F_filt.fastq.gz"))
filtRs <- file.path(path, "filtered", paste0(sample.names, "_R_filt.fastq.gz"))

Filtrar lecturas

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

save(out, file="out")

Función filterAndTrim()

La función filterAndTrim() del paquete dada2 se utiliza para filtrar y recortar las lecturas crudas en formato FASTQ, eliminando secuencias de baja calidad antes del análisis.

truncLen=c(290,200): recorta las lecturas forward a 290 bases y las reverse a 200 bases.

maxN=0: descarta cualquier lectura con nucleótidos ambiguos (N).

maxEE=c(2,2): elimina lecturas con un número esperado de errores mayor a 2 (tanto en forward como en reverse).

truncQ=2: recorta la lectura en el primer punto donde la calidad baja a Q=2 o menos.

rm.phix=TRUE: elimina contaminantes derivados del control PhiX.

compress=TRUE: guarda los archivos filtrados en formato comprimido.

multithread=FALSE: ejecuta el proceso en un solo hilo (en Linux/Mac puede ponerse TRUE para usar múltiples núcleos).

Se puede utilizar trimLeft para remover secuencias de primers.

2. Dereplicar, aprender errores y ejecutar DADA2

Lecturas Forward

errF <- learnErrors(filtFs, multithread=FALSE)
save(errF, file="errF")

Lecturas Reverse

errR <- learnErrors(filtRs, multithread=FALSE)
save(errR, file="errR")

Función learnErrors()

La función learnErrors() del paquete dada2 se utiliza para aprender los perfiles de error a partir de las lecturas filtradas.
El algoritmo de DADA2 modela explícitamente las probabilidades de error de cada posición de la secuencia, lo cual es esencial para distinguir entre verdaderas variantes biológicas y errores de secuenciación.

filtFs: conjunto de lecturas forward filtradas que se utilizan para estimar el modelo de errores.

multithread=FALSE: ejecuta el cálculo en un solo núcleo; en Linux/Mac puede ponerse en TRUE para acelerar el proceso.

Visualizar las tasas de error estimadas

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

Función plotErrors()

La función plotErrors() del paquete dada2 se utiliza para visualizar los perfiles de error aprendidos con la función learnErrors().

Descripción:
- Muestra la tasa observada de errores (puntos) en comparación con las tasas de error esperadas según el modelo (líneas).
- Cada gráfico corresponde a un tipo de sustitución de bases (A→C, A→G, etc.).
- Sirve como control de calidad para verificar que el modelo de errores se ajusta correctamente a los datos.

Parámetros comunes:
- err: objeto de errores obtenido con learnErrors().
- nominalQ=TRUE: agrega al gráfico las líneas correspondientes a los valores teóricos de calidad (Q-scores) para comparación.

Resultado esperado:
Un gráfico donde las curvas del modelo (líneas) se alinean de forma razonable con los puntos observados. Esto indica que la estimación de errores es adecuada y puede usarse para la inferencia de variantes de secuencia exactas (ASVs).

Dereplicar lecturas idénticas

derepFs <- derepFastq(filtFs, verbose=TRUE)
names(derepFs) <- sample.names
save(derepFs, file= "derepFs")

derepRs <- derepFastq(filtRs, verbose=TRUE)
names(derepRs) <- sample.names
save(derepRs, file= "derepRs")

Función derepFastq()

La función derepFastq() del paquete dada2 se utiliza para colapsar lecturas idénticas en una única secuencia representativa.
Este proceso se llama dereplicación y tiene como objetivo reducir la redundancia en los datos y acelerar el análisis posterior.

Ejecutar dada2

Lecturas Forward

dadaFs <- dada(derepFs, err=errF, pool = TRUE , multithread=TRUE)
save(dadaFs, file="dadaFs")

Lecturas Reverse

dadaRs <- dada(derepRs, err=errR, pool = TRUE , multithread=TRUE)
save(dadaRs, file="dadaRs")

Función dada()

La función dada() es el núcleo del paquete dada2, y se utiliza para realizar la inferencia de variantes de secuencia exactas (ASVs) a partir de las lecturas dereplicadas.
En lugar de agrupar lecturas en OTUs por un umbral de similitud, DADA2 utiliza un modelo estadístico para distinguir errores de secuenciación de verdaderas diferencias biológicas.

derepFs: lecturas forward previamente dereplicadas.

err=errF: modelo de errores aprendido con learnErrors(), usado para evaluar la probabilidad de que una diferencia entre secuencias sea un error o una variante real.

pool=TRUE: combina todas las muestras para aumentar la sensibilidad en la detección de variantes poco abundantes (si se deja en FALSE, cada muestra se procesa de manera independiente).

multithread=TRUE: permite usar múltiples núcleos para acelerar el cálculo.

Un objeto de clase dada almacena la información de la ASVs a partir de las lecturas dereplicadas. Dentro de cada elemento se encuentra:

$denoised: las secuencias inferidas como variantes reales (ASVs).

$clustering: asignación de cada lectura a su secuencia representativa.

$map: correspondencia entre lecturas originales y ASVs.

$birth_subs / $birth_q: detalles de la primera diferencia de calidad usada en la inferencia.

$err_in / $err_out: tasas de error utilizadas y corregidas en el proceso.

Inspeccionando el objeto dada:

dadaFs[[1]]
dadaRs[[1]]

mean(dadaFs[[1]]$quality)
mean(dadaRs[[1]]$quality)

3. Construir la tabla de secuencias

Unir las lecturas forward y reverse y construir la tabla de secuencias

mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs, verbose=TRUE)
save(mergers, file="mergers")

Inspeccionar el data.frame de fusiones de la primera muestra

head(mergers[[1]])

Función mergePairs()

La función mergePairs() del paquete dada2 se utiliza para unir las lecturas forward y reverse de una misma muestra, siempre que presenten una región solapada suficiente y consistente.
Este paso es crucial en análisis de secuencias de amplicones, ya que la región resultante es de mayor calidad y cubre toda la zona de interés.

  • sequence: la secuencia consenso resultante de la fusió o merge entre forward y reverse.
  • abundance: el número de lecturas que dieron lugar a esa secuencia fusionada.
  • forward: índice de la ASV en dadaFs que aportó la lectura forward.
  • reverse: índice de la ASV en dadaRs que aportó la lectura reverse.
  • nmatch: número de bases que coincidieron en la región solapada.
  • nmismatch: número de discrepancias detectadas en la región solapada.
  • nindel: número de inserciones/deleciones observadas.
  • prefer: cuál de las dos lecturas (forward o reverse) se usó para resolver un conflicto en la fusión.

Ahora podemos construir una tabla de variantes de secuencia de amplicón (ASV), una versión de mayor resolución de la tabla de OTUs producida por otros métodos.

seqtab <- makeSequenceTable(mergers)
save(seqtab, file="seqtab")

Función makeSequenceTable()

La función makeSequenceTable() del paquete dada2 se utiliza para construir la tabla de secuencias a partir de los objetos resultantes de mergePairs().
En esta tabla, cada fila corresponde a una muestra y cada columna a una ASV, de forma análoga a una matriz de conteos.

Verificar el número de variantes de secuencia

dim(seqtab)

Inspeccionar la distribución de longitudes de secuencia

table(nchar(getSequences(seqtab)))

seq_dis <- table(nchar(getSequences(seqtab)))

barplot(seq_dis,
        col = "lightblue",
        border = "darkblue",
        main = "Distribución de longitudes de ASVs",
        xlab = "Longitud de secuencia (bp)",
        ylab = "Número de ASVs")

Evaluación de longitudes de secuencia

Un paso de control de calidad importante consiste en verificar la longitud de las secuencias presentes en la tabla (seqtab).

Notas

  • Lo habitual es observar una sola longitud dominante, correspondiente a la región de amplicón objetivo (por ejemplo, ~250–300 bp para 16S, dependiendo de los primers usados).

  • Pueden aparecer unas pocas secuencias con longitudes atípicas (más cortas o más largas), que suelen ser artefactos y pueden eliminarse en pasos posteriores.

  • Las secuencias que son mucho más largas o más cortas de lo esperado pueden ser el resultado de unión inespecífica de los primers.

  • Puede eliminar las secuencias con longitudes no objetivo.

seqtab2 <- seqtab[,nchar(colnames(seqtab)) %in% seq(440,466)]
save(seqtab2, file="seqtab2")

Filtrado por longitud del amplicón

Después de construir la tabla de secuencias, es recomendable conservar únicamente las variantes cuya longitud corresponde al rango esperado del amplicón (definido por los primers utilizados).
Esto ayuda a eliminar posibles artefactos de secuenciación o ensamblaje.

colnames(seqtab): obtiene las secuencias (ASVs) como nombres de columna.

nchar(…): calcula la longitud de cada secuencia.

seq(440,466): define el rango esperado de longitudes (por ejemplo, de 440 a 466 bp).

%in%: selecciona solo aquellas secuencias cuya longitud cae dentro de ese rango.

El resultado (seqtab2) es una nueva tabla que incluye únicamente las secuencias con la longitud biológica esperada.

Verificar nuevamente la distribución

table(nchar(getSequences(seqtab2)))

seq_dis2 <- table(nchar(getSequences(seqtab2)))

barplot(seq_dis2,
        col = "lightblue",
        border = "darkblue",
        main = "Distribución de longitudes de ASVs filtradas",
        xlab = "Longitud de secuencia (bp)",
        ylab = "Número de ASVs")

Eliminar quimeras

seqtab.nochim <- removeBimeraDenovo(seqtab2, method="consensus", multithread=TRUE, verbose=TRUE)
save(seqtab.nochim,file="seqtab.nochim")
dim(seqtab.nochim)
sum(seqtab.nochim)/sum(seqtab2)

Función removeBimeraDenovo()

La función removeBimeraDenovo() del paquete dada2 se utiliza para detectar y eliminar secuencias quiméricas (chimeras).
Las quimeras son artefactos de PCR que resultan de la combinación de fragmentos de diferentes secuencias durante la amplificación y por lo tanto no representan variantes biológicas reales.

seqtab2: tabla de secuencias filtrada por longitud esperada del amplicón.

method=“consensus”: identifica quimeras combinando información de todas las muestras, aumentando la robustez de la detección.

multithread=TRUE: acelera el cálculo usando múltiples núcleos.

verbose=TRUE: muestra información detallada sobre el número de quimeras detectadas y eliminadas.

Inspeccionar la longitud media de las secuencias

dna <- DNAStringSet(getSequences(seqtab.nochim))
meansequencelength <- as.data.frame(table(nchar(getSequences(dna))))
meansequencelength$Var1 <- as.numeric(as.character(meansequencelength$Var1))
meansequencelength$total <- meansequencelength$Var1*meansequencelength$Freq
sum(meansequencelength$total)/sum(meansequencelength$Freq)

Función DNAStringSet()

La función DNAStringSet() del paquete Biostrings se utiliza para crear un objeto especializado que almacena conjuntos de secuencias de ADN.
Este tipo de objeto es mucho más eficiente que un simple vector de caracteres, ya que está optimizado para trabajar con datos biológicos y permite aplicar fácilmente funciones de manipulación y análisis de secuencias.

getSequences(seqtab.nochim): extrae todas las secuencias únicas (ASVs) de la tabla sin quimeras.

Seguimiento de las lecturas a lo largo del pipeline

Como verificación final de nuestro progreso, revisaremos el número de lecturas que pasaron por cada paso del pipeline:

getN <- function(x) sum(getUniques(x))
track <- cbind(out, sapply(dadaFs, getN), sapply(dadaRs, getN), sapply(mergers, getN),rowSums(seqtab.nochim))

Seguimiento del número de lecturas en cada paso

Para evaluar la eficiencia del pipeline de DADA2, es recomendable llevar un control del número de lecturas que se mantienen en cada etapa del procesamiento.

Si se procesa una sola muestra, elimine las llamadas a sapply: por ejemplo, reemplace sapply(dadaFs, getN) con getN(dadaFs)

colnames(track) <- c("input", "filtered", "denoisedF", "denoisedR", "merged", "nonchim")
rownames(track) <- sample.names
head(track)

Etiquetado de la tabla de seguimiento

Para mejorar la interpretación de la matriz track, se añaden nombres descriptivos a las columnas (etapas del pipeline) y a las filas (muestras).

4. Asignar taxonomía

taxa <- assignTaxonomy(seqtab.nochim, "ruta/a/silva_nr_v132_train_set.fa.gz", multithread=TRUE)

Función assignTaxonomy()

La función assignTaxonomy() del paquete dada2 se utiliza para asignar taxonomía a las variantes de secuencia exactas (ASVs) obtenidas después del filtrado, corrección de errores y eliminación de quimeras.

Descripción:
- Compara cada secuencia de la tabla (seqtab.nochim) contra una base de datos de referencia taxonómica (por ejemplo, SILVA, RDP o Greengenes).
- Utiliza un clasificador basado en la probabilidad Naïve Bayesian classifier, entrenado con secuencias anotadas.
- Devuelve un objeto donde cada ASV queda anotada con su taxonomía jerárquica:
Reino → Filo → Clase → Orden → Familia → Género

Parámetros principales:
- seqtab.nochim: tabla de secuencias sin quimeras sobre la cual se asignará la taxonomía.
- "ruta/a/silva_nr_v132_train_set.fa.gz": archivo de referencia previamente entrenado (en este ejemplo, SILVA v132).
- multithread=TRUE: permite usar múltiples núcleos para acelerar el proceso.

Resultado esperado:
- Una matriz o data.frame donde:
- Cada fila corresponde a un ASV.
- Cada columna representa un nivel taxonómico.
- Las celdas contienen la clasificación más probable para cada secuencia.

taxa <- addSpecies(taxa, "ruta/a/silva_species_assignment_v132.fa.gz")

Función addSpecies()

La función addSpecies() del paquete dada2 se utiliza para afinar la asignación taxonómica hasta el nivel de especie, siempre que la base de datos de referencia lo permita.

Descripción:
- Toma como entrada el objeto taxa previamente generado con assignTaxonomy().
- Compara cada ASV contra una base de datos que contiene secuencias de referencia a nivel de especie.
- Asigna la especie más probable cuando existe una coincidencia exacta o muy cercana.

Parámetros principales:
- taxa: objeto de taxonomía obtenido con assignTaxonomy().
- "ruta/a/silva_species_assignment_v132.fa.gz": archivo de referencia con asignaciones a nivel de especie (en este caso, de la base SILVA v132).

Resultado esperado:
- El objeto taxa se actualiza con una columna adicional llamada Especie, que contiene la anotación más probable para cada ASV.
- Si no hay coincidencia confiable, la celda queda vacía (NA).

Importancia:
La asignación hasta especie puede ser muy útil para estudios clínicos o ecológicos detallados, aunque debe interpretarse con precaución ya que depende de la calidad y cobertura de la base de datos de referencia.

save(taxa, file="taxa")

La base de datos para la asignación de taxonomía pueden descargarse desde este link.
Puede explorar las versión mas reciente de la base de datos SILVA aquí.

Por razones de tiempo la asignación taxonómica no se realizará durante el taller.

Preguntas a responder durante el taller

  1. ¿Qué entiende por ASV?
  2. ¿Qué pasos son necesarios para la generación de las ASVs?
  3. ¿A qué podría atribuirse un número elevado de secuencias quiméricas?
  4. ¿Cómo se asigna taxonomía a los ASVs?