Traitement avec SAMtools

Prérequis

Cette section prolonge la leçon 3.3 (SAM/BAM, FLAG, CIGAR, read groups). Vous devez distinguer une lecture primaire d'un alignement secondaire ou supplémentaire et reconnaître les drapeaux 0x4 (non aligné), 0x100 (secondaire), 0x400 (duplicat) et 0x800 (supplémentaire). Les formats de 1.2 et Biopython de 1.3 servent de rappel ; l'alignement est introduit en 2.1 et développé dans le module 3. Les commandes ciblent SAMtools 1.19+ et n'ont pas été exécutées ici : les sorties sont illustratives.

Introduction

Un fichier d'alignements n'est pas encore un résultat biologique. Pour appeler un variant, mesurer une couverture ou compter les lectures d'un gène, il faut un fichier dont les lectures sont ordonnées, dont les paires se répondent, dont les duplicats sont identifiés et dont les compteurs sont cohérents. Une duplication non traitée gonfle la profondeur locale, fausse la fréquence allélique apparente et peut faire passer un artefact de préparation de bibliothèque pour un vrai signal. À l'inverse, une duplication marquée puis supprimée trop tôt efface une information utile au diagnostic de la bibliothèque.

Le défi computationnel vient de l'échelle et des dépendances entre étapes. Un fichier BAM de séquençage whole-genome contient des centaines de millions d'enregistrements, soit plusieurs dizaines de gigaoctets. Les paires sont dispersées dans le fichier, leurs deux extrémités peuvent être séparées par des mégabases, et le marquage des duplicats doit comparer des lectures voisines sans charger tout le fichier en mémoire. Chaque étape impose donc un ordre précis : regrouper les noms pour reconstituer les paires, compléter les informations de mate, puis trier par coordonnées pour comparer les positions.

SAMtools est la boîte à outils de référence pour ces opérations. Il expose le format, les conversions, le tri, l'indexation, le marquage des duplicats, le filtrage et plusieurs résumés statistiques. Savoir enchaîner ses sous-commandes, c'est savoir transformer un alignement brut en un fichier auditable, reproductible et directement utilisable par BCFtools.

Théorie

Un conteneur, un index

Le format SAM (Sequence Alignment/Map) est un texte tabulaire : un en-tête commence par @HD, @SQ ou @RG, puis chaque ligne décrit une lecture alignée. Le format BAM (Binary Alignment/Map) contient les mêmes informations en binaire, compressées en blocs BGZF (Blocked GNU Zip Format). CRAM est plus compact, mais BAM suffit ici : il est indexable et largement accepté.

L'en-tête fixe l'ordre des références (@SQ), décrit les read groups (@RG) et déclare le tri (@HD SO). Un index .bai ou .csi associe des intervalles de coordonnées à des blocs du fichier. Il permet de lire une région sans parcourir tout le BAM et alimente idxstats. samtools index exige un tri par coordonnées ; -c produit un CSI pour les grands chromosomes.

Un BAM non trié reste lisible, mais n'est ni indexable ni exploitable par les outils qui attendent des coordonnées croissantes. Le tri structure le fichier.

Deux tris, deux contrats

Le tri par nom (name sort) regroupe les enregistrements qui partagent le même QNAME. Les deux extrémités d'une paire deviennent adjacentes, ce qui permet à fixmate de les lire ensemble. samtools sort -n utilise un ordre alphanumérique naturel : read9 précède read10 ; -N impose un ordre lexicographique. samtools collate est une alternative plus rapide lorsque seul le regroupement compte.

Le tri par coordonnées (coordinate sort) range les enregistrements selon l'ordre des @SQ, puis POS, puis le brin. C'est l'ordre attendu par markdup, index et les appels par région. Un fichier trié par nom ne peut pas être indexé : l'index suppose que toutes les lectures d'une région sont contiguës.

Le coût est celui d'un tri externe. Avec $n$ enregistrements et une mémoire de $M$, le tri produit $r = \lceil n/M \rceil$ runs, fusionnés en $\lceil \log_k r \rceil$ passes. Le temps est en $O(n \log n)$, le disque temporaire en $O(n)$. -T choisit le préfixe des runs, -m la mémoire par thread.

fixmate : rendre les paires utilisables

fixmate remplit les champs du mate : RNEXT, PNEXT, TLEN, les drapeaux MREVERSE et MUNMAP, et les tags MC (mate CIGAR) et MQ (mate mapping quality). Il vérifie l'orientation et peut retirer PROPER_PAIR lorsqu'il est manifestement faux. Son entrée doit être groupée par nom : il compare deux enregistrements consécutifs et n'utilise qu'un petit tampon.

L'option -m ajoute ms, le score du mate : la somme des qualités de base du mate, en ne comptant que les bases de qualité au moins 15. markdup s'en sert pour comparer les deux extrémités et conserver le meilleur couple. Sans -m, il s'arrête avec no ms score tag. L'option -r retire les lectures secondaires et non alignées, ce qui allège le fichier mais compromet l'audit.

fixmate refuse un fichier déclaré trié par coordonnées : les mates n'y sont pas adjacents, et l'outil ne peut pas les apparier sans recharger tout le fichier. Il faut donc trier par nom d'abord.

markdup : même extrémité, même orientation

samtools markdup identifie les duplicats dans un fichier trié par coordonnées et passé par fixmate -m. Il utilise les positions non rognées (unclipped) des extrémités 5', la référence, l'orientation et le caractère gauche ou droit de la lecture. Pour une paire, il combine les deux extrémités à l'aide de MC ; pour une lecture seule, il utilise sa position 5' et son orientation. La lecture de meilleur score de qualité, complétée par ms, est conservée ; les autres reçoivent 0x400.

Le marquage ne compare pas les séquences complètes : deux fragments différents aux mêmes extrémités sont marqués comme duplicats, et deux lectures décalées d'une base ne le sont pas, sauf si le rognage (clipping) aligne leurs extrémités. C'est une définition positionnelle rapide.

Par défaut, markdup marque sans supprimer : les lectures restent dans le BAM avec 0x400, ce qui préserve la traçabilité. -r les retire, -s imprime un résumé, -f l'écrit dans un fichier, -l borne la longueur attendue, -T choisit les temporaires et -S marque aussi les supplémentaires des duplicats. Les duplicats optiques sont traités avec -d, discuté plus bas.

Duplicats PCR et duplicats optiques

Un duplicat PCR (polymerase chain reaction) provient de la même molécule d'ADN amplifiée plusieurs fois pendant la préparation de la bibliothèque. Les copies partagent extrémités et orientation ; elles ne sont pas indépendantes. Un duplicat optique est un artefact de lecture : deux clusters distincts, très proches sur la même tuile, se contaminent optiquement. Les molécules peuvent différer, mais leurs coordonnées de lecture sont voisines.

markdup traite d'abord les duplicats positionnels, puis classe les optiques lorsque -d est fourni. La distance utilise les coordonnées x et y du nom Illumina : deux lectures sont optiques si $\max(|\Delta x|, |\Delta y|) \le D$. La documentation suggère environ 100 pour HiSeq et 2500 pour NovaSeq. Les duplicats optiques reçoivent dt:Z:SQ, les autres dt:Z:LB. Dans SAMtools 1.19, -d vaut 0 par défaut : sans cette option, aucun duplicat optique n'est identifié.

La distinction guide l'action : un taux élevé de duplicats PCR signale une bibliothèque peu complexe ou une amplification excessive, tandis qu'un taux élevé de duplicats optiques signale un problème de lecture. Le premier appelle plus de matériel ou moins de cycles, le second un réglage de l'instrument.

Lire les compteurs sans se tromper

samtools flagstat compte des enregistrements par drapeau. La première ligne donne le total QC-passed et QC-failed ; primary, secondary et supplementary le partitionnent. duplicates compte tous les 0x400, primary duplicates le sous-ensemble primaire, et mapped les enregistrements sans 0x4. Les lignes paired in sequencing, read1, read2, properly paired, with itself and mate mapped et singletons décrivent les paires.

Les pourcentages n'ont pas tous le même dénominateur : mapped est rapporté au total, properly paired et singletons aux paires. Un taux de 100 % de paires correctes ne signifie pas que toutes les lectures le sont, et flagstat compte des lectures : une paire correcte contribue deux enregistrements.

samtools idxstats lit l'index et affiche référence, longueur, segments alignés et non alignés ; la ligne * regroupe les lectures sans coordonnée. Les comptes portent sur des segments, secondaires et supplémentaires inclus, et une lecture peut être comptée plusieurs fois. Il sert à la répartition par référence, pas à estimer des molécules uniques. Sans index, il parcourt tout le fichier.

samtools stats produit des lignes SN : séquences, lectures alignées, paires correctes, duplicats, longueur et qualité moyennes, insert, erreurs. Avec -r reference.fa, il ajoute mésappariements et taux d'erreur. flagstat ne lit que les drapeaux ; stats descend aux bases.

samtools view filtre et convertit : -b, -h, -c, -q, -f (tous les bits), -F (au moins un bit). samtools quickcheck vérifie l'en-tête et le bloc de fin d'un BAM ou CRAM, pas chaque enregistrement. Silencieux signifie « structure plausible », pas « fichier correct ».

Files d'attente, mémoire et coût

Les trois étapes ont des profils mémoire différents. fixmate ne garde que deux enregistrements : son coût est négligeable, mais il exige l'adjacence des mates. Le tri est le consommateur principal : -m fixe la mémoire par thread et -@ le nombre de threads, si bien que la consommation approximative est leur produit, augmentée des tampons. Avec -@ 8 -m 1G, prévoyez plusieurs gigaoctets et de l'espace disque pour les runs temporaires.

markdup entretient une table de hachage des signatures et une file des lectures de la fenêtre courante. La mémoire croît avec la densité de duplicats et le nombre de lectures qui partagent une position ; un locus très couvert coûte donc plus cher. -l borne la longueur attendue et -T déplace les temporaires.

Le coût total est dominé par les entrées-sorties et les tris. Chaque passage réécrit le BAM, et la compression BGZF coûte du temps CPU. Augmenter -@ accélère compression et tri, mais sature mémoire et disques. Mesurez sur un échantillon.

Exemple travaillé

Prenons une bibliothèque de $L = 10^{7}$ fragments uniques et $n = 5 \times 10^{7}$ lectures séquencées, avec un tirage supposé uniforme avec remise. Le nombre attendu de fragments distincts observés est

La fraction de lectures dupliquées attendue vaut alors

soit environ 80 %. Ce modèle ignore l'amplification PCR et les biais de bibliothèque, mais montre l'essentiel : au-delà d'une certaine profondeur, séquencer davantage ajoute surtout des copies. Le taux observé doit donc être interprété avec la complexité de la bibliothèque, pas seulement avec la profondeur.

Pour le marquage, considérons deux paires sur chr1 : la première a son extrémité gauche en position 100 et son extrémité droite en position 300, en orientation FR ; la seconde a exactement les mêmes extrémités et la même orientation. Elles partagent la signature de markdup : l'une est conservée, l'autre marquée. Une troisième paire à 101 et 300 a une signature différente et n'est pas un duplicat, même si elle se chevauche largement. Si le rognage ramène sa position non rognée à 100, elle rejoint la signature et devient candidate.

Implémentation

Les commandes ciblent SAMtools 1.19+ et n'ont pas été exécutées ici ; les sorties sont illustratives. La chaîne part de demo.bam.

set -euo pipefail
 
samtools sort -n -o demo.namesort.bam demo.bam
samtools fixmate -m demo.namesort.bam demo.fixmate.bam
samtools sort -o demo.coordsort.bam demo.fixmate.bam
samtools markdup -s demo.coordsort.bam demo.markdup.bam
samtools index demo.markdup.bam
samtools flagstat demo.markdup.bam
samtools idxstats demo.markdup.bam

Le premier tri rapproche les mates ; fixmate -m remplit MC, MQ et ms ; le second tri range les positions ; markdup -s marque et affiche ses comptes ; index crée demo.markdup.bam.bai ; les deux dernières commandes résument. Si l'entrée est déjà triée par nom, le premier tri peut être omis, jamais fixmate.

Read group : une métadonnée à ne pas perdre

Le read group identifie l'échantillon, la bibliothèque et l'unité de séquençage ; attachez-le dès l'alignement. Avec BWA-MEM, on passe une chaîne @RG dont les tabulations sont échappées :

bwa mem -R '@RG\tID:demo\tSM:synthetic\tPL:ILLUMINA\tLB:lib1\tPU:unit1' \
  reference.fa reads_R1.fastq reads_R2.fastq > demo.sam
 
samtools sort -n -o demo.namesort.bam demo.sam
samtools fixmate -m demo.namesort.bam demo.fixmate.bam
samtools sort -o demo.coordsort.bam demo.fixmate.bam
samtools markdup -s demo.coordsort.bam demo.markdup.bam
samtools index demo.markdup.bam

ID doit être unique, SM désigne l'échantillon, LB la bibliothèque et PU l'unité de séquençage. Un ID dupliqué ou un SM absent complique l'appel de variants et l'interprétation des duplicats. Chaque commande SAMtools ajoute une ligne @PG avec sa version et sa ligne de commande : c'est une partie de la provenance.

Automatiser la chaîne en Python

La fonction suivante construit la liste ordonnée des commandes sans les exécuter : la chaîne reste testable avant tout lancement.

def markdup_pipeline(prefix):
    if not prefix or any(character.isspace() for character in prefix):
        raise ValueError("prefix must be non-empty and contain no whitespace")
    return [
        f"samtools sort -n -o {prefix}.namesort.bam {prefix}.bam",
        f"samtools fixmate -m {prefix}.namesort.bam {prefix}.fixmate.bam",
        f"samtools sort -o {prefix}.coordsort.bam {prefix}.fixmate.bam",
        f"samtools markdup -s {prefix}.coordsort.bam {prefix}.markdup.bam",
        f"samtools index {prefix}.markdup.bam",
        f"samtools flagstat {prefix}.markdup.bam",
        f"samtools idxstats {prefix}.markdup.bam",
    ]
 
 
for command in markdup_pipeline("demo"):
    print(command)

Sortie illustrative :

samtools sort -n -o demo.namesort.bam demo.bam
samtools fixmate -m demo.namesort.bam demo.fixmate.bam
samtools sort -o demo.coordsort.bam demo.fixmate.bam
samtools markdup -s demo.coordsort.bam demo.markdup.bam
samtools index demo.markdup.bam
samtools flagstat demo.markdup.bam
samtools idxstats demo.markdup.bam

Dans un pipeline réel, ajoutez -@ aux commandes qui le supportent, -m au tri et -T aux temporaires. Vérifiez chaque code de retour : un samtools interrompu peut laisser un BAM partiel trompeur.

Lire flagstat et idxstats

Un parseur de flagstat doit séparer les comptes QC-passed et QC-failed, retirer le pourcentage entre parenthèses et refuser une catégorie inconnue. Le code ci-dessous conserve les catégories réellement présentes, ce qui accepte les anciennes sorties dépourvues de primary.

import re
 
LABELS = {
    "primary": "primary",
    "secondary": "secondary",
    "supplementary": "supplementary",
    "duplicates": "duplicates",
    "primary duplicates": "primary_duplicates",
    "mapped": "mapped",
    "primary mapped": "primary_mapped",
    "paired in sequencing": "paired",
    "read1": "read1",
    "read2": "read2",
    "properly paired": "properly_paired",
    "with itself and mate mapped": "mate_mapped",
    "singletons": "singletons",
    "with mate mapped to a different chr": "different_chr",
    "with mate mapped to a different chr (mapQ>=5)": "different_chr_mapq5",
}
 
COUNT = re.compile(r"^(\d+)\s*\+\s*(\d+)\s+(.+)$")
PERCENT = re.compile(
    r"\s*\((N/A|\d+(?:\.\d+)?%?)\s*:\s*(N/A|\d+(?:\.\d+)?%?)\)\s*$"
)
 
 
def parse_flagstat(text):
    passed, failed, percentages = {}, {}, {}
    for line_number, raw_line in enumerate(text.splitlines(), start=1):
        line = raw_line.strip("\r\n")
        if not line.strip() or line.lstrip().startswith("#"):
            continue
        match = COUNT.match(line)
        if not match:
            raise ValueError(f"line {line_number}: invalid flagstat line")
        label = match.group(3).strip()
        percent = PERCENT.search(label)
        values = None
        if percent:
            label = label[:percent.start()].strip()
            values = (
                None if percent.group(1) == "N/A" else float(percent.group(1).rstrip("%")),
                None if percent.group(2) == "N/A" else float(percent.group(2).rstrip("%")),
            )
        key = "total" if label.startswith("in total") else LABELS.get(label)
        if key is None:
            raise ValueError(f"line {line_number}: unknown category {label!r}")
        if key in passed:
            raise ValueError(f"line {line_number}: duplicate category")
        passed[key] = int(match.group(1))
        failed[key] = int(match.group(2))
        if values is not None:
            percentages[key] = {"passed": values[0], "failed": values[1]}
    if "total" not in passed:
        raise ValueError("missing total line")
    return {"passed": passed, "failed": failed, "percentages": percentages}

Un parseur d'idxstats est plus simple : quatre colonnes séparées par des tabulations, des entiers non négatifs, une longueur strictement positive pour les références et nulle pour la ligne *.

def parse_idxstats(text):
    rows, seen = [], set()
    for line_number, raw_line in enumerate(text.splitlines(), start=1):
        line = raw_line.strip("\r\n")
        if not line.strip() or line.lstrip().startswith("#"):
            continue
        fields = line.split("\t")
        if len(fields) != 4:
            raise ValueError(f"line {line_number}: expected 4 columns")
        name = fields[0].strip()
        if not name or name in seen:
            raise ValueError(f"line {line_number}: invalid reference name")
        try:
            length, mapped, unmapped = (int(field) for field in fields[1:])
        except ValueError as exc:
            raise ValueError(f"line {line_number}: invalid number") from exc
        if min(length, mapped, unmapped) < 0:
            raise ValueError(f"line {line_number}: negative count")
        if name == "*" and length != 0:
            raise ValueError(f"line {line_number}: '*' must have length 0")
        if name != "*" and length <= 0:
            raise ValueError(f"line {line_number}: length must be positive")
        seen.add(name)
        rows.append({
            "name": name,
            "length": length,
            "mapped": mapped,
            "unmapped": unmapped,
        })
    if not rows:
        raise ValueError("no idxstats rows")
    return rows

Ces fonctions n'exécutent aucun binaire : elles testent le contrat de sortie. Un contrôle utile vérifie total = primary + secondary + supplementary, mapped <= total et properly_paired <= paired. Un écart entre les mapped d'idxstats et de flagstat mérite un avertissement, pas une erreur : alignements secondaires, supplémentaires et index peuvent l'expliquer.

Filtrer, vérifier, résumer

Quelques commandes de contrôle, toutes illustratives :

samtools quickcheck -v demo.markdup.bam
samtools view -c -f 0x2 demo.markdup.bam
samtools view -c -F 0x404 demo.markdup.bam
samtools view -b -q 30 -F 0x404 -o demo.filtered.bam demo.markdup.bam
samtools view -h demo.markdup.bam chr1:1-100
samtools stats -r reference.fa demo.markdup.bam > demo.stats.txt
grep '^SN' demo.stats.txt | head -n 12

quickcheck ne produit rien si tout va bien et renvoie un code non nul sinon. view -c -f 0x2 compte les paires correctes ; -F 0x404 exclut non alignés et duplicats. La requête par région exige l'index ; stats -r exige une référence indexée par samtools faidx.

Reproductibilité

Pour une chaîne reproductible, conservez la version (samtools --version), la commande complète, le read group, la référence et son index, les options de tri et de mémoire, et les journaux de markdup -s, flagstat, idxstats et stats. Les lignes @PG documentent déjà une partie de l'historique ; ne les supprimez pas avec --no-PG sans archiver les commandes.

Les comptes doivent être reproductibles ; l'identité octet pour octet du BAM ne l'est pas forcément. Compression, threads, ordre des enregistrements de même clé et temporaires changent les octets sans changer les alignements. Comparez statistiques et checksums d'entrée. Épinglez SAMtools 1.19+ et relisez le manuel : markdup et flagstat ont évolué.

Étude de cas

Le fichier datasets/ngs/demo.sam est un exemple synthétique. Son en-tête déclare @HD VN:1.6 SO:coordinate, une référence chr1 de 1000 pb et un read group ID:demo SM:synthetic PL:ILLUMINA. Il contient trois enregistrements : la paire pair1 (FLAG 99 et 147) alignée en 12M aux positions 10 et 40, avec un insert de 42 pb, et la lecture seule pair2 (FLAG 0) alignée en 12M à la position 100, de MAPQ 30, dont quatre bases ont une qualité très faible.

Question : la chaîne de marquage préserve-t-elle les comptes, et que disent les résumés sur cette mini-bibliothèque ? Les étapes sont les suivantes.

samtools sort -n -o demo.namesort.bam datasets/ngs/demo.sam
samtools fixmate -m demo.namesort.bam demo.fixmate.bam
samtools sort -o demo.coordsort.bam demo.fixmate.bam
samtools markdup -s demo.coordsort.bam demo.markdup.bam
samtools index demo.markdup.bam
samtools quickcheck -v demo.markdup.bam
samtools flagstat demo.markdup.bam
samtools idxstats demo.markdup.bam

Après fixmate -m, la table suivante (sortie illustrative, déduite des trois enregistrements) montre les tags ajoutés ; les qualités I valent Q40, donc ms vaut 480.

markdup -s examine trois lectures, n'en écrit aucune comme duplicat et signale zéro duplicat. Le résumé illustratif suivant ne montre que les lignes utiles :

READ: 3
WRITTEN: 3
EXAMINED: 3
PAIRED: 2
SINGLE: 1
DUPLICATE PRIMARY TOTAL: 0
DUPLICATE TOTAL: 0

flagstat doit retrouver trois lectures primaires, deux lectures appariées, deux correctement appariées et aucun duplicat. Sortie illustrative :

3 + 0 in total (QC-passed reads + QC-failed reads)
3 + 0 primary
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)

idxstats doit afficher une seule référence et trois segments alignés, puis la ligne * à zéro. Sortie illustrative :

chr1	1000	3	0
*	0	0	0

L'absence de duplicat ne prouve pas que markdup fonctionne : avec une seule paire à une position unique, il n'y a rien à marquer. Pour tester le marquage, dupliquez pair1 sous un nouveau QNAME avec des qualités plus faibles : les deux copies doivent porter 0x400, flagstat doit compter deux duplicats et markdup -s deux lectures dupliquées appariées.

pair2 illustre un autre point : une lecture seule de MAPQ 30 avec quatre bases de qualité 2 reste alignée mais contribue peu à un appel de variant ; les compteurs globaux masquent cette hétérogénéité.

Erreurs fréquentes

Quiz

Exercices

python -m pytest exercises/module-3/3.4-traitement-samtools -q

python -m pytest exercises/module-3/3.4-traitement-samtools -q

python -m pytest exercises/module-3/3.4-traitement-samtools/projet -q

Résumé

• Le BAM est un conteneur binaire compressé ; l'index .bai ou .csi permet les requêtes par région et idxstats, et exige un tri par coordonnées.
• sort -n rapproche les mates, fixmate -m ajoute MC, MQ et ms, puis sort range les positions.
• markdup exige coordonnées et tags de fixmate ; il marque 0x400 par défaut et ne supprime qu'avec -r.
• La signature d'un duplicat repose sur les extrémités 5' non rognées et l'orientation, pas sur l'identité des séquences.
• Les duplicats PCR viennent de l'amplification ; les optiques de clusters voisins, détectés avec -d et les coordonnées du nom.
• flagstat compte des enregistrements et sépare QC-passed et QC-failed ; ses pourcentages ont des dénominateurs différents.
• idxstats décrit la répartition par référence, stats descend aux bases et aux inserts, view filtre et quickcheck vérifie l'en-tête et le bloc de fin.

Références

• Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., Marth, G., Abecasis, G. & Durbin, R. (2009). The Sequence Alignment/Map format and SAMtools. Bioinformatics, 25(16), 2078–2079.
• Danecek, P., Bonfield, J. K., Liddle, J., Marshall, J., Ohan, V., Pollard, M. O., Whitwham, A., Keane, T., McCarthy, S. A., Davies, R. M. & Li, H. (2021). Twelve years of SAMtools and BCFtools. GigaScience, 10(2), giab008.
• Lander, E. S. & Waterman, M. S. (1988). Genomic mapping by fingerprinting random clones: a mathematical analysis. Genomics, 2(3), 231–239.
• SAMtools. Documentation des sous-commandes. https://www.htslib.org/doc/samtools.html
• SAMtools. SAM/BAM Format Specification Working Group. https://samtools.github.io/hts-specs/SAMv1.pdf