Accessibility Tools

Web Accessibility plugin by DJ-Extensions.com

- Le blog participatif de bioinformatique francophone depuis 2012 -

Snakemake aller plus loin avec la parallélisation

-

par

lelouar
dans

Bon­jour à tous, bien­ve­nue dans un nou­vel épi­sode de tuto­riels sur Sna­ke­make (épi­sode pré­cé­dent).

Aujourd'hui nous allons voir ensemble com­ment paral­lé­li­ser faci­le­ment par la don­née grâce à Sna­ke­make. L'idée géné­rale consiste à décou­per les fichiers bruts au début de notre pipe­line et de les ras­sem­bler après les étapes lourdes en cal­cul.

cluster_parallelisation

Nous allons éga­le­ment voir com­ment uti­li­ser le fichier de confi­gu­ra­tion au for­mat Json. Ce fichier est l'équivalent d'un dictionnaire/​tables de hachage et peut être uti­li­sé pour sto­cker les variables glo­bales ou les para­mètres uti­li­sés par les règles. Ain­si, il sera plus facile de géné­ra­li­ser votre work­flow et modi­fier les para­mètres des règles sans avoir à tou­cher au Sna­ke­file.

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" : "nar­row"
            },
            "IP_​B" : {
                "fastq" : "sample2",
                "fragment_​length" : "200",
                "peak_​type" : "broad"
            }
        }
    }
}

Pour l'utiliser dans un fichier Make­file il faut ajou­ter le mot clé :

1
config­file : "path/to/config.json"

On accède ensuite aux élé­ments comme si c'était un simple dico :

1
maVa­riable = config["maClé"]

Un seul mot clé pour paralléliser

Seule­ment un nou­veau mot clé (dyna­mic) et deux règles (décou­per et fusion­ner) sont néces­saire pour paral­lé­li­ser.

Je vais prendre pour exemple le work­flow du tuto pré­cé­dent. Ici l'étape limi­tante était bwa aln (rule bwa_​aln) car cette com­mande n'a pas d'option de paral­lé­li­sa­tion. Nous allons ain­si dépas­ser cette limite grâce aux deux règles sui­vantes :

1
2
3
4
5
6
7
8
9
rule split_​fastq :
input :
"{base}.fastq"
params :
lines_​per_​file =  config["lines_​per_​file"]
out­put :
dyna­mic("{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 :
dyna­mic("{base}.fastq.split.{split_id}.bam")
out­put :
"{base}.bam"
shell :
"sam­tools merge {out­put} {input}"

Ces deux règles seront à adap­ter si vous vou­lez les uti­li­ser dans le work­flow du tuto pré­cé­dent, ici j'ai juste vou­lu sim­pli­fier au maxi­mum pour mon­trer la puis­sance de cette fonc­tion­na­li­té.

Note : l'option –clus­ter per­met d'utiliser les com­mandes d'un ordon­nan­ceur (ex. –clus­ter 'qsub').

Aller plus loin dans l'automatisation

Le fichier configfile.json peut nous per­mettre d'automatiser la géné­ra­tion de fichiers cibles (ceux que l'on veut à la fin de notre work­flow). Dans l'exemple sui­vant, nous allons uti­li­ser le fichiers de confi­gu­ra­tion pré­sen­té 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 tar­get :
    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 work­flow avec la paral­lé­li­sa­tion :

Python
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
################################
####### map­ping avec BWA #######
################################
# Tous vos fichers fastq doivent être dans raw_data/xxx.fq.gz
#
# L'extension doit être exac­te­ment fq.gz
#
# Pour des don­nées pai­red-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
# out­put :
#     mapping/hg38/sampleX_c3-pe.sort.bam
#
# Dans le fichier config on elève le _​R1/​_​R2 :
# ex :
#   "data": {
#         "input" : {
#             "expX" : {
#                 "fastq" : "sam­pleX_c3-pe"
#             }
#         }
#   }
#############
# pro­ces­sus :
#     Index du génome de réfé­rence
#     Map­ping de deux jeux de don­nées en fonc­tion du fichier config (bwa)
#         Ali­gne­ment sur l'index du génome (bwa aln)
#         Trans­for­ma­tion des coor­don­nées de l'index en coor­don­nées
#             géno­miques (bwa sampe/​samse)
#     Ran­ge­ment des don­nées par coor­don­nées géno­mique (sam­tools sort)
#     Indexer le fichier ordon­née (sam­tools index)
#
################################
################################
 
# Pour récu­pé­rer la ver­sion des outils uti­li­sés (voir rule sort/​index)
import sub­pro­cess
 
 
# Ver­sion de bwa
bwa_​version = "0.7.10"
 
config­file : "path/to/config.json"
 
# on veri­fie que le nombre de ligne pour la découpe des fichier fastq
# est bien est un mul­tiple de 4
import sys
if config["lines_​per_​file"] %4 != 0 :
    sys.exit()
 
# Fichiers en sor­tie du pipe­line
rule tar­get :
    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 lec­tures (séquences)
rule bwa_​index :
    input :
        "{path_to_genome}.fa"
    out­put :
        "{path_to_genome}.fa.bwt"
    mes­sage :
        "Index du génome de réfé­rence : {wildcards.path_to_genome}"
    ver­sion : bwa_​version
    threads : 1
    shell :
        "bwa index {input}"
 
 
#décou­page du fichier fastq
rule split_​fastq :
    input :
        "{base}.fq.gz"
    out­put :
        dyna­mic("{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."
 
 
# Ali­gne­ment des lec­tures sur le génome de réfé­rence
rule bwa_​aln :
    input :
        bwt = lamb­da wild­cards : config["genome"][wild­cards.genome]["index"] + '.bwt',
        fas­ta = lamb­da wild­cards : config["genome"][wild­cards.genome]["index"],
        fq = "{exp}/raw_data/{samples}.fq.split.{split_id}"
    out­put :
        temp("{exp}/mapping/{genome}/{samples}.{split_id}.sai")
    ver­sion : bwa_​version
    threads : 20
    log :
        "{exp}/mapping/{genome}/{samples}.bwa_aln.log"
    shell :
        "bwa aln ‑t {threads} {input.fasta} {input.fq} >{out­put} 2>{log}"
 
 
# Si les don­nées sont appa­riées (pai­red-end), autre­ment dit si les frag­ments
# ont été séquen­cés dans les deux sens, cette règle sera uti­li­sée
rule bwa_​sampe_​to_​bam :
    input :
        fas­ta = "{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}"
    out­put :
        #fichier tem­po­raire
        temp("{exp}/mapping/{genome}/{samples}.{split_id}.bam")
    log :
        "{exp}/mapping/{genome}/{samples}.samse.{split_id}.log"
    threads : 2
    ver­sion : bwa_​version
    mes­sage :
        "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} | "
        "sam­tools view ‑Sbh — > {out­put}"
# Si on a séquen­cé les frag­ments dans un sens (Single read) cette règle sera uti­li­sée
rule bwa_​samse_​to_​bam :
    input :
        fas­ta = "{genome}.fa",
        sai = "{exp}/mapping/{genome}/{samples}.{split_id}.sai",
        fq = "{exp}/raw_data/{samples}.fq.split.{split_id}"
    out­put :
        #fichier tem­po­raire
        temp("{exp}/mapping/{genome}/{samples}.fq.split.{split_id}.bam")
    log :
        "{exp}/mapping/{genome}/{samples}.samse.{split_id}.log"
    threads : 2
    mes­sage :
        "SE : sai –> bam ({wildcards.genome}/{exp}/{wildcards.samples} : id : fq.split.{split_id})"
    shell :
        "bwa samse {input.fasta} {input.sai} {input.fq} "
        "2>{log} | "
        "sam­tools view ‑Sbh — > {out­put}"
 
 
# fusion de tout les fichiers bam
rule merge_​bam :
    input :
        dyna­mic("{base}.fq.split.{split_id}.bam")
    out­put :
        temp("{base}.bam")
    shell :
        "sam­tools merge {out­put} {input}"
 
 
# Pour ran­ger les don­nées par coor­don­nées géno­miques
rule sort_​bam :
    input :
        "{base}.bam"
    out­put :
        # fichier pro­té­gé en écri­ture
        pro­tec­ted("{base}.sort.bam")
    threads : 10
    log :
        "{base}.sort.log"
    mes­sage :
        "Ordonne le fichier {wildcards.base}.bam"
    ver­sion : # pour récu­pé­rer la ver­sion de l'outil avec une com­mande shell
        sub­pro­cess.getout­put(
            "sam­tools –ver­sion | "
            "head ‑1 | "
            "cut ‑d' ' ‑f2"
        )
    params :
        niveau_​compression = "9",
        memoire_​max_​par_​cpu = "1G"
    shell :
        "sam­tools sort "
        "-l {params.niveau_compression} "
        "-m {params.memoire_max_par_cpu} "
        "-@ {threads} "
        "-f {input} "
        "{out­put} "
        "2>{log}"
# Pour indexer le fichier trié (balise le fichier pour accé­der aux don­nées plus rapi­de­ment)
rule index_​bam :
    input :
        "{base}.sort.bam"
    out­put :
        "{base}.sort.bam.bai"
    prio­ri­ty : 5 # aug­mente la prio­ri­té de la règle (défaut = 1)
    ver­sion :
        sub­pro­cess.getout­put(
            "sam­tools –ver­sion | "
            "head ‑1 | "
            "cut ‑d' ' ‑f2"
        )
    threads : 1
    shell :
        "sam­tools index {input}"

Conclusion

Nous avons fran­chi un pas de plus dans l'utilisation de Sna­ke­make avec seule­ment un mot clé et deux nou­velles règles. Main­te­nant vous allez pou­voir faci­le­ment amé­lio­rer l'efficacité de vos work­flows en rédui­sant le temps de cal­cul des tâches les plus lourdes. Cepen­dant il faut gar­der en tête qu'une paral­lé­li­sa­tion à outrance n'est pas for­cé­ment une bonne stra­té­gie. Par exemple, si je décide de décou­per un fichier conte­nant 1000 lignes en 1000 fichiers d'une seule ligne et que j'ai seule­ment deux pauvres pro­ces­seurs à ma dis­po­si­tion, je risque une perte de temps et non un gain. C'est donc à vous de choi­sir judi­cieu­se­ment la stra­té­gie de paral­lé­li­sa­tion en fonc­tion de la machine uti­li­sée, de la taille des fichiers à décou­per, des logiciels/​scripts et du temps sup­plé­men­taire que vos deux nou­velles règles vont ajou­ter au work­flow.

Si vous avez des tâches très gour­mandes et un clus­ter de cal­cul à dis­po­si­tion sachez que le gain de temps peut être impres­sion­nant.

 

Je remer­cie Nol­wenn, NiGo­Pol et tadai­ma pour leur relec­ture et leur remarques.

Note : pour ce tuto­riel, je me suis ins­pi­ré de celui-ci : https://​github​.com/​i​n​o​d​b​/​s​n​a​k​e​m​a​k​e​-​p​a​r​a​l​l​e​l​-​bwa

Vous avez aimé ? Dites-le nous !

Moyenne : 0 /​ 5. Nb de votes : 0

Pas encore de vote pour cet article.

Partagez cet article

Lelouar

lelouar

Après un DUT génie biologique option bioinfo à Aurillac, j'ai poursuivi mes études avec une licence professionnelle (traitement des données génomiques) à Carcassonne puis un master (bioinformatique et biologie des systèmes) à Toulouse. J'ai ensuite eu une expérience d'un an à l'Institut de Génétique Humaine à Montpellier où je réalise actuellement une thèse.

Pour continuer la lecture :

Commentaires

Laisser un commentaire

Pour insérer du code dans vos commentaires, utilisez les balises <code> et <\code>.