Newer
Older
#!/usr/bin/env nextflow
/*
========================================================================================
nf-core/hic
========================================================================================
nf-core/hic Analysis Pipeline.
#### Homepage / Documentation
https://github.com/nf-core/hic
----------------------------------------------------------------------------------------
*/
/*TOOO
- outputs
- multiqc
- update version tools
- install + compile
*/
def helpMessage() {
log.info"""
=======================================================
,--./,-.
___ __ __ __ ___ /,-._.--~\'
|\\ | |__ __ / ` / \\ |__) |__ } {
| \\| | \\__, \\__/ | \\ |___ \\`-._,-`-,
`._,._,\'
nf-core/hic v${workflow.manifest.version}
=======================================================
This pipeline is a Nextflow version of the HiC-Pro pipeline for Hi-C data processing.
See https://github.com/nservant/HiC-Pro for details.
Usage:
The typical command for running the pipeline is as follows:
nextflow run nf-core/hic --reads '*_R{1,2}.fastq.gz' -profile conda
Mandatory arguments:
--reads Path to input data (must be surrounded with quotes)
--genome Name of iGenomes reference
-profile Configuration profile to use. Can use multiple (comma separated)
Available: conda, docker, singularity, awsbatch, test and more.
References If not specified in the configuration file or you wish to overwrite any of the references.
--bwt2_index Path to Bowtie2 index
--fasta Path to Fasta reference
--chromosome_size Path to chromosome size file
--restriction_fragments Path to restriction fragment file (bed)
--bwt2_opts_end2end Options for bowtie2 end-to-end mappinf (first mapping step)
--bwt2_opts_trimmed Options for bowtie2 mapping after ligation site trimming
--min_mapq Minimum mapping quality values to consider
--restriction_site Cutting motif(s) of restriction enzyme(s) (comma separated)
--ligation_site Ligation motifs to trim (comma separated)
--min_restriction_fragment_size Minimum size of restriction fragments to consider
--max_restriction_framgnet_size Maximum size of restriction fragmants to consider
--min_insert_size Minimum insert size of mapped reads to consider
--max_insert_size Maximum insert size of mapped reads to consider
--min_cis_dist Minimum intra-chromosomal distance to consider
--rm_singleton Remove singleton reads
--rm_multi Remove multi-mapped reads
--rm_dup Remove duplicates
--bin_size Bin size for contact maps (comma separated)
--ice_max_iter Maximum number of iteration for ICE normalization
--ice_filter_low_count_perc Percentage of low counts columns/rows to filter before ICE normalization
--ice_filter_high_count_perc Percentage of high counts columns/rows to filter before ICE normalization
--ice_eps Convergence criteria for ICE normalization
--outdir The output directory where the results will be saved
--email Set this parameter to your e-mail address to get a summary e-mail with details of the run sent to you when the workflow exits
-name Name for the pipeline run. If not specified, Nextflow will automatically generate a random mnemonic.
--awsqueue The AWSBatch JobQueue that needs to be set when running on AWSBatch
--awsregion The AWS Region for your AWS Batch job to run on
""".stripIndent()
}
/**********************************************************
* SET UP CONFIGURATION VARIABLES
*/
// Show help emssage
if (params.help){
helpMessage()
exit 0
}
// Check if genome exists in the config file
if (params.genomes && params.genome && !params.genomes.containsKey(params.genome)) {
exit 1, "The provided genome '${params.genome}' is not available in the iGenomes file. Currently the available genomes are ${params.genomes.keySet().join(", ")}"
}
// Reference index path configuration
// Define these here - after the profiles are loaded with the iGenomes paths
params.bwt2_index = params.genome ? params.genomes[ params.genome ].bowtie2 ?: false : false
params.fasta = params.genome ? params.genomes[ params.genome ].fasta ?: false : false
//params.chromosome_size = false
//params.restriction_fragments = false
// Has the run name been specified by the user?
// this has the bonus effect of catching both -name and --name
custom_runName = params.name
if( !(workflow.runName ==~ /[a-z]+_[a-z]+/) ){
custom_runName = workflow.runName
}
if( workflow.profile == 'awsbatch') {
// AWSBatch sanity checking
if (!params.awsqueue || !params.awsregion) exit 1, "Specify correct --awsqueue and --awsregion parameters on AWSBatch!"
if (!workflow.workDir.startsWith('s3') || !params.outdir.startsWith('s3')) exit 1, "Specify S3 URLs for workDir and outdir parameters on AWSBatch!"
// Check workDir/outdir paths to be S3 buckets if running on AWSBatch
// related: https://github.com/nextflow-io/nextflow/issues/813
if (!workflow.workDir.startsWith('s3:') || !params.outdir.startsWith('s3:')) exit 1, "Workdir or Outdir not on S3 - specify S3 Buckets for each to run on AWSBatch!"
}
// Stage config files
ch_multiqc_config = Channel.fromPath(params.multiqc_config)
ch_output_docs = Channel.fromPath("$baseDir/docs/output.md")
/**********************************************************
* SET UP CHANNELS
*/
/*
* input read files
*/
if (params.readPaths){
raw_reads = Channel.create()
raw_reads_2 = Channel.create()
Channel
.from( params.readPaths )
.map { row -> [ row[0], [file(row[1][0]), file(row[1][1])]] }
.separate( raw_reads, raw_reads_2 ) { a -> [tuple(a[0], a[1][0]), tuple(a[0], a[1][1])] }
}else{
raw_reads = Channel.create()
raw_reads_2 = Channel.create()
Channel
.fromFilePairs( params.reads )
.separate( raw_reads, raw_reads_2 ) { a -> [tuple(a[0], a[1][0]), tuple(a[0], a[1][1])] }
}
raw_reads = raw_reads.concat( raw_reads_2 )
// SPlit fastq files
// https://www.nextflow.io/docs/latest/operator.html#splitfastq
lastPath = params.bwt2_index.lastIndexOf(File.separator)
bwt2_dir = params.bwt2_index.substring(0,lastPath+1)
bwt2_base = params.bwt2_index.substring(lastPath+1)
Channel.fromPath( bwt2_dir , checkIfExists: true)
.ifEmpty { exit 1, "Genome index: Provided index not found: ${params.bwt2_index}" }
.into { bwt2_index_end2end; bwt2_index_trim }
lastPath = params.fasta.lastIndexOf(File.separator)
bwt2_base = params.fasta.substring(lastPath+1)
.ifEmpty { exit 1, "Genome index: Fasta file not found: ${params.fasta}" }
.set { fasta_for_index }
}
else {
exit 1, "No reference genome specified!"
}
//println (bwt2_dir)
//println (bwt2_base)
Channel.fromPath( params.chromosome_size , checkIfExists: true)
.set {chromosome_size}
.ifEmpty { exit 1, "Chromosome sizes: Fasta file not found: ${params.fasta}" }
.set { fasta_for_chromsize }
}
else {
exit 1, "No chromosome size specified!"
}
// Restriction fragments
if ( params.restriction_fragments ){
Channel.fromPath( params.restriction_fragments, checkIfExists: true )
.ifEmpty { exit 1, "Restriction fragments: Fasta file not found: ${params.fasta}" }
else {
exit 1, "No restriction fragments file specified!"
}
// Resolutions for contact maps
map_res = Channel.from( params.bins_size.tokenize(',') )
/**********************************************************
* SET UP LOGS
*/
// Header log info
log.info """=======================================================
,--./,-.
___ __ __ __ ___ /,-._.--~\'
|\\ | |__ __ / ` / \\ |__) |__ } {
| \\| | \\__, \\__/ | \\ |___ \\`-._,-`-,
`._,._,\'
nf-core/hic v${workflow.manifest.version}"
======================================================="""
def summary = [:]
summary['Pipeline Name'] = 'nf-core/hic'
summary['Pipeline Version'] = workflow.manifest.version
summary['Run Name'] = custom_runName ?: workflow.runName
summary['Reads'] = params.reads
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
summary['Max Memory'] = params.max_memory
summary['Max CPUs'] = params.max_cpus
summary['Max Time'] = params.max_time
summary['Output dir'] = params.outdir
summary['Working dir'] = workflow.workDir
summary['Container Engine'] = workflow.containerEngine
if(workflow.containerEngine) summary['Container'] = workflow.container
summary['Current home'] = "$HOME"
summary['Current user'] = "$USER"
summary['Current path'] = "$PWD"
summary['Working dir'] = workflow.workDir
summary['Output dir'] = params.outdir
summary['Script dir'] = workflow.projectDir
summary['Config Profile'] = workflow.profile
if(workflow.profile == 'awsbatch'){
summary['AWS Region'] = params.awsregion
summary['AWS Queue'] = params.awsqueue
}
if(params.email) summary['E-mail Address'] = params.email
log.info summary.collect { k,v -> "${k.padRight(15)}: $v" }.join("\n")
log.info "========================================="
def create_workflow_summary(summary) {
def yaml_file = workDir.resolve('workflow_summary_mqc.yaml')
yaml_file.text = """
id: 'nf-core-hic-summary'
description: " - this information is collected when the pipeline is started."
section_name: 'nf-core/hic Workflow Summary'
section_href: 'https://github.com/nf-core/hic'
plot_type: 'html'
data: |
<dl class=\"dl-horizontal\">
${summary.collect { k,v -> " <dt>$k</dt><dd><samp>${v ?: '<span style=\"color:#999999;\">N/A</a>'}</samp></dd>" }.join("\n")}
</dl>
""".stripIndent()
return yaml_file
}
/*
* Parse software version numbers
*/
process get_software_versions {
output:
file 'software_versions_mqc.yaml' into software_versions_yaml
script:
echo $workflow.manifest.version > v_pipeline.txt
echo $workflow.nextflow.version > v_nextflow.txt
bowtie2 --version > v_bowtie2.txt
python --version > v_python.txt
samtools --version > v_samtools.txt
scrape_software_versions.py > software_versions_mqc.yaml
"""
/****************************************************
* PRE-PROCESSING
*/
if(!params.bwt2_index && params.fasta){
process makeBowtie2Index {
tag "$bwt2_base"
//publishDir path: { params.saveReference ? "${params.outdir}/reference_genome" : params.outdir },
// saveAs: { params.saveReference ? it : null }, mode: 'copy'
file "bowtie2_index" into bwt2_index_end2end
file "bowtie2_index" into bwt2_index_trim
bwt2_base = fasta.toString() - ~/(\.fa)?(\.fasta)?(\.fas)?$/
mkdir bowtie2_index
bowtie2-build ${fasta} bowtie2_index/${bwt2_base}
"""
}
}
if(!params.chromosome_size && params.fasta){
process makeChromSize {
tag "$fasta"
//publishDir path: { params.saveReference ? "${params.outdir}/reference_genome" : params.outdir },
// saveAs: { params.saveReference ? it : null }, mode: 'copy'
input:
file fasta from fasta_for_chromsize
output:
file "*.size" into chromosome_size
script:
"""
samtools faidx ${fasta}
cut -f1,2 ${fasta}.fai > chrom.size
//publishDir path: { params.saveReference ? "${params.outdir}/reference_genome" : params.outdir },
// saveAs: { params.saveReference ? it : null }, mode: 'copy'
input:
file fasta from fasta_for_resfrag
output:
digest_genome.py -r ${params.restriction_site} -o restriction_fragments.bed ${fasta}
/****************************************************
* MAIN WORKFLOW
*/
* STEP 1 - Two-steps Reads Mapping
*/
process bowtie2_end_to_end {
tag "$prefix"
input:
set val(sample), file(reads) from raw_reads
file index from bwt2_index_end2end.collect()
output:
set val(prefix), file("${prefix}_unmap.fastq") into unmapped_end_to_end
set val(prefix), file("${prefix}.bam") into end_to_end_bam
script:
prefix = reads.toString() - ~/(\.fq)?(\.fastq)?(\.gz)?$/
def bwt2_opts = params.bwt2_opts_end2end
"""
bowtie2 --rg-id BMG --rg SM:${prefix} \\
${bwt2_opts} \\
-p ${task.cpus} \\
--un ${prefix}_unmap.fastq \\
-U ${reads} | samtools view -F 4 -bS - > ${prefix}.bam
"""
}
process trim_reads {
tag "$prefix"
input:
set val(prefix), file(reads) from unmapped_end_to_end
output:
set val(prefix), file("${prefix}_trimmed.fastq") into trimmed_reads
script:
"""
cutsite_trimming --fastq $reads \\
--cutsite ${params.ligation_site} \\
--out ${prefix}_trimmed.fastq
"""
}
process bowtie2_on_trimmed_reads {
tag "$prefix"
input:
set val(prefix), file(reads) from trimmed_reads
file index from bwt2_index_trim.collect()
output:
set val(prefix), file("${prefix}_trimmed.bam") into trimmed_bam
script:
prefix = reads.toString() - ~/(_trimmed)?(\.fq)?(\.fastq)?(\.gz)?$/
"""
bowtie2 --rg-id BMG --rg SM:${prefix} \\
${params.bwt2_opts_trimmed} \\
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
-U ${reads} | samtools view -bS - > ${prefix}_trimmed.bam
"""
}
process merge_mapping_steps{
tag "$bam1 + $bam2"
input:
set val(prefix), file(bam1), file(bam2) from end_to_end_bam.join( trimmed_bam )
output:
set val(sample), file("${prefix}_bwt2merged.bam") into bwt2_merged_bam
script:
sample = prefix.toString() - ~/(_R1)?(_R2)?(_val_1)?(_val_2)?$/
"""
samtools merge -@ ${task.cpus} \\
-f ${prefix}_bwt2merged.bam \\
${bam1} ${bam2}
samtools sort -@ ${task.cpus} -m 800M \\
-n -T /tmp/ \\
-o ${prefix}_bwt2merged.sorted.bam \\
${prefix}_bwt2merged.bam
mv ${prefix}_bwt2merged.sorted.bam ${prefix}_bwt2merged.bam
"""
process combine_mapped_files{
tag "$sample = $r1_prefix + $r2_prefix"
input:
set val(sample), file(aligned_bam) from bwt2_merged_bam.groupTuple()
output:
set val(sample), file("${sample}_bwt2pairs.bam") into paired_bam
script:
r1_bam = aligned_bam[0]
r1_prefix = r1_bam.toString() - ~/_bwt2merged.bam$/
r2_bam = aligned_bam[1]
r2_prefix = r2_bam.toString() - ~/_bwt2merged.bam$/
def opts = ""
opts = params.rm_singleton ? "${opts}" : "--single ${opts}"
opts = params.rm_multi ? "${opts}" : "--multi ${opts}"
mergeSAM.py -f ${r1_bam} -r ${r2_bam} -o ${sample}_bwt2pairs.bam -q ${params.min_mapq} ${opts}
* STEP2 - DETECT VALID PAIRS
*/
process get_valid_interaction{
tag "$sample"
input:
set val(sample), file(pe_bam) from paired_bam
val frag_file from res_frag_file
output:
set val(sample), file("*.validPairs") into valid_pairs
set val(sample), file("*.validPairs") into valid_pairs_4cool
def opts = ""
if (params.min_cis_dist != "") opts="${opts} -d ${params.min_cis_dist}"
if (params.min_insert_size != "") opts="${opts} -s ${params.min_insert_size}"
if (params.max_insert_size != "") opts="${opts} -l ${params.max_insert_size}"
if (params.min_restriction_fragment_size != "") opts="${opts} -t ${params.min_restriction_fragment_size}"
if (params.max_restriction_fragment_size != "") opts="${opts} -m ${params.max_restriction_fragment_size}"
mapped_2hic_fragments.py -f ${frag_file} -r ${pe_bam} ${opts}
"""
}
/*
* STEP3 - BUILD MATRIX
*/
set val(sample), file(vpairs), val(mres) from valid_pairs.combine(map_res)
file("*.matrix") into raw_maps
build_matrix --matrix-format upper --binsize ${mres} --chrsizes ${chrsize} --ifile ${vpairs} --oprefix ${sample}_${mres}
/*
* STEP 4 - NORMALIZE MATRIX
*/
process run_iced{
tag "$rmaps"
input:
file(rmaps) from raw_maps
output:
file("*iced.matrix") into iced_maps
script:
prefix = rmaps.toString() - ~/(\.matrix)?$/
"""
ice --filter_low_counts_perc ${params.ice_filer_low_count_perc} \
--results_filename ${prefix}_iced.matrix \
--filter_high_counts_perc ${params.ice_filer_high_count_perc} \
--max_iter ${params.ice_max_iter} --eps ${params.ice_eps} --remove-all-zeros-loci --output-bias 1 --verbose 1 ${rmaps}
"""
/*
* STEP 5 - COOLER FILE
process generate_cool{
tag "$sample"
input:
set val(sample), file(vpairs) from valid_pairs_4cool
output:
file("*mcool") into cool_maps
script:
"""
hicpro2higlass.sh -i $vpairs -r 5000 -c ${chrsize} -n
"""
}
*/
process multiqc {
publishDir "${params.outdir}/MultiQC", mode: 'copy'
input:
file multiqc_config from ch_multiqc_config
// TODO nf-core: Add in log files from your new processes for MultiQC to find!
file ('fastqc/*') from fastqc_results.collect().ifEmpty([])
file ('software_versions/*') from software_versions_yaml
file workflow_summary from create_workflow_summary(summary)
output:
file "*multiqc_report.html" into multiqc_report
file "*_data"
script:
rtitle = custom_runName ? "--title \"$custom_runName\"" : ''
rfilename = custom_runName ? "--filename " + custom_runName.replaceAll('\\W','_').replaceAll('_+','_') + "_multiqc_report" : ''
// TODO nf-core: Specify which MultiQC modules to use with -m for a faster run time
"""
multiqc -f $rtitle $rfilename --config $multiqc_config .
"""
}
process output_documentation {
publishDir "${params.outdir}/Documentation", mode: 'copy'
input:
file output_docs from ch_output_docs
output:
file "results_description.html"
script:
"""
markdown_to_html.r $output_docs results_description.html
"""
}
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
/*
* Completion e-mail notification
*/
workflow.onComplete {
// Set up the e-mail variables
def subject = "[nf-core/hic] Successful: $workflow.runName"
if(!workflow.success){
subject = "[nf-core/hic] FAILED: $workflow.runName"
}
def email_fields = [:]
email_fields['version'] = workflow.manifest.version
email_fields['runName'] = custom_runName ?: workflow.runName
email_fields['success'] = workflow.success
email_fields['dateComplete'] = workflow.complete
email_fields['duration'] = workflow.duration
email_fields['exitStatus'] = workflow.exitStatus
email_fields['errorMessage'] = (workflow.errorMessage ?: 'None')
email_fields['errorReport'] = (workflow.errorReport ?: 'None')
email_fields['commandLine'] = workflow.commandLine
email_fields['projectDir'] = workflow.projectDir
email_fields['summary'] = summary
email_fields['summary']['Date Started'] = workflow.start
email_fields['summary']['Date Completed'] = workflow.complete
email_fields['summary']['Pipeline script file path'] = workflow.scriptFile
email_fields['summary']['Pipeline script hash ID'] = workflow.scriptId
if(workflow.repository) email_fields['summary']['Pipeline repository Git URL'] = workflow.repository
if(workflow.commitId) email_fields['summary']['Pipeline repository Git Commit'] = workflow.commitId
if(workflow.revision) email_fields['summary']['Pipeline Git branch/tag'] = workflow.revision
email_fields['summary']['Nextflow Version'] = workflow.nextflow.version
email_fields['summary']['Nextflow Build'] = workflow.nextflow.build
email_fields['summary']['Nextflow Compile Timestamp'] = workflow.nextflow.timestamp
// Render the TXT template
def engine = new groovy.text.GStringTemplateEngine()
def tf = new File("$baseDir/assets/email_template.txt")
def txt_template = engine.createTemplate(tf).make(email_fields)
def email_txt = txt_template.toString()
// Render the HTML template
def hf = new File("$baseDir/assets/email_template.html")
def html_template = engine.createTemplate(hf).make(email_fields)
def email_html = html_template.toString()
// Render the sendmail template
def smail_fields = [ email: params.email, subject: subject, email_txt: email_txt, email_html: email_html, baseDir: "$baseDir" ]
def sf = new File("$baseDir/assets/sendmail_template.txt")
def sendmail_template = engine.createTemplate(sf).make(smail_fields)
def sendmail_html = sendmail_template.toString()
// Send the HTML e-mail
if (params.email) {
try {
if( params.plaintext_email ){ throw GroovyException('Send plaintext e-mail, not HTML') }
// Try to send HTML e-mail using sendmail
[ 'sendmail', '-t' ].execute() << sendmail_html
log.info "[nf-core/hic] Sent summary e-mail to $params.email (sendmail)"
} catch (all) {
// Catch failures and try with plaintext
[ 'mail', '-s', subject, params.email ].execute() << email_txt
log.info "[nf-core/hic] Sent summary e-mail to $params.email (mail)"
}
}
// Write summary e-mail HTML to a file
def output_d = new File( "${params.outdir}/Documentation/" )
if( !output_d.exists() ) {
output_d.mkdirs()
}
def output_hf = new File( output_d, "pipeline_report.html" )
output_hf.withWriter { w -> w << email_html }
def output_tf = new File( output_d, "pipeline_report.txt" )
output_tf.withWriter { w -> w << email_txt }
log.info "[nf-core/hic] Pipeline Complete"
}