Il y a donc un fichier d'entrée, qui est le résultat de la tâche /download/. Et il y a les deux fichiers de sortie. Enfin, pour faire le travail, j'ai opté pour un script Python cette fois. =snakemake= reconnaît le langage par l'extension =.py=.
Il y a donc un fichier d'entrée, qui est le résultat de la tâche /download/. Et il y a les deux fichiers de sortie, un pour les résultats et un pour les messages d'erreur. Enfin, pour faire le travail, j'ai opté pour un script Python cette fois. =snakemake= reconnaît le langage par l'extension =.py=.
Le contenu de ce script est presque un copier-coller d'un document computationnel du module 3, plus précisément du document que j'ai montré dans le parcours Emacs/Org-Mode:
Le contenu de ce script est presque un copier-coller d'un document computationnel du module 3, plus précisément du document que j'ai montré dans le parcours Emacs/Org-Mode:
#+begin_src python :exports both :tangle incidence_syndrome_grippal_snakemake/scripts/preprocess.py
#+begin_src python :exports both :tangle incidence_syndrome_grippal_snakemake/scripts/preprocess.py
Le workflow que je viens de montrer produit 7 fichiers. Ce n'est pas beaucoup. On peut les nommer à la main, un par un, sans difficulté. Dans la vraie vie, par exemple en bioinformatique, un workflow peut facilement gérer des centaines ou milliers de fichiers, par exemple un fichier par séquence d'acides aminés dans une étude de protéomique. Dans une telle situation, il faut définir un schéma pour nommer les fichiers de façon systématique, et introduire des boucles dans le workflow dont les itérations seront idéalement exécutées en parallèle. Je vais illustrer ceci avec une autre décomposition de l'analyse de l'incidence du syndrome grippal.
Le workflow que je viens de montrer produit 7 fichiers. Ce n'est pas beaucoup. On peut les nommer à la main, un par un, sans difficulté. Dans la vraie vie, par exemple en bioinformatique, un workflow peut facilement gérer des centaines ou milliers de fichiers, par exemple un fichier par séquence d'acides aminés dans une étude de protéomique. Dans une telle situation, il faut définir un schéma pour nommer les fichiers de façon systématique, et introduire des boucles dans le workflow dont les itérations seront idéalement exécutées en parallèle. Je vais illustrer ceci avec une variante de l'analyse de l'incidence du syndrome grippal. Elle utilise une forme plus détaillée des données brutes dans laquelle les incidence sont repertoriées par région plutôt que pour la France entière. Il faut donc répéter le calcul de l'incidence annuelle 13 fois, une fois pour chaque région. Pour simplifier un peu, le résultat principal de ce nouveau workflow sera un fichier qui contient, pour chaque région, l'année dans laquelle l'incidence était la plus élevée. Il n'y a donc pas d'histogramme.
Dans cette nouvelle décomposition en tâches, je vais calculer l'incidence annuelle pour chaque année séparément. Comme nous avons les données pour 33 ans, ceci fait 33 tâches à la place d'une seule dans la version que j'ai montrée avant. Pour ce calcul trivial, qui fait simplement la somme d'une cinquantaine d'entiers, ça n'a aucun sense. On va même perdre en efficacité, malgré le parallélisme. Mais il est facile d'imaginer un calcul plus compliqué à la place de cette simple sommme. Mon but ici est de montrer les techniques pour définir le workflow, qui serviront aussi pour mieux comprendre comment fonctionne =snakemake=.
Pour cette deuxième version, je crée un nouveau répertoire, et j'y fais une copie de tous les scripts, car la plupart ne nécessite pas de modification:
Pour cette deuxième version, je crée un nouveau répertoire, et j'y fais une copie du script de pré-traitement, qui reste identique:
#+begin_src sh :session *snakemake2* :results output :exports both
#+begin_src sh :session *snakemake2* :results output :exports both
Commençons en haut: j'ai mis la règle =all= au début pour pouvoir être paresseux à l'exécution: la simple commande =snakemake= déclenchera l'ensemble des calculs. Et =all=, c'est simplement l'histogramme des incidences annuelles ici. Tous les autres fichiers sont des résultats intermédiaires. Par simplicité, j'ai décidé de ne plus faire les plots de la suite chronologiques - vous l'avez vu assez de fois maintenant. Dans la règle =all=, je n'ai pas écrit le nom du fichier image de l'histgramme, mais à la place j'ai utilisé une référence indirecte, =rules.histogram.output=, ce qui veut dire tous les fichiers de sortie de la règle =histogram=. Ceci évite d'écrire le même nom deux fois, et potentiellement introduire des erreurs en le faisant.
Commençons en haut: j'ai mis la règle =all= au début pour pouvoir être paresseux à l'exécution: la simple commande =snakemake= déclenchera l'ensemble des calculs. Et =all=, c'est simplement le fichier qui résume les années du pic maximal pour chaque région.
Les deux règles suivantes, =download= et =preprocess=, sont les mêmes qu'avant, à un détail de notation près: dans =preprocess=, j'ai également utilisé la référence indirecte =rules.download.output= à la place du nom du fichier en question. Et j'ai rajouté le fichier du script comme fichier d'entrée pour que =snakemake= refasse le calcul après chaque modification du script.
Les trois règles finales sont la partie vraiment intéressante. Commençons par la dernière, =histogram=. Elle produit un fichier image, comme avant. Elle le fait en appelant un script en R, comme avant. Mais à l'entrée, elle réclame un fichier par an. L'expression =expand(...)= produit une liste de noms de fichier en remplaçant ={year}= dans le modèle donné par les éléments de la liste =range(1986, 2019)=, qui contient les entiers de =1986= à =2018=. Si cela vous rappele le langage Python, vous avez raison - le nom =snakemake= n'est pas une coïncidence, c'est écrit en Python !
La règles =histogram= réclame donc les fichiers =data/1986-incidence.txt=,
=data/1987-incidence.txt=, =data/1988-incidence.txt=, etc. Comme ces fichiers n'existent pas au départ, Une autre règle doit les produire. Et c'est la règle =annual_incidence=.qui le fait. En fait, elle s'applique à la production de tout fichier dont le nom à la forme =data/{year}-incidence.txt=. Quand =histogram= réclame le fichier =data/1986-incidence.txt=, =snakemake= trouve que la règle =annual_incidence= est applicable si on remplace ={year}= par =1986=. Il faut donc exécuter le script =scripts/annual_incidence.py= avec le fichier d'entrée =data/1986.csv=. Sauf que... ce fichier n'existe pas non plus. Pas grave, car la règle =extract_one_year= peut le produire! Il suffit d'appeler le script =scripts/extract_one_year.py= avec à l'entrée le fichier =rules.preprocess.output.data=, autrement dit le fichier =data/preprocessed-weekly-incidence.csv=, que fournit la règles =preprocess=.
La boucle qui fait le calcul pour chaque année est donc contenue dans la spécification des entrées de la règle =histogram=, qui en consomme le résultat. Et si vous regardez bien, c'est le principe de fonctionnement de =snakemake= partout: on demande un résultat, =snakemake= cherche la règle qui le produit, et applique cette règle après avoir assuré l'existence de ses entrée par exactement le même mécanisme. =snakemake= traite donc le workflow en commençant par la fin, les résultats, et en remontant vers les données d'entrée.
Il ne reste plus qu'à regarder les trois scripts réferencés dans les règles. Le premier, =scripts/extract_one_year.py=, lit le fichier des données pré-traitées et en extrait la part d'une année. Comme expliqué dans le module 3, la part de l'année N va du 1er août de l'année N-1 jusqu'au 31 juillet de l'année N et inclut ainsi le pic de l'année N qui se situe en janvier ou février. L'année est passée comme paramètre, défini dans la section =params= du =Snakefile= et récupéré en Python comme =snakemake.params=.
Un point important est la vérification de l'année. J'ai utilisé le nom suggestif =year= pour la partie variable des noms de fichier dans la règle =extract_one_year=, mais pour =snakemake=, ce n'est qu'un nom. Si je demande =snakemake data/foo.csv=, la même règle va être appliquée avec =foo= comme valeur de =year= ! Il faut donc que le script vérifie la validité du paramètre.
Dans la règle =download=, seul le nom du fichier de données a changé par rapport à avant. J'ai trouvé le nom du fichier "par région" sur le site Web du Réseau Sentinelles. C'est après qu'il y a le plus grand changement: la définition d'une variable =REGIONS=, qui est une liste des 13 régions administratives, dont les noms sont écrits exactement comme dans le fichier des données. On devrait récupérer cette liste du fichier de façon automatique, et je montrerai plus tard comment faire. Pour l'instant, je préfère copier la liste manuellement dans le =Snakefile= afin de ne pas introduire trop de nouveautés d'aun seul coup. La variable =REGIONS= est utilisée immédiatement après, pour définir les fichiers de sortie de la règle =split_by_region=. La fonction =expand= produit une liste des noms de fichier en insérant le nom de la région au bon endroit dans le modèle.
#+begin_src python :exports both :tangle incidence_syndrome_grippal_snakemake_parallele/scripts/extract_one_year.py
# Libraries used by this script:
import csv # for reading and writing CSV files
import os # for filename manipulation
# Read the CSV file into memory
with open(snakemake.input[0], "r") as csvfile:
data = []
csv_reader = csv.reader(csvfile)
for row in csv_reader:
data.append(row)
assert data[0] == ['week_starting', 'incidence']
del data[0]
# Get the year from the parameter object.
# Check that it is a valid year number, i.e. a four-digit integer.
year = snakemake.params.year
assert len(year) == 4
assert str(int(year)) == year
year = int(year)
# Check that we have data for that year.
# The dates are in ISO format, meaning that string
# comparison is equivalent to date comparison.
# There is thus no need to convert the date string
Le script =scripts/annual_incidence.py= fonctionnent d'après les mêmes principes. Il n'a pas besoin d'un paramètre pour indiquer l'année, car il n'en a pas besoin. Il lit un fichier CSV et fait la somme des nombres dans la deuxième colonne, c'est tout.
Le rôle de la règle =split_by_region= est de découper les données téléchargées en un fichier par région, afin de pouvoir traiter les régions en parallèle et avec les même scripts que nous avons déjà. Le script appliqué par la règle est assez simple:
#+begin_src python :exports both :tangle incidence_syndrome_grippal_snakemake_parallele/scripts/annual_incidence.py
#+begin_src python :exports both :tangle incidence_syndrome_grippal_par_region_snakemake/scripts/split-by-region.py
# Libraries used by this script:
import os
import csv # for reading CSV files
# Read the CSV file into memory
# Read the CSV file into memory
with open(snakemake.input[0], "r") as csvfile:
data = open(snakemake.input[0], 'rb').read()
data = []
# Decode the Latin-1 character set,
csv_reader = csv.reader(csvfile)
# remove white space at both ends,
for row in csv_reader:
# and split into lines.
data.append(row)
lines = data.decode('latin-1') \
assert data[0] == ['week_starting', 'incidence']
.strip() \
del data[0]
.split('\n')
# Compute total incidence
incidence = sum(int(row[1]) for row in data)
# Write the output file
with open(snakemake.output[0], "w") as output_file:
output_file.write(str(incidence))
output_file.write("\n")
#+end_src
Reste le script R qui fait l'histogramme. Rien à signaler, autre que peut-être la façon de lire tous les fichiers d'entrée avec une seule ligne de code avec la fonction =lapply=. A noter que =snakemake@input= est la liste des tous les fichiers d'entrée, y compris le nom du script lui-même, qu'il faut supprimer bien sûr.
#+begin_src R :exports both :tangle incidence_syndrome_grippal_snakemake_parallele/scripts/annual-incidence-histogram.R
# Read in the data. The last input file is the name of the script,
# so it needs to ne removed from the list before reading all the files.
files = snakemake@input
datafiles = files[-length(files)]
data = as.numeric(lapply(datafiles, function(fn) read.table(fn)[[1]]))
# Plot the histogram
# Separate header from data table
png(filename=snakemake@output[[1]])
comment = lines[0]
hist(data,
header = lines[1]
breaks=10,
table = [line.split(',') for line in lines[2:]]
xlab="Annual incidence",
ylab="Number of observations",
# Find all the regions mentioned in the table
main="")
regions = set(record[-1] for record in table)
dev.off()
# Write CSV files for each region
for region in regions:
filename = 'data/weekly-incidence-' + region + '.csv'
with open(filename, 'w') as output_file:
# The other scripts expect a comment in the first line,
# so write a minimal one to make them happy.
output_file.write('#\n')
output_file.write(header)
output_file.write('\n')
for record in table:
# Write only the records for right region
if record[-1] == region:
output_file.write(','.join(record))
output_file.write('\n')
#+end_src
#+end_src
Je pourrais maintenant taper =snakemake= et voir une longue liste de calculs défiler devant mes yeux. Je me retiens encore un peu pour illustrer ce que j'ai expliqué avant. En fait, je vais d'abord demander juste l'incidence annuelle de 2008:
Avant de continuer, voyons déjà ce que ça donne:
#+begin_src sh :session *snakemake2* :results output :exports both
#+begin_src sh :session *snakemake2* :results output :exports both
Les trois règles suivantes, =preprocess=, =plot=, et =annual_incidence= sont presques les mêmes qu'avant. Ce qui a changé, c'est la partie =-{region}= dans les noms des fichiers. Il faut interpréter le mot entre les accolades ("region") comme un nom de variable. La règle =preprocess=, par exemple, peut produire tout fichier qui a la forme "data/preprocessed-weekly-incidence-{region}.csv" si on lui donne le fichier "data/weekly-incidence-{region}.csv" avec la même valeur pour ={region}=. Etant donné les fichiers que nous avons obtenu par =split_by_region=, nous pouvons donc demander à snakemake le fichier "data/preprocessed-weekly-incidence-CORSE.csv", et snakemake va appliquer la règle =preprocess= au fichier d'entrée "data/weekly-incidence-CORSE.csv" que nous avons déjà. Faison-le:
rule annual_incidence:
input: data/2008.csv, scripts/annual_incidence.py
#+begin_src sh :session *snakemake2* :results output :exports both
On peut bien suivre l'exécution des tâches: d'abord =download=, puis =preprocess=, =extract_one_year=, et =annual_incidence=. Regardons le contenu de ce petit fichier:
#+begin_src sh :session *snakemake2* :results output :exports both
#+begin_src sh :session *snakemake2* :results output :exports both
cat data/2008-incidence.txt
ls data/preprocessed*
#+end_src
#+end_src
#+RESULTS:
#+RESULTS:
: 2975925
: data/preprocessed-weekly-incidence-CORSE.csv
Si maintenant je demande une autre année, seulement les tâches =extract_one_year= et =annual_incidence= devraient s'afficher, car les deux première sont déjà faites. Voyons:
Le script =preprocess.py= n'a d'ailleurs pas changé du tout. Un workflow permet donc de séparer la logistique de la gestion des données du code qui fait les calculs.
Le même mécanisme permet de demander l'incidence annuelle pour la Corse:
#+begin_src sh :session *snakemake2* :results output :exports both
#+begin_src sh :session *snakemake2* :results output :exports both
snakemake data/1992-incidence.txt
snakemake data/annual-incidence-CORSE.csv
#+end_src
#+end_src
#+RESULTS:
#+RESULTS:
...
@@ -1171,35 +1157,59 @@ Rules claiming more threads will be scaled down.
...
@@ -1171,35 +1157,59 @@ Rules claiming more threads will be scaled down.
Ça marche ! Je peux alors attaquer la totale, mais je vais supprimer l'affichage de tous les détails de l'exécution (option =-q=), pour éviter de voir les années défiler devant mes yeux!
Snakemake nous dit d'ailleurs explicitement quelle règle a été appliquée (=annual_incidence=), avec quel fichier d'entrée (=data/preprocessed-weekly-incidence-CORSE.csv=), et avec quel fichier de sortie (=data/annual-incidence-CORSE.csv=).
A la fin du workflow, il y a une nouvelle règle, =peak_years=, qui extrait l'année du pic maximal de chaque fichier d'incience annuelle, et produit un fichier résumant ces années par région. Sa seule particularité est la spécification des fichiers d'entrée, qui utilise la fonction =expand= exactement comme on l'a vu pour les fichiers résultats de la règle =split_by_region=. Le script Python associé est assez simple:
#+begin_src python :exports both :tangle incidence_syndrome_grippal_par_region_snakemake/scripts/peak-years.py
# Libraries used by this script:
import csv # for reading CSV files
import os # for path manipulations
with open(snakemake.output[0], 'w') as result_file:
for filename in snakemake.input:
region = '-'.join(os.path.splitext(filename)[0].split('-')[2:])
with open(filename, 'r') as csv_file:
csv_reader = csv.reader(csv_file)
csv_reader.__next__()
peak_year = None
peak_incidence = 0
for year, incidence in csv_reader:
incidence = int(incidence)
if incidence > peak_incidence:
peak_incidence = incidence
peak_year = year
result_file.write(region)
result_file.write(', ')
result_file.write(peak_year)
result_file.write('\n')
#+end_src
Dans ce workflow, nous avons donc introduit une boucle sur les régions en jouant avec les noms des fichiers. Chaque fichier du workflow précédent a été remplacé par une version "régionalisée", avec le suffix =-{region}= dans le nom. Ce qui déclenche la boucle, c'est la fonction =expand= dans notre =Snakefile=. Le grand avantage d'une telle boucle, par rapport à une boucle standard en Python ou R, est la parallélisation automatique. Sur une machine avec suffisamment de processeurs, toutes les 13 régions seront traitées en même temps. Mon ordinateur portable n'a qu'un processeur à deux coeurs, donc =snakemake= traite seulement deux régions à la foi.
Je vais maintenant lancer le calcul total - avec une petite précaution, l'option =-q= ("quiet") qui dit à snakemake d'être moins bavard:
#+begin_src sh :session *snakemake2* :results output :exports both
#+begin_src sh :session *snakemake2* :results output :exports both
En regardant bien le début du rapport que snakemake a fourni, on voit que =preprocess= et =annual_incidence= sont comptés 12 fois: une fois par région, moins la Corse que j'ai déjà traitée à la main. Une fois =all= et =peak_years=, ça a l'air bon. Et le résultat est là:
#+begin_src sh :session *snakemake2* :results output :exports both
cat data/peak-year-all-regions.txt
#+end_src
#+RESULTS:
#+begin_example
AUVERGNE-RHONE-ALPES, 2009
BOURGOGNE-FRANCHE-COMTE, 1986
BRETAGNE, 1996
CENTRE-VAL-DE-LOIRE, 1996
CORSE, 1989
GRAND EST, 2000
HAUTS-DE-FRANCE, 2013
ILE-DE-FRANCE, 1989
NORMANDIE, 1990
NOUVELLE-AQUITAINE, 1989
OCCITANIE, 2013
PAYS-DE-LA-LOIRE, 1989
PROVENCE-ALPES-COTE-D-AZUR, 1986
#+end_example
Un dernier détail à noter: la règle =plot= est bien dans mon =Snakefile=, mais elle n'a jamais été appliquée, et il n'y a aucun plot. C'est simplement parce que la règle =all= ne réclame que la production du fichier =data/peak-year-all-regions.txt=. J'aurais pu rajouter les plots, mais je ne l'ai pas fait. Ceci ne m'empêche pas de les demander explicitement:
#+begin_src sh :session *snakemake2* :results output :exports both
Enfin, je vais tenter de produire le dessin du graphe des tâches, m'attendant à une graphique un peu disproportionnée à cause du grand nombre de tâches:
Enfin, je vais tenter de produire le dessin du graphe des tâches, comme je l'ai fait avant pour un workflow nettement plus simple. Voyons...
#+begin_src sh :session *snakemake2* :results output :exports both
#+begin_src sh :session *snakemake2* :results output :exports both
snakemake -q --forceall --dag all | dot -Tpng > graph.png
snakemake -q --forceall --dag all | dot -Tpng > graph.png
