Mapping avec BWA-MEM

Prérequis

Cette section prolonge le contrôle qualité de la section 3.1 : vous devez savoir lire un rapport FastQC et interpréter des lectures filtrées. Les formats FASTA et FASTQ sont supposés connus depuis 1.2, ainsi que le principe d'un score d'alignement vu en 2.2 et 2.4. Les commandes sont données pour Bash ; sous Windows, utilisez WSL ou un conteneur, car les pipelines reproducibles du cours supposent un shell POSIX.

Introduction

Le séquençage d'un génome ne produit pas un génome : il produit des centaines de millions de fragments courts, lus indépendamment, dont on ignore la position. Le mapping, ou alignement sur une référence, consiste à retrouver pour chaque lecture l'origine la plus probable sur un génome connu. Cette étape transforme un tas de séquences en coordonnées. Elle conditionne tout ce qui suit : appel de variants, génotypage, comptage d'expression, détection de fusions, études d'association. Une lecture mal placée ne provoque pas d'erreur visible ; elle déplace simplement une observation, parfois vers un paralogue, parfois vers une région répétée.

Pour Escherichia coli K-12 MG1655, la référence est le chromosome de NC_000913.3, long d'environ 4,6 Mb. Une expérience à 30× avec des lectures de 150 pb contient environ $30 \times 4{,}6\times10^6 / 150 \approx 9{,}2\times10^5$ lectures. Chercher chaque lecture naïvement dans la référence coûterait de l'ordre de $10^{15}$ comparaisons de caractères. Un index ramène la recherche d'un seed à quelques dizaines d'opérations, mais introduit ses propres choix : longueur minimale des seeds, score de match, pénalité de gap, gestion des extrémités.

Le défi n'est donc pas seulement algorithmique. Il faut aussi conserver la provenance des lectures, distinguer un alignement unique d'un alignement ambigu, apparier correctement les deux fichiers d'une paire, et documenter la référence et les versions. Un BAM sans read group ni version d'index n'est pas reproductible, même si ses alignements sont excellents. Nous allons construire cette chaîne avec BWA-MEM, puis lire ses sorties avec SAMtools.

Théorie

L'index : transformer la référence en structure interrogeable

Une référence est un fichier FASTA, mais un aligneur ne le parcourt pas linéairement à chaque lecture. bwa index écrit cinq fichiers (.amb, .ann, .bwt, .pac, .sa) à partir d'un préfixe ; samtools faidx ajoute un .fai pour les accès aléatoires. Le .bwt contient la transformée de Burrows-Wheeler et .sa le tableau des suffixes, qui forment un index FM (Ferragina-Manzini) : chercher un motif de longueur $k$ coûte $O(k)$ étapes, et la construction est linéaire en la taille $L$ de la référence, avec une constante élevée.

L'index est indissociable du FASTA qui l'a produit. Modifier une base, renommer un contig ou changer de version déplace les coordonnées et invalide l'index. BWA retrouve les fichiers par leur préfixe ; un index BWA-MEM2 n'est pas lisible par BWA-MEM, et réciproquement. Notez la somme de contrôle du FASTA : c'est la seule preuve que l'index et la référence sont les mêmes.

Seed-and-extend : trouver, chaîner, étendre

BWA-MEM appartient aux aligneurs seed-and-extend. Une correspondance exacte maximale (maximal exact match, MEM) est un intervalle $[i,j]$ de la lecture et $[p,p+\ell)$ de la référence tels que $\ell=j-i+1$, $r_{i..j}=R_{p..p+\ell-1}$, non extensible ni à gauche ni à droite. Le programme extrait des MEM, privilégie les plus longs (SMEM), les place, puis regroupe les seeds compatibles en chaînes colinéaires : même contig, même orientation, ordre et espacement compatibles. Chaque chaîne est étendue par un alignement local à gaps affines, avec $A=1$, $B=4$, $O=6$, $E=1$. Le coût d'un gap de longueur $k$ est $g(k)=O+kE$.

La longueur minimale d'un seed est -k, 19 par défaut. Une lecture plus courte que 19 pb ne fournit aucun seed et reste non alignée, même identique à la référence. Une lecture de 40 pb avec une erreur au centre laisse deux MEM de 19 et 20 pb ; une lecture de 30 pb avec une erreur au milieu ne laisse que 14 et 15 pb, sous le seuil.

Prenons une référence de 40 pb, ACGTTACGTACGTGCAACGTTAGCTAGCTAGGCTAACGTA, et trois lectures de 40 pb. La lecture exacte donne 40M et un score de 40. Une substitution en position 5 donne 39 matches et un mismatch : $39-4=35$, avec un MEM de 35 pb. Deux substitutions, en positions 5 et 25, donnent $38-8=30$ et des MEM de 4, 19 et 15 pb : un seul seed atteint le seuil. La MAPQ ne se déduit pas de ce calcul.

Lectures appariées, insert size et flags

La plupart des protocoles Illumina produisent des paires : un fragment est lu à ses deux extrémités, R1 dans un sens et R2 dans l'autre. Les deux fichiers doivent être fournis dans l'ordre, R1 puis R2, et se correspondre position par position. BWA-MEM ne réapparie pas les noms : la $i$-ème lecture de R1 et la $i$-ème de R2 forment une paire. Des fichiers désynchronisés produisent des paires fausses, même si chaque lecture s'aligne.

Pour une paire en orientation FR (forward-reverse), l'insert size se définit comme la distance entre le début de la lecture la plus à gauche et la fin de la lecture la plus à droite :

Avec une lecture gauche de 150 pb démarrant en position 100 et une lecture droite de 150 pb démarrant en position 300, $I = 300 + 150 - 100 = 350$ pb. Le champ TLEN vaut $+350$ pour la lecture la plus à gauche et $-350$ pour l'autre. Un insert aberrant, une orientation inattendue ou un mate sur un autre contig produit une paire discordante ; BWA-MEM peut tenter un sauvetage (mate rescue) autour de la position attendue.

Les flags SAM portent cette information. 0x1 indique une lecture appariée, 0x2 une paire correctement appariée, 0x4 une lecture non alignée, 0x10 le brin reverse, 0x20 un mate reverse, 0x40 la première lecture et 0x80 la seconde. 99 = 64 + 32 + 2 + 1 et 147 = 128 + 16 + 2 + 1 forment la signature classique d'une paire FR correcte.

Read groups : donner une identité aux lectures

Un read group est un ensemble de métadonnées attaché aux lectures d'une même unité technique. Dans l'en-tête SAM, il s'écrit @RG suivi de champs séparés par des tabulations. Cinq champs structurent le pipeline :

SM relie plusieurs read groups au même échantillon ; LB permet le marquage des duplicats par bibliothèque ; PU distingue les flowcells, lanes et barcodes ; PL documente la technologie et influence certains modèles d'erreur. Sans ces champs, un variant caller ne peut pas distinguer deux échantillons séquencés dans le même run, et le marquage des duplicats mélange des bibliothèques. Le tag RG:Z:<ID> est répété sur chaque enregistrement aligné.

BWA-MEM accepte le read group avec -R. La chaîne contient des séquences \t que BWA convertit en tabulations ; dans le shell, on protège l'argument par des guillemets simples :

bwa mem -R '@RG\tID:lane1\tSM:ecoli_ctrl\tPL:ILLUMINA\tLB:lib1\tPU:fc1.lane1' ref.fa R1.fq.gz R2.fq.gz

Déclarez un read group par lane et par bibliothèque, de préférence au mapping ; samtools addreplacerg permet de les ajouter après coup. Vérifiez l'en-tête avec samtools view -H et les premiers enregistrements avec samtools view : un @RG présent dans l'en-tête mais absent des lectures signale un problème de pipeline.

MAPQ, alignements secondaires et supplémentaires

La qualité de mapping (mapping quality, MAPQ) est un score Phred qui résume la confiance dans la position choisie :

MAPQ 60 correspond à $P=10^{-6}$, MAPQ 30 à $P=10^{-3}$, MAPQ 0 à une confiance nulle. BWA-MEM plafonne à 60 et combine score du meilleur alignement, score des alternatives et abondance des seeds ; ce n'est pas une probabilité calibrée. Une MAPQ de 60 ne garantit pas la justesse biologique, seulement l'absence de concurrent sérieux.

Un read peut avoir un alignement primaire, des secondaires (0x100, positions alternatives du même read) et des supplémentaires (0x800, parties d'un read chimérique, reliées par SA). Ce ne sont pas des duplicats : un duplicat est une molécule indépendante de la même bibliothèque, marquée plus tard par samtools markdup.

Pour compter les fragments, filtrez les secondaires et les supplémentaires avec samtools view -F 0x900 ; -q 20 ne les supprime pas. L'option -M de BWA-MEM marque les petites parties d'un split alignement comme secondaires plutôt que supplémentaires, pour la compatibilité avec certains outils historiques.

Choisir le génome de référence

Pour le fil rouge, la référence est E. coli K-12 MG1655, accession NC_000913.3. L'accession inclut un numéro de version : deux versions successives de la même accession ne sont pas interchangeables, car des corrections de séquence déplacent les coordonnées. Téléchargez le FASTA depuis RefSeq, vérifiez l'identifiant et conservez la somme de contrôle. Le nom des contigs compte autant que la séquence : un BAM aligné sur NC_000913.3 ne pourra pas être utilisé avec un VCF qui référence chr, et inversement.

Pour E. coli, la référence est un chromosome unique, sans contigs alternatifs ni décoys. Pour un génome humain, il faut choisir entre l'assemblage RefSeq, l'assemblage Ensembl et les versions avec ou sans contigs ALT et decoy ; ce choix change les alignements dans les régions dupliquées. Fixez une référence, utilisez-la pour l'index, le mapping, l'appel de variants et l'annotation, et ne la remplacez jamais silencieusement. samtools faidx doit être exécuté sur le même FASTA que bwa index.

BWA-MEM, BWA-MEM2 et minimap2

BWA-MEM est la référence pour les lectures courtes Illumina, en particulier appariées. BWA-MEM2 réimplémente le même algorithme avec des instructions vectorielles : plus rapide, plus gourmand en mémoire, avec un index différent. Ses sorties sont conçues pour être équivalentes, mais validez-les sur vos données et ne mélangez pas les index. Minimap2 vise d'abord les lectures longues (Oxford Nanopore, PacBio) et les assemblages ; son préréglage -x sr existe pour les lectures courtes, avec des conventions de MAPQ et de clipping différentes.

Reproductibilité

Une commande de mapping n'est reproductible que si l'on fixe six éléments : version de l'aligneur, version de SAMtools, FASTA de référence avec sa somme de contrôle, index exact, paramètres et ordre des fichiers. Ajoutez le read group et le nom de l'échantillon. bwa mem écrit une ligne @PG avec la commande et la version ; conservez-la, mais ne considérez pas le BAM comme stable octet pour octet : compression, ordre des records et lignes @PG peuvent varier. Comparez les alignements triés, pas les sommes de contrôle du BAM.

Le nombre de threads ne devrait pas changer les alignements, mais il peut changer l'ordre de production des records et la taille des lots. L'option -K fixe la taille des lots et facilite la comparaison entre exécutions. samtools sort normalise l'ordre par coordonnées ; pour les positions identiques, l'ordre peut encore dépendre de l'implémentation. Un manifeste qui enregistre versions, chemins, sommes de contrôle et read groups vaut mieux qu'un long message dans un carnet.

Implémentation

Nous partons d'un FASTA de référence et de deux FASTQ appariés. La préparation tient en trois commandes : télécharger, vérifier, indexer.

set -euo pipefail
 
efetch -db nuccore -id NC_000913.3 -format fasta > refs/NC_000913.3.fa
sha256sum refs/NC_000913.3.fa > refs/NC_000913.3.fa.sha256
bwa index refs/NC_000913.3.fa
samtools faidx refs/NC_000913.3.fa

set -euo pipefail n'est pas décoratif : dans un pipeline, si bwa mem échoue et que samtools sort réussit sur une entrée vide, le code de sortie final peut être zéro. Sans pipefail, un BAM vide passe pour un succès.

La commande de production enchaîne le mapping et le tri. Le read group est passé à BWA, qui le recopie dans chaque enregistrement ; samtools sort écrit un BAM trié par coordonnées.

bwa mem -t 8 -R '@RG\tID:lane1\tSM:ecoli_ctrl\tPL:ILLUMINA\tLB:lib1\tPU:fc1.lane1' \
  refs/NC_000913.3.fa \
  data/ecoli_ctrl_R1.fastq.gz \
  data/ecoli_ctrl_R2.fastq.gz \
  | samtools sort -@ 8 -o ecoli_ctrl.bam -
samtools index ecoli_ctrl.bam
samtools quickcheck -v ecoli_ctrl.bam

Le tiret final de samtools sort lui indique de lire l'entrée standard. samtools index crée ecoli_ctrl.bam.bai, nécessaire aux accès par région. samtools quickcheck vérifie l'en-tête et la fin du fichier ; il ne valide pas la biologie, seulement l'intégrité du BAM.

Contrôlez ensuite le read group et les compteurs :

samtools view -H ecoli_ctrl.bam | grep '^@RG'
samtools view ecoli_ctrl.bam | head -n 2
samtools flagstat ecoli_ctrl.bam
samtools idxstats ecoli_ctrl.bam

La première commande doit afficher une ligne @RG avec les cinq champs. La deuxième doit montrer RG:Z:lane1 sur chaque lecture. flagstat résume les lectures totales, appariées, correctement appariées, non alignées et de MAPQ suffisante ; idxstats donne le nombre d'alignements par contig. Aucun de ces compteurs ne prouve que l'alignement est correct, mais leur cohérence détecte les erreurs de pipeline.

Les fonctions Python de l'exercice construisent la commande et valident les métadonnées sans exécuter BWA. La construction sous forme de liste évite les problèmes de quoting et permet de tester les arguments un par un.

def build_bwa_command(reference, r1, r2, sample, threads=8):
    read_group = f"@RG\\tID:{sample}\\tSM:{sample}\\tPL:ILLUMINA"
    return [
        "bwa", "mem",
        "-t", str(threads),
        "-R", read_group,
        reference, r1, r2,
    ]

Le read group contient ici les trois champs minimaux de la commande de cours. Pour un vrai projet, ajoutez LB et PU, et vérifiez que SM correspond au manifeste. Le mini-projet impose les cinq champs.

Étude de cas

Nous utilisons les lectures appariées synthétiques du dépôt, datasets/ngs/reads_R1.fastq et reads_R2.fastq, et la référence réelle E. coli K-12 MG1655 NC_000913.3.

Les fichiers contiennent deux paires : pair1/1 et pair1/2, puis pair2/1 et pair2/2. Chaque lecture mesure 12 pb, soit 48 pb au total. La seconde paire a quatre bases à Q2 à son extrémité 3'. Pour une de ces lectures, la probabilité d'au moins une erreur vaut $1-(1-10^{-4})^8(1-0{,}630957)^4 \approx 0{,}9815$, soit environ 98 %. Le jeu est donc très petit et localement de mauvaise qualité ; il sert à tester la mécanique, pas à conclure sur E. coli.

La question est simple : ces lectures peuvent-elles être mappées avec les paramètres par défaut ? Après indexation, lancez la commande avec le read group et sans modifier -k. BWA-MEM utilise alors une longueur minimale de seed de 19 pb. Aucune lecture de 12 pb ne peut fournir un MEM de 19 pb : les quatre lectures ressortent non alignées, avec le flag 0x4. Le résumé attendu de samtools flagstat ressemble à ceci.

4 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 mapped (0.00% : N/A)
4 + 0 unmapped
4 + 0 paired in sequencing
2 + 0 read1
2 + 0 read2

Cette sortie est illustrative : elle décrit le résultat attendu, elle n'a pas été produite par une exécution de BWA dans cette session. Le point important n'est pas le zéro, mais sa cause : le seuil de seed est un filtre silencieux. Une lecture trop courte est invisible, même parfaitement identique à la référence.

Pour vérifier l'effet du paramètre, relancez avec -k 8. Les lectures peuvent alors produire des seeds, et certains alignements apparaissent. Leur MAPQ est souvent faible, car des seeds de 8 pb sont abondants dans un génome de 4,6 Mb ; le nombre exact d'alignements dépend de la référence et ne doit pas être interprété comme une validation. La bonne conclusion est qu'il faut adapter la longueur des lectures au projet, pas baisser -k jusqu'à obtenir des alignements.

Le fichier datasets/ngs/demo.sam montre la structure attendue d'un BAM correct, avec un en-tête @RG et des flags de paire. Il est fourni comme sortie illustrative, pas comme résultat de BWA.

@HD	VN:1.6	SO:coordinate
@SQ	SN:chr1	LN:1000
@RG	ID:demo	SM:synthetic	PL:ILLUMINA
pair1	99	chr1	10	60	12M	=	40	42	ACGTACGTACGT	IIIIIIIIIIII	RG:Z:demo
pair1	147	chr1	40	60	12M	=	10	-42	TGCATGCATGCA	IIIIIIIIIIII	RG:Z:demo
pair2	0	chr1	100	30	12M	*	0	0	TTTTGGGGCCCC	IIIIIIII####	RG:Z:demo

Les flags 99 et 147 décrivent une paire correcte, et 0 une lecture non appariée. Le read group demo est répété sur chaque ligne. Ce fichier ne prouve pas que BWA produirait ces coordonnées ; il montre les champs à contrôler. Dans un vrai projet, on vérifierait aussi la distribution des insert sizes, la proportion de paires correctement appariées et la couverture par contig.

Erreurs fréquentes

Quiz

Exercices

python -m pytest exercises/module-3/3.2-mapping-bwa-mem -q

python -m pytest exercises/module-3/3.2-mapping-bwa-mem -q

python -m pytest exercises/module-3/3.2-mapping-bwa-mem/projet -q

Résumé

• bwa index construit un index FM lié au FASTA ; l'index et la référence doivent rester appariés et documentés par une somme de contrôle.
• BWA-MEM cherche des MEM, les chaîne en ordre colinéaire, puis étend chaque chaîne par un alignement local à gaps affines.
• -k 19 est un filtre silencieux : une lecture plus courte ne produit aucun seed et ressort non alignée.
• Les lectures appariées se fournissent dans l'ordre R1 puis R2 ; BWA apparie par position, pas par nom.
• Un read group complet contient ID, SM, PL, LB et PU ; SM relie les fichiers à l'échantillon biologique.
• MAPQ mesure la confiance dans la position ; secondaires (0x100) et supplémentaires (0x800) ne sont ni des duplicats ni des lectures non alignées.
• La reproductibilité exige versions, index, commande, ordre des fichiers et read groups ; comparez les alignements, pas les octets du BAM.

Références

• Li, H. & Durbin, R. (2009). Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics, 25(14), 1754–1760.
• Li, H. (2013). Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv:1303.3997.
• Danecek, P. et al. (2021). Twelve years of SAMtools and BCFtools. GigaScience, 10(2), giab008.
• Vasimuddin, M., Misra, S., Li, H. & Aluru, S. (2019). Efficient architecture-aware acceleration of BWA-MEM for multicore systems. IEEE IPDPS, 314–324.
• Li, H. (2018). Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics, 34(18), 3094–3100.
• BWA. Manual page and source code. https://bio-bwa.sourceforge.net/bwa.shtml
• SAMtools. Documentation. https://www.htslib.org/doc/samtools.html
• SAM/BAM. Sequence Alignment/Map Format Specification. https://samtools.github.io/hts-specs/SAMv1.pdf