SNP avec BCFtools

Prérequis

Cette section prolonge la section 3.4 : vous devez savoir trier un BAM, marquer les duplicats avec fixmate -m puis markdup, indexer le fichier et lire les champs FLAG, CIGAR et MAPQ. Les notions de profondeur et de qualité de base viennent de 1.2. Le pipeline suppose une référence FASTA indexée par samtools faidx et un BAM trié et indexé. Les exercices n'utilisent que Python standard.

Introduction

Un génome séquencé ne révèle pas directement ses variants. Les lectures alignées forment un faisceau de preuves indirectes : à chaque position, plusieurs lectures couvrent la référence, certaines portent une base différente, d'autres sont mal appariées ou proviennent d'un paralogue. Décider qu'une différence est un vrai variant — et non une erreur de séquençage, un artefact de PCR ou un défaut d'alignement — est le cœur de l'appel de variants. L'enjeu est considérable : une mutation dans gyrA peut conférer une résistance aux fluoroquinolones chez Escherichia coli K-12, un variant somatique peut initier une tumeur, et une substitution dans un gène de réparation de l'ADN peut expliquer une maladie héréditaire. La sortie brute d'un appelant reste une liste d'hypothèses à filtrer et à valider.

Le défi computationnel consiste à transformer un BAM en un VCF (Variant Call Format) défendable. Il faut agréger les lectures par position, estimer la probabilité de chaque génotype, tenir compte d'un a priori de mutation, puis décider si le site est variable. BCFtools fournit deux étages : bcftools mpileup calcule les vraisemblances de génotype à partir du pileup, et bcftools call applique le modèle statistique de Li (2011). Le résultat contient des sites multialléliques, des indels représentés de plusieurs façons et des artefacts liés aux répétitions. Un pipeline reproductible normalise, indexe et filtre ces enregistrements avec des seuils explicites.

Nous suivons la chaîne de l'alignement au VCF filtré, puis lisons les champs qui jugent un appel. L'étude de cas utilise datasets/ngs/demo.vcf, un fichier synthétique de trois variants. Les commandes BCFtools sont fournies avec leurs options ; les sorties d'outils non exécutés sont signalées comme illustratives.

Théorie

Du pileup au génotype

À une position donnée, mpileup empile les bases des lectures qui couvrent le site. Pour chaque lecture, il connaît la base observée, sa qualité et sa qualité de mapping. Il en déduit une vraisemblance pour chaque génotype possible : homozygote référence, hétérozygote, homozygote alternatif, et combinaisons pour les sites multialléliques.

bcftools call combine ces vraisemblances avec un a priori de mutation et choisit le génotype le plus probable. Le modèle suppose que les erreurs de lecture sont indépendantes et que la profondeur reflète l'échantillonnage. Les duplicats de PCR, les lectures chevauchantes d'un même fragment et les régions répétées violent cette hypothèse : deux lectures identiques ne valent pas deux preuves indépendantes.

Vraisemblance, QUAL et GQ

Pour un site et des données $D$, la vraisemblance d'un génotype $G$ est $P(D \mid G)$. La règle de Bayes donne la probabilité posterior :

où $P(G)$ résume l'a priori de mutation et d'hétérozygotie. L'appelant retient le génotype de posterior maximal. La qualité du site, notée QUAL, est une probabilité déguisée en score de Phred :

Un QUAL de 20 correspond à $10^{-2}$, soit environ 1 % de chance que le site soit en réalité invariant selon le modèle ; un QUAL de 30 correspond à 0,1 %. Ce n'est pas la probabilité que le variant soit vrai, ni un pourcentage de lectures alternatives. GQ répond à une autre question :

Un site peut avoir un QUAL élevé et un GQ médiocre si plusieurs génotypes restent plausibles. Les deux scores ne se remplacent pas.

Lire un VCF : GT, DP et AD

Un VCF mêle des champs de site et des champs de format répétés par échantillon.

Dans GT, l'allèle 0 désigne REF, 1 le premier ALT, 2 le second. La barre oblique / note un génotype non phasé ; la barre verticale | un génotype phasé. Ainsi 0/1 est hétérozygote, 1/1 homozygote pour le premier alternatif, 1/2 hétérozygote composite et ./. manquant. La fraction d'allèles alternatifs se calcule avec AD :

Pour un hétérozygote diploïde sans biais, on attend $f_{\text{alt}} \approx 0{,}5$. Une valeur de 0,2 sur 10 lectures n'est pas nécessairement un faux positif : l'échantillonnage binomial produit une variance non nulle. DP et AD ne sont pas interchangeables : DP compte toutes les lectures retenues, AD seulement celles attribuées à chaque allèle.

Multiallélices et normalisation

Un site peut porter plusieurs alternatives : ALT=C,A avec GT=1/2. Beaucoup d'outils d'annotation et de comparaison supposent un VCF biallélique. bcftools norm -m -any découpe chaque site multiallélique en plusieurs enregistrements, un par allèle alternatif. Le génotype 1/2 devient deux lignes cohérentes, annotables par allèle.

La normalisation traite aussi la représentation des indels. Un même événement peut s'écrire de plusieurs façons selon l'ancre choisie. Avec REF=GCA et ALT=GTA, le suffixe A est commun : la forme minimale est GC/GT, et la substitution est C>T à la position suivante. bcftools norm -f reference.fa aligne l'indel à gauche dans les répétitions, retire les bases communes et vérifie que REF correspond à la référence. Sans cette étape, deux fichiers décrivent le même indel à des positions différentes et les comparaisons échouent.

Après découpage, des doublons peuvent apparaître. bcftools norm -d exact supprime les enregistrements strictement identiques. La normalisation doit utiliser exactement la référence qui a servi à l'alignement : un contig absent, une version différente ou un REF incorrect produisent des erreurs ou des variants fantômes.

SNV, indels et faux positifs

Une SNV (single-nucleotide variant) remplace une base par une autre. Un indel ajoute ou retire une ou plusieurs bases. Un MNV remplace plusieurs bases contiguës ; certains outils le découpent en SNV, d'autres le conservent comme un bloc. Un variant complexe combine substitutions et changements de longueur. Les indels sont plus sensibles à l'ambiguïté d'alignement, aux répétitions en tandem et aux homopolymères, où la polymérase peut glisser.

Les faux positifs ont des causes récurrentes. Une lecture de paralogue se mappe sur le gène voisin ; les lectures alternatives semblent cohérentes. Un artefact de PCR produit de nombreux duplicats qui gonflent la profondeur sans apporter de preuve indépendante. Un biais de brin, où les lectures alternatives viennent toutes du même brin, signale souvent une erreur de séquençage ou de préparation. Les indels dans les homopolymères et les zones de faible complexité sont également fragiles. Enfin, une contamination ou un échantillon mal identifié peut créer des variants massifs et incohérents.

Un truth set est un ensemble de variants considérés comme vrais, obtenu par séquençage orthogonal. Pour le génome humain, les jeux de référence du consortium GIAB (Genome in a Bottle) servent à mesurer précision et rappel. Pour E. coli, simulez des lectures ou validez par Sanger un échantillon de sites. Un truth set sert à évaluer un pipeline, pas à filtrer les variants d'un nouvel échantillon ni à prouver qu'un variant est pathogène.

Filtrer sans confondre les échelles

Un filtre dur supprime les enregistrements ; un filtre souple les conserve avec une étiquette dans FILTER. Le filtre dur allège les fichiers, le filtre souple garde la trace des sites écartés. bcftools view -i inclut les sites qui satisfont une expression, -e exclut ceux qui la satisfont, et -f PASS ne garde que la valeur PASS.

bcftools view -f PASS -i 'QUAL>=20 && INFO/DP>=10 && INFO/MQ>=30' variants.vcf.gz

Les seuils sont inclusifs et se lisent avec les opérateurs &&, ||, ! et des parenthèses. Il faut citer l'expression pour que le shell ne l'interprète pas. Les seuils ne sont pas universels : QUAL>=20 et DP>=10 sont un point de départ pour un génome bactérien ou un exome, pas une règle biologique. Un seuil de profondeur trop bas laisse passer des génotypes non soutenus ; un seuil trop haut élimine les vrais hétérozygotes dans les régions difficiles. MQ>=30 écarte les lectures mal mappées, mais peut supprimer un vrai variant dans une région dupliquée. GQ>=20 contrôle la confiance du génotype, pas la qualité du site.

Le choix de INFO/DP ou de FORMAT/DP dépend de la question. Pour un échantillon unique, les deux sont souvent proches. Pour plusieurs échantillons, INFO/DP est la somme du site et FORMAT/DP la profondeur de chaque génotype ; filtrer sur la somme peut masquer un échantillon à couverture nulle.

Couverture, mapping et limites

La profondeur conditionne la détection. Si les lectures arrivent selon un processus de Poisson de moyenne $\lambda$, la probabilité de n'observer aucune lecture est $e^{-\lambda}$. Pour un hétérozygote, la profondeur attendue de l'allèle alternatif vaut $\lambda/2$, donc la probabilité de manquer totalement l'allèle est $e^{-\lambda/2}$ et la probabilité de le voir au moins une fois est

À 10 lectures au total, cette probabilité vaut environ 99,3 % ; à 3 lectures, 77,7 % ; à 2 lectures, 63,2 %. Ces chiffres supposent un échantillonnage parfait. En pratique, le biais de GC, la capture, les régions répétées et les erreurs de base réduisent la profondeur effective. La moyenne de couverture cache des trous locaux : un génome à 30× peut avoir des positions à 0×.

À l'inverse, une profondeur élevée n'est pas une garantie. Les lectures mal mappées s'accumulent dans les paralogues et les répétitions, et un seuil de profondeur maximal (-d) tronque le pileup. Les indels déplacent l'alignement et créent des piles locales. Le mapping de reads courts sur une référence divergente produit des variants faux. La qualité de mapping, la couverture des deux brins et l'uniformité sont donc aussi importantes que la profondeur moyenne.

Exemple travaillé

Prenons un site couvert par 20 lectures, dont 12 portent la référence et 8 l'alternative. La fraction alternative vaut $8/20 = 0{,}4$. Pour un hétérozygote diploïde, on attend 10 lectures alternatives ; l'écart de 2 lectures représente $\sqrt{20 \times 0{,}5 \times 0{,}5} \approx 2{,}24$ lectures, soit moins d'un écart-type binomial. Le génotype 0/1 reste parfaitement plausible. Avec 5 lectures, la puissance reste faible et un seuil DP>=10 écarte le site.

Ce calcul suppose des lectures indépendantes. Si les 8 lectures alternatives sont des duplicats d'un même fragment, la profondeur effective est bien inférieure à 20 et l'intervalle de confiance est faux. C'est pourquoi AD, DP et les indicateurs de duplicats se lisent ensemble.

Implémentation

Le pipeline ci-dessous utilise BCFtools 1.19 ou ultérieur et SAMtools 1.19 ou ultérieur. La référence doit être la même que celle du BAM, avec un index .fai ; le BAM doit être trié et indexé. Les commandes sont écrites pour Bash, avec set -euo pipefail : une erreur dans le pipe interrompt le script au lieu de produire un VCF partiel.

set -euo pipefail
 
REF="reference.fa"
BAM="aln.bam"
 
# 1. Vraisemblances de génotype, puis appel des sites variants.
bcftools mpileup -f "$REF" -q 20 -Q 20 -d 500 -Ou "$BAM" \
  | bcftools call -mv -Oz -o variants.vcf.gz
 
# 2. Index tabix du VCF compressé.
bcftools index -t variants.vcf.gz
 
# 3. Normalisation et découpage des sites multialléliques.
bcftools norm -f "$REF" -m -any -Oz -o variants.norm.vcf.gz variants.vcf.gz
bcftools index -t variants.norm.vcf.gz
 
# 4. Filtre dur : PASS, qualité, profondeur totale.
bcftools view -f PASS -i 'QUAL>=20 && INFO/DP>=10' \
  -Oz -o variants.filtered.vcf.gz variants.norm.vcf.gz
bcftools index -t variants.filtered.vcf.gz
 
# 5. Contrôle : nombre de variants et répartition.
bcftools stats variants.filtered.vcf.gz > variants.stats

-Ou écrit du BCF non compressé, adapté à un pipe ; -Oz écrit du VCF compressé par BGZF, lisible par zcat et indexable. -q 20 écarte les lectures de qualité de mapping inférieure à 20, -Q 20 les bases de qualité inférieure à 20, et -d 500 plafonne la profondeur à 500 lectures par site. Sans -d, la valeur par défaut peut tronquer un site très couvert ; avec une valeur trop basse, vous perdez des lectures utiles. -m active l'appel multiallélique et -v ne conserve que les sites variants.

La normalisation intervient avant le filtre pour que les expressions portent sur des enregistrements bialléliques. bcftools norm -f "$REF" -m -any aligne les indels à gauche, vérifie REF et découpe les multivariants. Si des doublons subsistent, ajoutez bcftools norm -d exact dans une seconde passe. L'indexation doit suivre chaque écriture : un VCF compressé sans .tbi ne peut pas être interrogé par région.

Pour inspecter les champs utiles, bcftools query évite d'ouvrir le VCF à la main :

bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\t%INFO/DP[\t%GT\t%DP\t%AD]\n' \
  variants.filtered.vcf.gz

L'annotation fonctionnelle se fait après le filtrage, avec VEP, SnpEff ou une base de variants. bcftools annotate -a annotations.vcf.gz -c ID transfère par exemple les identifiants d'une base de connaissances. Annotez un VCF déjà normalisé, sinon les correspondances d'indels et de multiallélices sont fragiles.

Étude de cas

Le fichier datasets/ngs/demo.vcf est un VCF synthétique de 1 000 pb, avec un contig chr1 et un échantillon S1. Il contient trois variants et ne définit que INFO/DP, FORMAT/GT et FORMAT/DP : ni AD, ni MQ, ni GQ. Question : quels variants passent les seuils QUAL>=20, DP>=10 et FILTER=PASS, et que peut-on réellement conclure ?

Un filtre Python équivalent, exécuté le 2026-11-01 sur ce fichier, donne n_total=3, n_pass=2 et n_filtered=1. Le variant 120 cumule trois motifs de rejet : FILTER=LowQual, QUAL=5 inférieur à 20 et DP=3 inférieur à 10. La commande BCFtools correspondante produirait les deux mêmes lignes ; sa sortie est illustrative, car BCFtools n'a pas été exécuté dans cet environnement :

bcftools view -f PASS -i 'QUAL>=20 && INFO/DP>=10' datasets/ngs/demo.vcf

chr1	50	.	A	G	35	PASS	DP=18	GT:DP	0/1:18
chr1	200	.	G	A	60	PASS	DP=25	GT:DP	1/1:25

Le variant 50 a une qualité moyenne, une profondeur suffisante et un génotype hétérozygote. Sans AD, on ne peut pas vérifier l'équilibre des allèles : un 0/1 avec 17 lectures référence et 1 lecture alternative serait suspect, mais le VCF ne permet pas de le savoir. Le variant 200 a une bonne qualité et une profondeur confortable ; 1/1 suggère un homozygote alternatif, mais l'absence de MQ empêche d'exclure un paralogue ou une région dupliquée. Le variant 120 ne doit pas être sauvé en abaissant seulement QUAL : sa profondeur de 3 ne permet pas un génotype fiable.

Enfin, demo.vcf est synthétique : les positions ne correspondent pas au génome NC_000913.3 et aucune conclusion biologique sur un gène n'est possible. Le fichier vérifie la mécanique de lecture et de filtrage, pas un pipeline sur données réelles.

Erreurs fréquentes

Quiz

Exercices

python -m pytest exercises/module-3/3.5-appel-variants-bcftools -q

python -m pytest exercises/module-3/3.5-appel-variants-bcftools -q

python -m pytest exercises/module-3/3.5-appel-variants-bcftools/projet -q

Résumé

• mpileup produit les vraisemblances, call -mv choisit les génotypes et garde les sites variants.
• -Ou transporte un BCF non compressé dans le pipe, -Oz écrit un VCF compressé et indexable.
• QUAL qualifie le site, GQ le génotype, MQ le mapping ; DP et AD décrivent la profondeur et l'équilibre des allèles.
• bcftools norm -f ref -m -any aligne les indels, vérifie REF et découpe les multiallélices ; -d exact retire les doublons.
• Les filtres QUAL, DP, MQ et GQ sont des points de départ à adapter à la couverture, à la ploïdie et aux régions difficiles.
• Les indels, les répétitions et les paralogues produisent des faux positifs que la profondeur seule ne révèle pas.
• Un truth set évalue un pipeline ; il ne remplace ni la validation orthogonale ni l'analyse des limites de couverture et de mapping.

Références

• Li, H. (2011). A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics, 27(21), 2987–2993.
• Danecek, P. et al. (2021). Twelve years of SAMtools and BCFtools. GigaScience, 10(2), giab008.
• Li, H. et al. (2009). The Sequence Alignment/Map format and SAMtools. Bioinformatics, 25(16), 2078–2079.
• Danecek, P. et al. (2011). The variant call format and VCFtools. Bioinformatics, 27(15), 2156–2158.
• DePristo, M. A. et al. (2011). A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nature Genetics, 43(5), 491–498.
• BCFtools. Documentation and manual pages. https://samtools.github.io/bcftools/bcftools.html
• SAMtools. Documentation and manual pages. https://www.htslib.org/doc/samtools.html