Astuce :
Snakemake pour les nuls (ou comment créer un pipeline facilement ?)

Bonjour à tous, et bienvenue dans le premier épisode d'une (longue ?) série de prise en main de l'outil dédié au pipeline : Snakemake.

Si vous ne connaissez pas encore cet outil, c'est que vous êtes sûrement passés à côté de cet article écrit par Nisaea. Alors, quel sera les bénéfices de retranscrire vos pipelines déjà tout prêt en Snakefile ?

Lisibilité du code, gestion des ressources et reproductibilité

Lorsque vous êtes sur le point de publier,  il va bien falloir expliquer aux futurs lecteurs comment vous avez obtenu les résultats. Cela permettra ainsi aux autres bioinformaticiens de partir de vos données brutes, et de retrouver les mêmes résultats que vous. On approche ici un point clé de la bioinformatique : la reproductibilité. Quels ont été les outils utilisés, à quelles versions, avec quels paramètres... Tous ces petits "trucs" qui permettent d'obtenir votre résultat. Un article scientifique avec des résultats majeurs découverts grâce à la bioinformatique DEVRAIT être associé à un pipeline/script permettant de reproduire les résultats à l'identique.

Ces dernières années ont vu le volume de données augmenter de façon quasi exponentielle. En parallèle, les ressources mises à disposition (telles que les clusters de calcul) n'ont cessé de croître en taille (nombre de cpu) et en puissance de calcul (vitesse des cpu), permettant ainsi de rester dans la course. Cependant, pour les exploiter au mieux, il est nécessaire de paralléliser nos tâches. En effet, beaucoup d'outils de bioinformatique ont des options permettant d'utiliser simultanément plusieurs cpu  (ex.: bwa, STAR, ...). D'autres, en revanche, n'ont pas été conçus pour être utilisés dans un environnement multi-cpu (awk, macs2, bedtools, ...). Dans ces cas précis, il faut soit paralléliser à la main (lancer plusieurs tâches en fond avec l'ajout d'un '&' à la fin de la commande), soit mettre les opérations en file d'attente et sous-exploiter la machine sur laquelle vous travaillez (gestion séquentielle des tâches). Nous verrons dans un prochain article qu'avec Snakemake on peut, très facilement, aller plus loin dans la parallélisation.

L'inventeur de Snakemake, Johannes Köster [1], a su tirer avantage de Python (syntaxe claire) et de GNU make (parallélisation, gestion des ressources,  système de règles). Avec l'ajout de nombreuses fonctionnalités comme le tracking de version/paramètres, la création de rapport html ou la visualisation des tâches et de leurs dépendances sous forme de graph (DAG), ce langage a tout le potentiel pour avoir un bel avenir.

Principes généraux

Tout comme GNU make, Snakemake fonctionne sur le principe des règles (rule). Une règle est un ensemble d'élément qui vont permettre de créer un fichier, c'est en quelque sorte une étape dans le pipeline.

Chaque règle a au moins un fichier de sortie (output) et un (ou plusieurs) fichier(s) d'entrée (input) (j'ai volontairement commencé par l'output, vision Snakemakienne, je pars de la fin pour aller au début).

La toute première règle du Snakefile doit définir les fichiers que l'on veut à la fin du traitement (fichier cible/target). L'ordre (et le choix) des règles est établi automatiquement à l'aide du nom des fichiers/dossiers cibles. On va ainsi remonter les fichiers de règle en règle jusqu'à trouver une règle avec un input plus récent que l'output. Ensuite, cette règle et toutes celles qui suivent vont être exécutées.

Snakemake fonctionne avec le nom des fichiers

Le plus simple serait de commencer par un exemple :

Vous l'aurez deviné, s'il existe data/raw/s1.fq et data/raw/s2.fq, plus récents que data/raw/s1.fq.gz et data/raw/s2.fq.gz alors la règle gzip va créer/remplacer les cibles. De plus, on peut paralléliser les opérations en passant en paramètre le nombre de processeurs à utiliser (option -j N, --jobs=N).

Chaque règle peut avoir plusieurs mots clés

Snakemake utilise des mots clés pour définir les règles.

Entrée/sortie

input = fichier(s) à utiliser pour créer la sortie
output = fichier(s) qui seront générés avec la règle

wildcards = Ce sont des variables qui sont définies avec le {} et une partie du nom/chemin des fichiers de sortie, cela permet de généraliser les règles. Avec cet exemple nous avons utilisé deux wildcards (genome et sample), on peut également utiliser des expressions régulières pour bien délimiter la définition d'un wildcards comme pour genome qui vaut "hg[0-9]+".

Utilisation de l'objet wildcards (par section) :

  • log : "mapping/{genome}/{sample}.sort.log"
  • shell/run/message :  "mapping/{wildcards.genome}/{wildcards.sample}.sort.log"
  • params : param1 = lambda wildcards : wildcards.genome + '.fa'

Explication :

En dehors de input/output/log on peut utiliser directement le nom de variable car l'objet wildcards n'est connu que de shell/run/message. C'est pour cela que l'on passe à la forme  "wildcards.variable" dans ses 3 dernières sections. En dehors de shell/run/message on peut utiliser "lambda wildcards : wildcards.variable". Ici la fonction lambda permet d'utiliser directement l'objet wildcards dans les sections.

Traitement des fichiers créés

temps = fichier temporaire supprimé après utilisation.
protected =  ce fichier sera protégé en écriture après création.

Commande pour générer le(s) fichier(s) sortie (run/shell)

shell = utiliser une commande UNIX pour créer le(s) fichier(s) sortie.
run = utiliser du code python pour construire le(s) fichier(s) sortie.

note: on doit choisir entre run et shell
note2 : avec run on peut utiliser la fonction shell()

Autres paramètres des règles

threads = nombre de processeurs maximum utilisables
message = message à afficher au début
log = fichier log
priority = permet de prioriser les règles (ex: quand plusieurs de règles utilisent le même input)
params* = paramètres des commandes utilisées
version* = version de la règle (ou de la commande utilisée)

*: ces deux mots clés permettent d'associer les paramètres et les versions des règles aux fichiers de sortie. Ceci permet de (re)générer un fichier de sortie si un paramètre ou une version de la règle a été modifiée. Il y a aussi du code tracking, pour suivre l'évolution du pipeline et des fichiers générés.

Un exemple plus concret :

Ici nous allons créer un pipeline d'analyse données génomique (ex : ChIP-seq) qui va aligner/mapper de petites séquences issues du séquençage haut débit (reads/lectures) sur un génome de référence.

Avec environ 130 lignes, ce code est déjà prêt à être utilisé de manière intensive sur pc, serveur ou cluster de calcul 🙂

A partir de ce Snakefile, on peut générer un graphe représentant la dépendance entre les règles avec la commande  :

 

rulegraph

Graphe représentant la dépendance entre les règles

On peut également générer un graphe plus complet qui prend aussi en compte les wildcards avec la commande :

 

dag

Graphe représentant le pipeline complet

 

Cet article touche à sa fin,  je vous donne rendez-vous au prochain épisode, d'ici là, vous pouvez poser des questions dans les commentaires (ils seront peut être utilisés pour les prochains articles).

Et si ça vous a plu, n'hésitez pas à partager et nous faire part de vos astuces, proposez-nous des règles ou même des Snakefiles dans les commentaires, si cela a du succès nous ouvrirons peut-être un dépôt git pour regrouper toutes ces contributions.

Un grand merci aux relecteurs Yoann M., Lroy et Kumquat[um] pour leur remarques et corrections 😉

[1] Johannes Köster, "Parallelization, Scalability, and Reproducibility in Next-Generation Sequencing Analysis", TU Dortmund 2014

  • À propos de
  • Après un DUT génie biologique option bioinfo à Aurillac, j'ai poursuivi mes études avec une licence professionnelle (traitement des données génomiques) à Carcassonne puis un master (bioinformatique et biologie des systèmes) à Toulouse. J'ai ensuite eu une expérience d'un an à l'Institut de Génétique Humaine à Montpellier où je réalise actuellement une thèse.

Un commentaire sur “Snakemake pour les nuls (ou comment créer un pipeline facilement ?)

  1. Bonjour,

    Comme très justement écrit en introduction, il serait intéressant d\'avoir accès aux fichiers tests présents dans cet exemple afin d\'en tester la reproductibilité.

    Il semble que le génome soit accessible ici : http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz

    Pour les fichiers raw de l\'exemple (sampleX_c3-pe_R1.fq.gz et sampleX_c3-pe_R2.fq.gz), cela semble plus complexe.
    Une idée où les trouver?

    Merci

Laisser un commentaire