Formaliser ses protocoles avec Snakemake

Vous avez dit protocole ?

workflow
Mon pro­chain pro­to­cole. Sans aucun doute un pas de géant pour la science.

Qui dit bio­lo­gie et bio­in­for­ma­tique dit pro­to­cole expé­ri­men­tal. C'est le cœur de la démarche scien­ti­fique, et un for­ma­lisme adap­té est la clef pour assu­rer la repro­duc­ti­bi­li­té des expé­riences, et ain­si garan­tir la vali­da­tion des décou­vertes par la com­mu­nau­té. En paillasse, les solu­tions pour for­ma­li­ser et conser­ver les pro­to­coles sont plu­tôt natu­rel­le­ment prag­ma­tiques et éprou­vées par les années. Le papier a encore de beaux jours devant lui. Cepen­dant, il est assez dif­fi­cile d'adapter le prin­cipe du cahier de manip au monde de l'informatique, car il est inef­fi­cace tel quel.

À la place, beau­coup ont adop­té les sys­tèmes de ges­tion de work­flow, qui per­mettent de retrans­crire de manière assez natu­relle la notion de repro­duc­ti­bi­li­té au sup­port infor­ma­tique. Beau­coup d'outils spé­cia­li­sés à la bio­in­for­ma­tique ont depuis vu le jour : Galaxy (dont j'avais déjà par­lé ici), Taver­na, Erga­tis, etc. Ils per­mettent de créer, exé­cu­ter et par­ta­ger des pro­to­coles entiers, mais la plu­part sont déve­lop­pés pour offrir une inter­face gra­phique à un uti­li­sa­teur novice en pro­gram­ma­tion, et sont peu voire pas adap­tés à l'utilisation en ligne de com­mande. De plus, l'utilisateur est sou­vent dépen­dant d'un ensemble d'outils pré-inté­grés au sys­tème, ce qui est sou­vent han­di­ca­pant.

Rendez-moi mon shell !

GNU
Hono­rable GNU com­pi­lant du libre avec l'ancestral Make

L'outil de réfé­rence pour créer des work­flows en ligne de com­mande, vers lequel la plu­part des gens se tournent encore aujourd'hui, est GNU Make. Il est bien connu, docu­men­té, éprou­vé, et per­met de défi­nir des work­flows rami­fiés. Un de ses points forts est le fait que si une par­tie du work­flow seule­ment doit être modi­fiée, il n'exécutera que les par­ties néces­saires. Étant don­nés les volumes des don­nées mani­pu­lées en bio­in­for­ma­tique, c'est un avan­tage de taille.

Même s'il est géné­ra­le­ment uti­li­sé pour la com­pi­la­tion de logi­ciels, il est très adap­té à la démarche scien­ti­fique. En effet, lorsque vous défi­nis­sez un pro­to­cole expé­ri­men­tal, il y a de grandes chances pour qu'il soit à peu de choses près direc­te­ment trans­po­sable en Make­file, les dif­fé­rentes étapes du pro­to­cole sont sou­vent assi­mi­lable à des règles pour GNU Make.

Mighty Python

Pytroll
Jeune et frin­gant Python trol­lant éhon­té­ment le véné­rable GNU. Ah, la fougue de la jeu­nesse, quelle imper­ti­nence !

Seule­ment voi­là : GNU Make est un dino­saure pous­sié­reux sacré­ment com­pli­qué à prendre en main, et long à maî­tri­ser. Heu­reu­se­ment, il existe des ini­tia­tives dont le but est de pro­po­ser des alter­na­tives plus modernes à cette hono­rable anti­qui­té. C'est le cas, par exemple, de SCONS, Ant, Waf, Rake, Aap, Nin­ja et Sna­ke­make. Même si je ne trai­te­rai ici que de ce der­nier, qui pro­pose une solu­tion simi­laire à GNU Make basé sur Python, il n'est pas exclu que d'autres alter­na­tives puissent être aus­si valables dans ce contexte.

Le prin­cipe est simi­laire : On défi­nit des étapes/​règles qui ont cha­cune des fichiers d'entrée, des fichiers de sor­tie, et des choses à exé­cu­ter pour pro­duire les fichiers de sor­tie à par­tir des fichiers d'entrée. On les écrit dans un fichier nom­mé Sna­ke­file (équi­valent d'un Make­file chez GNU Make) et on lance l'exécution avec la com­mande sna­ke­make. Il pro­pose en plus de retrans­crire le work­flow en graphe acy­clique orien­té (DAG) pour mieux le visua­li­ser, ce qui per­met notam­ment de véri­fier si le Sna­ke­file cor­res­pond bien au pro­to­cole expé­ri­men­tal dési­ré.

Il n'est d'ailleurs tel­le­ment pas sur­pre­nant d'utiliser Sna­ke­make dans ce contexte que le pre­mier exemple four­ni dans leur docu­men­ta­tion est un exemple de bio­in­for­ma­tique. Ça ne s'invente pas. En voi­ci la tra­duc­tion :

Cuf­flinks est un outil per­met­tant d'assembler des trans­crits, esti­mer leur abon­dance, et effec­tuer des ana­lyses d'expression dif­fé­ren­tielle sur des don­nées RNA-Seq. Cet exemple montre com­ment créer un work­flow typique pour Cuf­flinks avec sna­ke­make. Il sup­pose que des don­nées RNA-Seq ali­gnées sont four­nies sous for­mat BAM pour quatre échan­tillons 101–104. Pour chaque échan­tillon, les trans­crits sont assem­blés avec '''cuf­flinks''' (rule assem­bly). Les assem­blages sont fusion­nés en un fichier GTF avec '''cuff­merge''' (rule merge_​assemblies). Une com­pa­rai­son avec la piste GTF de hg19 est ensuite effec­tuée (rule compare_​assemblies). Enfin, les expres­sions dif­fé­ren­tielles sont cal­cu­lées sur les trans­crits trou­vés (rule dif­fexp).

TRACK = 'hg19.gtf' # should be fetched from the cufflinks page since special fields are expected
REF = 'hg19.fa'

CLASS1 = '101 102'.split()
CLASS2 = '103 104'.split()
SAMPLES = CLASS1 + CLASS2

rule all:
input: 'diffexp/isoform_exp.diff', 'assembly/comparison'

rule assembly:
input: 'mapped/{sample}.bam'
output: 'assembly/{sample}/transcripts.gtf', dir='assembly/{sample}'
threads: 4
shell: 'cufflinks --num-threads {threads} -o {output.dir} --frag-bias-correct {REF} {input}'

rule compose_merge:
input: expand('assembly/{sample}/transcripts.gtf', sample=SAMPLES)
output: txt='assembly/assemblies.txt'
run:
with open(output.txt, 'w') as out:
out.write('\n'.join(input))

rule merge_assemblies:
input: 'assembly/assemblies.txt'
output: 'assembly/merged/merged.gtf', dir='assembly/merged'
shell: 'cuffmerge -o {output.dir} -s {REF} {input}'

rule compare_assemblies:
input: 'assembly/merged/merged.gtf'
output: 'assembly/comparison/all.stats', dir='assembly/comparison'
shell: 'cuffcompare -o {output.dir}all -s {REF} -r {TRACK} {input}'

rule diffexp:
input: expand('mapped/{sample}.bam', sample=SAMPLES), gtf='assembly/merged/merged.gtf'
output: 'diffexp/gene_exp.diff', 'diffexp/isoform_exp.diff'
threads: 8
run:
classes = [','.join(expand('mapped/{sample}.bam', sample=cls)) for cls in (CLASS1, CLASS2)]
shell('cuffdiff --num-threads {threads} {gtf} {classes[0]} {classes[1]}')

On peut remar­quer qu'il s'agit ni plus ni moins d'un script python : il est tout à fait pos­sible d'utiliser tous les objets et fonc­tions dis­po­nibles dans le lan­gage, et impor­ter des modules en fonc­tion des besoins. Pour chaque règle, il est pos­sible de deman­der direc­te­ment l'exécution de code Python direc­te­ment inté­gré dans le Sna­ke­file lorsqu'il est pla­cé dans le bloc "run" (voir rule compose_​merge). Bref, si vous êtes fami­lier avec ce lan­gage, vous ne serez pas dépay­sé.

Retours d'expériences (pun intended)

Myron Fass
Jean-Clo­taire, indé­crot­table pytho­niste, découvre Sna­ke­make avec une joie non dis­si­mu­lée.

Bien enten­du, je ne vous en par­le­rais pas si je ne l'avais pas tes­té moi-même. Pour être tout à fait hon­nête, je pense qu'une solide expé­rience de déve­lop­pe­ment en Python n'est pas du luxe. Je trouve Sna­ke­make plus user-friend­ly que Make, mais il reste un sys­tème com­plexe qui deman­de­ra de la patience et quelques efforts à prendre en main. Cela dit, quand vous mani­pu­lez de grands volumes de don­nées, c'est agréable de pou­voir faci­le­ment choi­sir quoi réexé­cu­ter en fonc­tion des besoins.

Je l'utilise dans le cadre d'un pro­jet de recherche où les pro­to­coles ne sont pas pré-éta­blis, et dans ce contexte l'avantage d'avoir un sys­tème faci­le­ment édi­table et exé­cu­table en ligne de com­mande se fait évident. J'apprécie beau­coup la sou­plesse du sys­tème, il faut avouer que c'est autre chose que les scripts shell cra­cras ou les Make­files abs­cons. Sna­ke­make a encore du che­min à faire, il y a pas mal de fonc­tion­na­li­tés qui manquent un peu ou semblent être encore en chan­tier, mais il est tout à fait uti­li­sable en l'état, et je l'ai adop­té.

L'option per­met­tant de visua­li­ser le pipe­line sous forme de DAG est éga­le­ment très appré­ciable, car elle per­met de visua­li­ser l'état d'avancement de son exé­cu­tion. Les étapes déjà effec­tuées s'affichent en poin­tillés et celles qui seront exé­cu­tées au pro­chain appel de la com­mande sont entou­rées d'un trait plein. Voi­ci ce que ça donne pour un Sna­ke­make de mon pro­jet :

DAG Snakefile
Com­mande uti­li­sée (néces­site Gra­ph­viz):
sna­ke­make –dag | dot ‑Tpng ‑o dag.png

Je vous invite à visi­ter leur wiki si vous sou­hai­tez en savoir plus. Sna­ke­make est héber­gé sur Bit­bu­cket : https://​bit​bu​cket​.org/​j​o​h​a​n​n​e​s​k​o​e​s​t​e​r​/​s​n​a​k​e​m​ake

_​_​_​_​_​_​_​_​_​_​_​_​_​_​_​_​_​_​_​_​_​_​

Le pre­mier à me dire ce qui ne va pas dans mon superbe pro­to­cole décrit dans la pre­mière illus­tra­tion (à part les noms débiles) aura gagné le coco­tier.

En effet, il est incom­pa­tible avec Sna­ke­make. Pour­quoi ?

_​_​_​_​_​_​_​_​_​_​_​_​_​_​_​_​_​_​_​_​_​_​

Dis­clai­mer : Tous les conte­nus uti­li­sés pour les images finales de cet article sont soit dans le domaine public (1) soit sous licence Crea­tive Com­mons Attri­bu­tion-Share Alike 3.0 Unported (2,3,4), soit de ma pro­duc­tion.

 Mer­ci à Sam pour m'avoir fait décou­vrir Sna­ke­make et ain­si don­né l'idée géné­rale de cet article, et à mes relec­teurs ché­ris pour leur patience, leurs cor­rec­tion et leurs idées : Guillaume, Mica, Nal­lias et ZaZoOo (sisi miss, au moins pour la patience!)



Pour continuer la lecture :


Commentaires

6 réponses à “Formaliser ses protocoles avec Snakemake”

  1. Le petit sou­ci du pre­mier work­flow est son cycle (le retour cyclique sur "bar"), non ?
    Les ges­tion­naire de work­flow n'aiment pas ça !

    A par­tir de votre article j'ai jeté un coup d'oeil à qques outils ruby : pwrake semble être un bon can­di­dat, je n'ai pas encore tes­té, si qqn l'a déjà éprou­vé ?

    1. _
      /_'. _
      _ / '-.
      -.;),–'
      '--.--.
      / |/-/'._
      |/ |=|
      |_|
      |-| - cocotier |=| |_| |=|
      """`

      pwrake a l'air cool aus­si ! (mer­ci, j'étais lamen­ta­ble­ment pas­sée à côté) mais pas facile de trou­ver de la doc qui ne soit pas en japo­nais… pour info pour ceux qui seraient inté­res­sés : Agile paral­lel bio­in­for­ma­tics work­flow mana­ge­ment using Pwrake

      1. bah bra­vo, heu­reu­se­ment que j'ai mis une éti­quette, WP m'a sau­va­ge­ment rati­boi­sé mon coco­tier… :/​

        1. mer­ci pour la doc pwrake, j'y jet­te­rai un coup d'oeil !

  2. Un incon­vé­nient : sna­ke­make il fonc­tionne avec python 3…

    1. Faux incon­vé­nient 🙂 va bien fal­loir lais­ser python2 au musée un jour quand même ! Embrasse donc le pro­grès et la moder­ni­té cher ami ! ^^

      Nan sérieux, soyons hon­nêtes, j'ai com­pi­lé, ins­tal­lé et main­tiens les 2.7 et 3.3 en modules au labo (pasque faut pas comp­ter sur Cen­tOS pour du récent), c'est juste too easy. Ça prend juste un peu de temps à com­pi­ler, faut bien deman­der à com­pi­ler toutes les libs dont t'as besoin, mais sinon vrai­ment ça se fait tout seul ou presque. Et une fois que c'est fait c'est fait.

Laisser un commentaire