-
Notifications
You must be signed in to change notification settings - Fork 15
/
nf_bisulfite_PBAT_DirtyHarry
executable file
·374 lines (289 loc) · 20.8 KB
/
nf_bisulfite_PBAT_DirtyHarry
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
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
#!/usr/bin/env nextflow
nextflow.enable.dsl=2
// Last modified 24 02 2021
// This pipeline will first run paired-end alignments, with --unmapped specified, and then
// run single-end alignments afterwards
params.outdir = "."
params.genome = ""
params.pbat = '9' // default is 9, but this can be changed
params.verbose = false
params.single_end = false // default mode is auto-detect. NOTE: params are handed over automatically
params.fastqc_args = ''
params.fastq_screen_args = ''
params.trim_galore_args = ''
params.bismark_args = ''
params.deduplicate_bismark_args = ''
params.bismark_methylation_extractor_args = ''
params.bismark2bedGraph_args = ''
params.bismark2summary_args = ''
params.bismark2report_args= ''
params.multiqc_args = ''
params.help = false
// Show help message and exit
if (params.help){
helpMessage()
exit 0
}
params.list_genomes = false;
if (params.list_genomes){
println ("[WORKLFOW] List genomes selected")
}
if (params.verbose){
println ("[WORKFLOW] FASTQC ARGS: " + params.fastqc_args)
println ("[WORKFLOW] FASTQ SCREEN ARGS ARE: " + params.fastq_screen_args)
println ("[WORKFLOW] TRIM GALORE ARGS: " + params.trim_galore_args)
println ("[WORKFLOW] BISMARK ARGS ARE: " + params.bismark_args)
println ("[WORKFLOW] BISMARK DEDUPLICATION ARGS ARE: " + params.deduplicate_bismark_args)
println ("[WORKFLOW] BISMARK METHYLATION EXTRACTOR ARGS ARE: " + params.bismark_methylation_extractor_args)
}
include { makeFilesChannel; getFileBaseNames } from './nf_modules/files.mod.nf'
include { getGenome } from './nf_modules/genomes.mod.nf'
include { listGenomes } from './nf_modules/genomes.mod.nf'
if (params.list_genomes){
listGenomes() // this lists all available genomes, and exits
}
genome = getGenome(params.genome)
include { FASTQC } from './nf_modules/fastqc.mod.nf'
include { FASTQC as FASTQC2 } from './nf_modules/fastqc.mod.nf'
include { FASTQ_SCREEN } from './nf_modules/fastq_screen.mod.nf' params(bisulfite: '--bisulfite')
// Adding pbat flag for PBAT processing. For Trim Galore we pass on the number of bp that need to be clipped at the 5' end
include { TRIM_GALORE } from './nf_modules/trim_galore.mod.nf' params(pbat: params.pbat)
// Adding pbat flag for PBAT processing
include { BISMARK } from './nf_modules/bismark.mod.nf' params(genome: genome, pbat: true, unmapped: true)
include { BISMARK as BISMARK_UN1 } from './nf_modules/bismark.mod.nf' params(genome: genome, pbat: true, read_identity: "1")
include { BISMARK as BISMARK_UN2 } from './nf_modules/bismark.mod.nf' params(genome: genome, read_identity: "2")
include { BISMARK_DEDUPLICATION } from './nf_modules/bismark_deduplication.mod.nf'
include { BISMARK_DEDUPLICATION as BISMARK_DEDUP1} from './nf_modules/bismark_deduplication.mod.nf'
include { BISMARK_DEDUPLICATION as BISMARK_DEDUP2} from './nf_modules/bismark_deduplication.mod.nf'
// Adding pbat flag for PBAT processing
include { BISMARK_METHYLATION_EXTRACTOR } from './nf_modules/bismark_methylation_extractor.mod.nf' params(pbat: true)
include { BISMARK_METHYLATION_EXTRACTOR as METHXTRACT_UN1 } from './nf_modules/bismark_methylation_extractor.mod.nf'
include { BISMARK_METHYLATION_EXTRACTOR as METHXTRACT_UN2 } from './nf_modules/bismark_methylation_extractor.mod.nf'
include { BISMARK2BEDGRAPH } from './nf_modules/bismark2bedGraph.mod.nf' params(dirty_harry: true)
include { BISMARK2REPORT } from './nf_modules/bismark2report.mod.nf'
include { BISMARK2SUMMARY } from './nf_modules/bismark2summary.mod.nf'
include { MULTIQC } from './nf_modules/multiqc.mod.nf'
file_ch = makeFilesChannel(args)
workflow {
main:
FASTQC (file_ch, params.outdir, params.fastqc_args, params.verbose)
TRIM_GALORE (file_ch, params.outdir, params.trim_galore_args, params.verbose)
FASTQ_SCREEN (TRIM_GALORE.out.reads, params.outdir, params.fastq_screen_args, params.verbose)
FASTQC2 (TRIM_GALORE.out.reads, params.outdir, params.fastqc_args, params.verbose)
BISMARK (TRIM_GALORE.out.reads, params.outdir, params.bismark_args, params.verbose)
BISMARK_DEDUPLICATION (BISMARK.out.bam, params.outdir, params.deduplicate_bismark_args, params.verbose)
BISMARK_UN1 (BISMARK.out.unmapped1, params.outdir, params.bismark_args, params.verbose)
BISMARK_DEDUP1 (BISMARK_UN1.out.bam, params.outdir, params.deduplicate_bismark_args, params.verbose)
BISMARK_UN2 (BISMARK.out.unmapped2, params.outdir, params.bismark_args, params.verbose)
BISMARK_DEDUP2 (BISMARK_UN2.out.bam, params.outdir, params.deduplicate_bismark_args, params.verbose)
BISMARK_METHYLATION_EXTRACTOR (BISMARK_DEDUPLICATION.out.bam, params.outdir, params.bismark_methylation_extractor_args, params.verbose)
METHXTRACT_UN1 (BISMARK_DEDUP1.out.bam, params.outdir, params.bismark_methylation_extractor_args, params.verbose)
METHXTRACT_UN2 (BISMARK_DEDUP2.out.bam, params.outdir, params.bismark_methylation_extractor_args, params.verbose)
// BISMARK_METHYLATION_EXTRACTOR.out.context_files_CG.view()
// METHXTRACT_UN1.out.context_files_CG.view()
// METHXTRACT_UN2.out.context_files_CG.view()
// To merge all relevant CpG calls, we will need a join statement instead of mix
// left.join(right).view(); Joining occurs with the first element as the key
ch_bedGraph_unmapped = METHXTRACT_UN1.out.context_files_CG
.join(METHXTRACT_UN2.out.context_files_CG, by: [0])
// .view()
ch_bismark2bedGraph = BISMARK_METHYLATION_EXTRACTOR.out.context_files_CG
.join(ch_bedGraph_unmapped, by: [0])
.map { it -> [ it[0], [it[1], it[2], it[3]].flatten() ] }
//.view()
BISMARK2BEDGRAPH (ch_bismark2bedGraph, params.outdir, params.bismark2bedGraph_args, params.verbose)
// merging channels for Bismark reports
bismark_report_ch = BISMARK.out.bam.mix(
BISMARK.out.report,
BISMARK_UN1.out.bam.ifEmpty([]),
BISMARK_UN1.out.report.ifEmpty([]),
BISMARK_UN2.out.bam.ifEmpty([]),
BISMARK_UN2.out.report.ifEmpty([]),
BISMARK_DEDUPLICATION.out.report.ifEmpty([]),
BISMARK_DEDUP1.out.report.ifEmpty([]),
BISMARK_DEDUP2.out.report.ifEmpty([]),
BISMARK_METHYLATION_EXTRACTOR.out.report.ifEmpty([]),
BISMARK_METHYLATION_EXTRACTOR.out.mbias.ifEmpty([]),
METHXTRACT_UN1.out.report.ifEmpty([]),
METHXTRACT_UN1.out.mbias.ifEmpty([]),
METHXTRACT_UN2.out.report.ifEmpty([]),
METHXTRACT_UN2.out.mbias.ifEmpty([]),
).collect()
BISMARK2REPORT (bismark_report_ch, params.outdir, params.bismark2report_args, params.verbose)
// merging channels for Bismark summary report
// PAIRED-END only as it fails with a mixture of single-end and paired-end files
bismark_summary_report_ch = BISMARK.out.bam.mix(
BISMARK.out.report,
BISMARK_DEDUPLICATION.out.report.ifEmpty([]),
BISMARK_METHYLATION_EXTRACTOR.out.report.ifEmpty([]),
BISMARK_METHYLATION_EXTRACTOR.out.mbias.ifEmpty([]),
).collect()
// bismark_report_ch.subscribe { println "Got: $it" }
BISMARK2SUMMARY (bismark_summary_report_ch, params.outdir, params.bismark2summary_args, params.verbose)
// merging channels for MultiQC
multiqc_ch = FASTQC.out.report.mix(
TRIM_GALORE.out.report,
FASTQ_SCREEN.out.report.ifEmpty([]),
FASTQC2.out.report.ifEmpty([]),
BISMARK.out.report.ifEmpty([]),
BISMARK_UN1.out.report.ifEmpty([]),
BISMARK_UN2.out.report.ifEmpty([]),
BISMARK_DEDUPLICATION.out.report.ifEmpty([]),
BISMARK_DEDUP1.out.report.ifEmpty([]),
BISMARK_DEDUP2.out.report.ifEmpty([]),
BISMARK_METHYLATION_EXTRACTOR.out.report.ifEmpty([]),
BISMARK_METHYLATION_EXTRACTOR.out.mbias.ifEmpty([]),
METHXTRACT_UN1.out.report.ifEmpty([]),
METHXTRACT_UN1.out.mbias.ifEmpty([]),
METHXTRACT_UN2.out.report.ifEmpty([]),
METHXTRACT_UN2.out.mbias.ifEmpty([]),
).collect()
MULTIQC (multiqc_ch, params.outdir, params.multiqc_args, params.verbose)
}
// Since workflows with very long command lines tend to fail to get rendered at all, I was experimenting with a
// minimal execution summary report so we at least know what the working directory was...
workflow.onComplete {
def msg = """\
Pipeline execution summary
---------------------------
Jobname : ${workflow.runName}
Completed at: ${workflow.complete}
Duration : ${workflow.duration}
Success : ${workflow.success}
workDir : ${workflow.workDir}
exit status : ${workflow.exitStatus}
"""
.stripIndent()
sendMail(to: "${workflow.userName}@babraham.ac.uk", subject: 'Minimal pipeline execution report', body: msg)
}
def helpMessage() {
log.info"""
>>
SYNOPSIS:
PBAT - Dirty Harry: an abomination for paired-end data...
This workflow runs an entire Bisulfite-seq processing pipeline for PBAT (post-bisulfite adaptor tagging) samples
on FastQ files, first for paired-end end file, and then in single-end mode for unmapped files. This includes QC,
quality-/adapter trimming, contamination QC (post-trimming), alignments to a genome using Bismark (see below
deduplication, methylation extraction, coverage file generation, and finally generates aa aggregate MultiQC report.
The workflow is suitable for paired-end PBAT-type bisulfite sequencing experiments, by default assuming 9N oligos
for the pulldown reaction. The term ‘Dirty Harry’ was coined in ~2014 because it doesn’t seem to be the cleanest of methods…:
(1) PE alignment (--pbat) to start, while grabbing the unmapped R1 and R2 reads out the end using the option ‘--unmapped’.
Properly aligned PE reads should be methylation extracted using --no_overlap.
(2) unmapped R1 is then mapped SE (--pbat)
(3) unmapped R2 is then mapped SE (default mode, i.e. directional)
Single-end aligned R1 and R2 can then be methylation extracted normally as they should in theory map to different places
in the genome anyway so don’t require attention to overlapping reads. Note that paired-end and single-end processing is
PE/SE specific until the methylation extraction stage. [CpG]* methylation calls of PE as well as both SE aligments are then
fed to a single bismark2bedGraph process to generate a single, merged coverage file (ends in "_DH.cov.gz").
Here is a graphical representation of the workflow:
--- FastQC
--- Trim Galore
|
--- FastQ Screen
--- FastQC
--- Bismark (PE)
|
--- Bismark deduplication (PE)
|
--- Bismark methylation extraction (PE)
---unmapped R1
|
---Bismark deduplication (SE R1)
|
--- Bismark methylation extraction (SE R1)
---unmapped R2
|
--- Bismark deduplication (SE R2)
|
--- Bismark methylation extraction (SE R2)
---bismark2bedGraph (PE, SE R1, SE R2)
--- bismark2report*
--- bismark2summary*
--- MultiQC*
* These steps run only once ALL other jobs have completed.
By default all these steps are submitted as jobs to the Babraham stone compute cluster.
By default, the involved tools are run in the following way:
------------------------------------------------------------
FastQC: defaults ('-q')
FastQ Screen: '--bisulfite'
Trim Galore: adapter auto-detection; '--paired --clip_r1 9 --clip_r2 9'
Bismark: PE: '--unmapped --pbat', SE R1: '--pbat', SE R2: defaults
Bismark deduplication: defaults
Bismark methylation extraction: '--bedGraph --buffer 10G --parallel 4'
bismark2bedGraph 'dirty_harry'
To add additional parameters to any of the programs, consider supplying tool-specific arguments (see --toolname_args="..." below).
==============================================================================================================
USAGE:
nf_bisulfite_PBAT_DirtyHarry [options] --genome <genomeID> <input files>
Mandatory arguments:
====================
<input files> List of input files, e.g. '*fastq.gz' or '*fq.gz'. Files are automatically processed as
single-end (SE) or paired end (PE) files (if file pairs share the same base-name, and differ only
by a read number, e.g. 'base_name_R1.fastq.gz' and 'base_name_R2.fastq.gz' (or R3, R4). For
PE files, only Read 1 is run through FastQ Screen (as typically R1 and R2 produce nearly identical
contamination profiles). To run PE files in single-end mode, please see '--single_end' below.
--genome [str] Genome build ID to be used for the alignment, e.g. GRCh38 (latest human genome) or GRCm38
(latest mouse genome build). To list all available genomes, see '--list_genomes' below.
Tool-specific options:
======================
For all following options, please note that the format: ="your options" needs to be strictly adhered to in order to work correctly.
--fastqc_args="[str]" This option can take any number of options that are compatible with FastQC to modify its default
behaviour. For more detailed information on available options please refer to the FastQC documentation,
or run 'fastqc --help' on the command line. As an example, to run FastQC without grouping of bases when
reads are >50bp and use a specific file with non-default adapter sequences, use:
' --fastqc_args="--nogroup --adapters ./non_default_adapter_file.txt" '. [Default: None]
--fastq_screen_args="[str]" This option can take any number of options that are compatible with FastQ Screen to modify its
default behaviour. For more detailed information on available options please refer to the FastQ Screen
documentation, or run 'fastq_screen --help' on the command line. For instance, to process a bisulfite
converted library with fairly relaxed parameters, you could use:
' --fastq_screen_args="--bisulfite --score_min L,0,-0.6" '. [Default: None]
--trim_galore_args="[str]" This option can take any number of options that are compatible with Trim Galore to modify its
default trimming behaviour. For more detailed information on available options please refer
to the Trim Galore User Guide, or run 'trim_galore --help' on the command line. As an example, to trim
off the first 10bp from the 5' of R1 and 5bp of R2, use:
' --trim_galore_args="--clip_r1 10 --clip_r2 5" '. [Default: None]
--bismark_args="[str]" This option can take any number of options that are compatible with Bismark to modify its
default mapping behaviour. For more detailed information on available options please refer
to the Bismark User Guide, or run 'bismark --help' on the comamnd line. As an example, to run somewhat
relaxed alignments, use: ' --bismark_args="--score_min L,0,-0.4" '. [Default: None]
--bismark_deduplication_args="[str]" This option can take any number of options that are compatible with deduplicate_bismark
to modify its default deduplication behaviour. For more detailed information on available options please
refer to the Bismark User Guide, or run 'deduplicate_bismark --help'. As an example, to take UMI sequences
into account, use: ' --bismark_deduplication_args="--barcode" '. [Default: None]
--bismark_methylation_extractor_args="[str]" This option can take any number of options that are compatible with bismark_methylation_extractor
to modify its default exctraction behaviour. For more detailed information on available options please refer
to the Bismark User Guide, or run 'bismark_methylation_extractor --help' on the command line. As an example,
to run the extraction and coverage file generation for a genome containing thousands of scaffolds
or contigs, use: ' --bismark_methylation_extractor_args="--gazillion" '. [Default: None]
--bismark2bedGraph_args="[str]" This option can take any number of options that are compatible with bismark2bedGraph to modify its
default behaviour. For more detailed information on available options please refer to the Bismark
User Guide, or run 'bismark2bedGraph --help' on the command line. As an example, to run the coverage file
generation for a genome containing thousands of scaffolds or contigs, use:
' --bismark2bedGraph_args="--gazillion" '. [Default: None]
Other options:
==============
--outdir [str] Path to the output directory. [Default: current working directory]
--pbat [int] Number of nucleotides used for the PBAT pulldown reaction. This number of bases will be removed
from the 5'-end of both R1 and R2 prior to mapping. Default: 9.
--list_genomes List all genome builds that are currently available to choose from. To see this list
of available genomes with more detailed information about paths and indexes, run
the command as '--list_genomes --verbose'
--single_end Force files of a read pair to be treated as single-end files. [Default: auto-detect]
--verbose More verbose status messages. [Default: OFF]
--help Displays this help message and exits.
Workflow options:
=================
Please note the single '-' hyphen for the following options!
-resume If a pipeline workflow has been interrupted or stopped (e.g. by accidentally closing a laptop),
this option will attempt to resume the workflow at the point it got interrupted by using
Nextflow's caching mechanism. This may save a lot of time.
-bg Sends the entire workflow into the background, thus disconnecting it from the terminal session.
This option launches a daemon process (which will keep running on the headnode) that watches over
your workflow, and submits new jobs to the SLURM queue as required. Use this option for big pipeline
jobs, or whenever you do not want to watch the status progress yourself. Upon completion, the
pipeline will send you an email with the job details. This option is HIGHLY RECOMMENDED!
-process.executor=local Temporarily changes where the workflow is executed to the 'local' machine. See also Nextflow config
file for more details. [Default: slurm]
<<
""".stripIndent()
}