Bonjour à tous, bienvenue dans un nouvel épisode de tutoriels sur Snakemake (épisode précédent).
Aujourd'hui nous allons voir ensemble comment paralléliser facilement par la donnée grâce à Snakemake. L'idée générale consiste à découper les fichiers bruts au début de notre pipeline et de les rassembler après les étapes lourdes en calcul.
Nous allons également voir comment utiliser le fichier de configuration au format Json. Ce fichier est l'équivalent d'un dictionnaire/tables de hachage et peut être utilisé pour stocker les variables globales ou les paramètres utilisés par les règles. Ainsi, il sera plus facile de généraliser votre workflow et modifier les paramètres des règles sans avoir à toucher au Snakefile.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 |
{ "lines_per_file": "10000", "genome": { "hg38" : { "index" : "path/to/genome.fasta", "gtf" : "path/to/gtf", "chr_len" : "path/to/genome.fasta.fai" } }, "data": { "input" : { "Input_X" : { "fastq" : "sample0", "fragment_length" : "200" } }, "IP" : { "IP_A" : { "fastq" : "sample1", "fragment_length" : "200", "peak_type" : "narrow" }, "IP_B" : { "fastq" : "sample2", "fragment_length" : "200", "peak_type" : "broad" } } } } |
Pour l'utiliser dans un fichier Makefile il faut ajouter le mot clé :
|
1 |
configfile: "path/to/config.json" |
On accède ensuite aux éléments comme si c'était un simple dico :
|
1 |
maVariable = config["maClé"] |
Un seul mot clé pour paralléliser
Seulement un nouveau mot clé (dynamic) et deux règles (découper et fusionner) sont nécessaire pour paralléliser.
Je vais prendre pour exemple le workflow du tuto précédent. Ici l'étape limitante était bwa aln (rule bwa_aln) car cette commande n'a pas d'option de parallélisation. Nous allons ainsi dépasser cette limite grâce aux deux règles suivantes :
|
1 2 3 4 5 6 7 8 9 |
rule split_fastq: input: "{base}.fastq" params: lines_per_file = config["lines_per_file"] output: dynamic("{base}.fastq.split.{split_id}") shell: "split -d -l {params.lines_per_file} {input} {base}.fastq.split." |
|
1 2 3 4 5 6 7 |
rule merge_bam: input: dynamic("{base}.fastq.split.{split_id}.bam") output: "{base}.bam" shell: "samtools merge {output} {input}" |
Ces deux règles seront à adapter si vous voulez les utiliser dans le workflow du tuto précédent, ici j'ai juste voulu simplifier au maximum pour montrer la puissance de cette fonctionnalité.
Note : l'option --cluster permet d'utiliser les commandes d'un ordonnanceur (ex. --cluster 'qsub').
Aller plus loin dans l'automatisation
Le fichier configfile.json peut nous permettre d'automatiser la génération de fichiers cibles (ceux que l'on veut à la fin de notre workflow). Dans l'exemple suivant, nous allons utiliser le fichiers de configuration présenté plus tôt pour générer les fichiers cibles. Note : {exp} et {samples} vont de pairs.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 |
rule target: input: [ expand( "{exp}/mapping/{genome}/{samples}.sorted.{bam}", exp=EXP, genome=config["genome"].keys(), samples=config["data"][cond][EXP]["fastq"], bam = ["bam", "bam.bai"] ) for cond in config["data"].keys() for EXP in config["data"][cond].keys() ], |
Exemple du workflow avec la parallélisation :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 |
################################ ####### mapping avec BWA ####### ################################ # Tous vos fichers fastq doivent être dans raw_data/xxx.fq.gz # # L'extension doit être exactement fq.gz # # Pour des données paired-end (2 fichiers : _R1/_R2): # ex : # input: # hg38.fa # expX/raw_data/sampleX_c3-pe_R1.fq.gz # expX/raw_data/sampleX_c3-pe_R2.fq.gz # output: # mapping/hg38/sampleX_c3-pe.sort.bam # # Dans le fichier config on elève le _R1/_R2: # ex: # "data": { # "input" : { # "expX" : { # "fastq" : "sampleX_c3-pe" # } # } # } ############# # processus : # Index du génome de référence # Mapping de deux jeux de données en fonction du fichier config (bwa) # Alignement sur l'index du génome (bwa aln) # Transformation des coordonnées de l'index en coordonnées # génomiques (bwa sampe/samse) # Rangement des données par coordonnées génomique (samtools sort) # Indexer le fichier ordonnée (samtools index) # ################################ ################################ # Pour récupérer la version des outils utilisés (voir rule sort/index) import subprocess # Version de bwa bwa_version = "0.7.10" configfile: "path/to/config.json" # on verifie que le nombre de ligne pour la découpe des fichier fastq # est bien est un multiple de 4 import sys if config["lines_per_file"] %4 != 0: sys.exit() # Fichiers en sortie du pipeline rule target: input: [ expand( "{exp}/mapping/{genome}/{samples}.sort.{bam}", exp=EXP, genome=config["genome"].keys(), samples=config["data"][cond][EXP]["fastq"], bam = ["bam", "bam.bai"] ) for cond in config["data"].keys() for EXP in config["data"][cond].keys() ] # Index du génome de référence avant l'alignement des lectures (séquences) rule bwa_index: input: "{path_to_genome}.fa" output: "{path_to_genome}.fa.bwt" message: "Index du génome de référence : {wildcards.path_to_genome}" version: bwa_version threads: 1 shell: "bwa index {input}" #découpage du fichier fastq rule split_fastq: input: "{base}.fq.gz" output: dynamic("{base}.fq.split.{split_id}") params: lines_per_file = config["lines_per_file"] shell: "zcat {input} | split -d -l {params.lines_per_file} - {base}.fq.split." # Alignement des lectures sur le génome de référence rule bwa_aln: input: bwt = lambda wildcards: config["genome"][wildcards.genome]["index"] + '.bwt', fasta = lambda wildcards: config["genome"][wildcards.genome]["index"], fq = "{exp}/raw_data/{samples}.fq.split.{split_id}" output: temp("{exp}/mapping/{genome}/{samples}.{split_id}.sai") version: bwa_version threads: 20 log: "{exp}/mapping/{genome}/{samples}.bwa_aln.log" shell: "bwa aln -t {threads} {input.fasta} {input.fq} >{output} 2>{log}" # Si les données sont appariées (paired-end), autrement dit si les fragments # ont été séquencés dans les deux sens, cette règle sera utilisée rule bwa_sampe_to_bam: input: fasta = "{genome}.fa", sai1 = "{exp}/mapping/{genome}/{samples}_R1.{split_id}.sai", fq1 = "{exp}/raw_data/{samples}_R1.fq.split.{split_id}", sai2 = "{exp}/mapping/{genome}/{samples}_R2.{split_id}.sai", fq2 = "{exp}/raw_data/{samples}_R2.fq.split.{split_id}" output: #fichier temporaire temp("{exp}/mapping/{genome}/{samples}.{split_id}.bam") log: "{exp}/mapping/{genome}/{samples}.samse.{split_id}.log" threads: 2 version: bwa_version message: "PE : sai --> bam ({wildcards.genome}/{exp}/{wildcards.samples}, id: {split_id})" shell: "bwa sampe " "{input.fasta} " "{input.sai1} {input.sai2} " "{input.fq1} {input.fq2} " "2>{log} | " "samtools view -Sbh - > {output}" # Si on a séquencé les fragments dans un sens (Single read) cette règle sera utilisée rule bwa_samse_to_bam: input: fasta = "{genome}.fa", sai = "{exp}/mapping/{genome}/{samples}.{split_id}.sai", fq = "{exp}/raw_data/{samples}.fq.split.{split_id}" output: #fichier temporaire temp("{exp}/mapping/{genome}/{samples}.fq.split.{split_id}.bam") log: "{exp}/mapping/{genome}/{samples}.samse.{split_id}.log" threads: 2 message: "SE : sai --> bam ({wildcards.genome}/{exp}/{wildcards.samples} : id: fq.split.{split_id})" shell: "bwa samse {input.fasta} {input.sai} {input.fq} " "2>{log} | " "samtools view -Sbh - > {output}" # fusion de tout les fichiers bam rule merge_bam: input: dynamic("{base}.fq.split.{split_id}.bam") output: temp("{base}.bam") shell: "samtools merge {output} {input}" # Pour ranger les données par coordonnées génomiques rule sort_bam: input: "{base}.bam" output: # fichier protégé en écriture protected("{base}.sort.bam") threads: 10 log: "{base}.sort.log" message: "Ordonne le fichier {wildcards.base}.bam" version: # pour récupérer la version de l'outil avec une commande shell subprocess.getoutput( "samtools --version | " "head -1 | " "cut -d' ' -f2" ) params: niveau_compression = "9", memoire_max_par_cpu = "1G" shell: "samtools sort " "-l {params.niveau_compression} " "-m {params.memoire_max_par_cpu} " "-@ {threads} " "-f {input} " "{output} " "2>{log}" # Pour indexer le fichier trié (balise le fichier pour accéder aux données plus rapidement) rule index_bam: input: "{base}.sort.bam" output: "{base}.sort.bam.bai" priority: 5 # augmente la priorité de la règle (défaut = 1) version: subprocess.getoutput( "samtools --version | " "head -1 | " "cut -d' ' -f2" ) threads: 1 shell: "samtools index {input}" |
Conclusion
Nous avons franchi un pas de plus dans l'utilisation de Snakemake avec seulement un mot clé et deux nouvelles règles. Maintenant vous allez pouvoir facilement améliorer l'efficacité de vos workflows en réduisant le temps de calcul des tâches les plus lourdes. Cependant il faut garder en tête qu'une parallélisation à outrance n'est pas forcément une bonne stratégie. Par exemple, si je décide de découper un fichier contenant 1000 lignes en 1000 fichiers d'une seule ligne et que j'ai seulement deux pauvres processeurs à ma disposition, je risque une perte de temps et non un gain. C'est donc à vous de choisir judicieusement la stratégie de parallélisation en fonction de la machine utilisée, de la taille des fichiers à découper, des logiciels/scripts et du temps supplémentaire que vos deux nouvelles règles vont ajouter au workflow.
Si vous avez des tâches très gourmandes et un cluster de calcul à disposition sachez que le gain de temps peut être impressionnant.
Je remercie Nolwenn, NiGoPol et tadaima pour leur relecture et leur remarques.
Note : pour ce tutoriel, je me suis inspiré de celui-ci : https://github.com/inodb/snakemake-parallel-bwa
