Pipeline de variants

Prérequis

Cette section assemble les briques du module 3 : contrôle qualité avec 3.1, mapping avec 3.2, traitement SAM/BAM avec 3.4 et appel de variants avec 3.5. Elle suppose aussi le cadrage de la section 7.1 et les principes de reproductibilité de 0.1. Les exercices n'utilisent que la bibliothèque standard de Python ; les commandes d'outils sont fournies à titre illustratif.

Introduction

Un variant intéressant ne sort jamais d'un seul programme. Il faut contrôler les lectures, les aligner, corriger les groupes de lectures, marquer les duplicats, appeler les sites variables, normaliser les indels, filtrer, puis rassembler des statistiques. Chaque étape produit un fichier que la suivante consomme. Si une étape est relancée sans que les suivantes le sachent, le résultat final mélange deux versions des données. Si une étape échoue silencieusement, le rapport peut décrire un génome qui n'a jamais été analysé.

L'enjeu biologique est direct. Chez Escherichia coli K-12 MG1655, une substitution dans gyrA peut modifier la sensibilité aux fluoroquinolones ; dans un projet clinique, un variant somatique de faible fréquence peut changer la décision thérapeutique. La question n'est pas seulement « quel variant appeler ? », mais « ce variant survit-il à une chaîne de traitement complète, versionnée et reproductible ? ». Un appel isolé, sans contrôle de couverture ni de qualité de mapping, n'est pas une preuve.

L'enjeu computationnel est l'orchestration. Un script monolithique fonctionne sur un échantillon, puis devient difficile à reprendre, à paralléliser et à auditer. Les gestionnaires de workflow transforment les commandes en graphe de dépendances : chaque règle déclare ses entrées, ses sorties et sa commande. Le moteur calcule l'ordre, parallélise les branches indépendantes, conserve les sorties valides et signale les fichiers manquants. Snakemake et Nextflow diffèrent par la syntaxe, pas par la logique de graphe.

Nous construisons un pipeline QC, mapping, traitement BAM, appel, filtrage et rapport, puis nous lisons journaux et checksums. L'étude de cas utilise datasets/ngs/demo.sam, un alignement synthétique de trois lectures, sans prétendre décrire un variant réel.

Théorie

Du script au graphe

Un pipeline se représente par un graphe orienté $G=(V,E)$. Chaque sommet $v \in V$ est une règle ou une tâche ; chaque arc $(u,v) \in E$ signifie que la sortie de $u$ est une entrée de $v$. L'arc va du producteur vers le consommateur. Un pipeline valide est un graphe acyclique orienté, ou DAG (directed acyclic graph) : si un cycle existait, aucune règle du cycle ne pourrait démarrer, puisque chacune attendrait une sortie produite plus tard.

Un ordre topologique est une permutation des sommets telle que $u$ précède $v$ pour tout arc $(u,v)$. Il n'est pas unique : deux branches indépendantes peuvent avancer dans n'importe quel ordre. Le calcul supprime répétitivement les sommets de degré entrant nul ; avec une file et des listes d'adjacence, la complexité est $O(|V|+|E|)$. Un tri alphabétique des sommets prêts rend le résultat déterministe.

Une chaîne QC, mapping, marquage, appel, filtrage, rapport a des degrés entrants 0, 1, 1, 1, 1 et 1 : l'ordre est unique. Si le rapport attend aussi les statistiques de mapping, deux branches avancent en parallèle et l'ordre n'est plus unique.

Un DAG décrit les dépendances de fichiers, pas les dépendances scientifiques. Deux règles peuvent être correctement ordonnées et produire un résultat faux si la référence, les paramètres ou le modèle sont inadaptés. Le graphe garantit l'enchaînement, pas la validité biologique.

Anatomie d'une règle Snakemake

Snakemake étend Python avec des règles déclaratives. Une règle nomme ses entrées, ses sorties, son journal, ses ressources et la commande shell qui produit les sorties. Le moteur part de la cible demandée, remonte les dépendances jusqu'aux fichiers d'entrée, puis exécute les jobs dans l'ordre requis.

rule marquage_duplicats:
    input:
        "aln/{sample}.positionsorted.bam"
    output:
        "aln/{sample}.dedup.bam"
    log:
        "logs/markdup/{sample}.log"
    threads: 4
    shell:
        "samtools markdup -s {input} {output} 2> {log}"

La sortie est le contrat. Un fichier absent déclenche la règle ; un fichier présent et plus récent que ses entrées peut être conservé. Les entrées nommées rendent le shell lisible : {input.bam} désigne un fichier précis, {input} la liste complète. Les accolades sont interprétées par Snakemake, pas par le shell.

log ne redirige rien à lui seul. Il déclare un fichier que la commande doit écrire, par exemple avec 2> {log}. Un journal déclaré mais jamais écrit peut faire échouer le job ou laisser un fichier vide. threads réserve des cœurs et alimente {threads} ; resources déclare la mémoire ou la durée sans unité automatique.

La règle all est une cible, pas un job. Elle déclare souvent les entrées finales ; le moteur vérifie qu'elles existent sans exécuter de commande. Un journal peut donc être complet sans ligne rule all.

Wildcards, expansion et collisions

Une wildcard est une variable de nom de fichier encadrée d'accolades. Dans qc/{sample}_R1.fastq, {sample} capture la partie du nom qui correspond au motif. Snakemake instancie une règle par valeur nécessaire : deux échantillons produisent deux jobs, pas une boucle implicite. expand("results/{sample}.vcf.gz", sample=SAMPLES) construit explicitement une liste de cibles.

Les wildcards ne sont pas des expressions régulières. Le moteur compare les motifs de sortie et d'entrée, puis vérifie qu'une valeur unique satisfait toutes les contraintes. Si {sample} peut capturer A_1 ou A, deux instanciations sont possibles ; un fichier supplémentaire ou une contrainte lève l'ambiguïté. wildcard_constraints: sample="[^/]+" évite qu'un slash entre dans un nom.

Les collisions de sorties sont plus graves qu'une ambiguïté. Deux règles qui écrivent le même chemin se contredisent : le moteur peut en choisir une ou produire un fichier partiel. Les underscores, points et tirets dans les noms d'échantillons causent de mauvaises captures. Fixez une convention, testez snakemake -n --reason et vérifiez que chaque sortie a un producteur unique.

Nextflow exprime la même idée avec des canaux (channels). Un canal transporte des valeurs et des fichiers ; un processus est instancié pour chaque élément. Channel.fromFilePairs("data/*_{R1,R2}.fastq") crée des tuples échantillon, lecture 1, lecture 2. Snakemake raisonne en noms de fichiers et en motifs, Nextflow en flux de données et en processus.

Reprise, cache et checksums

Un moteur de workflow ne rejoue pas tout à chaque appel. Snakemake compare les dates de modification et d'autres déclencheurs ; Nextflow calcule une empreinte de tâche à partir des entrées, du script, des paramètres et de l'environnement, puis réutilise le résultat stocké dans work/. Cette reprise économise des heures de calcul si les dépendances sont correctement déclarées.

Une dépendance cachée casse la reprise. Si une règle lit un fichier de configuration sans le déclarer en input ou en params, le moteur ne voit pas le changement et conserve une sortie périmée. Si une commande écrit un fichier annexe non déclaré, ce fichier peut manquer sur une autre machine. Tout ce qui influence la sortie doit être une entrée, un paramètre ou une version d'environnement.

Les checksums complètent les dates. Un manifeste au format sha256sum associe un condensat de 64 caractères hexadécimaux à un chemin :

9f2c...  data/demo.sam

La vérification recalcule le condensat et le compare à la valeur attendue. Elle détecte une modification, une troncature ou un fichier remplacé, mais ne prouve ni l'authenticité, ni la provenance, ni la validité biologique. Un fichier faux peut avoir un checksum cohérent.

snakemake --rerun-incomplete nettoie les sorties incomplètes ; nextflow run -resume réutilise les tâches terminées dont l'empreinte correspond. Aucun des deux ne répare un résultat scientifiquement faux ni ne reprend un processus au milieu d'une commande. Ils restaurent un état de fichiers, pas une intention.

Versions, environnements et conteneurs

Un pipeline reproductible fixe trois couches : le code, les dépendances logicielles et les données. Le code se versionne avec Git ; les dépendances se figent dans un environnement conda, un verrou ou une image ; les données s'identifient par une accession, une date et un checksum.

Snakemake accepte conda: "workflow/envs/ngs.yaml" par règle et --use-conda pour créer les environnements, ainsi que container: pour exécuter une règle dans une image. Nextflow utilise container dans chaque processus avec Docker, Podman ou Singularity. Une image doit être épinglée par un tag stable et, mieux, par un digest : image@sha256:... ne change pas, contrairement à latest.

Les versions de BWA-MEM 0.7.17, SAMtools 1.19+ et BCFtools 1.19+ doivent apparaître dans le journal ou un fichier de versions. samtools markdup a remplacé rmdup et exige une chaîne précise ; bcftools norm a évolué dans sa gestion des multiallélices. Mélanger des binaires système non épinglés et des conteneurs rend le pipeline irreproductible, même à Snakefile identique.

Parallélisation, erreurs et journalisation

Le moteur parallélise les branches indépendantes du DAG. Avec $p$ travailleurs, le gain brut s'écrit $S_p = T_1/T_p$ et l'efficacité $E_p = S_p/p$. Un pipeline dominé par une étape non parallélisable gagne peu à multiplier les cœurs ; les ressources déclarées évitent que dix jobs de mapping se disputent la mémoire.

Une erreur doit interrompre la branche concernée, pas produire un fichier partiel. En shell, set -euo pipefail fait échouer un pipe si une commande intermédiaire échoue. Sans cette option, bwa mem | samtools sort peut renvoyer le code de samtools sort et masquer un échec de bwa. Snakemake signale le job en échec ; --keep-going poursuit les branches indépendantes et --retries relance un job. Nextflow propose errorStrategy 'retry' et maxRetries.

Un bon journal contient l'horodatage, le nom de la règle, le statut et l'erreur. Il ne contient ni secret ni chemin absolu inutile. Séparez stdout et stderr : un rapport fastp sur stdout n'est pas un journal d'erreur. Un fichier par échantillon et par règle facilite la reprise ciblée.

Exemple travaillé

Considérons six étapes : QC, mapping, marquage des duplicats, appel, filtrage, rapport. Les dépendances sont QC vers mapping, mapping vers marquage, marquage vers appel, appel vers filtrage, filtrage vers rapport. Les degrés entrants sont 0, 1, 1, 1, 1, 1. L'ordre topologique est unique : QC, mapping, marquage, appel, filtrage, rapport. Si le rapport dépend aussi des statistiques de mapping, la branche statistiques part de mapping et le rapport attend les deux branches ; l'ordre n'est plus unique, mais le rapport reste dernier.

Sur datasets/ngs/demo.sam, trois alignements de 12 pb couvrent un contig de 1 000 pb. La profondeur moyenne vaut $(3 \times 12)/1000 = 0{,}036$, soit 0,036 lecture par base. Ce calcul, réalisé sur un fichier synthétique, montre que le jeu ne permet aucun appel biologique sérieux : il sert à vérifier la mécanique des fichiers et des dépendances, pas à découvrir un variant.

Implémentation

Le Snakefile ci-dessous décrit la chaîne complète. Les commandes lourdes ne sont pas exécutées dans cet environnement ; les sorties sont étiquetées « illustratives » et les journaux « illustratifs ». Les versions cibles sont Snakemake 8.x, BWA-MEM 0.7.17, SAMtools 1.19+ et BCFtools 1.19+.

configfile: "config/config.yaml"
 
SAMPLES = config["samples"]
REF = "reference/ecoli.fa"
 
shell.prefix("set -euo pipefail; ")
 
# Cible finale : un rapport HTML par échantillon.
rule all:
    input:
        expand("results/{sample}.report.html", sample=SAMPLES)
 
# Index FASTA et BWA partagés par tous les jobs.
rule index_reference:
    input:
        REF
    output:
        "reference/ecoli.fa.fai",
        "reference/ecoli.fa.amb",
        "reference/ecoli.fa.ann",
        "reference/ecoli.fa.bwt",
        "reference/ecoli.fa.pac",
        "reference/ecoli.fa.sa"
    log:
        "logs/index_reference.log"
    conda:
        "workflow/envs/ngs.yaml"
    shell:
        "samtools faidx {input} && bwa index {input} > {log} 2>&1"
 
# Contrôle qualité, trimming et rapports fastp.
rule fastp:
    input:
        r1="data/{sample}_R1.fastq",
        r2="data/{sample}_R2.fastq"
    output:
        r1="qc/{sample}_R1.trim.fastq",
        r2="qc/{sample}_R2.trim.fastq",
        json="qc/{sample}.fastp.json",
        html="qc/{sample}.fastp.html"
    log:
        "logs/fastp/{sample}.log"
    threads: 4
    conda:
        "workflow/envs/ngs.yaml"
    shell:
        "fastp -i {input.r1} -I {input.r2} -o {output.r1} -O {output.r2} "
        "--detect_adapter_for_pe --qualified_quality_phred 20 "
        "--length_required 50 --thread {threads} "
        "--json {output.json} --html {output.html} > /dev/null 2> {log}"
 
# Mapping BWA-MEM puis tri par nom, requis par fixmate.
rule bwa_mem:
    input:
        r1="qc/{sample}_R1.trim.fastq",
        r2="qc/{sample}_R2.trim.fastq",
        fai="reference/ecoli.fa.fai",
        bwt="reference/ecoli.fa.bwt"
    output:
        "aln/{sample}.namesorted.bam"
    log:
        "logs/bwa/{sample}.log"
    threads: 8
    conda:
        "workflow/envs/ngs.yaml"
    shell:
        "bwa mem -t {threads} -R '@RG\\tID:{sample}\\tSM:{sample}\\tPL:ILLUMINA' "
        "reference/ecoli.fa {input.r1} {input.r2} 2> {log} "
        "| samtools sort -n -o {output} -"
 
# Chaîne imposée : fixmate -m, tri par coordonnées, markdup, index.
rule fixmate:
    input:
        "aln/{sample}.namesorted.bam"
    output:
        "aln/{sample}.fixmate.bam"
    log:
        "logs/fixmate/{sample}.log"
    conda:
        "workflow/envs/ngs.yaml"
    shell:
        "samtools fixmate -m {input} {output} 2> {log}"
 
rule sort_fixmate:
    input:
        "aln/{sample}.fixmate.bam"
    output:
        "aln/{sample}.positionsorted.bam"
    log:
        "logs/sort/{sample}.log"
    threads: 4
    conda:
        "workflow/envs/ngs.yaml"
    shell:
        "samtools sort -@ {threads} -o {output} {input} 2> {log}"
 
rule markdup:
    input:
        "aln/{sample}.positionsorted.bam"
    output:
        "aln/{sample}.dedup.bam"
    log:
        "logs/markdup/{sample}.log"
    conda:
        "workflow/envs/ngs.yaml"
    shell:
        "samtools markdup -s {input} {output} 2> {log}"
 
rule index_bam:
    input:
        "aln/{sample}.dedup.bam"
    output:
        "aln/{sample}.dedup.bam.bai"
    log:
        "logs/index/{sample}.log"
    conda:
        "workflow/envs/ngs.yaml"
    shell:
        "samtools index {input} 2> {log}"
 
rule flagstat:
    input:
        bam="aln/{sample}.dedup.bam",
        bai="aln/{sample}.dedup.bam.bai"
    output:
        "qc/{sample}.flagstat"
    log:
        "logs/flagstat/{sample}.log"
    conda:
        "workflow/envs/ngs.yaml"
    shell:
        "samtools flagstat {input.bam} > {output} 2> {log}"
 
# Appel brut : mpileup alimente call dans un pipe BCF.
rule call:
    input:
        bam="aln/{sample}.dedup.bam",
        bai="aln/{sample}.dedup.bam.bai",
        fai="reference/ecoli.fa.fai",
        bwt="reference/ecoli.fa.bwt"
    output:
        "variants/{sample}.raw.vcf.gz"
    log:
        "logs/call/{sample}.log"
    conda:
        "workflow/envs/ngs.yaml"
    shell:
        "bcftools mpileup -f reference/ecoli.fa -q 20 -Q 20 -d 500 -Ou "
        "{input.bam} 2> {log} | bcftools call -mv -Oz -o {output}"
 
rule normalize:
    input:
        "variants/{sample}.raw.vcf.gz"
    output:
        "variants/{sample}.norm.vcf.gz"
    log:
        "logs/normalize/{sample}.log"
    conda:
        "workflow/envs/ngs.yaml"
    shell:
        "bcftools norm -f reference/ecoli.fa -m -any -Oz -o {output} {input} 2> {log}"
 
rule filter:
    input:
        "variants/{sample}.norm.vcf.gz"
    output:
        "variants/{sample}.filtered.vcf.gz"
    log:
        "logs/filter/{sample}.log"
    conda:
        "workflow/envs/ngs.yaml"
    shell:
        "bcftools view -f PASS -i 'QUAL>=20 && INFO/DP>=10' "
        "-Oz -o {output} {input} 2> {log}"
 
# Agrégation des sorties : le rapport ne recalcule pas les variants.
rule report:
    input:
        vcf="variants/{sample}.filtered.vcf.gz",
        stats="qc/{sample}.flagstat",
        qc="qc/{sample}.fastp.json"
    output:
        "results/{sample}.report.html"
    log:
        "logs/report/{sample}.log"
    container:
        "docker://python:3.12.7-slim"
    shell:
        "python workflow/scripts/report.py --vcf {input.vcf} "
        "--flagstat {input.stats} --fastp {input.qc} > {output} 2> {log}"

Les règles de mapping suivent la chaîne imposée : sort -n, puis fixmate -m, puis sort, puis markdup, puis index. fixmate doit voir les lectures groupées par nom pour renseigner les champs de mate ; markdup doit voir un BAM trié par coordonnées. Inverser ces deux tris produit des duplicats mal marqués ou une erreur. flagstat dépend à la fois du BAM et de son index : l'index n'est pas nécessaire à flagstat, mais le déclarer rend la dépendance explicite et évite de lire un BAM non indexé dans une autre règle.

Les commandes bcftools mpileup et bcftools call sont dans un pipe. -Ou produit un BCF non compressé adapté au pipe, -Oz un VCF compressé. bcftools norm -m -any découpe les multiallélices avant le filtre, pour que les expressions portent sur des enregistrements bialléliques. Le filtre QUAL>=20 et INFO/DP>=10 est un point de départ, pas une règle universelle.

Le rapport est une règle à part, avec un conteneur Python minimal. Le script report.py lit le VCF, le flagstat et le JSON fastp, puis écrit un HTML. Il ne recalcule pas les variants : il agrège des fichiers déjà produits. Cette séparation permet de rejouer le rapport sans relancer le mapping.

L'exécution se prépare avec un essai à blanc, puis avec des cœurs et des environnements explicites :

set -euo pipefail
 
snakemake -n --reason --printshellcmds
snakemake --cores 8 --use-conda --rerun-incomplete --latency-wait 30
snakemake -n --dag | dot -Tsvg > dag.svg

-n ne lance rien ; --reason explique pourquoi chaque job serait exécuté. --latency-wait laisse au système de fichiers le temps de voir les sorties d'un job avant de conclure qu'elles manquent. --dag produit un graphe Graphviz ; il faut dot pour le convertir en SVG.

Nextflow exprime le même enchaînement avec des processus et des canaux. L'aperçu suivant ne remplace pas le Snakefile : il montre le modèle DSL2, où chaque processus est instancié par les éléments du canal.

nextflow.enable.dsl=2
 
process FASTP {
    tag "$sample"
    cpus 4
    input:
    tuple val(sample), path(r1), path(r2)
    output:
    tuple val(sample), path("${sample}_R1.trim.fastq"), path("${sample}_R2.trim.fastq")
    script:
    """
    fastp -i $r1 -I $r2 -o ${sample}_R1.trim.fastq -O ${sample}_R2.trim.fastq \\
      --detect_adapter_for_pe --thread $task.cpus
    """
}
 
workflow {
    reads = Channel.fromFilePairs("data/*_{R1,R2}.fastq")
    FASTP(reads)
}

Avec Nextflow, -resume réutilise les tâches dont l'empreinte n'a pas changé, -with-trace écrit un tableau de suivi et -with-report un rapport HTML. Le cache vit dans work/ ; le supprimer revient à repartir de zéro. Snakemake conserve son état dans .snakemake/ et propose --report pour un rapport d'exécution. Dans les deux cas, le journal et les checksums restent nécessaires pour vérifier ce qui a réellement été produit.

Étude de cas

Le dossier datasets/ngs/ contient quatre fichiers synthétiques : demo.sam, demo.vcf, reads_R1.fastq et reads_R2.fastq. demo.sam décrit un contig chr1 de 1 000 pb, un échantillon synthetic et trois alignements de 12 pb ; demo.vcf contient trois enregistrements indépendants. Les lectures FASTQ font 12 pb, une longueur trop courte pour un mapping BWA-MEM réaliste. Nous utilisons donc demo.sam comme substitut de la sortie de mapping et nous ne prétendons pas appeler de nouveaux variants.

Question : que peut-on vérifier, et que ne peut-on pas conclure, lorsqu'on fait passer ce jeu dans un pipeline QC, mapping, traitement BAM, appel, filtrage et rapport ?

L'arbre de sorties attendu, illustratif, ressemble à ceci :

qc/demo_R1.trim.fastq
qc/demo_R2.trim.fastq
qc/demo.fastp.json
aln/demo.namesorted.bam
aln/demo.fixmate.bam
aln/demo.positionsorted.bam
aln/demo.dedup.bam
aln/demo.dedup.bam.bai
qc/demo.flagstat
variants/demo.raw.vcf.gz
variants/demo.norm.vcf.gz
variants/demo.filtered.vcf.gz
results/demo.report.html

Avec demo.sam, les étapes QC et mapping sont remplacées par une règle d'entrée qui copie ou convertit le SAM déjà aligné. La chaîne restante, fixmate, tri, markdup, index, appel, normalisation, filtrage et rapport, conserve les mêmes dépendances. Le rapport doit signaler que la couverture moyenne est de 0,036× et que le nombre de lectures est de trois ; il ne doit pas présenter de variant comme biologique.

Un journal synthétique de démonstration peut ressembler à ceci. Il est étiqueté illustratif : les commandes n'ont pas été exécutées sur un cluster.

[2026-11-01T10:00:00] INFO rule fixmate completed
[2026-11-01T10:00:01] INFO rule sort_fixmate completed
[2026-11-01T10:00:02] INFO rule markdup completed
[2026-11-01T10:00:03] INFO rule index_bam completed
[2026-11-01T10:00:04] INFO rule call started
[2026-11-01T10:00:05] ERROR rule call failed
[2026-11-01T10:00:06] INFO rule normalize skipped

Le résumé attendu est : deux règles terminées avant l'échec, une règle démarrée puis échouée, une règle ignorée et aucune règle report. La reprise doit commencer par call, puis relancer normalize, filter et report ; les étapes amont peuvent être réutilisées si leurs entrées n'ont pas changé. Un checksum de demo.sam identique avant et après confirme que l'entrée n'a pas bougé, pas que l'échec vient de l'appel.

demo.vcf ne doit pas être attribué à ce pipeline : ses trois variants synthétiques n'ont ni reads correspondantes ni profondeur crédible. La conclusion est méthodologique : le pipeline peut être ordonné, journalisé et repris ; aucune conclusion génétique n'est possible sur ces données.

Erreurs fréquentes

Quiz

Exercices

python -m pytest exercises/module-7/7.2-pipeline-variants -q

python -m pytest exercises/module-7/7.2-pipeline-variants -q

python -m pytest exercises/module-7/7.2-pipeline-variants/projet -q

Résumé

• Un pipeline est un DAG : chaque règle déclare ses entrées, ses sorties et sa commande ; l'ordre topologique se calcule en $O(|V|+|E|)$.
• Snakemake raisonne en fichiers et en wildcards, Nextflow en processus et en canaux ; les deux construisent un graphe à partir des dépendances déclarées.
• La chaîne BAM correcte est sort -n, fixmate -m, sort, markdup, index, puis flagstat et l'appel de variants.
• La reprise réutilise des sorties valides, mais ne répare pas une erreur scientifique ; tout ce qui influence un résultat doit être déclaré.
• Les checksums SHA-256 détectent une modification de contenu ; ils ne prouvent ni la provenance ni la validité biologique.
• Les versions doivent être épinglées dans un environnement ou une image, et les journaux doivent contenir le nom de la règle, le statut et l'erreur.
• Sur les données synthétiques de datasets/ngs/, le pipeline vérifie la mécanique et la reprise ; aucun variant réel n'est découvert.

Références

• Mölder, F. et al. (2021). Sustainable data analysis with Snakemake. F1000Research, 10, 33.
• Di Tommaso, P. et al. (2017). Nextflow enables reproducible computational workflows. Nature Biotechnology, 35(4), 316–319.
• Li, H. & Durbin, R. (2009). Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics, 25(14), 1754–1760.
• Li, H. et al. (2009). The Sequence Alignment/Map format and SAMtools. Bioinformatics, 25(16), 2078–2079.
• Danecek, P. et al. (2021). Twelve years of SAMtools and BCFtools. GigaScience, 10(2), giab008.
• Snakemake. Documentation officielle. https://snakemake.readthedocs.io/
• Nextflow. Documentation officielle. https://www.nextflow.io/docs/latest/
• BWA. Manual and source. https://bio-bwa.sourceforge.net/