¿Por qué familiarizarse con samtools?

En las dos entregas anteriores, completamos el alineamiento con BWA-MEM y STAR/HISAT2. Ahora debemos manipular los archivos BAM resultantes. GATK, IGV y scripts personalizados de Python — todos estos componentes dependen subyacentemente de los comandos de samtools. samtools es en sí mismo una herramienta de línea de comandos (CLI) y la interfaz pública de la biblioteca C htslib.

Dedicar solo una hora a aprenderlo facilitará enormemente cinco años de trabajo con pipelines bioinformáticos. No hay prisa.

SAM · BAM · CRAM — La relación entre los tres formatos

Estos tres formatos almacenan la misma información de maneras distintas.

La conversión es libre. Con samtools se puede ir en cualquiera de las tres direcciones.

Cinco comandos iniciales son suficientes

En la práctica diaria, solo se necesitan cinco comandos.

1. samtools view — Abramos un archivo

samtools view SAMPLE01.bam | head -3
samtools view -h SAMPLE01.bam | head -20   # incluir cabecera
samtools view -c SAMPLE01.bam              # número total de lecturas
samtools view -c -f 4 SAMPLE01.bam         # número de lecturas no alineadas
samtools view -c -F 256 SAMPLE01.bam       # número de lecturas excluyendo alineamientos secundarios

• -h: Imprime también la cabecera.
• -c: Cuenta únicamente el número de registros.
• -f N: Solo los reads con esta bandera activada (por ejemplo, 4 = unmapped).
• -F N: Excluye los reads con esta bandera activada (por ejemplo, 256 = secondary).

2. samtools sort — Etapa obligatoria de ordenación por coordenadas

samtools sort -@ 4 -o SAMPLE01.sorted.bam SAMPLE01.bam

BWA/STAR genera SAM/BAM en el orden original de las lecturas. GATK, IGV y bcftools requieren BAMs ordenados por coordenadas. sort realiza ese ordenamiento.

• -@ 4: 4 hilos.
• -o out.bam: salida.
• Si se requiere ordenamiento por nombre, use -n (caso específico: reintegración de pares, HTSeq).

3. samtools index — La clave para el acceso aleatorio

samtools index SAMPLE01.sorted.bam

Se crea el archivo .bai, y posteriormente IGV/bcftools/samtools view puede realizar consultas inmediatas en un segmento específico.

samtools view SAMPLE01.sorted.bam chr17:41196312-41277500   # solo la región de BRCA1

Sin índice, este comando escanea todo el archivo y puede tardar varios minutos; con él, termina en milisegundos. En los flujos de trabajo de Snakemake, es habitual que a sort le siga inmediatamente index.

4. samtools flagstat — Primer panel del flujo

samtools flagstat SAMPLE01.sorted.bam

Ejemplo de salida (resumen):

120,000,000 + 0 in total
115,200,000 + 0 mapped (96.00% : N/A)
119,800,000 + 0 paired in sequencing
114,900,000 + 0 properly paired (95.91% : N/A)
2,400,000 + 0 duplicates

Con solo revisar estas cinco líneas cada día, se obtiene un diagnóstico de salud del pipeline.

• mapped ratio: Si es inferior al 90%, sugiere discordancia con el genoma de referencia o contaminación.
• properly paired ratio: Si es inferior al 85%, indica problemas en el tamaño del inserto de la librería o residuos de adaptadores.
• duplicates: Si supera el 20%, implica una reducción en la complejidad de la librería.

5. samtools idxstats — Distribución de lecturas por cromosoma

samtools idxstats SAMPLE01.sorted.bam

Indica cuántas lecturas se han asignado por cada secuencia (cromosoma y contig). Con este resultado, se puede determinar el sexo (proporción de lecturas del cromosoma Y) o detectar anomalías mitocondriales de inmediato.

CRAM: ¿por qué la práctica profesional está migrando a este formato?

BAM almacena las secuencias para cada lectura individual. Sin embargo, la mayoría de las lecturas son casi idénticas al genoma de referencia. ¿No sería posible reducir significativamente el tamaño si solo se almacenan las diferencias con respecto a la referencia? Esta idea es la base de CRAM.

• CRAM codifica en binario únicamente las diferencias de cada lectura con respecto a la referencia.
• Un archivo BAM de secuenciación del genoma completo (WGS) humano de 100 GB se reduce a entre 30 y 50 GB al almacenarse como CRAM.
• Requiere la referencia: Para abrir un archivo CRAM, es necesario disponer del archivo FASTA de la referencia original.

Conversión BAM ↔ CRAM

# BAM → CRAM (es obligatorio indicar la referencia)
samtools view -T hg38.fa -C -o SAMPLE01.cram SAMPLE01.sorted.bam
samtools index SAMPLE01.cram   # crear .crai

# CRAM → BAM (conversión inversa)
samtools view -T hg38.fa -b -o SAMPLE01.bam SAMPLE01.cram

-T hg38.fa: Especificación de referencia. Debe ser la misma referencia utilizada durante el alineamiento. El resto del uso de CRAM es idéntico al de BAM.

Tres trampas prácticas de CRAM

1. Imposibilidad de reproducir sin la referencia: Si se conserva solo el archivo CRAM y se elimina la referencia, no será posible recuperar las lecturas originales. Es obligatorio archivar también la referencia.
2. Incompatibilidad de versiones de referencia: Si se alineó con hg19 pero se intenta abrir con hg38, las secuencias en cada posición no coincidirán, provocando errores. El estándar original indica registrar el SHA1 de la referencia en el encabezado del archivo CRAM.
3. Algunas herramientas antiguas no admiten CRAM: La mayoría de los nuevos pipelines sí lo admiten. Sin embargo, se debe seguir el principio de almacenar solo en CRAM cuando la fuente de la referencia sea clara.

Flujo práctico del pipeline — En una sola página

Se resume en una sola página el orden de uso de samtools desde BWA-MEM hasta la entrada en GATK.

# 1) Alineamiento y ordenación por coordenadas (streaming)
bwa-mem2 mem -t 16 -R "$RG" hg38.fa R1.fq.gz R2.fq.gz \
  | samtools sort -@ 4 -o S1.bam -

# 2) Índice
samtools index S1.bam

# 3) Panel del pipeline
samtools flagstat S1.bam > S1.flagstat.txt
samtools idxstats S1.bam > S1.idxstats.tsv

# 4) Conversión a CRAM para archivo
samtools view -T hg38.fa -C -o S1.cram S1.bam
samtools index S1.cram

# 5) Entrada al pipeline: pasar el BAM ordenado a MarkDuplicates (S07)

Incluso con 30 muestras, estos cinco pasos constituyen la estructura básica.

Cálculo manual — ¿Por qué se reduce tanto el tamaño de CRAM?

Si consideramos una cobertura WGS humana de 30x con 300 millones de lecturas, cada una de 150 pb, solo la secuencia codificada en BAM ocupa 45 GB. Sin embargo, ¿cuántas bases difieren de la referencia por lectura, en promedio?

• Polimorfismo del genoma humano individual: aproximadamente 4.000.000 de SNPs (0,13 % de los 3.000 millones de pb del genoma).
• En una lectura de 150 pb, estadísticamente hay 0,2 SNPs.
• Es decir, hay menos de un sitio diferente a la referencia por lectura.

CRAM almacena solo este "menos de uno", reduciendo la información de secuencia a menos del 1 % del original. El resto de la reducción (90 %) proviene de la codificación del encabezado y de la compresión con pérdida de las cadenas de calidad (opcional). La elegancia del diseño de CRAM radica en que permite ahorrar espacio de almacenamiento mientras se pueden restaurar las lecturas.

Movamos las manos en Colab

!apt-get install -y samtools >/dev/null
!wget -q https://sra-pub-src-1.s3.amazonaws.com/SRR1039508/SRR1039508_1.fastq.gz -O r1.fq.gz

# Familiarizarse con los comandos usando un BAM pequeño ya ordenado
!wget -q https://github.com/samtools/samtools/raw/master/examples/toy.bam
!samtools view -h toy.bam | head
!samtools flagstat toy.bam
!samtools index toy.bam
!samtools view toy.bam ref:100-500

Con este pequeño ejemplo, podrás familiarizarte con view, flagstat, index y consultas por intervalo en menos de 5 minutos.

Cinco comandos de una sola línea muy útiles

Estos son los comandos que más se consultan en los registros del pipeline.

# Extraer solo una región génica concreta a un BAM nuevo
samtools view -h -b in.bam chr17:41196312-41277500 > BRCA1.bam

# Extraer solo las lecturas no alineadas a FASTQ (para comprobar contaminación o especies cruzadas)
samtools view -b -f 4 in.bam | samtools fastq - > unmapped.fq

# Contar solo las lecturas marcadas como duplicadas
samtools view -c -f 1024 in.bam

# Cobertura por grupo de lecturas (resumen por muestra)
samtools view -H in.bam | grep '^@RG'

# Combinar dos BAM
samtools merge -@ 4 merged.bam S1.bam S2.bam S3.bam

Estos cinco comandos se utilizan todas las semanas.

Mapeo de CS

• Compresión basada en referencia: almacena solo las diferencias respecto a la referencia. Consulte el episodio de DryBench «Compresión basada en referencia» (reference-based-compression).
• BGZF (Block Gzip Format): BAM almacena bloques gzip separados, lo que permite el acceso aleatorio. El índice contiene los desplazamientos de esos bloques.
• CLI por streaming: conecta tuberías mediante - (stdin/stdout) de samtools, una aplicación directa de la filosofía Unix.

Conclusión

Con estos cinco comandos puede manejarse el 95 % del trabajo habitual con BAM. En el siguiente episodio (S06) veremos el panorama general de GATK4 Best Practices y en S07 MarkDuplicates recibirá el BAM ordenado que hemos preparado. Samtools seguirá apareciendo antes y después de numerosas etapas del pipeline GATK.

Si desea profundizar

• Danecek P. et al. (2021), Twelve years of SAMtools and BCFtools. GigaScience 10:giab008. Síntesis de doce años de desarrollo y decisiones de diseño de samtools.
• Hsi-Yang Fritz M. et al. (2011), Efficient storage of high throughput DNA sequencing data using reference-based compression. Genome Research 21:734. Artículo original de CRAM.
• Documentación oficial de htslib (https://www.htslib.org/): especificación primaria de la API y la CLI.
• EMBL-EBI Training — Working with BAM files, tutorial gratuito con numerosos casos de uso de cada comando.
• Broad Institute BroadE — Working with SAM files, vídeo con ejemplos reales de inspección de lecturas mediante samtools view.

Quien haya abierto un BAM y probado flagstat y view comprenderá con mucha más confianza las opciones del siguiente episodio sobre GATK. Dedique cinco minutos a practicarlas en Colab.