Didacticiel :
Rendre un pipeline Snakemake à l'épreuve des plateformes

"Trans-Alaska Pipeline" by Ted LaBar

Pour avoir été client des articles ("Snakemake pour les nuls""Formaliser ses protocoles avec Snakemake" et "Snakemake, aller plus loin avec la parallélisation") de mon prédécesseur lelouar, j'ai décidé d'apporter ma pierre à l'édifice et de continuer cette série sur Snakemake. Je vais ici vous parler de généralisation de pipeline pour l'utilisation intensive au sein d'une plateforme par exemple.

Pourquoi rendre mon pipeline générique?

C'est une très bonne question Jamy ! Les articles précédents se basent sur le fait que les fichiers attendus en fin de pipeline (règle target) sont tous écrits "en dur". Pas très pratique si vous devez utiliser ce pipeline pour des projets différents. Cela oblige à retoucher le code ce qui est pénible et source d'erreur. De plus le fait de devoir fixer le nombre de threads pour une règle donnée peut être assez limitant si le nombre de cœurs disponible au lancement est variable (par exemple si on est plusieurs à utiliser un serveur de calcul). Je vais donc aborder certaines techniques applicables pour rendre un pipeline plus adapté à différentes conditions d'utilisation.

Un Snakefile reste avant tout un script Python3

Il ne faut pas perdre de vue que tout code python écrit dans le Snakefile est exécuté normalement (a quelques petites exceptions près). On peut donc très bien englober les règles dans des conditions. Tout l’intérêt est de jouer avec le fichier de configuration afin de lancer une règle précise en fonction de ce qui est requis par l'utilisateur. Pour illustrer mon propos tout au long de cet article, je vais me baser sur le pipeline des articles précédents en l'adaptant à ma problématique.

Reprenons avec le fichier de configuration

Snakemake propose l'utilisation de fichiers de configuration au format JSON mais également au format YAML. C'est ce format que je vais utiliser.

Voici un petit pipeline simple réalisant l'alignement de fastq paired-end avec bwa-aln ou bien bwa-mem avec un fichier de configuration très simple :

Et le fichier de configuration :

A ce stade vous me direz qu'il n'y a pas grande différence avec un pipeline classique vu précédemment et vous aurez raison. Nous allons maintenant nous occuper de rendre tout ça un peu plus dynamique. Le but étant de reproduire plus ou moins le comportement des options de lancement d'un bon vieux script Python (argparse ce module qui vous veut du bien). Les valeurs des variables du fichier de configuration sont modifiables avec la commande de lancement Snakemake via l'option --config comme dans l'exemple ci-dessous :

snakemake --cores 30 --config genome=hg38.fa bwa=aln sampleDirectory=DNAseq

Il peut être intéressant de laisser des valeurs en dur dans le fichier de configuration afin d'avoir un comportement par défaut du pipeline. Par contre, attribuer une valeur à une variable non présente au sein du fichier de configuration (par exemple --config gemome=hg38.fa) la rajoute pour cette session de lancement. Vous voyez venir le problème ? Si vous faites une faute dans la variable que vous voulez modifier, non seulement une nouvelle variable sera créée mais votre pipeline se lancera sans soucis en conservant la valeur par défaut dans le fichier de configuration. Mais pas de panique ! Snakemake propose de valider un fichier de configuration selon un schéma YAML via le code suivant :

Pour le petit fichier de configuration utilisé ici, le schéma pourrait ressembler à celui ci-dessous :

La ligne finale spécifie que seules les variables présentes dans ce schéma doivent être présentes dans le fichier de configuration (ou via la ligne de commande snakemake). Pour aller plus loin dans l'écriture des schémas je vous invite à lire la documentation en ligne (bon courage !)

Rendons tout ça un peu plus général !

Maintenant que nous avons un pipeline bien défini avec fichier de configuration et validation YAML, occupons-nous de rendre ce pipeline applicable au plus grand nombre de cas possible. Vous avez sans doute remarqué que les fichiers attendus en fin de pipeline sont écrit en dur dans la règle target ce qui oblige d'aller modifier des chemins et des noms de fichiers dans le code du Snakefile ce qui n'est pas très agréable, vous en conviendrez. La solution est d'utiliser des fonctions d'input afin d'établir automatiquement la liste des fichiers attendus en analysant par exemple le contenu d'un dossier (donné en option) contenant nos échantillons à aligner. Cela permettrait de ne changer uniquement le dossier d'entrée via la commande --config et ne plus toucher au code du pipeline.

Prenons l'exemple suivant : Joe travaille au sein d'une plateforme de séquençage qui enchaîne les projets. Comment améliorer le pipeline du dessus afin de pouvoir l'utiliser pour tous les projets paired-end ? Dans un premier temps on peut détecter automatiquement les fichiers d'entrée et de sortie comme expliqué au paragraphe précédent mais on peut également adapter le nombre de threads alloués à chaque règle en fonction de la taille du projet. Dans cet exemple la plateforme exige que le pipeline fasse tourner le maximum d'échantillons en parallèle sur un serveur de 32 cœurs et que les fichiers issus d'un séquençage Illumina soient formatés sous la forme nom_S[numéro_individu]_L00[numéro_lane_de_séquençage]_R[1,2]_001.fastq.gz

Exemple pour un projet X de séquençage du rat en paired-end:

  • Rattus_projetX_S1_L001_R1_001.fastq.gz
  • Rattus_projetX_S1_L001_R2_001.fastq.gz

Pour commencer, grâce au fichier de config nous pouvons récupérer la liste des fichiers présents dans le dossier contenant ces deux fastq :

On peut créer une fonction retournant la liste des fichiers attendus en fonction du contenu du dossier d'input en utilisant les informations récupérées dans LIST_NAME :

Et donc la règle target s'écrit désormais comme cela :

Et hop ! Vous avez une fonction d'input automatique ! Bien évidemment ce code ne fonctionne que pour notre exemple précis, si les fastq ne sont pas formatés comme demandé, le pipeline ne trouvera pas vos échantillons. A vous d'adapter ce code aux spécificités de votre labo.

La fonction lambda

Si utiliser des fonctions Python classiques fonctionne parfaitement dans Snakemake, utiliser les wildcards de ce dernier est encore plus efficace. Pour continuer sur la lancée de la partie précédente, on peut écrire deux petites fonctions qui prendront en entrée la valeur de la wildcards définie par Snakemake afin d'avoir une continuité dans l'automatisation :

Ces deux fonctions sélectionnent les fichiers fastq R1 et R2 en fonction d'une wildcards donnée. Cette wildcards va prendre automatiquement comme valeur les noms présents dans LIST_NAME car c'est à partir de la règle target que Snakemake sait quels fichiers doivent être créés et donc quelle valeur attribuer aux wildcards.

Gérer l'attribution du nombre de threads

Étrangement, il n'existe pas de façon simple pour récupérer la valeur du nombre de cœurs disponible spécifiée au lancement (--cores ou -j). Mais l'idée est de récupérer cette valeur afin d'adapter le nombre de threads pour chaque règle en fonction du nombre d'échantillons à analyser.

J'ai donc utilisé sys.argv (la liste contenant les arguments de la commande de lancement d'un programme python) afin de récupérer cette valeur mais le comportement de snakemake étant assez particulier, la récupération doit se faire en deux temps. En effet, Snakemake exécute le code une première fois (votre commande de lancement) et sys.argv contient toute vos options de lancement. Ensuite le code s'exécute une seconde fois et Snakemake modifie les valeurs de sys.argv par ses propres valeurs. Oui c'est tordu ! Essayez de print sys.argv a l'exécution d'un pipeline vous allez voir.

Je vous donne donc la fonction toute faite ce qui vous évitera chercher deux heures comme moi :

Cette fonction vous permettra de récupérer le nombre de cœurs max disponible. Vous pourrez ensuite par exemple décider de lancer le plus d'analyse possible en parallèle quel que soit le nombre d'échantillons. Snakemake permet de mettre une variable au nombre de threads par règle au lieu d'un nombre fixe. Il est donc intéressant de créer une petite conditions qui adaptera ce nombre au nombre d'échantillons.

Au final

Avec toutes ces modifications, le pipeline est désormais capable de réaliser l'alignement sur bwa-mem ou aln au choix, en fonction d'un dossier d'input et d'optimiser tout seul le nombre de threads en fonction du nombre d'échantillons à traiter. Le pipeline a été testé sur une machine locale sous Ubuntu 19.04 avec 16Gb de Ram et un i5-8250U.

Je vous redonne le code total du pipeline optimisé :

Ce pipeline peut évidemment être encore optimisé notamment via l'ajout de messages personnalisés explicatif à l'exécution des règles par exemple. Il est également possible d'imaginer une adaptation à l'exécution sur cluster.

Je remercie les relecteurs de cet article Akira, Lins et lelouar pour leur aide et leur avis pendant la rédaction de cet article.

Source photo : "Trans-Alaska Pipeline"

  • À propos de
  • Ingénieur d'études en analyse de données de séquençage haut-débit au sein de la plateforme de séquençage Montpellier GenomiX

Laisser un commentaire