Análisis de Metagenoma usando SqueezeMeta

Autor/a

OMICs Analysis

Introducción

ADVERTENCIA ⚠️: En este tutorial vamos a echar a andar el flujo de trabajo SqueezeMeta, ten en cuenta que para este punto debes tener suficiente espacio de memoria (470Gb) en tu computadora para generar la base de datos que utiliza esta herramienta. También ten en cuenta tener mínimo 8Gb de RAM.

Los siguiente análisis se realizaron en un equipo Gigabyte Technology Co., Ltd. H67MA-USB3-B3 procesador Intel® Core™ i7-2600 CPU @ 3.40GHz × 8 con 32Gb de RAM. Con un sistema operativo Ubuntu 22.04.5 LTS.

Para este tutorial vamos a trabajar con datos del estudio “Effects of heat stress on 16S rDNA, metagenome and metabolome in Holstein cows at different growth stages”. Para fines educativos solo usaremos un subconjunto de los datos de las vacas Holstein en tres etapas de crecimiento, en condiciones normales y con estrés térmico.

SRA ID Condición Etapa de crecimiento
SRR19746212 Periodo Normal Novillas en crecimiento
SRR19746215 Estrés Térmico Novillas en crecimiento
SRR19746202 Periodo Normal Novillas
SRR19746205 Estrés Térmico Novillas
SRR19746208 Periodo Normal Vacas Lactantes
SRR19746219 Estrés Térmico Vacas Lactantes

Preparación de datos

Para descargar los datos usaremos RSA Toolkit y parallel-fastq-dump, a continuación el ejemplo para uno. Por otra parte vamos a generar nuestro directorio de trabajo.

mkdir 2025_Demo_sqm
cd 2025_Demo_sqm
mkdir data
cd data

Este mismo proceso se realizó para todos los SRA ID.

prefetch SRR19746212

parallel-fastq-dump \
--sra-id SRR19746212 \
-t 7 \
--outdir ./ \
--split-files \
--gzip

Pre procesamiento y control de calidad

Trimmomatic es una herramienta flexible y eficiente diseñada específicamente para el preprocesamiento de datos de secuenciación de próxima generación (NGS) generados por plataformas Illumina. Se utiliza principalmente para realizar tareas de recorte (trimming) y filtrado de lecturas en formato FASTQ, tanto en datos paired-end como single-end. Sus funciones clave incluyen:

  • Remoción de adaptadores: Elimina secuencias de adaptadores Illumina que pueden contaminar las lecturas, lo que es esencial para evitar problemas en alineamientos posteriores.

  • Recorte de bases de baja calidad: Corta bases al inicio o final de las lecturas que no cumplen con un umbral de calidad (por ejemplo, usando Phred scores), mejorando la calidad general de los datos.

  • Filtrado de lecturas: Descarta lecturas cortas, de baja calidad o con patrones no deseados, como secuencias de baja complejidad.

  • Otras operaciones: Puede realizar cropping (recorte fijo de longitud) y sliding window para evaluar calidad en ventanas móviles.

Esto es particularmente útil en flujos de trabajo de metagenómica, transcriptómica o genómica con datos Illumina, ya que ayuda a limpiar los datos crudos antes de análisis downstream como ensamblaje o mapeo, reduciendo errores y mejorando la precisión.

Cutadapt, por su parte, es una herramienta versátil enfocada en la remoción de adaptadores y secuencias no deseadas de lecturas de secuenciación de alto rendimiento, con un énfasis en datos Illumina. Sus usos principales en secuencias Illumina son:

  • Remoción de adaptadores: Busca y elimina adaptadores (como los universales de Illumina) en cualquier posición de las lecturas, ya sea al 3’ o 5’, o incluso en el medio, lo que es común en datos paired-end o single-end contaminados.

  • Recorte basado en calidad: Puede trimmed bases de baja calidad desde los extremos de las lecturas, similar a Trimmomatic, pero con opciones más flexibles para filtrado.

  • Filtrado y modificación de lecturas: Permite descartar lecturas demasiado cortas, con mismatches excesivos o que no cumplan criterios específicos, y soporta operaciones como demultiplexing (separación por barcodes).

  • Otras funcionalidades: Maneja poly-A tails, secuencias repetitivas o primers, y es eficiente para grandes volúmenes de datos.

Cutadapt

Cutadapt es util para limpiar lecturas paired-end de Illumina en archivos FASTQ:

  • -a y -A: Elimina adaptadores específicos de Illumina (forward: GATCGGAAGAGCACACGTCTGAACTCCAGTCA; reverse: AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT).

  • --trim-n: Recorta bases “N” (indeterminadas) de los extremos.

  • --match-read-wildcards: Permite comodines en las lecturas para coincidencias flexibles.

conda activate cutadapt

mkdir 1_QC
cd 2025_Demo_sqm
cutadapt \
-a GATCGGAAGAGCACACGTCTGAACTCCAGTCA \
-A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT \
--trim-n \
--match-read-wildcards -j 7 \
-o "1_QC/SRR19746202_1.fastq.gz" \
-p "1_QC/SRR19746202_2.fastq.gz" \
data/SRR19746202_1.fastq.gz \
data/SRR19746202_2.fastq.gz > SRR19746202_cutadapt.log 2>&1

Trimmomatic

El siguiente código lo que hace aplicar filtros:

  • LEADING:15: Elimina bases iniciales con calidad <15.

  • TRAILING:15: Elimina bases finales con calidad <15.

  • SLIDINGWINDOW:10:15: Recorta si el promedio de calidad en ventana de 10 bases es <15.

  • AVGQUAL:25: Descarta lecturas con calidad promedio <25.

  • MINLEN:50: Descarta lecturas <50 bases.

java -jar bin/trimmomatic.jar PE \
-threads 7 \
-phred33 \
"1_QC/SRR19746202_1.fastq.gz" "1_QC/SRR19746202_2.fastq.gz"\
"1_QC/SRR19746202_1_trimmed.fastq.gz" "SRR19746202_1_discard.fastq" \
"1_QC/SRR19746202_2_trimmed.fastq.gz" "SRR19746202_2_discard.fastq" \
LEADING:15 TRAILING:15 SLIDINGWINDOW:10:15 AVGQUAL:25 MINLEN:50 > "SRR19746202_trimmomatic.log" 2>&1

Quitar contaminación del Hospedero

sudo apt install bowtie2

Creación de base de datos del hospedero (Bos taurus)

Descarga de su genoma

cd data

wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/002/263/795/GCF_002263795.3_ARS-UCD2.0/GCF_002263795.3_ARS-UCD2.0_genomic.fna.gz
mv GCF_002263795.3_ARS-UCD2.0_genomic.fna.gz bos_taurus_host.fna.gz
mkdir db_host

Indexación para creación de base de datos

bowtie2-build --seed 123 bos_taurus_host.fna.gz db_host

Alineamiento contra hospedero

Alineamiento de fastq contra la base de datos

bowtie2 \
--very-sensitive-local \
-p 7 --seed 123 \
-x ./data/db_host/db_host \
-1 data/SRR19746202_1.fastq.gz \
-2 data/SRR19746202_2.fastq.gz \
-S 1_QC/host_filter/SRR19746202.sam

⚠️El archivo .sam es muy pesado, hasta 20Gb, es indispensable comprimirlo y eliminarlo si deseas.

samtools view \
-@ 7 \
-b -f 12 -F 256 \
1_QC/host_filter/SRR19746202.sam > 1_QC/host_filter/SRR19746202.bam
samtools sort \
-n 1_QC/host_filter/SRR19746202.bam \
-o 1_QC/host_filter/SRR19746202.sorted.bam

Creación de archivo FASTQ interleaved (pareado)

samtools bam2fq \
1_QC/host_filter/SRR19746202.sorted.bam > 1_QC/host_filter/SRR19746202.fastq

Automatización

Para hacerlo con todas las muestras automatizamos el proceso, y de paso comprimir el archivo .fastq, creamos un script llamado host_filter.sh:

host_filter.sh 
#!/usr/bin/env bash
set -euo pipefail

# host_filter.sh
# Pipeline simple: bowtie2 -> samtools (filtrar pares no mapeados) -> convertir a FASTQ comprimido (.fastq.gz)
# Uso: ./host_filter.sh                # usa la lista por defecto
#       ./host_filter.sh 19746205 19746208   # procesa muestras dadas (acepta con o sin prefijo SRR)

samples=(19746202 19746205 19746208 19746212 19746215 19746219)
if [ "$#" -gt 0 ]; then
  samples=("$@")
fi

INDEX="./data/db_host/db_host"
OUTDIR="1_QC/host_filter"
THREADS=7
SEED=123

mkdir -p "$OUTDIR"

command -v bowtie2 >/dev/null 2>&1 || { echo "ERROR: bowtie2 no está en PATH" >&2; exit 1; }
command -v samtools >/dev/null 2>&1 || { echo "ERROR: samtools no está en PATH" >&2; exit 1; }

for s in "${samples[@]}"; do
  SAMPLE="$s"
  if [[ "$SAMPLE" != SRR* ]]; then
    SAMPLE="SRR${SAMPLE}"
  fi

  R1="data/${SAMPLE}_1.fastq.gz"
  R2="data/${SAMPLE}_2.fastq.gz"
  SAM="$OUTDIR/${SAMPLE}.sam"
  LOG="$OUTDIR/${SAMPLE}.log"
  BAM="$OUTDIR/${SAMPLE}.bam"
  SORTED="$OUTDIR/${SAMPLE}.sorted.bam"
  FQ_GZ="$OUTDIR/${SAMPLE}.fastq.gz"

  echo "=== Procesando ${SAMPLE} ===" | tee -a "$LOG"

  if [ ! -f "$R1" ]; then
    echo "ERROR: $R1 no existe. Saltando ${SAMPLE}." | tee -a "$LOG"
    continue
  fi
  if [ ! -f "$R2" ]; then
    echo "ERROR: $R2 no existe. Saltando ${SAMPLE}." | tee -a "$LOG"
    continue
  fi

  # 1) Bowtie2 -> SAM (stdout+stderr al log)
  echo "[1] bowtie2 -> SAM" | tee -a "$LOG"
  bowtie2 --very-sensitive-local -p "$THREADS" --seed "$SEED" \
    -x "$INDEX" -1 "$R1" -2 "$R2" -S "$SAM" >>"$LOG" 2>&1

  # 2) Filtrar lecturas donde ambos pares NO se alinearon al host (-f 12) y convertir a BAM
  echo "[2] samtools view -b -f 12 -F 256 -> BAM" | tee -a "$LOG"
  samtools view -@ "$THREADS" -b -f 12 -F 256 "$SAM" -o "$BAM" 2>>"$LOG"

  # eliminar SAM para ahorrar espacio
  rm -f "$SAM"

  # 3) Ordenar por nombre
  echo "[3] samtools sort -n -> $SORTED" | tee -a "$LOG"
  samtools sort -n -@ "$THREADS" -o "$SORTED" "$BAM" 2>>"$LOG"
  rm -f "$BAM"

  # 4) Convertir a FASTQ comprimido (.fastq.gz)
  echo "[4] samtools bam2fq -> gzip -> $FQ_GZ" | tee -a "$LOG"
  samtools bam2fq "$SORTED" | gzip -c > "$FQ_GZ" 2>>"$LOG"

  # limpiar intermedios
  rm -f "$SORTED"

  echo "=== Hecho ${SAMPLE}. Salida: $FQ_GZ . Log: $LOG ===" | tee -a "$LOG"
  echo
done

echo "Pipeline finalizado. Archivos en: $OUTDIR"
chmod +x host_filter.sh

./host_filter

⚠️PEROOO, necesitamos las lecturas pareada para SqueezeMeta, es por ello que utilizamos el siguiente código de Python para separarlas. A este archivo le llamamos sort_fastq.py:

sort_fastq.py: 
#!/usr/bin/env python3
import gzip
import sys

if len(sys.argv) != 4:
    print("Uso: split_interleaved_fastq.py IN.fastq.gz OUT_1.fastq.gz OUT_2.fastq.gz", file=sys.stderr)
    sys.exit(1)

inp, out1, out2 = sys.argv[1], sys.argv[2], sys.argv[3]

with gzip.open(inp, "rt") as inf, gzip.open(out1, "wt") as o1, gzip.open(out2, "wt") as o2:
    buf = []
    for i, line in enumerate(inf, 1):
        buf.append(line)
        if i % 8 == 0:
            # primer par (4 líneas) -> out1; segundo par (4 líneas) -> out2
            o1.writelines(buf[:4])
            o2.writelines(buf[4:])
            buf = []
    # por si queda algo al final (archivo truncado o similar)
    if buf:
        if len(buf) <= 4:
            o1.writelines(buf)
        else:
            o1.writelines(buf[:4])
            o2.writelines(buf[4:])
chmod +x sort_fastq.py

python3 sort_fastq.py \
1_QC/host_filter/SRR19746202.fastq.gz \
1_QC/host_filter/SRR19746202_1.fastq.gz \
1_QC/host_filter/SRR19746202_2.fastq.gz

✅Extra. Con este código puedes ver en cada paso la cantidad de lecturas que hay en cada archivo FASTQ:

echo "R1 reads: $(( $(zcat 1_QC/host_filter/SRR19746202_1.fastq.gz | wc -l) / 4 ))" 
echo "R2 reads: $(( $(zcat 1_QC/host_filter/SRR19746202_2.fastq.gz | wc -l) / 4 ))"

SqueezeMeta

Para iniciar con el flujo de trabajo de SqueezeMeta necesitagenerar un archivo de table-separate-value con el nombre de los archivos, lo nombraremos samples.samples

samples.samples: 
SRR19746202 SRR19746202_1.fastq.gz  pair1
SRR19746202 SRR19746202_2.fastq.gz  pair2
SRR19746205 SRR19746205_1.fastq.gz  pair1
SRR19746205 SRR19746205_2.fastq.gz  pair2
SRR19746208 SRR19746208_1.fastq.gz  pair1
SRR19746208 SRR19746208_2.fastq.gz  pair2
SRR19746212 SRR19746212_1.fastq.gz  pair1
SRR19746212 SRR19746212_2.fastq.gz  pair2
SRR19746215 SRR19746215_1.fastq.gz  pair1
SRR19746215 SRR19746215_2.fastq.gz  pair2
SRR19746219 SRR19746219_1.fastq.gz  pair1
SRR19746219 SRR19746219_2.fastq.gz  pair2
conda activate SqueezeMeta

SqueezeMeta.pl \
-m sequential \
-s 1_QC/host_filter/samples.samples \
-f 1_QC/host_filter \
-t 7 \
--lowmem

conda deactivate