diff --git a/module6/ressources/incidence_syndrome_grippal_snakemake/Snakefile b/module6/ressources/incidence_syndrome_grippal_snakemake/Snakefile index 59feb8a3f4abbec69dcd3b50bd5aab6bedac9211..40bd400bdfb707e5435e171e565b2fd8ea3f362b 100644 --- a/module6/ressources/incidence_syndrome_grippal_snakemake/Snakefile +++ b/module6/ressources/incidence_syndrome_grippal_snakemake/Snakefile @@ -9,7 +9,7 @@ rule preprocess: "data/weekly-incidence.csv" output: data="data/preprocessed-weekly-incidence.csv", - errorlog="data/errors-from-preprocessing.csv" + errorlog="data/errors-from-preprocessing.txt" script: "scripts/preprocess.py" diff --git a/module6/ressources/incidence_syndrome_grippal_snakemake/scripts/preprocess.py b/module6/ressources/incidence_syndrome_grippal_snakemake/scripts/preprocess.py index 056b006083b98a06b97cbe66c58d34f14cb1e0ca..c3aef7bb28f38f5b7b0e19adbd366cfc2c32c3a2 100644 --- a/module6/ressources/incidence_syndrome_grippal_snakemake/scripts/preprocess.py +++ b/module6/ressources/incidence_syndrome_grippal_snakemake/scripts/preprocess.py @@ -64,5 +64,6 @@ with open(snakemake.output.errorlog, "a") as errorlog: # 2. the incidence estimate for that week with open(snakemake.output.data, "w") as csvfile: csv_writer = csv.writer(csvfile) + csv_writer.writerow(["week_starting", "incidence"]) for row in converted_data: csv_writer.writerow(row) diff --git a/module6/ressources/snakemake_tutorial_fr.org b/module6/ressources/snakemake_tutorial_fr.org index bbe23731eabb16ca9e3457c776b88523dd7d16b8..240f363d00c7d5d66755f7a650f05e092a3015c7 100644 --- a/module6/ressources/snakemake_tutorial_fr.org +++ b/module6/ressources/snakemake_tutorial_fr.org @@ -8,6 +8,497 @@ * Installer snakemake TODO * L'analyse de l'incidence du syndrome grippal revisitée -Nous allons reprendre l'exemple du module 3, l'analyse de l'incidence du syndrome grippal. -** 1ère étape: le téléchargement des données +Je vais reprendre l'exemple du module 3, l'analyse de l'incidence du syndrome grippal, et je vais refaire exactement la même analyse sous forme d'un workflow par =snakemake=. Ceci veut dire que pour l'instant, nous quittons le monde des documents computationnels que nous vous avons montré dans les modules 2 et 3, pour passer dans l'univers de la ligne de commande. Il y a des bonnes raisons pour cela, que je vous donnerai plus tard. Et vous verrez aussi le retour des documents computationnels à la fin de ce tutoriel, même si ce sera dans un rôle moins central. + +Un workflow est composé de tâches dont chacun correspond à un bout du calcul total. Une tâche est typiquement l'exécution d'une commande ou d'un script. Les tâches communiques entre eux par des fichiers - au moins dans la vision de =snakemake= (et d'autres descendants de =make=). Pour faire le lien avec la programmation dans un langage comme Python ou R, une tâche est l'appel à une fonction, et les paramètres et les valeurs de retour sont stockés dans des fichiers. + +Il y a beaucoup de liberté dans la décomposition d'un calcul en tâches d'un workflow. Souvent les critères sont plutôt techniques que scientifiques: une tâche peut alors correspondre à l'exécution d'un logiciel, ou à une étape du calcul qui est fait sur un ordinateur précis. Pour mon analyse je propose la décomposition suivante, qui est assez arbitraire : + 1. Téléchargement des données du site du Réseau Sentinelles + 2. Pré-traitement: extraction des données utilisées, vérifications + 3. Visualisation: génération des plots + 4. Calcul des incidences annuelles + 5. Calcul de l'histogramme des incidences annuelles + +Pour faire les calculs, je vais recycler le code du module 3, sans les commenter de nouveau ici. +** Préparation +Un workflow finit par utiliser beaucoup de fichiers, donc je commence par la création d'un répertoire qui contient tout: +#+begin_src sh :session *snakemake* :results output :exports both +mkdir incidence_syndrome_grippal_snakemake +cd incidence_syndrome_grippal_snakemake +#+end_src + +#+RESULTS: + +En plus, je crée un répertoire pour les fichiers de données, +#+begin_src sh :session *snakemake* :results output :exports both +mkdir data +#+end_src + +#+RESULTS: + +et un autre pour les scripts en Python ou R: +#+begin_src sh :session *snakemake* :results output :exports both +mkdir scripts +#+end_src + +#+RESULTS: + +** 1ère tâche: le téléchargement des données +Pour télécharger un fichier, inutile d'écrire du code: l'utilitaire =wget= fait ce qu'il faut. La ligne de commande +#+begin_src sh :session *snakemake* :results output :exports both +wget -O data/weekly-incidence.csv http://www.sentiweb.fr/datasets/incidence-PAY-3.csv +#+end_src + +#+RESULTS: +: --2019-08-28 16:20:33-- http://www.sentiweb.fr/datasets/incidence-PAY-3.csv +: Resolving www.sentiweb.fr (www.sentiweb.fr)... 134.157.220.17 +: Connecting to www.sentiweb.fr (www.sentiweb.fr)|134.157.220.17|:80... connected. +: HTTP request sent, awaiting response... 200 OK +: Length: unspecified [text/csv] +: Saving to: 'data/weekly-incidence.csv' +: ] 0 --.-KB/s data/weekly-inciden [ <=> ] 79.88K --.-KB/s in 0.05s +: +: 2019-08-28 16:20:34 (1.48 MB/s) - 'data/weekly-incidence.csv' saved [81800] + +fait ce qu'il faut, et dépose les données dans le fichier =data/weekly-incidence.csv=. Je le supprime parce que je veux faire le téléchargement dans mon workflow! +#+begin_src sh :session *snakemake* ::results output :exports both +rm data/weekly-incidence.csv +#+end_src + +#+RESULTS: + +Je vais commencer la rédaction du =Snakefile=, le fichier qui déinit mon workflow: +#+begin_src :exports both :tangle incidence_syndrome_grippal_snakemake/Snakefile +rule download: + output: + "data/weekly-incidence.csv" + shell: + "wget -O {output} http://www.sentiweb.fr/datasets/incidence-PAY-3.csv" +#+end_src +Un =Snakefile= consiste de /règles/ qui définissent les tâches. Chaque règle a un nom, ici j'ai choisi /download/. Une règle liste aussi les fichiers d'entrée (aucun dans ce cas) et de sortie (notre fichier de données). Enfin, il faut dire ce qui est à faire pour exécuter la tâche, ce qui est ici la commande =wget=. + +Pour exécuter cette tâche, il y a deux façons de faire: on peut demander à =snakemake= d'exécuter la règle =download=: +#+begin_src sh :session *snakemake* ::results output :exports both +snakemake download +#+end_src + +#+RESULTS: +| Building | DAG | of | jobs... | | | | | | | | +| Using | shell: | /bin/bash | | | | | | | | | +| Provided | cores: | 1 | | | | | | | | | +| Rules | claiming | more | threads | will | be | scaled | down. | | | | +| Job | counts: | | | | | | | | | | +| | count | jobs | | | | | | | | | +| | 1 | download | | | | | | | | | +| | 1 | | | | | | | | | | +| [Wed | Aug | 28 | 16:38:03 | 2019] | | | | | | | +| rule | download: | | | | | | | | | | +| output: | data/weekly-incidence.csv | | | | | | | | | | +| jobid: | 0 | | | | | | | | | | +| --2019-08-28 | 16:38:03-- | http://www.sentiweb.fr/datasets/incidence-PAY-3.csv | | | | | | | | | +| Resolving | www.sentiweb.fr | (www.sentiweb.fr)... | 134.157.220.17 | | | | | | | | +| Connecting | to | www.sentiweb.fr | (www.sentiweb.fr) | 134.157.220.17 | :80... | connected. | | | | | +| HTTP | request | sent, | awaiting | response... | 200 | OK | | | | | +| Length: | unspecified | [text/csv] | | | | | | | | | +| Saving | to: | 'data/weekly-incidence.csv' | | | | | | | | | +| ] | 0 | --.-KB/s | data/weekly-inciden | [ | <=> | ] | 79.88K | --.-KB/s | in | 0.02s | +| 2019-08-28 | 16:38:03 | (3.71 | MB/s) | 0 | 'data/weekly-incidence.csv' | saved | [81800] | | | | +| [Wed | Aug | 28 | 16:38:03 | 2019] | | | | | | | +| Finished | job | 0 | | | | | | | | | +| ) | done | | | | | | | | | | +| Complete | log: | /home/hinsen/projects/RR_MOOC/repos-session02/mooc-rr-ressources/module6/ressources/incidence_syndrome_grippal_snakemake/.snakemake/log/2019-08-28T163803.322303.snakemake.log | | | | | | | | | + +ou encore on peut lui demander de produire le résultat souhaité: +#+begin_src sh :session *snakemake* ::results output :exports both +snakemake data/weekly-incidence.csv +#+end_src + +#+RESULTS: +| Building | DAG | of | jobs... | +| Nothing | to | be | done. | +| Complete | log: | /home/hinsen/projects/RR_MOOC/repos-session02/mooc-rr-ressources/module6/ressources/incidence_syndrome_grippal_snakemake/.snakemake/log/2019-08-28T163822.540694.snakemake.log | | + +En regardant bien ce que =snakemake= dit au deuxième tour, il s'est rendu compte qu'il n'y a rien à faire, parce que le fichier souhaité existe déjà. Voici un premier avantage important d'un workflow: une tâche n'est exécutée que s'il est nécessaire. Quand une tâche met deux heures à exécuter, c'est appréciable. + +** 2ème tâche: le pré-traitement des données +La deuxième tâche est le pré-traitement: en partant du fichier téléchargé du Réseau Sentinelle, il faut extraire juste les éléments nécessaires, et il faut vérifier s'il y a des données manquantes ou des erreurs. Dans un document computationnel, j'avais procédé pas par pas, en inspectant les résultats à chaque étape. Dans mon workflow, le pré-traitement devient une seule tâche, exécutée en bloc. + +Il faut donc bien réfléchir à ce qu'on attend comme résultat. En fait, il faut deux fichiers de sortie: un qui contient les données qui seront analysées par la suite, et un autre qui contient les éventuels messages d'erreur. Avec ça, la deuxième règle s'écrit assez vite: +#+begin_src :exports both :tangle incidence_syndrome_grippal_snakemake/Snakefile +rule preprocess: + input: + "data/weekly-incidence.csv" + output: + data="data/preprocessed-weekly-incidence.csv", + errorlog="data/errors-from-preprocessing.txt" + script: + "scripts/preprocess.py" +#+end_src +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=. + +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 +# Libraries used by this script: +import datetime # for date conversion +import csv # for writing output to a CSV file + +# Read the CSV file into memory +data = open(snakemake.input[0], 'rb').read() +# Decode the Latin-1 character set, +# remove white space at both ends, +# and split into lines. +lines = data.decode('latin-1') \ + .strip() \ + .split('\n') +# Discard the first line, which contains a comment +data_lines = lines[1:] +# Split each line into columns +table = [line.split(',') for line in data_lines] + +# Remove records with missing data and write +# the removed records to a separate file for inspection. +with open(snakemake.output.errorlog, "w") as errorlog: + valid_table = [] + for row in table: + missing = any([column == '' for column in row]) + if missing: + errorlog.write("Missing data in record\n") + errorlog.write(str(row)) + errorlog.write("\n") + else: + valid_table.append(row) + +# Extract the two relevant columns, "week" and "inc" +week = [row[0] for row in valid_table] +assert week[0] == 'week' +del week[0] +inc = [row[2] for row in valid_table] +assert inc[0] == 'inc' +del inc[0] +data = list(zip(week, inc)) + +# Check for obviously out-of-range values +with open(snakemake.output.errorlog, "a") as errorlog: + for week, inc in data: + if len(week) != 6 or not week.isdigit(): + errorlog.write("Suspect value in column 'week': {week}\n") + if not inc.isdigit(): + errorlog.write("Suspect value in column 'inc': {inc}\n") + +# Convert year/week by date of the corresponding Monday, +# then sort by increasing date +converted_data = \ + [(datetime.datetime.strptime(year_and_week + ":1" , '%G%V:%u').date(), inc) + for year_and_week, inc in data] +converted_data.sort(key = lambda record: record[0]) + +# Check that consecutive dates are seven days apart +with open(snakemake.output.errorlog, "a") as errorlog: + dates = [date for date, _ in converted_data] + for date1, date2 in zip(dates[:-1], dates[1:]): + if date2-date1 != datetime.timedelta(weeks=1): + errorlog.write(f"{date2-date1} between {date1} and {date2}\n") + +# Write data to a CSV file with two columns: +# 1. the date of the Monday of each week, in ISO format +# 2. the incidence estimate for that week +with open(snakemake.output.data, "w") as csvfile: + csv_writer = csv.writer(csvfile) + csv_writer.writerow(["week_starting", "incidence"]) + for row in converted_data: + csv_writer.writerow(row) +#+end_src + +Ce qui saute aux yeux d'abord, c'est =snakemake.input[0]= comme nom de fichier. Le nom =snakemake= semble venir de nulle part: il n'est ni défini dans le script, ni importé d'un module. En fait, c'est bien =snakemake= qui définit ce nom dans l'interprète Python avant de lancer le script. Il permet d'accèder aux définitions du =Snakefile=, et notamment aux noms des fichiers. + +Sinon, il y a deux modifications par rapport au code du module 3. Premièrement, les messages d'erreurs sont écrits dans un fichier. Deuxièmement, les données finales sont écrites également dans un fichier, en utilisant le format CSV. + +Pour appliquer le pré-traitement, demandons à =snakemake=: +#+begin_src sh :session *snakemake* :results output :exports both +snakemake preprocess +#+end_src + +#+RESULTS: +#+begin_example +Building DAG of jobs... +Using shell: /bin/bash +Provided cores: 1 +Rules claiming more threads will be scaled down. +Job counts: + count jobs + 1 preprocess + 1 + +[Wed Aug 28 18:38:33 2019] +rule preprocess: + input: data/weekly-incidence.csv + output: data/preprocessed-weekly-incidence.csv, data/errors-from-preprocessing.txt + jobid: 0 + +[Wed Aug 28 18:38:34 2019] +Finished job 0. +) done +Complete log: /home/hinsen/projects/RR_MOOC/repos-session02/mooc-rr-ressources/module6/ressources/incidence_syndrome_grippal_snakemake/.snakemake/log/2019-08-28T183833.807546.snakemake.log +#+end_example + +Voyons s'il y a eu des problèmes: +#+begin_src sh :session *snakemake* :results output :exports both +cat data/errors-from-preprocessing.txt +#+end_src + +#+RESULTS: +: Missing data in record +: ['198919', '3', '0', '', '', '0', '', '', 'FR', 'France'] +: 14 days, 0:00:00 between 1989-05-01 and 1989-05-15 + +En effet, on avait vu dans le module 3 qu'il y a un point manquant dans ce jeu de données. Quant aux données, je vais afficher juste le début: +#+begin_src sh :session *snakemake* :results output :exports both +head -10 data/preprocessed-weekly-incidence.csv +#+end_src + +#+RESULTS: +#+begin_example +week_starting,incidence +1984-10-29,68422 +1984-11-05,135223 +1984-11-12,87330 +1984-11-19,72029 +1984-11-26,78620 +1984-12-03,101073 +1984-12-10,123680 +1984-12-17,101726 +1984-12-24,84830 +#+end_example + +Ça a l'air pas mal! +** 3ème tâche: préparer les plots +La règle pour faire les plots ne présente plus aucune surprise: +#+begin_src :exports both :tangle incidence_syndrome_grippal_snakemake/Snakefile +rule plot: + input: + "data/preprocessed-weekly-incidence.csv" + output: + "data/weekly-incidence-plot.png", + "data/weekly-incidence-plot-last-years.png" + script: + "scripts/incidence-plots.R" +#+end_src +Il y a les données pré-traitées à l'entrée, et deux fichiers image à la sortie, créées par un script, cette fois en langage R: +#+begin_src R :exports both :tangle incidence_syndrome_grippal_snakemake/scripts/incidence-plots.R +# Read in the data and convert the dates +data = read.csv(snakemake@input[[1]]) +data$week_starting <- as.Date(data$week_starting) + +# Plot the complete incidence dataset +png(filename=snakemake@output[[1]]) +plot(data, type="l", xlab="Date", ylab="Weekly incidence") +dev.off() + +# Zoom on the last four years +png(filename=snakemake@output[[2]]) +plot(tail(data, 4*52), type="l", xlab="Date", ylab="Weekly incidence") +dev.off() +#+end_src +Comme pour le script Python de l'étape précedente, l'accès aux noms des fichiers se fait par le nom =snakemake= qui est créé par... =snakemake=. + +Passons à l'exécution: +#+begin_src sh :session *snakemake* :results output :exports both +snakemake plot +#+end_src + +#+RESULTS: +#+begin_example +Building DAG of jobs... +Using shell: /bin/bash +Provided cores: 1 +Rules claiming more threads will be scaled down. +Job counts: + count jobs + 1 plot + 1 + +[Wed Aug 28 19:43:42 2019] +rule plot: + input: data/preprocessed-weekly-incidence.csv + output: data/weekly-incidence-plot.png, data/weekly-incidence-plot-last-years.png + jobid: 0 + +During startup - Warning messages: +1: Setting LC_COLLATE failed, using "C" +2: Setting LC_TIME failed, using "C" +3: Setting LC_MESSAGES failed, using "C" +4: Setting LC_MONETARY failed, using "C" +Warning message: +Y-%m-%d", tz = "GMT") : + unknown timezone 'zone/tz/2019b.1.0/zoneinfo/Europe/Paris' +null device + 1 +null device + 1 +[Wed Aug 28 19:43:43 2019] +Finished job 0. +) done +Complete log: /home/hinsen/projects/RR_MOOC/repos-session02/mooc-rr-ressources/module6/ressources/incidence_syndrome_grippal_snakemake/.snakemake/log/2019-08-28T194342.688492.snakemake.log +#+end_example + +Voici les deux plots: + +[[file:incidence_syndrome_grippal_snakemake/data/weekly-incidence-plot.png]] + +[[file:incidence_syndrome_grippal_snakemake/data/weekly-incidence-plot-last-years.png]] + +** 4ème tâche: calculer l'incidence annuelle +Écrire les règles pour =snakemake= devient vite une routine: +#+begin_src :exports both :tangle incidence_syndrome_grippal_snakemake/Snakefile +rule annual_incidence: + input: + "data/preprocessed-weekly-incidence.csv" + output: + "data/annual-incidence.csv" + script: + "scripts/annual-incidence.R" +#+end_src +Et le script en langage R ressemble fortement au code du module 3: +#+begin_src R :exports both :tangle incidence_syndrome_grippal_snakemake/scripts/annual-incidence.R +# Read in the data and convert the dates +data = read.csv(snakemake@input[[1]]) +names(data) <- c("date", "incidence") +data$date <- as.Date(data$date) + +# A function that extracts the peak for year N +yearly_peak = function(year) { + start = paste0(year-1,"-08-01") + end = paste0(year,"-08-01") + records = data$date > start & data$date <= end + sum(data$incidence[records]) + } + +# The years for which we have the full peak +years <- 1986:2018 + +# Make a new data frame for the annual incidences +annual_data = data.frame(year = years, + incidence = sapply(years, yearly_peak)) +# write output file +write.csv(annual_data, + file=snakemake@output[[1]], + row.names=FALSE) +#+end_src + +Allons-y! +#+begin_src sh :session *snakemake* :results output :exports both +snakemake annual_incidence +#+end_src + +#+RESULTS: +#+begin_example +Building DAG of jobs... +Using shell: /bin/bash +Provided cores: 1 +Rules claiming more threads will be scaled down. +Job counts: + count jobs + 1 annual_incidence + 1 + +[Wed Aug 28 19:51:26 2019] +rule annual_incidence: + input: data/preprocessed-weekly-incidence.csv + output: data/annual-incidence.csv + jobid: 0 + +During startup - Warning messages: +1: Setting LC_COLLATE failed, using "C" +2: Setting LC_TIME failed, using "C" +3: Setting LC_MESSAGES failed, using "C" +4: Setting LC_MONETARY failed, using "C" +Warning message: +Y-%m-%d", tz = "GMT") : + unknown timezone 'zone/tz/2019b.1.0/zoneinfo/Europe/Paris' +[Wed Aug 28 19:51:26 2019] +Finished job 0. +) done +Complete log: /home/hinsen/projects/RR_MOOC/repos-session02/mooc-rr-ressources/module6/ressources/incidence_syndrome_grippal_snakemake/.snakemake/log/2019-08-28T195126.602879.snakemake.log +#+end_example + +Voyons le début du résultat: +#+begin_src sh :session *snakemake* :results output :exports both +head -10 data/annual-incidence.csv +#+end_src + +#+RESULTS: +#+begin_example +"year","incidence" +1986,5100540 +1987,2861556 +1988,2766142 +1989,5460155 +1990,5233987 +1991,1660832 +1992,2576347 +1993,2703708 +1994,3515735 +#+end_example + +** 5ème tâche: l'histogramme +Et pour finir, encore un petit script en R: +#+begin_src :exports both :tangle incidence_syndrome_grippal_snakemake/Snakefile +rule histogram: + input: + "data/annual-incidence.csv" + output: + "data/annual-incidence-histogram.png" + script: + "scripts/annual-incidence-histogram.R" +#+end_src + +#+begin_src R :exports both :tangle incidence_syndrome_grippal_snakemake/scripts/annual-incidence-histogram.R +# Read in the data and convert the dates +data = read.csv(snakemake@input[[1]]) + +# Plot the histogram +png(filename=snakemake@output[[1]]) +hist(data$incidence, + breaks=10, + xlab="Annual incidence", + ylab="Number of observations", + main="") +dev.off() +#+end_src + +#+begin_src sh :session *snakemake* :results output :exports both +snakemake histogram +#+end_src + +#+RESULTS: +#+begin_example +Building DAG of jobs... +Using shell: /bin/bash +Provided cores: 1 +Rules claiming more threads will be scaled down. +Job counts: + count jobs + 1 histogram + 1 + +[Wed Aug 28 19:54:24 2019] +rule histogram: + input: data/annual-incidence.csv + output: data/annual-incidence-histogram.png + jobid: 0 + +During startup - Warning messages: +1: Setting LC_COLLATE failed, using "C" +2: Setting LC_TIME failed, using "C" +3: Setting LC_MESSAGES failed, using "C" +4: Setting LC_MONETARY failed, using "C" +null device + 1 +[Wed Aug 28 19:54:24 2019] +Finished job 0. +) done +Complete log: /home/hinsen/projects/RR_MOOC/repos-session02/mooc-rr-ressources/module6/ressources/incidence_syndrome_grippal_snakemake/.snakemake/log/2019-08-28T195424.640999.snakemake.log +#+end_example + +[[file:incidence_syndrome_grippal_snakemake/data/annual-incidence-histogram.png]]