Control de calidad y recorte de primers con Cutadapt

Importante

Antes de empezar: asegúrate de estar en el directorio del proyecto con el ambiente activado:

cd ~/metabarcoding-code
conda activate metabarcoding

Objetivo

En esta clase aprenderemos a realizar control de calidad (QA/QC) en archivos FASTQ usando herramientas bioinformáticas.

Pasos principales

Conteo de lecturas
Control de calidad con FastQC y MultiQC
Revisión de resultados (datos crudos)
Recorte de primers con Cutadapt
Revisión de resultados (datos recortados)

1. Conteo de lecturas

Cada lectura ocupa 4 líneas en un archivo FASTQ. Para contar lecturas:

for f in raw_fastqs/*.fastq; do printf "%s\t%s\n" "$(basename "$f")" "$(expr $(wc -l < "$f") / 4)"; done > final_data/conteo_lecturas.tsv

Visualizar en terminal

Ordena las muestras por número de lecturas (de mayor a menor) para detectar rápidamente cuáles tienen pocas lecturas:

column -t final_data/conteo_lecturas.tsv | sort -k2 -rn | head -20

Visualización opcional en R

Si quieres una gráfica para explorar el conteo por muestra, abre R en la terminal:

Dentro de R, corre:

library(tidyverse)

conteo <- read_tsv("final_data/conteo_lecturas.tsv",
                   col_names = c("archivo", "lecturas"))

conteo <- conteo %>%
  mutate(
    muestra = str_remove(archivo, "_R[12]_001\\.fastq.*"),
    sentido = str_extract(archivo, "R[12]")
  )

png("final_data/conteo_lecturas.png", width = 900, height = 700)

ggplot(conteo, aes(x = reorder(muestra, lecturas), y = lecturas, fill = sentido)) +
  geom_col(position = "dodge") +
  geom_hline(yintercept = 10000, linetype = "dashed", color = "red") +
  annotate("text", x = 1, y = 10000, label = "10,000 reads",
           hjust = -0.1, vjust = -0.5, color = "red", size = 3.5) +
  coord_flip() +
  scale_y_continuous(labels = scales::comma) +
  labs(
    title = "Lecturas por muestra (crudas)",
    x = NULL, y = "Numero de lecturas", fill = "Sentido"
  ) +
  theme_minimal(base_size = 11)

dev.off()
q()

La línea roja marca 10,000 lecturas — muestras por debajo merecen revisión antes de continuar.

2. Control de calidad con FastQC y MultiQC

Crear directorios y correr fastqc

Primero creamos el directorio de salida. final_data/ no existe aún en el repositorio — es donde se guardarán todos los resultados del pipeline:

mkdir -p final_data/fastqc final_data/multiqc final_data/fastqc_trimmed final_data/multiqc_trimmed

fastqc -q -t 4 -o final_data/fastqc raw_fastqs/*.fastq

Verificar conteo de reportes

ls -1 final_data/fastqc/*_fastqc.html | wc -l && echo "FastQC OK"

Ejecutar MultiQC (resume todos los FastQC)

multiqc final_data/fastqc -o final_data/multiqc

Verificar reporte MultiQC

ls final_data/multiqc/multiqc_report.html && echo "MultiQC OK"

3. Revisión de resultados (datos crudos)

Descargar resultados a tu computadora

Desde tu terminal local, usa rsync para descargar toda la carpeta final_data/. Solo baja archivos nuevos o modificados:

rsync -avz usrXX@200.23.162.240:/home/usrXX/metabarcoding-code/final_data/ ~/Downloads/final_data/

Abre el reporte de MultiQC en tu navegador:

open ~/Downloads/final_data/multiqc/multiqc_report.html

Tips de revisión en MultiQC:
- Per base sequence quality: define zonas de truncado (Q≥25–30).
- Adapter/Overrepresented/Primer content: confirma necesidad de recorte.
- Length distribution: válida rango del amplicón esperado.

4. Recorte de primers con Cutadapt

bash scripts/02b-cutadapt.sh

Notas del script

Corre siempre desde la raíz del proyecto: El script usa rutas relativas (raw_fastqs/, for_dada2/, etc.), por lo que debe ejecutarse desde ~/metabarcoding-code. Si lo corres desde otra carpeta, no encontrará los archivos.
Los directorios de salida se crean automáticamente: El script hace mkdir -p de for_dada2/ y cutadapt_reports/ antes de escribir, así que no necesitas crearlos manualmente.
Convención de nombres requerida: El script espera archivos con el patrón *_R1_001.fastq y *_R2_001.fastq en raw_fastqs/. Si tus archivos tienen otro formato de nombre, hay que ajustar el script.

Obtendrás

FASTQ recortados: for_dada2/
Reportes por muestra: cutadapt_reports/*_cutadapt.txt
Resumen: cutadapt_reports/overall_report.txt

¿Por qué hay que remover los primers?

Los primers son secuencias artificiales que se añaden durante la PCR para amplificar la región de interés. No forman parte del amplicón biológico y deben eliminarse antes de correr DADA2 porque:

DADA2 modela los errores de secuenciación base por base. Si los primers están presentes, introduce bases con un perfil de error diferente al resto de la lectura, lo que confunde el modelo.
Las secuencias de primers son idénticas en todas las muestras — si no se remueven, DADA2 las trata como variantes biológicas reales.
La longitud del amplicón real queda definida solo después del recorte.

¿Cómo funciona la búsqueda de primers?

En nuestras lecturas, cada read contiene el amplicón flanqueado por ambos primers. Como las lecturas (150 bp) son más largas que el inserto (~102 bp), el secuenciador lee más allá del amplicón y entra en el primer del otro extremo:

R1: [primer F 17bp]──amplicón 65bp──[primer R-rc 20bp]──adaptador...
R2: [primer R 20bp]──amplicón 65bp──[primer F-rc 17bp]──adaptador...

Cutadapt necesita recortar ambos primers de cada lectura — el del inicio (5’) y el del final (3’). Para esto usamos linked adapters, que buscan ambos primers simultáneamente y devuelven solo la secuencia entre ellos:

R1: -g “^FWD...RC_REV” — busca primer F al inicio y reverse-complement del primer R al final
R2: -G “^REV...RC_FWD” — busca primer R al inicio y reverse-complement del primer F al final

La sintaxis ^PRIMER1...PRIMER2 le dice a Cutadapt: “encuentra PRIMER1 anclado al inicio (^), encuentra PRIMER2 en cualquier posición después, y devuelve solo lo que hay entre ambos”.

Tasa de error (-e): Permite hasta un 10% de mismatches al buscar el primer. Esto es necesario porque los primers a veces tienen pequeñas variaciones entre muestras o errores de síntesis.

¿Qué pasa con lecturas donde no se encontró el primer?

Nuestro script usa --discard-untrimmed, que descarta las lecturas donde no se encontraron ambos primers. En metabarcoding esto es lo correcto: si una lectura no tiene los primers esperados, no corresponde al amplicón objetivo y es contaminación o un artefacto.

Lo que hace el script `02b-cutadapt.sh`

Define rutas, primers y parámetros configurables.
Verifica que cutadapt esté instalado y que existan los archivos de entrada.
Localiza todos los pares *_R1_001.fastq / *_R2_001.fastq en raw_fastqs/.
Para cada muestra: recorta primers, guarda FASTQs en for_dada2/ y un reporte individual en cutadapt_reports/.
Resume las estadísticas de todas las muestras en cutadapt_reports/overall_report.txt y overall_summary.tsv.

Visualiza el contenido del script:

cat scripts/02b-cutadapt.sh

Parámetros principales

Parámetro	Descripción	Valor por defecto
`-j`	Hilos en paralelo	4
`-g “^FWD...RC_REV”`	Linked adapter R1: primer F al inicio + primer R-rc al final	—
`-G “^REV...RC_FWD”`	Linked adapter R2: primer R al inicio + primer F-rc al final	—
`-e`	Tasa máxima de error al buscar el primer	0.10 (10%)
`-m`	Longitud mínima de la lectura tras el recorte	40 bp
`--discard-untrimmed`	Descarta lecturas sin ambos primers	Activado
`-o` / `-p`	Archivos de salida R1 y R2 recortados	—

Interpretación del reporte

Al terminar, revisa el resumen:

cat cutadapt_reports/overall_report.txt

El reporte por muestra tiene tres métricas clave:

Métrica	Qué indica	Valor esperado
`Read 1 with adapter`	% de lecturas R1 donde se encontró el primer forward	> 90%
`Read 2 with adapter`	% de lecturas R2 donde se encontró el primer reverse	> 90%
`Pairs written (passing filters)`	% de pares que pasaron el recorte y el filtro de longitud mínima (`-m`)	> 60%

Diferencia importante entre “reads with adapter” y “pairs written”:

”Reads with adapter” alto (~95%) significa que el primer se encontró correctamente — los primers están bien definidos.
”Pairs written” más bajo significa que después de recortar el primer, muchas lecturas quedaron por debajo de la longitud mínima (-m 40 por defecto) y fueron descartadas. Esto es normal en amplicones cortos donde el solapamiento entre R1 y R2 es grande.

Señales de alerta:

Patrón	Interpretación
`Read 1/2 with adapter` < 80%	El primer no se encuentra — verificar secuencia del primer en el script
`Read 1 with adapter` alto pero `Read 2` bajo (o viceversa)	Posible inversión de orientación en R1/R2
`Pairs written` < 30%	Muchas lecturas demasiado cortas tras el recorte — revisar `-m` o longitud del amplicón
`Read 1/2 with adapter` ≈ 0%	La muestra probablemente no contiene el amplicón objetivo — posible blanco, contaminación cruzada o muestra incorrecta

5. Revisión de resultados (datos recortados)

fastqc -q -t 4 -o final_data/fastqc_trimmed for_dada2/*.fastq
multiqc final_data/fastqc_trimmed -o final_data/multiqc_trimmed

Descarga los resultados actualizados (desde tu terminal local):

rsync -avz usrXX@200.23.162.240:/home/usrXX/metabarcoding-code/final_data/ ~/Downloads/final_data/
open ~/Downloads/final_data/multiqc_trimmed/multiqc_report.html

Objetivo

Pasos principales

1. Conteo de lecturas

Visualizar en terminal

Visualización opcional en R

2. Control de calidad con FastQC y MultiQC

Crear directorios y correr fastqc

Verificar conteo de reportes

Ejecutar MultiQC (resume todos los FastQC)

Verificar reporte MultiQC

3. Revisión de resultados (datos crudos)

Descargar resultados a tu computadora

4. Recorte de primers con Cutadapt

Obtendrás

¿Por qué hay que remover los primers?

¿Cómo funciona la búsqueda de primers?

¿Qué pasa con lecturas donde no se encontró el primer?

Lo que hace el script 02b-cutadapt.sh

Parámetros principales

Interpretación del reporte

5. Revisión de resultados (datos recortados)

Lo que hace el script `02b-cutadapt.sh`