¿Por qué sigue siendo BWA-MEM?

Han pasado más de 12 años desde que Heng Li publicara BWA-MEM en 2013. Desde entonces, han surgido Minimap2, Bowtie2, HISAT2 y STAR, mientras que DeepVariant incorpora directamente aprendizaje profundo. Sin embargo, las tuberías de genómica clínica, las mejores prácticas de GATK4, gnomAD, UK Biobank y All of Us siguen utilizando BWA-MEM para el alineamiento. ¿Por qué?

La razón es sencilla: BWA-MEM alinea correctamente todas las lecturas que pueden alinearse bien y admite honestamente aquellas que no pueden alinearse. Esta honestidad satisface los supuestos estadísticos de las etapas posteriores de la tubería. En esta entrada, exploraremos de dónde proviene esta estabilidad gracias a sus decisiones de diseño y qué aspectos debemos tener en cuenta en la práctica.

Conexión con episodios anteriores — Reutilización del nivel Micro

Dentro de BWA-MEM, se combinan los conceptos que aprendimos en el nivel M:

• M13～M15: Indexación previa del genoma de referencia mediante la construcción de una matriz sufijada (SA) → BWT → índice FM.
• M10～M11: División de las lecturas en semillas k-mer y búsqueda rápida de posiciones candidatas utilizando el índice FM (principio de BLAST).
• M06～M08: Alineamiento local con penalización de gap afín alrededor de las posiciones candidatas (variante de Smith-Waterman).

Es decir, BWA-MEM es una implementación práctica que ensambla los tres algoritmos que ya hemos derivado. Adoptar esta perspectiva facilita enormemente la comprensión de las opciones disponibles.

Indexación — realizar una vez y reutilizar

Se indexa el genoma de referencia (3.1 GB para hg38 en formato fa) con el siguiente comando:

bwa index -a bwtsw hg38.fa

Esta única operación genera archivos de índice de 5 a 6 GB (.bwt, .pac, .ann, .amb, .sa). Para el genoma humano, esto toma aproximadamente 30 a 60 minutos. Una vez creados, se reutilizan en todos los proyectos, por lo que es habitual almacenarlos en una carpeta compartida del servidor del equipo.

Si utiliza BWA-MEM2, se requiere un índice separado.

bwa-mem2 index hg38.fa

BWA-MEM2 produce exactamente los mismos resultados que el BWA-MEM original, pero es de 1,7 a 3 veces más rápido gracias al uso de instrucciones vectoriales (AVX-512). La mayoría de las bibliotecas en la nube recientes priorizan BWA-MEM2.

Ejecución básica — Una muestra

Este es el comando práctico más sencillo.

bwa-mem2 mem -t 16 \
  -R "@RG\tID:S1_lane1\tSM:SAMPLE01\tLB:lib1\tPL:ILLUMINA\tPU:unit1" \
  hg38.fa clean_R1.fq.gz clean_R2.fq.gz \
  | samtools sort -@ 4 -o SAMPLE01.bam -
samtools index SAMPLE01.bam

Dividámoslo en tres partes.

• bwa-mem2 mem: El propio ordenamiento. -t es el número de hilos, y -R es la cabecera del grupo de lectura (read group).
• samtools sort: Convierte el flujo SAM en un archivo BAM ordenado por coordenadas.
• samtools index: Genera el índice .bai (requerido por IGV y GATK).

La práctica habitual es conectar estos pasos con tuberías (|) para evitar escribir SAM en disco. Mantener un archivo SAM de 100 GB en el disco es un desperdicio de almacenamiento y de E/S.

Grupo de lectura — la etiqueta del flujo de trabajo

¿Qué es esa cadena de la opción -R? La etiqueta @RG es metadatos que indican en qué experimento provino cada lectura.

Importante: GATK MarkDuplicates determina las duplicaciones a nivel de biblioteca usando LB. HaplotypeCaller llama variantes por muestra utilizando SM. Si se asignan los grupos de lectura de forma descuidada, el resto del flujo de trabajo producirá silenciosamente resultados incorrectos. La mitad de los incidentes en flujos de trabajo clínicos provienen de errores en los grupos de lectura.

Etiquetas SAM/BAM — ¿dónde verificar qué?

Abramos una línea del resultado del ordenamiento.

SRR..._1  99  chr1  100501  60  150M  =  100601  250  ACGT...  IIII...  NM:i:2  MD:Z:75A20T53  AS:i:145  XS:i:87

Orden de los campos (de izquierda a derecha): QNAME · FLAG · RNAME · POS · MAPQ · CIGAR · RNEXT · PNEXT · TLEN · SEQ · QUAL · Etiquetas.

En la práctica, los campos que se consultan constantemente son tres.

1. MAPQ (calidad de mapeo)

• Entero 0～60. La lógica de -10 log10(probabilidad de posición incorrecta) es similar a Phred.
• Protocolo BWA-MEM: 60 = alineación única óptima, 0 = alineado con la misma calidad en múltiples posiciones (multi-mapper).
• GATK BQSR/HaplotypeCaller utiliza por defecto solo MAPQ de 20 o superior. Las lecturas con MAPQ bajo son descartadas silenciosamente en las etapas posteriores.

2. CIGAR

• Representación en cadena de texto de cómo se alineó la lectura con la referencia.
• M(coincidencia/discordancia común), I(inserción en la lectura), D(deleción en la lectura), S(soft-clip, permanece en la lectura pero se excluye del alineamiento), H(hard-clip, eliminado de la lectura).
• Ejemplo: 100M50S significa que solo los primeros 100 bp están alineados y los 50 bp posteriores fueron recortados por restos de adaptadores.

3. FLAG

• Máscara de bits de 12 bits. Contiene en un solo entero si la lectura es pareada, si es la primera, si no está mapeada, si es un duplicado, etc.
• Valores frecuentes: 4=unmapped, 1024=duplicado de PCR, 256=secundario, 2048=suplementario.
• Consejo: En broadinstitute.github.io/picard/explain-flags.html puede introducir el entero para obtener su significado.

Si puede leer estos tres campos con fluidez, habrá resuelto el 80% de la depuración de pipelines.

Alineamiento suplementario / quimérico — La trampa de la era de las lecturas largas

Entre los resultados especiales de BWA-MEM, lo que suele causar errores frecuentes en la práctica es la bandera 2048 (alineamiento suplementario).

• Una sola lectura se divide y se alinea en dos lugares distintos de la referencia. La parte frontal en chr1 y la parte posterior en chr7.
• BWA determina esto como un "alineamiento quimérico" y genera dos líneas para la lectura. Una es la primaria y la otra es la suplementaria (bandera 2048).
• Si se cuentan estas dos líneas como lecturas independientes, la cobertura aparecerá duplicada o se perderán señales de variantes estructurales (genes de fusión).

Reglas prácticas: antes de calcular la cobertura, filtre solo las lecturas primarias (samtools view -F 2304, excluyendo secondary+supplementary). Por el contrario, la detección de SV depende fundamentalmente de estas lecturas supplementary. Recuerde que el filtro debe adaptarse al uso.

Solo hay que conocer realmente unas pocas opciones

El manual de BWA-MEM enumera más de 30 opciones, pero en la práctica solo se utilizan tres.

• -t N: Hilos. Igual al número de núcleos de CPU.
• -M: Marcar lecturas divididas como secondary (compatibilidad con versiones antiguas de GATK). Las versiones más recientes de GATK4 no lo requieren, pero debe mantenerse si el equipo utiliza pipelines antiguos.
• -K 100000000: Tamaño de lote fijo. Garantiza que los resultados sean idénticos al cambiar el número de hilos (reproducibilidad).

La precisión del alineamiento ya está extremadamente bien ajustada en los valores predeterminados, por lo que no es necesario modificarla.

Orquestación de múltiples muestras — Un solo script Bash

¿Cómo ejecutar 30 muestras? Antes de aprender Snakemake o Nextflow, se comienza con este bucle en Bash.

#!/bin/bash
set -euo pipefail
REF=/data/ref/hg38.fa
OUT=/data/aligned

for R1 in fastq/*_R1.fq.gz; do
    SAMPLE=$(basename "$R1" _R1.fq.gz)
    R2=fastq/${SAMPLE}_R2.fq.gz
    RG="@RG\tID:${SAMPLE}\tSM:${SAMPLE}\tLB:${SAMPLE}_lib\tPL:ILLUMINA"

    echo "[$(date +%T)] Aligning $SAMPLE"
    bwa-mem2 mem -t 16 -K 100000000 -R "$RG" "$REF" "$R1" "$R2" \
        | samtools sort -@ 4 -o "$OUT/${SAMPLE}.bam" -
    samtools index "$OUT/${SAMPLE}.bam"
done

Tres cosas que aprenderás de este script.

1. set -euo pipefail: Si uno falla, se detiene inmediatamente. Las tuberías (pipelines) desconfían más que nada del fallo silencioso.
2. Extracción automática del nombre de la muestra desde el nombre del archivo.
3. Configuración automática del grupo de lectura por muestra.

El siguiente paso es migrar desde este esqueleto a GNU Parallel o Slurm.

Cálculo manual — La intuición detrás de MAPQ

Hagamos una idea aproximada de cómo BWA-MEM asigna MAPQ. Si los puntajes de alineación de dos posiciones candidatas son AS=145 y XS=87 (ejemplo real de etiqueta SAM), la diferencia entre el mejor y el segundo mejor es de 58 puntos. BWA convierte esta diferencia en una probabilidad y asigna aproximadamente un MAPQ de 60. Por el contrario, si es AS=145, XS=143, las dos posiciones están casi empatadas, por lo que MAPQ cae drásticamente a unos 3～5. Esto ocurre en regiones repetitivas o sitios de genes similares.

• Implicación: Una lectura con MAPQ 60 está estadísticamente asignada con confianza.
• Implicación: MAPQ 0 significa "se alinea igual de bien en varios lugares, por lo que no se puede determinar dónde está". Es útil como información, pero generalmente se excluye en la llamada de variantes.

Mapeo CS — Búsqueda de semillas + DP local

• Búsqueda de índice FM: La parte donde se encuentran rápidamente los k-mer semilla de la lectura en la referencia. Es exactamente igual a la sección "FM-index básico" de DryBench (fm-index-fundamentals).
• DP local: La etapa de expansión de semillas. Ajustes industriales de Smith-Waterman (M07) y gap afín (M08).
• Procesamiento por lotes: -K 100000000 fija el tamaño del lote para la reproducibilidad. Es una convención de CS para mantener la unicidad de los resultados durante el procesamiento en paralelo.

Conclusión

BWA-MEM ha sido el estándar durante 12 años. Hoy hemos confirmado que su razón es la etiquetado estable, MAPQ y manejo quimérico más que el refinamiento del algoritmo. En el próximo episodio (S04), veremos qué otras decisiones toman STAR/HISAT2 para lecturas con empalmes, como en RNA-seq. Después, en S05 organizaremos BAM/CRAM con samtools y entraremos en la tubería GATK4 desde S06.

Si quieres profundizar más

• Li H. (2013), Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv:1303.3997. El artículo original de BWA-MEM — los ajustes prácticos de seed-and-extend se detallan aquí.
• Vasimuddin M. et al. (2019), Efficient architecture-aware acceleration of BWA-MEM for multicore systems. Artículo sobre la vectorización SIMD de BWA-MEM.
• Broad Institute BroadE YouTube — Read Alignment with BWA Clase de 30 minutos. Subtítulos completos. Se presentan casos de errores en los grupos de lecturas (read groups) desde la perspectiva del pipeline GATK4.
• Blog personal de Heng Li: https://lh3.github.io/. Se publican frecuentemente anécdotas sobre la decisión de las opciones de BWA-MEM por parte del autor.
• Documentación oficial de SAMtools (https://www.htslib.org/doc/): Fuente original de las especificaciones SAM/BAM/CRAM.

Solo al intentar alinear una muestra manualmente se obtiene una verdadera comprensión de las etiquetas. Si selecciona un FASTQ pequeño en Galaxy EU y realiza el proceso desde BWA-MEM → samtools sort → visualización en IGV, el siguiente episodio le resultará mucho más sencillo.