Bon­jour à tous, et bien­ve­nue dans le pre­mier épi­sode d'une (longue ?) série de prise en main de l'outil dédié au pipe­line : Sna­ke­make.

Si vous ne connais­sez pas encore cet outil, c'est que vous êtes sûre­ment pas­sés à côté de cet article écrit par Nisaea. Alors, quel sera les béné­fices de retrans­crire vos pipe­lines déjà tout prêt en Sna­ke­file ?

Lisibilité du code, gestion des ressources et reproductibilité

Lorsque vous êtes sur le point de publier,  il va bien fal­loir expli­quer aux futurs lec­teurs com­ment vous avez obte­nu les résul­tats. Cela per­met­tra ain­si aux autres bio­in­for­ma­ti­ciens de par­tir de vos don­nées brutes, et de retrou­ver les mêmes résul­tats que vous. On approche ici un point clé de la bio­in­for­ma­tique : la repro­duc­ti­bi­li­té. Quels ont été les outils uti­li­sés, à quelles ver­sions, avec quels para­mètres… Tous ces petits "trucs" qui per­mettent d'obtenir votre résul­tat. Un article scien­ti­fique avec des résul­tats majeurs décou­verts grâce à la bio­in­for­ma­tique DEVRAIT être asso­cié à un pipeline/​script per­met­tant de repro­duire les résul­tats à l'identique.

Ces der­nières années ont vu le volume de don­nées aug­men­ter de façon qua­si expo­nen­tielle. En paral­lèle, les res­sources mises à dis­po­si­tion (telles que les clus­ters de cal­cul) n'ont ces­sé de croître en taille (nombre de cpu) et en puis­sance de cal­cul (vitesse des cpu), per­met­tant ain­si de res­ter dans la course. Cepen­dant, pour les exploi­ter au mieux, il est néces­saire de paral­lé­li­ser nos tâches. En effet, beau­coup d'outils de bio­in­for­ma­tique ont des options per­met­tant d'utiliser simul­ta­né­ment plu­sieurs cpu  (ex.: bwa, STAR, …). D'autres, en revanche, n'ont pas été conçus pour être uti­li­sés dans un envi­ron­ne­ment mul­ti-cpu (awk, macs2, bed­tools, …). Dans ces cas pré­cis, il faut soit paral­lé­li­ser à la main (lan­cer plu­sieurs tâches en fond avec l'ajout d'un '&' à la fin de la com­mande), soit mettre les opé­ra­tions en file d'attente et sous-exploi­ter la machine sur laquelle vous tra­vaillez (ges­tion séquen­tielle des tâches). Nous ver­rons dans un pro­chain article qu'avec Sna­ke­make on peut, très faci­le­ment, aller plus loin dans la paral­lé­li­sa­tion.

L'inventeur de Sna­ke­make, Johannes Kös­ter [1], a su tirer avan­tage de Python (syn­taxe claire) et de GNU make (paral­lé­li­sa­tion, ges­tion des res­sources,  sys­tème de règles). Avec l'ajout de nom­breuses fonc­tion­na­li­tés comme le tra­cking de version/​paramètres, la créa­tion de rap­port html ou la visua­li­sa­tion des tâches et de leurs dépen­dances sous forme de graph (DAG), ce lan­gage a tout le poten­tiel pour avoir un bel ave­nir.

Principes généraux

Tout comme GNU make, Sna­ke­make fonc­tionne sur le prin­cipe des règles (rule). Une règle est un ensemble d'élément qui vont per­mettre de créer un fichier, c'est en quelque sorte une étape dans le pipe­line.

Chaque règle a au moins un fichier de sor­tie (out­put) et un (ou plu­sieurs) fichier(s) d'entrée (input) (j'ai volon­tai­re­ment com­men­cé par l'output, vision Sna­ke­ma­kienne, je pars de la fin pour aller au début).

La toute pre­mière règle du Sna­ke­file doit défi­nir les fichiers que l'on veut à la fin du trai­te­ment (fichier cible/​target). L'ordre (et le choix) des règles est éta­bli auto­ma­ti­que­ment à l'aide du nom des fichiers/​dossiers cibles. On va ain­si remon­ter les fichiers de règle en règle jusqu'à trou­ver une règle avec un input plus récent que l'output. Ensuite, cette règle et toutes celles qui suivent vont être exé­cu­tées.

Snakemake fonctionne avec le nom des fichiers

Le plus simple serait de com­men­cer par un exemple :

Vous l'aurez devi­né, 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 paral­lé­li­ser les opé­ra­tions en pas­sant en para­mètre le nombre de pro­ces­seurs à uti­li­ser (option -j N, --jobs=N).

Chaque règle peut avoir plusieurs mots clés

Sna­ke­make uti­lise des mots clés pour défi­nir les règles.

Entrée/​sortie

input = fichier(s) à uti­li­ser pour créer la sor­tie
out­put = fichier(s) qui seront géné­rés avec la règle

wild­cards = Ce sont des variables qui sont défi­nies avec le {} et une par­tie du nom/​chemin des fichiers de sor­tie, cela per­met de géné­ra­li­ser les règles. Avec cet exemple nous avons uti­li­sé deux wild­cards (genome et sample), on peut éga­le­ment uti­li­ser des expres­sions régu­lières pour bien déli­mi­ter la défi­ni­tion d'un wild­cards comme pour genome qui vaut "hg[0-9]+".

Uti­li­sa­tion de l'objet wild­cards (par sec­tion) :

• log : "mapping/​{genome}/​{sample}.sort.log"
• shell/​run/​mes­sage :  "mapping/​{wild­cards.genome}/​{wild­cards.sample}.sort.log"
• params : param1 = lamb­da wild­cards : wildcards.genome + '.fa'

Expli­ca­tion :

En dehors de input/​output/​log on peut uti­li­ser direc­te­ment le nom de variable car l'objet wild­cards n'est connu que de shell/​run/​message. C'est pour cela que l'on passe à la forme  "wildcards.variable" dans ses 3 der­nières sec­tions. En dehors de shell/​run/​message on peut uti­li­ser "lamb­da wild­cards : wildcards.variable". Ici la fonc­tion lamb­da per­met d'utiliser direc­te­ment l'objet wild­cards dans les sec­tions.

Traitement des fichiers créés

temps = fichier tem­po­raire sup­pri­mé après uti­li­sa­tion.
pro­tec­ted =  ce fichier sera pro­té­gé en écri­ture après créa­tion.

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

shell = uti­li­ser une com­mande UNIX pour créer le(s) fichier(s) sor­tie.
run = uti­li­ser du code python pour construire le(s) fichier(s) sor­tie.

note : on doit choi­sir entre run et shell
note2 : avec run on peut uti­li­ser la fonc­tion shell()

Autres paramètres des règles

threads = nombre de pro­ces­seurs maxi­mum uti­li­sables
mes­sage = mes­sage à affi­cher au début
log = fichier log
prio­ri­ty = per­met de prio­ri­ser les règles (ex : quand plu­sieurs de règles uti­lisent le même input)
params* = para­mètres des com­mandes uti­li­sées
ver­sion* = ver­sion de la règle (ou de la com­mande uti­li­sée)

*: ces deux mots clés per­mettent d'associer les para­mètres et les ver­sions des règles aux fichiers de sor­tie. Ceci per­met de (re)générer un fichier de sor­tie si un para­mètre ou une ver­sion de la règle a été modi­fiée. Il y a aus­si du code tra­cking, pour suivre l'évolution du pipe­line et des fichiers géné­rés.

Un exemple plus concret :

Ici nous allons créer un pipe­line d'analyse don­nées géno­mique (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 envi­ron 130 lignes, ce code est déjà prêt à être uti­li­sé de manière inten­sive sur pc, ser­veur ou clus­ter de cal­cul 🙂

A par­tir de ce Sna­ke­file, on peut géné­rer un graphe repré­sen­tant la dépen­dance entre les règles avec la com­mande  :

 

On peut éga­le­ment géné­rer un graphe plus com­plet qui prend aus­si en compte les wild­cards avec la com­mande :

 

 

Cet article touche à sa fin,  je vous donne ren­dez-vous au pro­chain épi­sode, d'ici là, vous pou­vez poser des ques­tions dans les com­men­taires (ils seront peut être uti­li­sés pour les pro­chains articles).

Et si ça vous a plu, n'hésitez pas à par­ta­ger et nous faire part de vos astuces, pro­po­sez-nous des règles ou même des Sna­ke­files dans les com­men­taires, si cela a du suc­cès nous ouvri­rons peut-être un dépôt git pour regrou­per toutes ces contri­bu­tions.

Un grand mer­ci aux relec­teurs Yoann M., Lroy et Kumquat[um] pour leur remarques et cor­rec­tions 😉

[1] Johannes Kös­ter, "Paral­le­li­za­tion, Sca­la­bi­li­ty, and Repro­du­ci­bi­li­ty in Next-Gene­ra­tion Sequen­cing Ana­ly­sis", TU Dort­mund 2014