Formats SAM et BAM

Prérequis

Cette section s'appuie sur les fichiers FASTQ de la section 1.2 et sur les notions d'alignement présentées en 3.2. Vous devez savoir lire un enregistrement de séquençage, distinguer une lecture de son brin complémentaire et interpréter un score d'alignement. Côté programmation, les listes, les dictionnaires, les chaînes et les conversions d'entiers en Python sont supposés acquis. Aucun outil externe n'est nécessaire pour les exercices : tout se calcule à partir du texte SAM.

Introduction

Une expérience de séquençage produit des millions de lectures courtes. Pour savoir à quelle position du génome chacune correspond, on aligne ces lectures sur une séquence de référence. L'aligneur doit conserver bien plus que la seule coordonnée : le brin, les éventuelles insertions et délétions, la qualité de l'emplacement, l'appariement des lectures et le statut de chaque lecture dans un éventuel doublon. Ces informations sont condensées dans le format SAM (Sequence Alignment/Map), devenu le langage commun de la génomique.

L'enjeu biologique est direct. Une erreur de position décale tous les variants d'une base ; un drapeau mal interprété fait compter deux fois la même molécule ; une couverture mal calculée donne une fausse impression de profondeur et fausse l'appel de variants. Les duplicats techniques, les alignements secondaires et les régions épissées ont chacun une signature dans le format. Savoir les lire, c'est pouvoir contrôler un pipeline au lieu de subir ses sorties.

Le défi computationnel vient de l'échelle et de la richesse du format. Un fichier SAM humain complet pèse plusieurs centaines de gigaoctets ; BAM le compresse, CRAM le réduit encore en s'appuyant sur la référence. Un index permet d'accéder à une région sans parcourir le fichier entier, mais il impose un tri préalable. Enfin, le format est tabulaire et tolérant en apparence : une ligne mal découpée, une qualité de mauvaise longueur ou un champ manquant produisent souvent un résultat silencieusement faux plutôt qu'une erreur immédiate. Cette section apprend à lire, calculer et vérifier.

Théorie

Un format texte, onze colonnes

SAM est un format textuel : chaque alignement occupe une ligne, et les champs sont séparés par des tabulations. Les onze premiers champs sont obligatoires et ont une signification fixe. Les tags optionnels viennent ensuite, sous la forme TAG:TYPE:VALEUR, par exemple RG:Z:demo pour rattacher la lecture à un groupe de lectures.

Une ligne SAM n'est donc pas une ligne de texte libre. Le nombre de colonnes, l'ordre et le type de chaque champ comptent. Un simple split() sur les espaces peut fusionner ou décaler des champs ; nous y reviendrons.

L'en-tête

Les lignes d'en-tête commencent par @ et précèdent les alignements. @HD décrit la version du format et l'ordre de tri, souvent VN:1.6 et SO:coordinate. @SQ déclare chaque séquence de référence avec son nom (SN) et sa longueur (LN) : c'est ce qui donne un sens aux coordonnées de RNAME et POS. @RG décrit un groupe de lectures, avec un identifiant (ID), un échantillon (SM) et une plateforme (PL). @PG trace le programme qui a produit le fichier, et @CO contient un commentaire libre.

L'en-tête n'est pas décoratif. Sans @SQ, on ne connaît ni la longueur de la référence ni l'ensemble des chromosomes attendus ; sans @RG, on ne peut pas rattacher les lectures à un échantillon ni à une librairie, ce qui est indispensable pour l'appel de variants. samtools view -H fichier.bam affiche l'en-tête seul, et samtools view -h l'inclut avec les alignements.

FLAG : un entier, douze propriétés

Le champ FLAG est un entier compris entre 0 et 4095. Chaque bit représente une propriété indépendante : la valeur se décode avec des puissances de deux, pas avec une table de 4 096 entrées. La formule est simple : si le bit $i$ vaut 1, la propriété correspondante est vraie, et un entier est la somme des bits actifs.

Un FLAG de 99 vaut $64 + 32 + 2 + 1$ : la lecture est appariée, correctement appariée, son mate est sur le brin reverse et c'est la première de la paire. Un FLAG de 147 vaut $128 + 16 + 2 + 1$ : c'est la seconde lecture, sur le brin reverse, dans une paire correcte. Le bit 0x4 signale une lecture non alignée ; le bit 0x8, un mate non aligné. Une lecture primaire n'est ni secondaire (0x100) ni supplémentaire (0x800). Un secondaire est un alignement alternatif de la même lecture ; un supplémentaire correspond à une partie chimérique ou divisée. Un duplicat (0x400) est marqué par un outil dédié, pas par l'aligneur.

CIGAR : décrire un alignement avec des opérations

Le champ CIGAR (Compact Idiosyncratic Gapped Alignment Report) décrit la suite d'opérations qui transforme la lecture en alignement sur la référence. Chaque opération est une longueur suivie d'une lettre. Les neuf opérations standard se répartissent selon ce qu'elles consomment.

Deux formules résument l'essentiel. L'emprise sur la référence est $L_{\mathrm{ref}} = M + D + N + = + X$, tandis que la longueur de lecture décrite par le CIGAR est $L_{\mathrm{lecture}} = M + I + S + = + X$. Pour 5S90M2I10M3D5M, la référence consommée vaut $90 + 10 + 3 + 5 = 108$ bases et la lecture décrite vaut $5 + 90 + 2 + 10 + 5 = 112$ bases. Le soft clip conserve les bases dans SEQ mais ne les aligne pas ; le hard clip les retire complètement de SEQ. N sert surtout aux alignements épissés d'ARN, tandis que D correspond à une délétion dans la molécule alignée. Le coût de lecture d'un CIGAR est linéaire en nombre d'opérations, pas en longueur des opérations.

POS, MAPQ et TLEN

POS est la position 1-based de la première base de référence alignée. C'est une convention du format SAM, alors que Python indexe à partir de 0 et que le format BED utilise des intervalles 0-based semi-ouverts. Une conversion oubliée décale tout d'une base. Pour un alignement 12M à POS 10, les bases de référence couvertes vont de 10 à 21 inclus.

MAPQ mesure la confiance dans l'emplacement, sur une échelle inspirée de Phred : $Q = -10\log_{10}P$, où $P$ est la probabilité que l'emplacement soit faux. Une valeur de 60 correspond à une probabilité d'erreur d'environ $10^{-6}$. La valeur 255 signale une qualité inconnue, et 0 un alignement ambigu ou multi-positionné. MAPQ ne mesure pas la qualité des bases : c'est le rôle de QUAL.

TLEN (template length) donne la longueur du fragment observé entre les deux lectures d'une paire, du début le plus à gauche à la fin la plus à droite, bornes incluses. Sa valeur est positive pour la lecture la plus à gauche et négative pour l'autre ; 0 signale une information indisponible. Pour pair1, les deux alignements 12M commencent à 10 et 40 : la fin droite vaut $40 + 12 - 1 = 51$, donc $TLEN = 51 - 10 + 1 = 42$.

SEQ et QUAL

SEQ contient la séquence de la lecture, et QUAL les qualités Phred+33 de même longueur. Le caractère I correspond à $Q = 40$ et # à $Q = 2$ ; l'encodage Phred+33 signifie que la valeur ASCII est décalée de 33. Une qualité absente est notée *. Lorsqu'une lecture est alignée sur le brin reverse (0x10), SAM stocke SEQ sous forme complémentaire inverse de la lecture d'origine et inverse QUAL en conséquence. Il ne faut donc pas appliquer un second complément inverse sans vérifier le drapeau.

La longueur de SEQ doit correspondre à la somme des opérations M, I, S, =, X du CIGAR, et QUAL doit avoir la même longueur que SEQ. Ces deux égalités sont les vérifications les plus rentables d'un parseur : elles détectent immédiatement un décalage de colonnes ou une troncature de ligne. Les tags optionnels, comme NM:i:1 ou RG:Z:demo, ajoutent des informations sans modifier les onze champs obligatoires.

BAM, CRAM et index

BAM est la version binaire compressée de SAM : mêmes champs, mêmes en-têtes, mais un encodage compact et un flux BGZF. CRAM va plus loin en codant les lectures par rapport à la référence : le fichier est souvent plus petit, au prix d'une dépendance à la bonne version de la référence. BAM conserve toutes les informations SAM sans perte ; CRAM est sans perte par défaut, mais certaines options de compression des qualités peuvent être avec perte. Vérifiez donc les paramètres avant d'utiliser un CRAM comme archive.

Un index, .bai pour BAM ou .crai pour CRAM, associe des intervalles de coordonnées à des blocs du fichier. Il permet à samtools view fichier.bam chr1:10-51 de ne lire que la région demandée. L'index exige un fichier trié par coordonnées : samtools sort produit cet ordre, puis samtools index construit l'index. Un BAM non trié peut être lu séquentiellement, mais pas indexé correctement. samtools idxstats résume le nombre de lectures par référence, et samtools quickcheck -v vérifie l'intégrité de l'en-tête et de la fin de fichier.

Exemple travaillé

Prenons la première ligne d'alignement de datasets/ngs/demo.sam, un fichier synthétique (sortie illustrative) :

pair1	99	chr1	10	60	12M	=	40	42	ACGTACGTACGT	IIIIIIIIIIII	RG:Z:demo

Le nom est pair1, le FLAG 99 se décode en appariée, correctement appariée, mate reverse, première de la paire. La référence est chr1, la position 10, la qualité d'emplacement 60 et le CIGAR 12M. Le mate est sur la même référence (=), à la position 40, et TLEN vaut 42. La lecture couvre donc les positions 10 à 21, et son mate les positions 40 à 51. Les douze qualités I valent 40, et le tag RG:Z:demo rattache la lecture au groupe demo. Aucun bit 0x4, 0x100, 0x400 ou 0x800 n'est actif : la lecture est alignée, primaire et non dupliquée.

Calculer une couverture

La profondeur à une position $j$ de la référence est le nombre de lectures dont une base alignée couvre cette position. Toutes les bases ne comptent pas : M, = et X ajoutent de la profondeur, tandis que D et N avancent sur la référence sans ajouter de lecture, et I, S, H, P n'y contribuent pas. Pour un ensemble d'alignements, la profondeur moyenne sur une référence de longueur $L$ vaut $\bar{d} = \frac{1}{L}\sum_{j=1}^{L} d_j$.

Sur demo.sam, les trois alignements 12M couvrent 36 bases de référence au total, sans chevauchement. La profondeur moyenne est donc $36 / 1000 = 0,036$, et la fraction de bases couvertes vaut également $36 / 1000 = 3,6\,\%$. Une profondeur moyenne de 0,036× est dérisoire pour appeler des variants : elle sert ici à vérifier le calcul, pas à produire une conclusion biologique. L'algorithme parcourt chaque CIGAR et met à jour un tableau de profondeur ; son coût est $O(L + B)$, où $B$ est le nombre de bases alignées, avec une mémoire $O(L)$.

Implémentation

Les fonctions suivantes n'utilisent que la bibliothèque standard. Elles supposent qu'un enregistrement est un dictionnaire avec pos et cigar, ce qui suffit pour la couverture ; le parseur complet du mini-projet ajoutera les autres champs.

import re
 
FLAG_BITS = {
    "paired": 0x1, "proper_pair": 0x2, "unmapped": 0x4,
    "mate_unmapped": 0x8, "reverse": 0x10, "mate_reverse": 0x20,
    "first_in_pair": 0x40, "second_in_pair": 0x80,
    "secondary": 0x100, "qc_fail": 0x200,
    "duplicate": 0x400, "supplementary": 0x800,
}
REFERENCE_CONSUMING = frozenset("MDN=X")
 
 
def decode_flag(flag):
    if flag < 0:
        raise ValueError("flag must be non-negative")
    return {name: bool(flag & bit) for name, bit in FLAG_BITS.items()}
 
 
def parse_cigar(cigar):
    if cigar == "*":
        return []
    if not re.fullmatch(r"(?:\d+[MIDNSHP=X])+", cigar):
        raise ValueError(f"invalid CIGAR: {cigar!r}")
    ops = [(int(n), op) for n, op in re.findall(r"(\d+)([MIDNSHP=X])", cigar)]
    if any(n <= 0 for n, _ in ops):
        raise ValueError(f"invalid CIGAR: {cigar!r}")
    return ops
 
 
def reference_span(cigar):
    return sum(n for n, op in parse_cigar(cigar) if op in REFERENCE_CONSUMING)
 
 
def coverage_from_cigar(records, length):
    if length <= 0:
        raise ValueError("length must be positive")
    depth = [0] * length
    for record in records:
        pos, cigar = record["pos"], record["cigar"]
        ops = parse_cigar(cigar)
        if not ops:
            continue
        if pos < 1:
            raise ValueError("pos must be 1-based")
        cursor = pos - 1
        span = reference_span(cigar)
        if cursor + span > length:
            raise ValueError("alignment outside reference")
        for n, op in ops:
            if op in "M=X":
                for base in range(cursor, cursor + n):
                    depth[base] += 1
                cursor += n
            elif op in "DN":
                cursor += n
    return depth

Trois points sont critiques. D'abord, parse_cigar refuse les opérations inconnues et les longueurs nulles : un CIGAR silencieusement mal interprété fausse toute la couverture. Ensuite, le curseur de référence avance pour D et N mais aucune profondeur n'est ajoutée, alors que S et I ne déplacent pas ce curseur. Enfin, les positions sont converties une seule fois en indice 0-based avec pos - 1 ; mélanger les conventions dans la boucle est une source classique de décalage.

Côté ligne de commande, la chaîne minimale est la suivante. Elle convertit un SAM en BAM, trie, indexe, puis extrait une région. Le -b demande une sortie BAM ; -h conserve l'en-tête dans une sortie SAM.

samtools view -b datasets/ngs/demo.sam > demo.bam
samtools sort -o demo.sorted.bam demo.bam
samtools index demo.sorted.bam
samtools view -h demo.sorted.bam chr1:10-51
samtools quickcheck -v demo.sorted.bam
samtools flagstat demo.sorted.bam

Pour revenir au texte, samtools view -h demo.sorted.bam > demo.sam. Pour produire un CRAM, samtools view -C -T reference.fasta demo.sorted.bam > demo.cram, puis samtools index demo.cram. La validation ne se limite pas à quickcheck : samtools view détecte les lignes mal formées, samtools flagstat résume les drapeaux, et samtools idxstats vérifie la cohérence entre index et références. Un fichier qui s'ouvre n'est pas forcément un fichier correct.

Étude de cas

Le fichier datasets/ngs/demo.sam est un jeu synthétique de trois alignements sur une référence de 1 000 bases, assez petit pour vérifier chaque valeur à la main.

La question est double : le fichier respecte-t-il les contraintes minimales du format, et que disent ses comptages sur la profondeur et l'appariement ? Nous procédons en quatre étapes : lire l'en-tête, décoder les trois FLAG, vérifier les longueurs de SEQ et QUAL, puis calculer la couverture et les compteurs de type flagstat.

L'en-tête déclare VN:1.6, SO:coordinate et une seule référence chr1 de longueur 1 000. Les deux lectures pair1 forment une paire correcte : FLAG 99 et 147, TLEN 42 et −42, positions 10 et 40. La lecture pair2 a un FLAG 0 : elle n'est pas appariée, mais elle est alignée en position 100 avec un CIGAR 12M. Les qualités de pair2 se terminent par ####, soit quatre bases à $Q = 2$, un signal de fin de lecture dégradée.

Voici la sortie illustrative attendue de samtools flagstat pour ce fichier, reconstituée à partir des drapeaux et des spécifications ; elle n'a pas été produite par SAMtools dans cet environnement.

3 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
3 + 0 mapped (100.00% : N/A)
3 + 0 primary mapped (100.00% : N/A)
2 + 0 paired in sequencing
1 + 0 read1
1 + 0 read2
2 + 0 properly paired (100.00% : N/A)
2 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

L'interprétation est cohérente avec le contenu : trois lectures, toutes alignées et primaires, aucune secondaire, supplémentaire ou dupliquée. Deux lectures sont appariées et correctement appariées, une seule est read1 et une seule read2, et aucune ne finit seule avec un mate non aligné. La couverture de chr1 vaut 36 bases alignées, soit une profondeur moyenne de 0,036× et 3,6 % de bases couvertes. Cette profondeur est trop faible pour un appel de variants ; le fichier sert à valider le parseur et les comptages, pas à illustrer un vrai projet de séquençage. Le tag RG:Z:demo est présent sur chaque lecture, ce qui permet de rattacher les alignements à l'échantillon synthetic déclaré dans @RG.

Erreurs fréquentes

Quiz

Exercices

python -m pytest exercises/module-3/3.3-formats-sam-bam -q

python -m pytest exercises/module-3/3.3-formats-sam-bam -q

python -m pytest exercises/module-3/3.3-formats-sam-bam/projet -q

Résumé

• SAM stocke un alignement par ligne, avec onze champs obligatoires séparés par des tabulations.
• L'en-tête @SQ définit les références et leurs longueurs ; @RG rattache les lectures à un échantillon.
• FLAG est un champ de bits : paires, brin, non-alignement, duplicat, secondaire et supplémentaire se lisent par masques.
• CIGAR distingue ce qui consomme la lecture de ce qui consomme la référence ; S et H ne sont pas équivalents.
• POS est 1-based ; MAPQ mesure la confiance dans l'emplacement, pas la qualité des bases.
• BAM et CRAM sont des encodages compressés de SAM ; un index exige un tri par coordonnées.
• La couverture se calcule en parcourant les CIGAR et en n'ajoutant de la profondeur que pour les bases réellement alignées.
• Une validation explicite des champs, des longueurs et des positions vaut mieux qu'un parseur qui échoue au hasard.

Références

• Li, H., Handsaker, B., Wysoker, A. et al. (2009). The Sequence Alignment/Map format and SAMtools. Bioinformatics, 25(16), 2078–2079.
• Danecek, P., Bonfield, J. K., Liddle, J. et al. (2021). Twelve years of SAMtools and BCFtools. GigaScience, 10(2), giab008.
• SAM/BAM Specification. Sequence Alignment/Map Format Specification. https://samtools.github.io/hts-specs/SAMv1.pdf
• CRAM Specification. CRAM format specification. https://samtools.github.io/hts-specs/CRAMv3.pdf
• SAMtools. Documentation officielle : https://www.hts.org/
• HTSlib. Bibliothèque de référence pour SAM, BAM et CRAM : https://www.htslib.org/