main.nf

#!/usr/bin/env nextflow
/*
========================================================================================
                         nf-core/atacseq
========================================================================================
 nf-core/atacseq Analysis Pipeline.
 #### Homepage / Documentation
 https://github.com/nf-core/atacseq
----------------------------------------------------------------------------------------
*/

def helpMessage() {
    log.info nfcoreHeader()
    log.info"""
    Usage:

    The typical command for running the pipeline is as follows:

      nextflow run nf-core/atacseq --input design.csv --genome GRCh37 -profile docker

    Mandatory arguments:
      --input [file]                  Comma-separated file containing information about the samples in the experiment (see docs/usage.md) (Default: './design.csv')
      --fasta [file]                  Path to Fasta reference. Not mandatory when using reference in iGenomes config via --genome (Default: false)
      --gtf [file]                    Path to GTF file. Not mandatory when using reference in iGenomes config via --genome (Default: false)
      -profile [str]                  Configuration profile to use. Can use multiple (comma separated)
                                      Available: conda, docker, singularity, awsbatch, test

    Generic
      --single_end [bool]             Specifies that the input is single-end reads (Default: false)
      --seq_center [str]              Sequencing center information to be added to read group of BAM files (Default: false)
      --fragment_size [int]           Estimated fragment size used to extend single-end reads (Default: 0)
      --fingerprint_bins [int]        Number of genomic bins to use when calculating fingerprint plot (Default: 500000)

    References                        If not specified in the configuration file or you wish to overwrite any of the references
      --genome [str]                  Name of iGenomes reference (Default: false)
      --bwa_index [file]              Full path to directory containing BWA index including base name i.e. /path/to/index/genome.fa (Default: false)
      --gene_bed [file]               Path to BED file containing gene intervals (Default: false)
      --tss_bed [file]                Path to BED file containing transcription start sites (Default: false)
      --macs_gsize [str]              Effective genome size parameter required by MACS2. If using iGenomes config, values have only been provided when --genome is set as GRCh37, GRCm38, hg19, mm10, BDGP6 and WBcel235 (Default: false)
      --blacklist [file]              Path to blacklist regions (.BED format), used for filtering alignments (Default: false)
      --mito_name [str]               Name of Mitochondrial chomosome in genome fasta (e.g. chrM). Reads aligning to this contig are filtered out (Default: false)
      --save_reference [bool]         If generated by the pipeline save the BWA index in the results directory (Default: false)

    Trimming
      --clip_r1 [int]                 Instructs Trim Galore to remove bp from the 5' end of read 1 (or single-end reads) (Default: 0)
      --clip_r2 [int]                 Instructs Trim Galore to remove bp from the 5' end of read 2 (paired-end reads only) (Default: 0)
      --three_prime_clip_r1 [int]     Instructs Trim Galore to remove bp from the 3' end of read 1 AFTER adapter/quality trimming has been performed (Default: 0)
      --three_prime_clip_r2 [int]     Instructs Trim Galore to remove bp from the 3' end of read 2 AFTER adapter/quality trimming has been performed (Default: 0)
      --trim_nextseq [int]            Instructs Trim Galore to apply the --nextseq=X option, to trim based on quality after removing poly-G tails (Default: 0)
      --skip_trimming [bool]          Skip the adapter trimming step (Default: false)
      --save_trimmed [bool]           Save the trimmed FastQ files in the results directory (Default: false)

    Alignments
      --keep_mito [bool]              Reads mapping to mitochondrial contig are not filtered from alignments (Default: false)
      --keep_dups [bool]              Duplicate reads are not filtered from alignments (Default: false)
      --keep_multi_map [bool]         Reads mapping to multiple locations are not filtered from alignments (Default: false)
      --skip_merge_replicates [bool]  Do not perform alignment merging and downstream analysis by merging replicates i.e. only do this by merging resequenced libraries (Default: false)
      --save_align_intermeds [bool]   Save the intermediate BAM files from the alignment step - not done by default (Default: false)

    Peaks
      --narrow_peak [bool]            Run MACS2 in narrowPeak mode (Default: false)
      --broad_cutoff [float]          Specifies broad cutoff value for MACS2. Only used when --narrow_peak isnt specified (Default: 0.1)
      --min_reps_consensus [int]      Number of biological replicates required from a given condition for a peak to contribute to a consensus peak (Default: 1)
      --save_macs_pileup [bool]       Instruct MACS2 to create bedGraph files normalised to signal per million reads (Default: false)
      --skip_peak_qc [bool]           Skip MACS2 peak QC plot generation (Default: false)
      --skip_peak_annotation [bool]   Skip annotation of MACS2 and consensus peaks with HOMER (Default: false)
      --skip_consensus_peaks [bool]   Skip consensus peak generation (Default: false)

    Differential analysis
      --deseq2_vst [bool]             Use vst transformation instead of rlog with DESeq2 (Default: false)
      --skip_diff_analysis [bool]     Skip differential accessibility analysis (Default: false)

    QC
      --skip_fastqc [bool]            Skip FastQC (Default: false)
      --skip_picard_metrics [bool]    Skip Picard CollectMultipleMetrics (Default: false)
      --skip_preseq [bool]            Skip Preseq (Default: false)
      --skip_plot_profile [bool]      Skip deepTools plotProfile (Default: false)
      --skip_plot_fingerprint [bool]  Skip deepTools plotFingerprint (Default: false)
      --skip_ataqv [bool]             Skip Ataqv (Default: false)
      --skip_igv [bool]               Skip IGV (Default: false)
      --skip_multiqc [bool]           Skip MultiQC (Default: false)

    Other
      --outdir [file]                 The output directory where the results will be saved (Default: './results')
      --publish_dir_mode [str]        Mode for publishing results in the output directory. Available: symlink, rellink, link, copy, copyNoFollow, move (Default: copy)
      --email [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 (Default: false)
      --email_on_fail [email]         Same as --email, except only send mail if the workflow is not successful (Default: false)
      --max_multiqc_email_size [str]  Threshold size for MultiQC report to be attached in notification email. If file generated by pipeline exceeds the threshold, it will not be attached (Default: 25MB)
      -name [str]                     Name for the pipeline run. If not specified, Nextflow will automatically generate a random mnemonic (Default: false)

    AWSBatch
      --awsqueue [str]                The AWSBatch JobQueue that needs to be set when running on AWSBatch
      --awsregion [str]               The AWS Region for your AWS Batch job to run on
      --awscli [str]                  Path to the AWS CLI tool
    """.stripIndent()
}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
/* --                                                                     -- */
/* --                SET UP CONFIGURATION VARIABLES                       -- */
/* --                                                                     -- */
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

// Show help message
if (params.help) {
    helpMessage()
    exit 0
}

// 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
}

////////////////////////////////////////////////////
/* --         DEFAULT PARAMETER VALUES         -- */
////////////////////////////////////////////////////

// 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(", ")}"
}

// Configurable variables
params.fasta = params.genome ? params.genomes[ params.genome ].fasta ?: false : false
params.bwa_index = params.genome ? params.genomes[ params.genome ].bwa ?: false : false
params.gtf = params.genome ? params.genomes[ params.genome ].gtf ?: false : false
params.gene_bed = params.genome ? params.genomes[ params.genome ].bed12 ?: false : false
params.mito_name = params.genome ? params.genomes[ params.genome ].mito_name ?: false : false
params.macs_gsize = params.genome ? params.genomes[ params.genome ].macs_gsize ?: false : false
params.blacklist = params.genome ? params.genomes[ params.genome ].blacklist ?: false : false
params.anno_readme = params.genome ? params.genomes[ params.genome ].readme ?: false : false

// Global variables
def PEAK_TYPE = params.narrow_peak ? 'narrowPeak' : 'broadPeak'

////////////////////////////////////////////////////
/* --          CONFIG FILES                    -- */
////////////////////////////////////////////////////

// Pipeline config
ch_multiqc_config = file("$baseDir/assets/multiqc_config.yaml", checkIfExists: true)
ch_multiqc_custom_config = params.multiqc_config ? Channel.fromPath(params.multiqc_config, checkIfExists: true) : Channel.empty()
ch_output_docs = file("$baseDir/docs/output.md", checkIfExists: true)
ch_output_docs_images = file("$baseDir/docs/images/", checkIfExists: true)

// JSON files required by BAMTools for alignment filtering
if (params.single_end) {
    ch_bamtools_filter_config = file(params.bamtools_filter_se_config, checkIfExists: true)
} else {
    ch_bamtools_filter_config = file(params.bamtools_filter_pe_config, checkIfExists: true)
}

// Header files for MultiQC
ch_mlib_peak_count_header = file("$baseDir/assets/multiqc/mlib_peak_count_header.txt", checkIfExists: true)
ch_mlib_frip_score_header = file("$baseDir/assets/multiqc/mlib_frip_score_header.txt", checkIfExists: true)
ch_mlib_peak_annotation_header = file("$baseDir/assets/multiqc/mlib_peak_annotation_header.txt", checkIfExists: true)
ch_mlib_deseq2_pca_header = file("$baseDir/assets/multiqc/mlib_deseq2_pca_header.txt", checkIfExists: true)
ch_mlib_deseq2_clustering_header = file("$baseDir/assets/multiqc/mlib_deseq2_clustering_header.txt", checkIfExists: true)

ch_mrep_peak_count_header = file("$baseDir/assets/multiqc/mrep_peak_count_header.txt", checkIfExists: true)
ch_mrep_frip_score_header = file("$baseDir/assets/multiqc/mrep_frip_score_header.txt", checkIfExists: true)
ch_mrep_peak_annotation_header = file("$baseDir/assets/multiqc/mrep_peak_annotation_header.txt", checkIfExists: true)
ch_mrep_deseq2_pca_header = file("$baseDir/assets/multiqc/mrep_deseq2_pca_header.txt", checkIfExists: true)
ch_mrep_deseq2_clustering_header = file("$baseDir/assets/multiqc/mrep_deseq2_clustering_header.txt", checkIfExists: true)

////////////////////////////////////////////////////
/* --          VALIDATE INPUTS                 -- */
////////////////////////////////////////////////////

if (params.input)     { ch_input = file(params.input, checkIfExists: true) } else { exit 1, 'Samples design file not specified!' }
if (params.gtf)       { ch_gtf = file(params.gtf, checkIfExists: true) } else { exit 1, 'GTF annotation file not specified!' }
if (params.gene_bed)  { ch_gene_bed = file(params.gene_bed, checkIfExists: true) }
if (params.tss_bed)   { ch_tss_bed = file(params.tss_bed, checkIfExists: true) }
if (params.blacklist) { ch_blacklist = Channel.fromPath(params.blacklist, checkIfExists: true) } else { ch_blacklist = Channel.empty() }

if (params.fasta) {
    lastPath = params.fasta.lastIndexOf(File.separator)
    bwa_base = params.fasta.substring(lastPath+1)
    ch_fasta = file(params.fasta, checkIfExists: true)
} else {
    exit 1, 'Fasta file not specified!'
}

if (params.bwa_index) {
    lastPath = params.bwa_index.lastIndexOf(File.separator)
    bwa_dir =  params.bwa_index.substring(0,lastPath+1)
    bwa_base = params.bwa_index.substring(lastPath+1)
    Channel
        .fromPath(bwa_dir, checkIfExists: true)
        .set { ch_bwa_index }
}

// Save AWS IGenomes file containing annotation version
if (params.anno_readme && file(params.anno_readme).exists()) {
    file("${params.outdir}/genome/").mkdirs()
    file(params.anno_readme).copyTo("${params.outdir}/genome/")
}

////////////////////////////////////////////////////
/* --                   AWS                    -- */
////////////////////////////////////////////////////

if (workflow.profile.contains('awsbatch')) {
    // AWSBatch sanity checking
    if (!params.awsqueue || !params.awsregion) exit 1, 'Specify correct --awsqueue and --awsregion parameters on AWSBatch!'
    // Check outdir paths to be S3 buckets if running on AWSBatch
    // related: https://github.com/nextflow-io/nextflow/issues/813
    if (!params.outdir.startsWith('s3:')) exit 1, 'Outdir not on S3 - specify S3 Bucket to run on AWSBatch!'
    // Prevent trace files to be stored on S3 since S3 does not support rolling files.
    if (params.tracedir.startsWith('s3:')) exit 1, 'Specify a local tracedir or run without trace! S3 cannot be used for tracefiles.'
}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
/* --                                                                     -- */
/* --                       HEADER LOG INFO                               -- */
/* --                                                                     -- */
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

// Header log info
log.info nfcoreHeader()
def summary = [:]
summary['Run Name']               = custom_runName ?: workflow.runName
summary['Data Type']              = params.single_end ? 'Single-End' : 'Paired-End'
summary['Design File']            = params.input
summary['Genome']                 = params.genome ?: 'Not supplied'
summary['Fasta File']             = params.fasta
summary['GTF File']               = params.gtf
if (params.gene_bed)              summary['Gene BED File'] = params.gene_bed
if (params.tss_bed)               summary['TSS BED File'] = params.tss_bed
if (params.bwa_index)             summary['BWA Index'] = params.bwa_index
if (params.blacklist)             summary['Blacklist BED'] = params.blacklist
if (params.mito_name)             summary['Mitochondrial Contig'] = params.mito_name
summary['MACS2 Genome Size']      = params.macs_gsize ?: 'Not supplied'
summary['Min Consensus Reps']     = params.min_reps_consensus
if (params.macs_gsize)            summary['MACS2 Narrow Peaks'] = params.narrow_peak ? 'Yes' : 'No'
if (!params.narrow_peak)          summary['MACS2 Broad Cutoff'] = params.broad_cutoff
if (params.skip_trimming) {
    summary['Trimming Step']      = 'Skipped'
} else {
    summary['Trim R1']            = "$params.clip_r1 bp"
    summary['Trim R2']            = "$params.clip_r2 bp"
    summary["Trim 3' R1"]         = "$params.three_prime_clip_r1 bp"
    summary["Trim 3' R2"]         = "$params.three_prime_clip_r2 bp"
    summary['NextSeq Trim']       = "$params.trim_nextseq bp"
}
if (params.seq_center)            summary['Sequencing Center'] = params.seq_center
if (params.single_end)            summary['Fragment Size'] = "$params.fragment_size bp"
summary['Fingerprint Bins']       = params.fingerprint_bins
if (params.keep_mito)             summary['Keep Mitochondrial'] = 'Yes'
if (params.keep_dups)             summary['Keep Duplicates'] = 'Yes'
if (params.keep_multi_map)        summary['Keep Multi-mapped'] = 'Yes'
summary['Save Genome Index']      = params.save_reference ? 'Yes' : 'No'
if (params.save_trimmed)          summary['Save Trimmed'] = 'Yes'
if (params.save_align_intermeds)  summary['Save Intermeds'] =  'Yes'
if (params.save_macs_pileup)      summary['Save MACS2 Pileup'] = 'Yes'
if (params.skip_merge_replicates) summary['Skip Merge Replicates'] = 'Yes'
if (params.skip_peak_qc)          summary['Skip MACS2 Peak QC'] = 'Yes'
if (params.skip_peak_annotation)  summary['Skip Peak Annotation'] = 'Yes'
if (params.skip_consensus_peaks)  summary['Skip Consensus Peaks'] = 'Yes'
if (params.deseq2_vst)            summary['Use DESeq2 vst Transform'] = 'Yes'
if (params.skip_diff_analysis)    summary['Skip Differential Analysis'] = 'Yes'
if (params.skip_fastqc)           summary['Skip FastQC'] = 'Yes'
if (params.skip_picard_metrics)   summary['Skip Picard Metrics'] = 'Yes'
if (params.skip_preseq)           summary['Skip Preseq'] = 'Yes'
if (params.skip_plot_profile)     summary['Skip plotProfile'] = 'Yes'
if (params.skip_plot_fingerprint) summary['Skip plotFingerprint'] = 'Yes'
if (params.skip_ataqv)            summary['Skip Ataqv'] = 'Yes'
if (params.skip_igv)              summary['Skip IGV'] = 'Yes'
if (params.skip_multiqc)          summary['Skip MultiQC'] = 'Yes'
summary['Max Resources']          = "$params.max_memory memory, $params.max_cpus cpus, $params.max_time time per job"
if (workflow.containerEngine)     summary['Container'] = "$workflow.containerEngine - $workflow.container"
summary['Output Dir']             = params.outdir
summary['Launch Dir']             = workflow.launchDir
summary['Working Dir']            = workflow.workDir
summary['Script Dir']             = workflow.projectDir
summary['User']                   = workflow.userName
if (workflow.profile.contains('awsbatch')) {
    summary['AWS Region']         = params.awsregion
    summary['AWS Queue']          = params.awsqueue
    summary['AWS CLI']            = params.awscli
}
summary['Config Profile']         = workflow.profile
if (params.config_profile_description) summary['Config Description'] = params.config_profile_description
if (params.config_profile_contact)     summary['Config Contact']     = params.config_profile_contact
if (params.config_profile_url)         summary['Config URL']         = params.config_profile_url
if (params.email || params.email_on_fail) {
    summary['E-mail Address']     = params.email
    summary['E-mail on failure']  = params.email_on_fail
    summary['MultiQC Max Size']   = params.max_multiqc_email_size
}
log.info summary.collect { k,v -> "${k.padRight(21)}: $v" }.join('\n')
log.info "-\033[2m--------------------------------------------------\033[0m-"

// Check the hostnames against configured profiles
checkHostname()

// Show a big warning message if we're not running MACS
if (!params.macs_gsize) {
    def warnstring = params.genome ? "supported for '${params.genome}'" : 'supplied'
    log.warn "=================================================================\n" +
             "  WARNING! MACS genome size parameter not $warnstring.\n" +
             "  Peak calling, annotation and differential analysis will be skipped.\n" +
             "  Please specify value for '--macs_gsize' to run these steps.\n" +
             "======================================================================="
}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
/* --                                                                     -- */
/* --                     PARSE DESIGN FILE                               -- */
/* --                                                                     -- */
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

/*
 * PREPROCESSING: Reformat design file and check validitiy
 */
process CHECK_DESIGN {
    tag "$design"
    publishDir "${params.outdir}/pipeline_info", mode: params.publish_dir_mode

    input:
    path design from ch_input

    output:
    path '*.csv' into ch_design_reads_csv

    script:  // This script is bundled with the pipeline, in nf-core/atacseq/bin/
    """
    check_design.py $design design_reads.csv
    """
}

/*
 * Create channels for input fastq files
 */
if (params.single_end) {
    ch_design_reads_csv
        .splitCsv(header:true, sep:',')
        .map { row -> [ row.sample_id, [ file(row.fastq_1, checkIfExists: true) ] ] }
        .into { ch_raw_reads_fastqc;
                ch_raw_reads_trimgalore;
                design_replicates_exist;
                design_multiple_samples }
} else {
    ch_design_reads_csv
        .splitCsv(header:true, sep:',')
        .map { row -> [ row.sample_id, [ file(row.fastq_1, checkIfExists: true), file(row.fastq_2, checkIfExists: true) ] ] }
        .into { ch_raw_reads_fastqc;
                ch_raw_reads_trimgalore;
                design_replicates_exist;
                design_multiple_samples }
}

// Boolean value for replicates existing in design
replicatesExist = design_replicates_exist
                      .map { it -> it[0].split('_')[-2].replaceAll('R','').toInteger() }
                      .flatten()
                      .max()
                      .val > 1

// Boolean value for multiple groups existing in design
multipleGroups = design_multiple_samples
                     .map { it -> it[0].split('_')[0..-3].join('_') }
                     .flatten()
                     .unique()
                     .count()
                     .val > 1

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
/* --                                                                     -- */
/* --                     PREPARE ANNOTATION FILES                        -- */
/* --                                                                     -- */
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

/*
 * PREPROCESSING: Build BWA index
 */
if (!params.bwa_index) {
    process BWA_INDEX {
        tag "$fasta"
        label 'process_high'
        publishDir path: { params.save_reference ? "${params.outdir}/genome" : params.outdir },
            saveAs: { params.save_reference ? it : null }, mode: params.publish_dir_mode

        input:
        path fasta from ch_fasta

        output:
        path 'BWAIndex' into ch_bwa_index

        script:
        """
        bwa index -a bwtsw $fasta
        mkdir BWAIndex && mv ${fasta}* BWAIndex
        """
    }
}

/*
 * PREPROCESSING: Generate gene BED file
 */
if (!params.gene_bed) {
    process MAKE_GENE_BED {
        tag "$gtf"
        label 'process_low'
        publishDir "${params.outdir}/genome", mode: params.publish_dir_mode

        input:
        path gtf from ch_gtf

        output:
        path '*.bed' into ch_gene_bed

        script: // This script is bundled with the pipeline, in nf-core/atacseq/bin/
        """
        gtf2bed $gtf > ${gtf.baseName}.bed
        """
    }
}

/*
 * PREPROCESSING: Generate TSS BED file
 */
if (!params.tss_bed) {
    process MAKE_TSS_BED {
        tag "$bed"
        publishDir "${params.outdir}/genome", mode: params.publish_dir_mode

        input:
        path bed from ch_gene_bed

        output:
        path '*.bed' into ch_tss_bed

        script:
        """
        cat $bed | awk -v FS='\t' -v OFS='\t' '{ if(\$6=="+") \$3=\$2+1; else \$2=\$3-1; print \$1, \$2, \$3, \$4, \$5, \$6;}' > ${bed.baseName}.tss.bed
        """
    }
}

/*
 * PREPROCESSING: Prepare genome intervals for filtering
 */
process MAKE_GENOME_FILTER {
    tag "$fasta"
    publishDir "${params.outdir}/genome", mode: params.publish_dir_mode

    input:
    path fasta from ch_fasta
    path blacklist from ch_blacklist.ifEmpty([])

    output:
    path "$fasta"                                      // FASTA FILE FOR IGV
    path '*.fai'                                       // FAI INDEX FOR REFERENCE GENOME
    path '*.bed' into ch_genome_filter_regions         // BED FILE WITHOUT BLACKLIST REGIONS & MITOCHONDRIAL CONTIG FOR FILTERING
    path '*.txt' into ch_genome_autosomes              // TEXT FILE CONTAINING LISTING OF AUTOSOMAL CHROMOSOMES FOR ATAQV
    path '*.sizes' into ch_genome_sizes_mlib_bigwig,   // CHROMOSOME SIZES FILE FOR BEDTOOLS
                        ch_genome_sizes_mrep_bigwig

    script:
    blacklist_filter = params.blacklist ? "sortBed -i $blacklist -g ${fasta}.sizes | complementBed -i stdin -g ${fasta}.sizes" : "awk '{print \$1, '0' , \$2}' OFS='\t' ${fasta}.sizes"
    name_filter = params.mito_name ? "| awk '\$1 !~ /${params.mito_name}/ {print \$0}'": ''
    mito_filter = params.keep_mito ? '' : name_filter
    """
    samtools faidx $fasta
    get_autosomes.py ${fasta}.fai ${fasta}.autosomes.txt
    cut -f 1,2 ${fasta}.fai > ${fasta}.sizes
    $blacklist_filter $mito_filter > ${fasta}.include_regions.bed
    """
}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
/* --                                                                     -- */
/* --                        FASTQ QC                                     -- */
/* --                                                                     -- */
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

/*
 * STEP 1: FastQC
 */
process FASTQC {
    tag "$name"
    label 'process_medium'
    publishDir "${params.outdir}/fastqc", mode: params.publish_dir_mode,
        saveAs: { filename ->
                      filename.endsWith('.zip') ? "zips/$filename" : filename
                }

    when:
    !params.skip_fastqc

    input:
    tuple val(name), path(reads) from ch_raw_reads_fastqc

    output:
    path '*.{zip,html}' into ch_fastqc_reports_mqc

    script:
    // Added soft-links to original fastqs for consistent naming in MultiQC
    if (params.single_end) {
        """
        [ ! -f  ${name}.fastq.gz ] && ln -s $reads ${name}.fastq.gz
        fastqc -q -t $task.cpus ${name}.fastq.gz
        """
    } else {
        """
        [ ! -f  ${name}_1.fastq.gz ] && ln -s ${reads[0]} ${name}_1.fastq.gz
        [ ! -f  ${name}_2.fastq.gz ] && ln -s ${reads[1]} ${name}_2.fastq.gz
        fastqc -q -t $task.cpus ${name}_1.fastq.gz
        fastqc -q -t $task.cpus ${name}_2.fastq.gz
        """
    }
}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
/* --                                                                     -- */
/* --                        ADAPTER TRIMMING                             -- */
/* --                                                                     -- */
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

/*
 * STEP 2: Trim Galore!
 */
if (params.skip_trimming) {
    ch_trimmed_reads = ch_raw_reads_trimgalore
    ch_trimgalore_results_mqc = Channel.empty()
    ch_trimgalore_fastqc_reports_mqc = Channel.empty()
} else {
    process TRIMGALORE {
        tag "$name"
        label 'process_high'
        publishDir "${params.outdir}/trim_galore", mode: params.publish_dir_mode,
            saveAs: { filename ->
                          if (filename.endsWith('.html')) "fastqc/$filename"
                          else if (filename.endsWith('.zip')) "fastqc/zips/$filename"
                          else if (filename.endsWith('trimming_report.txt')) "logs/$filename"
                          else params.save_trimmed ? filename : null
                    }

        input:
        tuple val(name), path(reads) from ch_raw_reads_trimgalore

        output:
        tuple val(name), path('*.fq.gz') into ch_trimmed_reads
        path '*.txt' into ch_trimgalore_results_mqc
        path '*.{zip,html}' into ch_trimgalore_fastqc_reports_mqc

        script:
        // Calculate number of --cores for TrimGalore based on value of task.cpus
        // See: https://github.com/FelixKrueger/TrimGalore/blob/master/Changelog.md#version-060-release-on-1-mar-2019
        // See: https://github.com/nf-core/atacseq/pull/65
        def cores = 1
        if (task.cpus) {
            cores = (task.cpus as int) - 4
            if (params.single_end) cores = (task.cpus as int) - 3
            if (cores < 1) cores = 1
            if (cores > 4) cores = 4
        }

        // Added soft-links to original fastqs for consistent naming in MultiQC
        c_r1 = params.clip_r1 > 0 ? "--clip_r1 ${params.clip_r1}" : ''
        c_r2 = params.clip_r2 > 0 ? "--clip_r2 ${params.clip_r2}" : ''
        tpc_r1 = params.three_prime_clip_r1 > 0 ? "--three_prime_clip_r1 ${params.three_prime_clip_r1}" : ''
        tpc_r2 = params.three_prime_clip_r2 > 0 ? "--three_prime_clip_r2 ${params.three_prime_clip_r2}" : ''
        nextseq = params.trim_nextseq > 0 ? "--nextseq ${params.trim_nextseq}" : ''

        // Added soft-links to original fastqs for consistent naming in MultiQC
        if (params.single_end) {
            """
            [ ! -f  ${name}.fastq.gz ] && ln -s $reads ${name}.fastq.gz
            trim_galore --cores $cores --fastqc --gzip $c_r1 $tpc_r1 $nextseq ${name}.fastq.gz
            """
        } else {
            """
            [ ! -f  ${name}_1.fastq.gz ] && ln -s ${reads[0]} ${name}_1.fastq.gz
            [ ! -f  ${name}_2.fastq.gz ] && ln -s ${reads[1]} ${name}_2.fastq.gz
            trim_galore --cores $cores --paired --fastqc --gzip $c_r1 $c_r2 $tpc_r1 $tpc_r2 $nextseq ${name}_1.fastq.gz ${name}_2.fastq.gz
            """
        }
    }
}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
/* --                                                                     -- */
/* --                        ALIGN                                        -- */
/* --                                                                     -- */
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

/*
 * STEP 3.1: Map read(s) with bwa mem
 */
process BWA_MEM {
    tag "$name"
    label 'process_high'

    input:
    tuple val(name), path(reads) from ch_trimmed_reads
    path index from ch_bwa_index.collect()

    output:
    tuple val(name), path('*.bam') into ch_bwa_bam

    script:
    prefix = "${name}.Lb"
    rg = "\'@RG\\tID:${name}\\tSM:${name.split('_')[0..-2].join('_')}\\tPL:ILLUMINA\\tLB:${name}\\tPU:1\'"
    if (params.seq_center) {
        rg = "\'@RG\\tID:${name}\\tSM:${name.split('_')[0..-2].join('_')}\\tPL:ILLUMINA\\tLB:${name}\\tPU:1\\tCN:${params.seq_center}\'"
    }
    """
    bwa mem \\
        -t $task.cpus \\
        -M \\
        -R $rg \\
        ${index}/${bwa_base} \\
        $reads \\
        | samtools view -@ $task.cpus -b -h -F 0x0100 -O BAM -o ${prefix}.bam -
    """
}

/*
 * STEP 3.2: Convert BAM to coordinate sorted BAM
 */
process SORT_BAM {
    tag "$name"
    label 'process_medium'
    if (params.save_align_intermeds) {
        publishDir path: "${params.outdir}/bwa/library", mode: params.publish_dir_mode,
            saveAs: { filename ->
                          if (filename.endsWith('.flagstat')) "samtools_stats/$filename"
                          else if (filename.endsWith('.idxstats')) "samtools_stats/$filename"
                          else if (filename.endsWith('.stats')) "samtools_stats/$filename"
                          else filename
                    }
    }

    input:
    tuple val(name), path(bam) from ch_bwa_bam

    output:
    tuple val(name), path('*.sorted.{bam,bam.bai}') into ch_sort_bam_merge
    path '*.{flagstat,idxstats,stats}' into ch_sort_bam_flagstat_mqc

    script:
    prefix = "${name}.Lb"
    """
    samtools sort -@ $task.cpus -o ${prefix}.sorted.bam -T $name $bam
    samtools index ${prefix}.sorted.bam
    samtools flagstat ${prefix}.sorted.bam > ${prefix}.sorted.bam.flagstat
    samtools idxstats ${prefix}.sorted.bam > ${prefix}.sorted.bam.idxstats
    samtools stats ${prefix}.sorted.bam > ${prefix}.sorted.bam.stats
    """
}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
/* --                                                                     -- */
/* --                    MERGE LIBRARY BAM                                -- */
/* --                                                                     -- */
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

/*
 * STEP 4.1: Merge BAM files for all libraries from same sample replicate
 */
ch_sort_bam_merge
    .map { it -> [ it[0].split('_')[0..-2].join('_'), it[1] ] }
    .groupTuple(by: [0])
    .map { it ->  [ it[0], it[1].flatten() ] }
    .set { ch_sort_bam_merge }

process MERGED_LIB_BAM {
    tag "$name"
    label 'process_medium'
    publishDir "${params.outdir}/bwa/mergedLibrary", mode: params.publish_dir_mode,
        saveAs: { filename ->
                      if (filename.endsWith('.flagstat')) "samtools_stats/$filename"
                      else if (filename.endsWith('.idxstats')) "samtools_stats/$filename"
                      else if (filename.endsWith('.stats')) "samtools_stats/$filename"
                      else if (filename.endsWith('.metrics.txt')) "picard_metrics/$filename"
                      else params.save_align_intermeds ? filename : null
                }

    input:
    tuple val(name), path(bams) from ch_sort_bam_merge

    output:
    tuple val(name), path("*${prefix}.sorted.{bam,bam.bai}") into ch_mlib_bam_filter,
                                                                  ch_mlib_bam_preseq,
                                                                  ch_mlib_bam_ataqv
    path '*.{flagstat,idxstats,stats}' into ch_mlib_bam_stats_mqc
    path '*.txt' into ch_mlib_bam_metrics_mqc

    script:
    prefix = "${name}.mLb.mkD"
    bam_files = bams.findAll { it.toString().endsWith('.bam') }.sort()
    def avail_mem = 3
    if (!task.memory) {
        log.info '[Picard MarkDuplicates] Available memory not known - defaulting to 3GB. Specify process memory requirements to change this.'
    } else {
        avail_mem = task.memory.toGiga()
    }
    if (bam_files.size() > 1) {
        """
        picard -Xmx${avail_mem}g MergeSamFiles \\
            ${'INPUT='+bam_files.join(' INPUT=')} \\
            OUTPUT=${name}.sorted.bam \\
            SORT_ORDER=coordinate \\
            VALIDATION_STRINGENCY=LENIENT \\
            TMP_DIR=tmp
        samtools index ${name}.sorted.bam

        picard -Xmx${avail_mem}g MarkDuplicates \\
            INPUT=${name}.sorted.bam \\
            OUTPUT=${prefix}.sorted.bam \\
            ASSUME_SORTED=true \\
            REMOVE_DUPLICATES=false \\
            METRICS_FILE=${prefix}.MarkDuplicates.metrics.txt \\
            VALIDATION_STRINGENCY=LENIENT \\
            TMP_DIR=tmp

        samtools index ${prefix}.sorted.bam
        samtools idxstats ${prefix}.sorted.bam > ${prefix}.sorted.bam.idxstats
        samtools flagstat ${prefix}.sorted.bam > ${prefix}.sorted.bam.flagstat
        samtools stats ${prefix}.sorted.bam > ${prefix}.sorted.bam.stats
        """
    } else {
      """
      picard -Xmx${avail_mem}g MarkDuplicates \\
          INPUT=${bam_files[0]} \\
          OUTPUT=${prefix}.sorted.bam \\
          ASSUME_SORTED=true \\
          REMOVE_DUPLICATES=false \\
          METRICS_FILE=${prefix}.MarkDuplicates.metrics.txt \\
          VALIDATION_STRINGENCY=LENIENT \\
          TMP_DIR=tmp

      samtools index ${prefix}.sorted.bam
      samtools idxstats ${prefix}.sorted.bam > ${prefix}.sorted.bam.idxstats
      samtools flagstat ${prefix}.sorted.bam > ${prefix}.sorted.bam.flagstat
      samtools stats ${prefix}.sorted.bam > ${prefix}.sorted.bam.stats
      """
    }
}

/*
 * STEP 4.2: Filter BAM file at merged library-level
 */
process MERGED_LIB_BAM_FILTER {
    tag "$name"
    label 'process_medium'
    publishDir path: "${params.outdir}/bwa/mergedLibrary", mode: params.publish_dir_mode,
        saveAs: { filename ->
                      if (params.single_end || params.save_align_intermeds) {
                          if (filename.endsWith('.flagstat')) "samtools_stats/$filename"
                          else if (filename.endsWith('.idxstats')) "samtools_stats/$filename"
                          else if (filename.endsWith('.stats')) "samtools_stats/$filename"
                          else if (filename.endsWith('.sorted.bam')) filename
                          else if (filename.endsWith('.sorted.bam.bai')) filename
                          else null
                      }
                }

    input:
    tuple val(name), path(bam) from ch_mlib_bam_filter
    path bed from ch_genome_filter_regions.collect()
    path bamtools_filter_config from ch_bamtools_filter_config

    output:
    tuple val(name), path('*.{bam,bam.bai}') into ch_mlib_filter_bam
    tuple val(name), path('*.flagstat') into ch_mlib_filter_bam_flagstat
    path '*.{idxstats,stats}' into ch_mlib_filter_bam_stats_mqc

    script:
    prefix = params.single_end ? "${name}.mLb.clN" : "${name}.mLb.flT"
    filter_params = params.single_end ? '-F 0x004' : '-F 0x004 -F 0x0008 -f 0x001'
    dup_params = params.keep_dups ? '' : '-F 0x0400'
    multimap_params = params.keep_multi_map ? '' : '-q 1'
    blacklist_params = params.blacklist ? "-L $bed" : ''
    name_sort_bam = params.single_end ? '' : "samtools sort -n -@ $task.cpus -o ${prefix}.bam -T $prefix ${prefix}.sorted.bam"
    """
    samtools view \\
        $filter_params \\
        $dup_params \\
        $multimap_params \\
        $blacklist_params \\
        -b ${bam[0]} \\
        | bamtools filter \\
            -out ${prefix}.sorted.bam \\
            -script $bamtools_filter_config

    samtools index ${prefix}.sorted.bam
    samtools flagstat ${prefix}.sorted.bam > ${prefix}.sorted.bam.flagstat
    samtools idxstats ${prefix}.sorted.bam > ${prefix}.sorted.bam.idxstats
    samtools stats ${prefix}.sorted.bam > ${prefix}.sorted.bam.stats

    $name_sort_bam
    """
}

/*
 * STEP 4.3: Remove orphan reads from paired-end BAM file
 */
if (params.single_end) {
    ch_mlib_filter_bam
        .into { ch_mlib_rm_orphan_bam_metrics;
                ch_mlib_rm_orphan_bam_bigwig;
                ch_mlib_rm_orphan_bam_macs;
                ch_mlib_rm_orphan_bam_plotfingerprint;
                ch_mlib_rm_orphan_bam_mrep;
                ch_mlib_name_bam_mlib_counts;
                ch_mlib_name_bam_mrep_counts }

    ch_mlib_filter_bam_flagstat
        .into { ch_mlib_rm_orphan_flagstat_bigwig;
                ch_mlib_rm_orphan_flagstat_macs;
                ch_mlib_rm_orphan_flagstat_mqc }

    ch_mlib_filter_bam_stats_mqc
        .set { ch_mlib_rm_orphan_stats_mqc }
} else {
    process MERGED_LIB_BAM_REMOVE_ORPHAN {
        tag "$name"
        label 'process_medium'
        publishDir path: "${params.outdir}/bwa/mergedLibrary", mode: params.publish_dir_mode,
            saveAs: { filename ->
                          if (filename.endsWith('.flagstat')) "samtools_stats/$filename"
                          else if (filename.endsWith('.idxstats')) "samtools_stats/$filename"
                          else if (filename.endsWith('.stats')) "samtools_stats/$filename"
                          else if (filename.endsWith('.sorted.bam')) filename
                          else if (filename.endsWith('.sorted.bam.bai')) filename
                          else null
                    }

        input:
        tuple val(name), path(bam) from ch_mlib_filter_bam

        output:
        tuple val(name), path('*.sorted.{bam,bam.bai}') into ch_mlib_rm_orphan_bam_metrics,
                                                             ch_mlib_rm_orphan_bam_bigwig,
                                                             ch_mlib_rm_orphan_bam_macs,
                                                             ch_mlib_rm_orphan_bam_plotfingerprint,
                                                             ch_mlib_rm_orphan_bam_mrep
        tuple val(name), path("${prefix}.bam") into ch_mlib_name_bam_mlib_counts,
                                                    ch_mlib_name_bam_mrep_counts
        tuple val(name), path('*.flagstat') into ch_mlib_rm_orphan_flagstat_bigwig,
                                                 ch_mlib_rm_orphan_flagstat_macs,
                                                 ch_mlib_rm_orphan_flagstat_mqc
        path '*.{idxstats,stats}' into ch_mlib_rm_orphan_stats_mqc

        script: // This script is bundled with the pipeline, in nf-core/atacseq/bin/
        prefix = "${name}.mLb.clN"
        """
        bampe_rm_orphan.py ${bam[0]} ${prefix}.bam --only_fr_pairs

        samtools sort -@ $task.cpus -o ${prefix}.sorted.bam -T $prefix ${prefix}.bam
        samtools index ${prefix}.sorted.bam
        samtools flagstat ${prefix}.sorted.bam > ${prefix}.sorted.bam.flagstat
        samtools idxstats ${prefix}.sorted.bam > ${prefix}.sorted.bam.idxstats
        samtools stats ${prefix}.sorted.bam > ${prefix}.sorted.bam.stats
        """
    }
}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
/* --                                                                     -- */
/* --                 MERGE LIBRARY BAM POST-ANALYSIS                     -- */
/* --                                                                     -- */
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

/*
 * STEP 5.1: Preseq analysis after merging libraries and before filtering
 */
process MERGED_LIB_PRESEQ {
    tag "$name"
    label 'process_low'
    label 'error_ignore'
    publishDir "${params.outdir}/bwa/mergedLibrary/preseq", mode: params.publish_dir_mode

    when:
    !params.skip_preseq

    input:
    tuple val(name), path(bam) from ch_mlib_bam_preseq

    output:
    path '*.ccurve.txt' into ch_mlib_preseq_mqc
    path '*.log'

    script:
    prefix = "${name}.mLb.mkD"
    pe = params.single_end ? '' : '-pe'
    """
    preseq lc_extrap \\
        -output ${prefix}.ccurve.txt \\
        -verbose \\
        -bam \\
        $pe \\
        -seed 1 \\
        ${bam[0]}
    cp .command.err ${prefix}.command.log
    """
}

/*
 * STEP 5.2: Picard CollectMultipleMetrics after merging libraries and filtering
 */
process MERGED_LIB_PICARD_METRICS {
    tag "$name"
    label 'process_medium'
    publishDir path: "${params.outdir}/bwa/mergedLibrary", mode: params.publish_dir_mode,
        saveAs: { filename ->
                      if (filename.endsWith('_metrics')) "picard_metrics/$filename"
                      else if (filename.endsWith('.pdf')) "picard_metrics/pdf/$filename"
                      else null
                }

    when:
    !params.skip_picard_metrics

    input:
    tuple val(name), path(bam) from ch_mlib_rm_orphan_bam_metrics
    path fasta from ch_fasta

    output:
    path '*_metrics' into ch_mlib_collectmetrics_mqc
    path '*.pdf'

    script:
    prefix = "${name}.mLb.clN"
    def avail_mem = 3
    if (!task.memory) {
        log.info '[Picard MarkDuplicates] Available memory not known - defaulting to 3GB. Specify process memory requirements to change this.'
    } else {
        avail_mem = task.memory.toGiga()
    }
    """
    picard -Xmx${avail_mem}g CollectMultipleMetrics \\
        INPUT=${bam[0]} \\
        OUTPUT=${prefix}.CollectMultipleMetrics \\
        REFERENCE_SEQUENCE=$fasta \\
        VALIDATION_STRINGENCY=LENIENT \\
        TMP_DIR=tmp
    """
}

/*
 * STEP 5.3: Read depth normalised bigWig
 */
process MERGED_LIB_BIGWIG {
    tag "$name"
    label 'process_medium'
    publishDir "${params.outdir}/bwa/mergedLibrary/bigwig", mode: params.publish_dir_mode,
        saveAs: { filename ->
                      if (filename.endsWith('scale_factor.txt')) "scale/$filename"
                      else if (filename.endsWith('.bigWig')) filename
                      else null
                }

    input:
    tuple val(name), path(bam), path(flagstat) from ch_mlib_rm_orphan_bam_bigwig.join(ch_mlib_rm_orphan_flagstat_bigwig, by: [0])
    path sizes from ch_genome_sizes_mlib_bigwig.collect()

    output:
    tuple val(name), path('*.bigWig') into ch_mlib_bigwig_plotprofile
    path '*igv.txt' into ch_mlib_bigwig_igv
    path '*scale_factor.txt'

    script:
    prefix = "${name}.mLb.clN"
    pe_fragment = params.single_end ? '' : '-pc'
    extend = (params.single_end && params.fragment_size > 0) ? "-fs ${params.fragment_size}" : ''
    """
    SCALE_FACTOR=\$(grep 'mapped (' $flagstat | awk '{print 1000000/\$1}')
    echo \$SCALE_FACTOR > ${prefix}.scale_factor.txt
    genomeCoverageBed -ibam ${bam[0]} -bg -scale \$SCALE_FACTOR $pe_fragment $extend | sort -T '.' -k1,1 -k2,2n >  ${prefix}.bedGraph

    bedGraphToBigWig ${prefix}.bedGraph $sizes ${prefix}.bigWig

    find * -type f -name "*.bigWig" -exec echo -e "bwa/mergedLibrary/bigwig/"{}"\\t0,0,178" \\; > ${prefix}.bigWig.igv.txt
    """
}

/*
 * STEP 5.4: Generate gene body coverage plot with deepTools plotProfile and plotHeatmap
 */
process MERGED_LIB_PLOTPROFILE {
    tag "$name"
    label 'process_high'
    publishDir "${params.outdir}/bwa/mergedLibrary/deepTools/plotProfile", mode: params.publish_dir_mode

    when:
    !params.skip_plot_profile

    input:
    tuple val(name), path(bigwig) from ch_mlib_bigwig_plotprofile
    path bed from ch_gene_bed

    output:
    path '*.plotProfile.tab' into ch_mlib_plotprofile_mqc
    path '*.{gz,pdf,mat.tab}'

    script:
    prefix = "${name}.mLb.clN"
    """
    computeMatrix scale-regions \\
        --regionsFileName $bed \\
        --scoreFileName $bigwig \\
        --outFileName ${prefix}.computeMatrix.mat.gz \\
        --outFileNameMatrix ${prefix}.computeMatrix.vals.mat.tab \\
        --regionBodyLength 1000 \\
        --beforeRegionStartLength 3000 \\
        --afterRegionStartLength 3000 \\
        --skipZeros \\
        --smartLabels \\
        --numberOfProcessors $task.cpus

    plotProfile --matrixFile ${prefix}.computeMatrix.mat.gz \\
        --outFileName ${prefix}.plotProfile.pdf \\
        --outFileNameData ${prefix}.plotProfile.tab

    plotHeatmap --matrixFile ${prefix}.computeMatrix.mat.gz \\
        --outFileName ${prefix}.plotHeatmap.pdf \\
        --outFileNameMatrix ${prefix}.plotHeatmap.mat.tab
    """
}

/*
 * STEP 5.5: deepTools plotFingerprint
 */
process MERGED_LIB_PLOTFINGERPRINT {
    tag "$name"
    label 'process_high'
    publishDir "${params.outdir}/bwa/mergedLibrary/deepTools/plotFingerprint", mode: params.publish_dir_mode

    when:
    !params.skip_plot_fingerprint

    input:
    tuple val(name), path(bam) from ch_mlib_rm_orphan_bam_plotfingerprint

    output:
    path '*.raw.txt' into ch_mlib_plotfingerprint_mqc
    path '*.{txt,pdf}'

    script:
    prefix = "${name}.mLb.clN"
    extend = (params.single_end && params.fragment_size > 0) ? "--extendReads ${params.fragment_size}" : ''
    """
    plotFingerprint \\
        --bamfiles ${bam[0]} \\
        --plotFile ${prefix}.plotFingerprint.pdf \\
        $extend \\
        --labels $prefix \\
        --outRawCounts ${prefix}.plotFingerprint.raw.txt \\
        --outQualityMetrics ${prefix}.plotFingerprint.qcmetrics.txt \\
        --skipZeros \\
        --numberOfProcessors $task.cpus \\
        --numberOfSamples $params.fingerprint_bins
    """
}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
/* --                                                                     -- */
/* --                 MERGE LIBRARY PEAK ANALYSIS                         -- */
/* --                                                                     -- */
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

/*
 * STEP 6.1: Call peaks with MACS2 and calculate FRiP score
 */
process MERGED_LIB_MACS2 {
    tag "$name"
    label 'process_medium'
    publishDir "${params.outdir}/bwa/mergedLibrary/macs/${PEAK_TYPE}", mode: params.publish_dir_mode,
        saveAs: { filename ->
                      if (filename.endsWith('.tsv')) "qc/$filename"
                      else if (filename.endsWith('.igv.txt')) null
                      else filename
                }

    when:
    params.macs_gsize

    input:
    tuple val(name), path(bam), path(flagstat) from ch_mlib_rm_orphan_bam_macs.join(ch_mlib_rm_orphan_flagstat_macs, by: [0])
    path mlib_peak_count_header from ch_mlib_peak_count_header
    path mlib_frip_score_header from ch_mlib_frip_score_header

    output:
    tuple val(name), path("*$PEAK_TYPE") into ch_mlib_macs_homer,
                                              ch_mlib_macs_qc,
                                              ch_mlib_macs_consensus,
                                              ch_mlib_macs_ataqv
    path '*igv.txt' into ch_mlib_macs_igv
    path '*_mqc.tsv' into ch_mlib_macs_mqc
    path '*.{bed,xls,gappedPeak,bdg}'

    script:
    prefix = "${name}.mLb.clN"
    broad = params.narrow_peak ? '' : "--broad --broad-cutoff ${params.broad_cutoff}"
    format = params.single_end ? 'BAM' : 'BAMPE'
    pileup = params.save_macs_pileup ? '-B --SPMR' : ''
    """
    macs2 callpeak \\
        -t ${bam[0]} \\
        $broad \\
        -f $format \\
        -g $params.macs_gsize \\
        -n $prefix \\
        $pileup \\
        --keep-dup all \\
        --nomodel

    cat ${prefix}_peaks.${PEAK_TYPE} | wc -l | awk -v OFS='\t' '{ print "${name}", \$1 }' | cat $mlib_peak_count_header - > ${prefix}_peaks.count_mqc.tsv

    READS_IN_PEAKS=\$(intersectBed -a ${bam[0]} -b ${prefix}_peaks.${PEAK_TYPE} -bed -c -f 0.20 | awk -F '\t' '{sum += \$NF} END {print sum}')
    grep 'mapped (' $flagstat | awk -v a="\$READS_IN_PEAKS" -v OFS='\t' '{print "${name}", a/\$1}' | cat $mlib_frip_score_header - > ${prefix}_peaks.FRiP_mqc.tsv

    find * -type f -name "*.${PEAK_TYPE}" -exec echo -e "bwa/mergedLibrary/macs/${PEAK_TYPE}/"{}"\\t0,0,178" \\; > ${prefix}_peaks.igv.txt
    """
}

/*
 * STEP 6.2: Annotate peaks with HOMER
 */
process MERGED_LIB_MACS2_ANNOTATE {
    tag "$name"
    label 'process_medium'
    publishDir "${params.outdir}/bwa/mergedLibrary/macs/${PEAK_TYPE}", mode: params.publish_dir_mode

    when:
    params.macs_gsize && !params.skip_peak_annotation

    input:
    tuple val(name), path(peak) from ch_mlib_macs_homer
    path fasta from ch_fasta
    path gtf from ch_gtf

    output:
    path '*.txt' into ch_mlib_macs_annotate

    script:
    prefix = "${name}.mLb.clN"
    """
    annotatePeaks.pl \\
        $peak \\
        $fasta \\
        -gid \\
        -gtf $gtf \\
        -cpu $task.cpus \\
        > ${prefix}_peaks.annotatePeaks.txt
    """
}

/*
 * STEP 6.3: Aggregated QC plots for peaks, FRiP and peak-to-gene annotation
 */
process MERGED_LIB_MACS2_QC {
    label 'process_medium'
    publishDir "${params.outdir}/bwa/mergedLibrary/macs/${PEAK_TYPE}/qc", mode: params.publish_dir_mode

    when:
    params.macs_gsize && !params.skip_peak_annotation && !params.skip_peak_qc

    input:
    path peaks from ch_mlib_macs_qc.collect{ it[1] }
    path annos from ch_mlib_macs_annotate.collect()
    path mlib_peak_annotation_header from ch_mlib_peak_annotation_header

    output:
    path '*.tsv' into ch_mlib_peak_qc_mqc
    path '*.{txt,pdf}'

    script:  // This script is bundled with the pipeline, in nf-core/atacseq/bin/
    suffix = 'mLb.clN'
    """
    plot_macs_qc.r \\
        -i ${peaks.join(',')} \\
        -s ${peaks.join(',').replaceAll(".${suffix}_peaks.${PEAK_TYPE}","")} \\
        -o ./ \\
        -p macs_peak.${suffix}

    plot_homer_annotatepeaks.r \\
        -i ${annos.join(',')} \\
        -s ${annos.join(',').replaceAll(".${suffix}_peaks.annotatePeaks.txt","")} \\
        -o ./ \\
        -p macs_annotatePeaks.${suffix}

    cat $mlib_peak_annotation_header macs_annotatePeaks.${suffix}.summary.txt > macs_annotatePeaks.${suffix}.summary_mqc.tsv
    """
}

/*
 * STEP 6.4: Consensus peaks across samples, create boolean filtering file, SAF file for featureCounts and UpSetR plot for intersection
 */
process MERGED_LIB_CONSENSUS {
    label 'process_long'
    publishDir "${params.outdir}/bwa/mergedLibrary/macs/${PEAK_TYPE}/consensus", mode: params.publish_dir_mode,
        saveAs: { filename ->
                      if (filename.endsWith('.igv.txt')) null
                      else filename
                }

    when:
    params.macs_gsize && (replicatesExist || multipleGroups) && !params.skip_consensus_peaks

    input:
    path peaks from ch_mlib_macs_consensus.collect{ it[1] }

    output:
    path '*.bed' into ch_mlib_macs_consensus_bed
    path '*.saf' into ch_mlib_macs_consensus_saf
    path '*.boolean.txt' into ch_mlib_macs_consensus_bool
    path '*igv.txt' into ch_mlib_macs_consensus_igv
    path '*.intersect.{txt,plot.pdf}'

    script: // scripts are bundled with the pipeline, in nf-core/atacseq/bin/
    suffix = 'mLb.clN'
    prefix = "consensus_peaks.${suffix}"
    mergecols = params.narrow_peak ? (2..10).join(',') : (2..9).join(',')
    collapsecols = params.narrow_peak ? (['collapse']*9).join(',') : (['collapse']*8).join(',')
    expandparam = params.narrow_peak ? '--is_narrow_peak' : ''
    """
    sort -T '.' -k1,1 -k2,2n ${peaks.collect{it.toString()}.sort().join(' ')} \\
        | mergeBed -c $mergecols -o $collapsecols > ${prefix}.txt

    macs2_merged_expand.py \\
        ${prefix}.txt \\
        ${peaks.collect{it.toString()}.sort().join(',').replaceAll("_peaks.${PEAK_TYPE}","")} \\
        ${prefix}.boolean.txt \\
        --min_replicates $params.min_reps_consensus \\
        $expandparam

    awk -v FS='\t' -v OFS='\t' 'FNR > 1 { print \$1, \$2, \$3, \$4, "0", "+" }' ${prefix}.boolean.txt > ${prefix}.bed

    echo -e "GeneID\tChr\tStart\tEnd\tStrand" > ${prefix}.saf
    awk -v FS='\t' -v OFS='\t' 'FNR > 1 { print \$4, \$1, \$2, \$3,  "+" }' ${prefix}.boolean.txt >> ${prefix}.saf

    sed -i 's/.${suffix}//g' ${prefix}.boolean.intersect.txt
    plot_peak_intersect.r -i ${prefix}.boolean.intersect.txt -o ${prefix}.boolean.intersect.plot.pdf

    find * -type f -name "${prefix}.bed" -exec echo -e "bwa/mergedLibrary/macs/${PEAK_TYPE}/consensus/"{}"\\t0,0,0" \\; > ${prefix}.bed.igv.txt
    """
}

/*
 * STEP 6.5: Annotate consensus peaks with HOMER, and add annotation to boolean output file
 */
process MERGED_LIB_CONSENSUS_ANNOTATE {
    label 'process_medium'
    publishDir "${params.outdir}/bwa/mergedLibrary/macs/${PEAK_TYPE}/consensus", mode: params.publish_dir_mode

    when:
    params.macs_gsize && (replicatesExist || multipleGroups) && !params.skip_consensus_peaks && !params.skip_peak_annotation

    input:
    path bed from ch_mlib_macs_consensus_bed
    path bool from ch_mlib_macs_consensus_bool
    path fasta from ch_fasta
    path gtf from ch_gtf

    output:
    path '*.annotatePeaks.txt'

    script:
    prefix = 'consensus_peaks.mLb.clN'
    """
    annotatePeaks.pl \\
        $bed \\
        $fasta \\
        -gid \\
        -gtf $gtf \\
        -cpu $task.cpus \\
        > ${prefix}.annotatePeaks.txt

    cut -f2- ${prefix}.annotatePeaks.txt | awk 'NR==1; NR > 1 {print \$0 | "sort -T '.' -k1,1 -k2,2n"}' | cut -f6- > tmp.txt
    paste $bool tmp.txt > ${prefix}.boolean.annotatePeaks.txt
    """
}

/*
 * STEP 6.6: Count reads in consensus peaks with featureCounts
 */
process MERGED_LIB_CONSENSUS_COUNTS {
    label 'process_medium'
    publishDir "${params.outdir}/bwa/mergedLibrary/macs/${PEAK_TYPE}/consensus", mode: params.publish_dir_mode

    when:
    params.macs_gsize && (replicatesExist || multipleGroups) && !params.skip_consensus_peaks

    input:
    path bams from ch_mlib_name_bam_mlib_counts.collect{ it[1] }
    path saf from ch_mlib_macs_consensus_saf.collect()

    output:
    path '*featureCounts.txt' into ch_mlib_macs_consensus_counts
    path '*featureCounts.txt.summary' into ch_mlib_macs_consensus_counts_mqc

    script:
    prefix = 'consensus_peaks.mLb.clN'
    bam_files = bams.findAll { it.toString().endsWith('.bam') }.sort()
    pe_params = params.single_end ? '' : '-p --donotsort'
    """
    featureCounts \\
        -F SAF \\
        -O \\
        --fracOverlap 0.2 \\
        -T $task.cpus \\
        $pe_params \\
        -a $saf \\
        -o ${prefix}.featureCounts.txt \\
        ${bam_files.join(' ')}
    """
}

/*
 * STEP 6.7: Differential analysis with DESeq2
 */
process MERGED_LIB_CONSENSUS_DESEQ2 {
    label 'process_medium'
    publishDir "${params.outdir}/bwa/mergedLibrary/macs/${PEAK_TYPE}/consensus/deseq2", mode: params.publish_dir_mode,
        saveAs: { filename ->
                      if (filename.endsWith('.igv.txt')) null
                      else filename
                }

    when:
    params.macs_gsize && replicatesExist && multipleGroups && !params.skip_consensus_peaks && !params.skip_diff_analysis

    input:
    path counts from ch_mlib_macs_consensus_counts
    path mlib_deseq2_pca_header from ch_mlib_deseq2_pca_header
    path mlib_deseq2_clustering_header from ch_mlib_deseq2_clustering_header

    output:
    path '*.tsv' into ch_mlib_macs_consensus_deseq_mqc
    path '*igv.txt' into ch_mlib_macs_consensus_deseq_comp_igv
    path '*.{RData,results.txt,pdf,log}'
    path 'sizeFactors'
    path '*vs*/*.{pdf,txt}'
    path '*vs*/*.bed'

    script:
    prefix = 'consensus_peaks.mLb.clN'
    bam_ext = params.single_end ? '.mLb.clN.sorted.bam' : '.mLb.clN.bam'
    vst = params.deseq2_vst ? '--vst TRUE' : ''
    """
    featurecounts_deseq2.r \\
        --featurecount_file $counts \\
        --bam_suffix '$bam_ext' \\
        --outdir ./ \\
        --outprefix $prefix \\
        --outsuffix .mLb.clN \\
        --cores $task.cpus \\
        $vst \\

    cat $mlib_deseq2_pca_header ${prefix}.pca.vals.txt > ${prefix}.pca.vals_mqc.tsv
    cat $mlib_deseq2_clustering_header ${prefix}.sample.dists.txt > ${prefix}.sample.dists_mqc.tsv

    find * -type f -name "*.FDR0.05.results.bed" -exec echo -e "bwa/mergedLibrary/macs/${PEAK_TYPE}/consensus/deseq2/"{}"\\t255,0,0" \\; > ${prefix}.igv.txt
    """
}

/*
 * STEP 6.8: Run ataqv on BAM file and corresponding peaks
 */
process MERGED_LIB_ATAQV {
    tag "$name"
    label 'process_medium'
    publishDir "${params.outdir}/bwa/mergedLibrary/ataqv/${PEAK_TYPE}", mode: params.publish_dir_mode

    when:
    !params.skip_ataqv

    input:
    tuple val(name), path(bam), path(peak) from ch_mlib_bam_ataqv.join(ch_mlib_macs_ataqv, by: [0])
    path autosomes from ch_genome_autosomes.collect()
    path tss_bed from ch_tss_bed

    output:
    path '*.json' into ch_mlib_ataqv

    script:
    suffix = 'mLb.clN'
    peak_param = params.macs_gsize ? "--peak-file ${peak}" : ''
    mito_param = params.mito_name ? "--mitochondrial-reference-name ${params.mito_name}" : ''
    """
    ataqv \\
        --threads $task.cpus \\
        $peak_param  \\
        --tss-file $tss_bed \\
        --metrics-file  ${name}.${suffix}.ataqv.json \\
        --name $name \\
        --ignore-read-groups \\
        --autosomal-reference-file $autosomes \\
        $mito_param \\
        NA \\
        ${bam[0]}
    """
}

/*
 * STEP 6.9: Run ataqv mkarv on all JSON files to render web app
 */
process MERGED_LIB_ATAQV_MKARV {
    label 'process_medium'
    publishDir "${params.outdir}/bwa/mergedLibrary/ataqv/${PEAK_TYPE}", mode: params.publish_dir_mode

    when:
    !params.skip_ataqv

    input:
    path json from ch_mlib_ataqv.collect()

    output:
    path 'html'

    script:
    """
    mkarv \\
        --concurrency $task.cpus \\
        --force \\
        ./html/ \\
        ${json.join(' ')}
    """
}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
/* --                                                                     -- */
/* --                    MERGE REPLICATE BAM                              -- */
/* --                                                                     -- */
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

/*
 * STEP 7: Merge library BAM files across all replicates
 */
ch_mlib_rm_orphan_bam_mrep
    .map { it -> [ it[0].split('_')[0..-2].join('_'), it[1] ] }
    .groupTuple(by: [0])
    .map { it ->  [ it[0], it[1].flatten() ] }
    .set { ch_mlib_rm_orphan_bam_mrep }

process MERGED_REP_BAM {
    tag "$name"
    label 'process_medium'
    publishDir "${params.outdir}/bwa/mergedReplicate", mode: params.publish_dir_mode,
        saveAs: { filename ->
                      if (filename.endsWith('.flagstat')) "samtools_stats/$filename"
                      else if (filename.endsWith('.idxstats')) "samtools_stats/$filename"
                      else if (filename.endsWith('.stats')) "samtools_stats/$filename"
                      else if (filename.endsWith('.metrics.txt')) "picard_metrics/$filename"
                      else filename
                }

    input:
    tuple val(name), path(bams) from ch_mlib_rm_orphan_bam_mrep

    output:
    tuple val(name), path("*${prefix}.sorted.{bam,bam.bai}") into ch_mrep_bam_bigwig,
                                                                  ch_mrep_bam_macs
    tuple val(name), path('*.flagstat') into ch_mrep_bam_flagstat_bigwig,
                                             ch_mrep_bam_flagstat_macs,
                                             ch_mrep_bam_flagstat_mqc
    path '*.{idxstats,stats}' into ch_mrep_bam_stats_mqc
    path '*.txt' into ch_mrep_bam_metrics_mqc

    when:
    !params.skip_merge_replicates && replicatesExist

    script:
    prefix = "${name}.mRp.clN"
    bam_files = bams.findAll { it.toString().endsWith('.bam') }.sort()
    def avail_mem = 3
    if (!task.memory) {
        log.info '[Picard MarkDuplicates] Available memory not known - defaulting to 3GB. Specify process memory requirements to change this.'
    } else {
        avail_mem = task.memory.toGiga()
    }
    if (bam_files.size() > 1) {
        """
        picard -Xmx${avail_mem}g MergeSamFiles \\
            ${'INPUT='+bam_files.join(' INPUT=')} \\
            OUTPUT=${name}.sorted.bam \\
            SORT_ORDER=coordinate \\
            VALIDATION_STRINGENCY=LENIENT \\
            TMP_DIR=tmp
        samtools index ${name}.sorted.bam

        picard -Xmx${avail_mem}g MarkDuplicates \\
            INPUT=${name}.sorted.bam \\
            OUTPUT=${prefix}.sorted.bam \\
            ASSUME_SORTED=true \\
            REMOVE_DUPLICATES=true \\
            METRICS_FILE=${prefix}.MarkDuplicates.metrics.txt \\
            VALIDATION_STRINGENCY=LENIENT \\
            TMP_DIR=tmp

        samtools index ${prefix}.sorted.bam
        samtools flagstat ${prefix}.sorted.bam > ${prefix}.sorted.bam.flagstat
        samtools idxstats ${prefix}.sorted.bam > ${prefix}.sorted.bam.idxstats
        samtools stats ${prefix}.sorted.bam > ${prefix}.sorted.bam.stats
        """
    } else {
      """
      ln -s ${bams[0]} ${prefix}.sorted.bam
      ln -s ${bams[1]} ${prefix}.sorted.bam.bai
      touch ${prefix}.MarkDuplicates.metrics.txt
      samtools flagstat ${prefix}.sorted.bam > ${prefix}.sorted.bam.flagstat
      samtools idxstats ${prefix}.sorted.bam > ${prefix}.sorted.bam.idxstats
      samtools stats ${prefix}.sorted.bam > ${prefix}.sorted.bam.stats
      """
    }
}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
/* --                                                                     -- */
/* --              MERGE REPLICATE BAM POST-ANALYSIS                      -- */
/* --                                                                     -- */
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

/*
 * STEP 8.1: Read depth normalised bigWig
 */
process MERGED_REP_BIGWIG {
    tag "$name"
    label 'process_medium'
    publishDir "${params.outdir}/bwa/mergedReplicate/bigwig", mode: params.publish_dir_mode,
        saveAs: { filename ->
                      if (filename.endsWith('scale_factor.txt')) "scale/$filename"
                      else if (filename.endsWith('.bigWig')) filename
                      else null
                }

    when:
    !params.skip_merge_replicates && replicatesExist

    input:
    tuple val(name), path(bam), path(flagstat) from ch_mrep_bam_bigwig.join(ch_mrep_bam_flagstat_bigwig, by: [0])
    path sizes from ch_genome_sizes_mrep_bigwig.collect()

    output:
    tuple val(name), path('*.bigWig') into ch_mrep_bigwig
    path '*igv.txt' into ch_mrep_bigwig_igv
    path '*scale_factor.txt'

    script:
    prefix = "${name}.mRp.clN"
    pe_fragment = params.single_end ? '' : '-pc'
    extend = (params.single_end && params.fragment_size > 0) ? "-fs ${params.fragment_size}" : ''
    """
    SCALE_FACTOR=\$(grep 'mapped (' $flagstat | awk '{print 1000000/\$1}')
    echo \$SCALE_FACTOR > ${prefix}.scale_factor.txt
    genomeCoverageBed -ibam ${bam[0]} -bg -scale \$SCALE_FACTOR $pe_fragment $extend | sort -T '.' -k1,1 -k2,2n >  ${prefix}.bedGraph

    bedGraphToBigWig ${prefix}.bedGraph $sizes ${prefix}.bigWig

    find * -type f -name "*.bigWig" -exec echo -e "bwa/mergedReplicate/bigwig/"{}"\\t0,0,178" \\; > ${prefix}.bigWig.igv.txt
    """
}

/*
 * STEP 8.2: Call peaks with MACS2 and calculate FRiP score
 */
process MERGED_REP_MACS2 {
    tag "$name"
    label 'process_medium'
    publishDir "${params.outdir}/bwa/mergedReplicate/macs/${PEAK_TYPE}", mode: params.publish_dir_mode,
        saveAs: { filename ->
                      if (filename.endsWith('.tsv')) "qc/$filename"
                      else if (filename.endsWith('.igv.txt')) null
                      else filename
                }

    when:
    !params.skip_merge_replicates && replicatesExist && params.macs_gsize

    input:
    tuple val(name), path(bam), path(flagstat) from ch_mrep_bam_macs.join(ch_mrep_bam_flagstat_macs, by: [0])
    path mrep_peak_count_header from ch_mrep_peak_count_header
    path mrep_frip_score_header from ch_mrep_frip_score_header

    output:
    tuple val(name), path("*$PEAK_TYPE") into ch_mrep_macs_homer,
                                              ch_mrep_macs_qc,
                                              ch_mrep_macs_consensus
    path '*igv.txt' into ch_mrep_macs_igv
    path '*_mqc.tsv' into ch_mrep_macs_mqc
    path '*.{bed,xls,gappedPeak,bdg}'

    script:
    prefix = "${name}.mRp.clN"
    broad = params.narrow_peak ? '' : "--broad --broad-cutoff ${params.broad_cutoff}"
    format = params.single_end ? 'BAM' : 'BAMPE'
    pileup = params.save_macs_pileup ? '-B --SPMR' : ''
    """
    macs2 callpeak \\
        -t ${bam[0]} \\
        $broad \\
        -f $format \\
        -g $params.macs_gsize \\
        -n $prefix \\
        $pileup \\
        --keep-dup all \\
        --nomodel

    cat ${prefix}_peaks.${PEAK_TYPE} | wc -l | awk -v OFS='\t' '{ print "${name}", \$1 }' | cat $mrep_peak_count_header - > ${prefix}_peaks.count_mqc.tsv

    READS_IN_PEAKS=\$(intersectBed -a ${bam[0]} -b ${prefix}_peaks.${PEAK_TYPE} -bed -c -f 0.20 | awk -F '\t' '{sum += \$NF} END {print sum}')
    grep 'mapped (' $flagstat | awk -v a="\$READS_IN_PEAKS" -v OFS='\t' '{print "${name}", a/\$1}' | cat $mrep_frip_score_header - > ${prefix}_peaks.FRiP_mqc.tsv

    find * -type f -name "*.${PEAK_TYPE}" -exec echo -e "bwa/mergedReplicate/macs/${PEAK_TYPE}/"{}"\\t0,0,178" \\; > ${prefix}_peaks.igv.txt
    """
}

/*
 * STEP 8.3: Annotate peaks with HOMER
 */
process MERGED_REP_MACS2_ANNOTATE {
    tag "$name"
    label 'process_medium'
    publishDir "${params.outdir}/bwa/mergedReplicate/macs/${PEAK_TYPE}", mode: params.publish_dir_mode

    when:
    !params.skip_merge_replicates && replicatesExist && params.macs_gsize && !params.skip_peak_annotation

    input:
    tuple val(name), path(peak) from ch_mrep_macs_homer
    path fasta from ch_fasta
    path gtf from ch_gtf

    output:
    path '*.txt' into ch_mrep_macs_annotate

    script:
    prefix = "${name}.mRp.clN"
    """
    annotatePeaks.pl \\
        $peak \\
        $fasta \\
        -gid \\
        -gtf $gtf \\
        -cpu $task.cpus \\
        > ${prefix}_peaks.annotatePeaks.txt
    """
}

/*
 * STEP 8.4: Aggregated QC plots for peaks, FRiP and peak-to-gene annotation
 */
process MERGED_REP_MACS2_QC {
    label 'process_medium'
    publishDir "${params.outdir}/bwa/mergedReplicate/macs/${PEAK_TYPE}/qc", mode: params.publish_dir_mode

    when:
    !params.skip_merge_replicates && replicatesExist && params.macs_gsize && !params.skip_peak_qc && !params.skip_peak_annotation

    input:
    path peaks from ch_mrep_macs_qc.collect{ it[1] }
    path annos from ch_mrep_macs_annotate.collect()
    path mrep_peak_annotation_header from ch_mrep_peak_annotation_header

    output:
    path '*.tsv' into ch_mrep_peak_qc_mqc
    path '*.{txt,pdf}'

    script:  // This script is bundled with the pipeline, in nf-core/atacseq/bin/
    suffix = 'mRp.clN'
    """
    plot_macs_qc.r \\
        -i ${peaks.join(',')} \\
        -s ${peaks.join(',').replaceAll(".${suffix}_peaks.${PEAK_TYPE}","")} \\
        -o ./ \\
        -p macs_peak.${suffix}

    plot_homer_annotatepeaks.r \\
        -i ${annos.join(',')} \\
        -s ${annos.join(',').replaceAll(".${suffix}_peaks.annotatePeaks.txt","")} \\
        -o ./ \\
        -p macs_annotatePeaks.${suffix}

    cat $mrep_peak_annotation_header macs_annotatePeaks.${suffix}.summary.txt > macs_annotatePeaks.${suffix}.summary_mqc.tsv
    """
}

/*
 * STEP 8.5: Consensus peaks across samples, create boolean filtering file, SAF file for featureCounts and UpSetR plot for intersection
 */
process MERGED_REP_CONSENSUS {
    label 'process_long'
    publishDir "${params.outdir}/bwa/mergedReplicate/macs/${PEAK_TYPE}/consensus", mode: params.publish_dir_mode,
        saveAs: { filename ->
                      if (filename.endsWith('.igv.txt')) null
                      else filename
                }

    when:
    !params.skip_merge_replicates && replicatesExist && params.macs_gsize && multipleGroups && !params.skip_consensus_peaks

    input:
    path peaks from ch_mrep_macs_consensus.collect{ it[1] }

    output:
    path '*.bed' into ch_mrep_macs_consensus_bed
    path '*.saf' into ch_mrep_macs_consensus_saf
    path '*.boolean.txt' into ch_mrep_macs_consensus_bool
    path '*igv.txt' into ch_mrep_macs_consensus_igv
    path '*.intersect.{txt,plot.pdf}'

    script: // scripts are bundled with the pipeline, in nf-core/atacseq/bin/
    suffix = 'mRp.clN'
    prefix = "consensus_peaks.${suffix}"
    mergecols = params.narrow_peak ? (2..10).join(',') : (2..9).join(',')
    collapsecols = params.narrow_peak ? (['collapse']*9).join(',') : (['collapse']*8).join(',')
    expandparam = params.narrow_peak ? '--is_narrow_peak' : ''
    """
    sort -T '.' -k1,1 -k2,2n ${peaks.collect{it.toString()}.sort().join(' ')} \\
        | mergeBed -c $mergecols -o $collapsecols > ${prefix}.txt

    macs2_merged_expand.py \\
        ${prefix}.txt \\
        ${peaks.collect{it.toString()}.sort().join(',').replaceAll("_peaks.${PEAK_TYPE}","")} \\
        ${prefix}.boolean.txt \\
        $expandparam

    awk -v FS='\t' -v OFS='\t' 'FNR > 1 { print \$1, \$2, \$3, \$4, "0", "+" }' ${prefix}.boolean.txt > ${prefix}.bed

    echo -e "GeneID\tChr\tStart\tEnd\tStrand" > ${prefix}.saf
    awk -v FS='\t' -v OFS='\t' 'FNR > 1 { print \$4, \$1, \$2, \$3,  "+" }' ${prefix}.boolean.txt >> ${prefix}.saf

    sed -i 's/.${suffix}//g' ${prefix}.boolean.intersect.txt
    plot_peak_intersect.r -i ${prefix}.boolean.intersect.txt -o ${prefix}.boolean.intersect.plot.pdf

    find * -type f -name "${prefix}.bed" -exec echo -e "bwa/mergedReplicate/macs/${PEAK_TYPE}/consensus/"{}"\\t0,0,0" \\; > ${prefix}.bed.igv.txt
    """
}

/*
 * STEP 8.6: Annotate consensus peaks with HOMER, and add annotation to boolean output file
 */
process MERGED_REP_CONSENSUS_ANNOTATE {
    label 'process_medium'
    publishDir "${params.outdir}/bwa/mergedReplicate/macs/${PEAK_TYPE}/consensus", mode: params.publish_dir_mode

    when:
    !params.skip_merge_replicates && replicatesExist && params.macs_gsize && multipleGroups && !params.skip_consensus_peaks && !params.skip_peak_annotation

    input:
    path bed from ch_mrep_macs_consensus_bed
    path bool from ch_mrep_macs_consensus_bool
    path fasta from ch_fasta
    path gtf from ch_gtf

    output:
    path '*.annotatePeaks.txt'

    script:
    prefix = 'consensus_peaks.mRp.clN'
    """
    annotatePeaks.pl \\
        $bed \\
        $fasta \\
        -gid \\
        -gtf $gtf \\
        -cpu $task.cpus \\
        > ${prefix}.annotatePeaks.txt

    cut -f2- ${prefix}.annotatePeaks.txt | awk 'NR==1; NR > 1 {print \$0 | "sort -T '.' -k1,1 -k2,2n"}' | cut -f6- > tmp.txt
    paste $bool tmp.txt > ${prefix}.boolean.annotatePeaks.txt
    """
}

/*
 * STEP 8.7: Count reads in consensus peaks with featureCounts
 */
process MERGED_REP_CONSENSUS_COUNTS {
    label 'process_medium'
    publishDir "${params.outdir}/bwa/mergedReplicate/macs/${PEAK_TYPE}/consensus", mode: params.publish_dir_mode

    when:
    !params.skip_merge_replicates && replicatesExist && params.macs_gsize && multipleGroups && !params.skip_consensus_peaks

    input:
    path bams from ch_mlib_name_bam_mrep_counts.collect{ it[1] }
    path saf from ch_mrep_macs_consensus_saf.collect()

    output:
    path '*featureCounts.txt' into ch_mrep_macs_consensus_counts
    path '*featureCounts.txt.summary' into ch_mrep_macs_consensus_counts_mqc

    script:
    prefix = 'consensus_peaks.mRp.clN'
    bam_files = bams.findAll { it.toString().endsWith('.bam') }.sort()
    pe_params = params.single_end ? '' : '-p --donotsort'
    """
    featureCounts \\
        -F SAF \\
        -O \\
        --fracOverlap 0.2 \\
        -T $task.cpus \\
        $pe_params \\
        -a $saf \\
        -o ${prefix}.featureCounts.txt \\
        ${bam_files.join(' ')}
    """
}

/*
 * STEP 8.8: Differential analysis with DESeq2
 */
process MERGED_REP_CONSENSUS_DESEQ2 {
    label 'process_medium'
    publishDir "${params.outdir}/bwa/mergedReplicate/macs/${PEAK_TYPE}/consensus/deseq2", mode: params.publish_dir_mode,
        saveAs: { filename ->
                      if (filename.endsWith('.igv.txt')) null
                      else filename
                }

    when:
    !params.skip_merge_replicates && replicatesExist && params.macs_gsize && multipleGroups && !params.skip_consensus_peaks && !params.skip_diff_analysis

    input:
    path counts from ch_mrep_macs_consensus_counts
    path mrep_deseq2_pca_header from ch_mrep_deseq2_pca_header
    path mrep_deseq2_clustering_header from ch_mrep_deseq2_clustering_header

    output:
    path '*.tsv' into ch_mrep_macs_consensus_deseq_mqc
    path '*igv.txt' into ch_mrep_macs_consensus_deseq_comp_igv
    path '*.{RData,results.txt,pdf,log}'
    path 'sizeFactors'
    path '*vs*/*.{pdf,txt}'
    path '*vs*/*.bed'

    script:
    prefix = 'consensus_peaks.mRp.clN'
    bam_ext = params.single_end ? '.mLb.clN.sorted.bam' : '.mLb.clN.bam'
    vst = params.deseq2_vst ? '--vst TRUE' : ''
    """
    featurecounts_deseq2.r \\
        --featurecount_file $counts \\
        --bam_suffix '$bam_ext' \\
        --outdir ./ \\
        --outprefix $prefix \\
        --outsuffix .mRp.clN \\
        --cores $task.cpus \\
        $vst

    cat $mrep_deseq2_pca_header ${prefix}.pca.vals.txt > ${prefix}.pca.vals_mqc.tsv
    cat $mrep_deseq2_clustering_header ${prefix}.sample.dists.txt > ${prefix}.sample.dists_mqc.tsv

    find * -type f -name "*.FDR0.05.results.bed" -exec echo -e "bwa/mergedReplicate/macs/${PEAK_TYPE}/consensus/deseq2/"{}"\\t255,0,0" \\; > ${prefix}.igv.txt
    """
}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
/* --                                                                     -- */
/* --                             IGV                                     -- */
/* --                                                                     -- */
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

/*
 * STEP 9: Create IGV session file
 */
process IGV {
    publishDir "${params.outdir}/igv/${PEAK_TYPE}", mode: params.publish_dir_mode

    when:
    !params.skip_igv

    input:
    path fasta from ch_fasta

    path bigwigs from ch_mlib_bigwig_igv.collect().ifEmpty([])
    path peaks from ch_mlib_macs_igv.collect().ifEmpty([])
    path consensus_peaks from ch_mlib_macs_consensus_igv.collect().ifEmpty([])
    path differential_peaks from ch_mlib_macs_consensus_deseq_comp_igv.collect().ifEmpty([])

    path rbigwigs from ch_mrep_bigwig_igv.collect().ifEmpty([])
    path rpeaks from ch_mrep_macs_igv.collect().ifEmpty([])
    path rconsensus_peaks from ch_mrep_macs_consensus_igv.collect().ifEmpty([])
    path rdifferential_peaks from ch_mrep_macs_consensus_deseq_comp_igv.collect().ifEmpty([])

    output:
    path '*.{txt,xml}'

    script: // scripts are bundled with the pipeline, in nf-core/atacseq/bin/
    """
    cat *.txt > igv_files.txt
    igv_files_to_session.py igv_session.xml igv_files.txt ../../genome/${fasta.getName()} --path_prefix '../../'
    """
}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
/* --                                                                     -- */
/* --                          MULTIQC                                    -- */
/* --                                                                     -- */
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

/*
 * Parse software version numbers
 */
process get_software_versions {
    publishDir "${params.outdir}/pipeline_info", mode: params.publish_dir_mode,
        saveAs: { filename ->
                      if (filename.indexOf('.csv') > 0) filename
                      else null
                }

    output:
    path 'software_versions_mqc.yaml' into ch_software_versions_mqc
    path 'software_versions.csv'

    script:
    """
    echo $workflow.manifest.version > v_pipeline.txt
    echo $workflow.nextflow.version > v_nextflow.txt
    fastqc --version > v_fastqc.txt
    trim_galore --version > v_trim_galore.txt
    echo \$(bwa 2>&1) > v_bwa.txt
    samtools --version > v_samtools.txt
    bedtools --version > v_bedtools.txt
    echo \$(bamtools --version 2>&1) > v_bamtools.txt
    echo \$(plotFingerprint --version 2>&1) > v_deeptools.txt || true
    picard MarkDuplicates --version &> v_picard.txt  || true
    echo \$(R --version 2>&1) > v_R.txt
    python -c "import pysam; print(pysam.__version__)" > v_pysam.txt
    echo \$(macs2 --version 2>&1) > v_macs2.txt
    touch v_homer.txt
    echo \$(ataqv --version 2>&1) > v_ataqv.txt
    echo \$(featureCounts -v 2>&1) > v_featurecounts.txt
    preseq &> v_preseq.txt
    multiqc --version > v_multiqc.txt
    scrape_software_versions.py &> software_versions_mqc.yaml
    """
}

Channel.from(summary.collect{ [it.key, it.value] })
    .map { k,v -> "<dt>$k</dt><dd><samp>${v ?: '<span style=\"color:#999999;\">N/A</a>'}</samp></dd>" }
    .reduce { a, b -> return [a, b].join("\n            ") }
    .map { x -> """
    id: 'nf-core-atacseq-summary'
    description: " - this information is collected when the pipeline is started."
    section_name: 'nf-core/atacseq Workflow Summary'
    section_href: 'https://github.com/nf-core/atacseq'
    plot_type: 'html'
    data: |
        <dl class=\"dl-horizontal\">
            $x
        </dl>
    """.stripIndent() }
    .set { ch_workflow_summary }

/*
 * STEP 10: MultiQC
 */
process MULTIQC {
    publishDir "${params.outdir}/multiqc/${PEAK_TYPE}", mode: params.publish_dir_mode

    when:
    !params.skip_multiqc

    input:
    path (multiqc_config) from ch_multiqc_config
    path (mqc_custom_config) from ch_multiqc_custom_config.collect().ifEmpty([])

    path ('software_versions/*') from ch_software_versions_mqc.collect()
    path workflow_summary from ch_workflow_summary.collectFile(name: 'workflow_summary_mqc.yaml')

    path ('fastqc/*') from ch_fastqc_reports_mqc.collect().ifEmpty([])
    path ('trimgalore/*') from ch_trimgalore_results_mqc.collect().ifEmpty([])
    path ('trimgalore/fastqc/*') from ch_trimgalore_fastqc_reports_mqc.collect().ifEmpty([])

    path ('alignment/library/*') from ch_sort_bam_flagstat_mqc.collect()

    path ('alignment/mergedLibrary/*') from ch_mlib_bam_stats_mqc.collect()
    path ('alignment/mergedLibrary/*') from ch_mlib_rm_orphan_flagstat_mqc.collect{it[1]}
    path ('alignment/mergedLibrary/*') from ch_mlib_rm_orphan_stats_mqc.collect()
    path ('alignment/mergedLibrary/picard_metrics/*') from ch_mlib_bam_metrics_mqc.collect()
    path ('alignment/mergedLibrary/picard_metrics/*') from ch_mlib_collectmetrics_mqc.collect()
    path ('macs/mergedLibrary/*') from ch_mlib_macs_mqc.collect().ifEmpty([])
    path ('macs/mergedLibrary/*') from ch_mlib_peak_qc_mqc.collect().ifEmpty([])
    path ('macs/mergedLibrary/consensus/*') from ch_mlib_macs_consensus_counts_mqc.collect().ifEmpty([])
    path ('macs/mergedLibrary/consensus/*') from ch_mlib_macs_consensus_deseq_mqc.collect().ifEmpty([])
    path ('preseq/*') from ch_mlib_preseq_mqc.collect().ifEmpty([])
    path ('deeptools/*') from ch_mlib_plotprofile_mqc.collect().ifEmpty([])
    path ('deeptools/*') from ch_mlib_plotfingerprint_mqc.collect().ifEmpty([])

    path ('alignment/mergedReplicate/*') from ch_mrep_bam_flagstat_mqc.collect{it[1]}.ifEmpty([])
    path ('alignment/mergedReplicate/*') from ch_mrep_bam_stats_mqc.collect().ifEmpty([])
    path ('alignment/mergedReplicate/*') from ch_mrep_bam_metrics_mqc.collect().ifEmpty([])
    path ('macs/mergedReplicate/*') from ch_mrep_macs_mqc.collect().ifEmpty([])
    path ('macs/mergedReplicate/*') from ch_mrep_peak_qc_mqc.collect().ifEmpty([])
    path ('macs/mergedReplicate/consensus/*') from ch_mrep_macs_consensus_counts_mqc.collect().ifEmpty([])
    path ('macs/mergedReplicate/consensus/*') from ch_mrep_macs_consensus_deseq_mqc.collect().ifEmpty([])

    output:
    path '*multiqc_report.html' into ch_multiqc_report
    path '*_data'

    script:
    rtitle = custom_runName ? "--title \"$custom_runName\"" : ''
    rfilename = custom_runName ? "--filename " + custom_runName.replaceAll('\\W','_').replaceAll('_+','_') + "_multiqc_report" : ''
    custom_config_file = params.multiqc_config ? "--config $mqc_custom_config" : ''
    """
    multiqc . -f $rtitle $rfilename $custom_config_file
    """
}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
/* --                                                                     -- */
/* --                       REPORTS/DOCUMENTATION                         -- */
/* --                                                                     -- */
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

/*
 * STEP 11: Output description HTML
 */
process output_documentation {
    publishDir "${params.outdir}/pipeline_info", mode: params.publish_dir_mode

    input:
    path output_docs from ch_output_docs
    path images from ch_output_docs_images

    output:
    path 'results_description.html'

    script:
    """
    markdown_to_html.py $output_docs -o results_description.html
    """
}

/*
 * Completion e-mail notification
 */
workflow.onComplete {

    // Set up the e-mail variables
    def subject = "[nf-core/atacseq] Successful: $workflow.runName"
    if (!workflow.success) {
        subject = "[nf-core/atacseq] 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

    // On success try attach the multiqc report
    def mqc_report = null
    try {
        if (workflow.success) {
            mqc_report = ch_multiqc_report.getVal()
            if (mqc_report.getClass() == ArrayList) {
                log.warn "[nf-core/atacseq] Found multiple reports from process 'multiqc', will use only one"
                mqc_report = mqc_report[0]
            }
        }
    } catch (all) {
        log.warn "[nf-core/atacseq] Could not attach MultiQC report to summary email"
    }

    // Check if we are only sending emails on failure
    email_address = params.email
    if (!params.email && params.email_on_fail && !workflow.success) {
        email_address = params.email_on_fail
    }

    // 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: email_address, subject: subject, email_txt: email_txt, email_html: email_html, baseDir: "$baseDir", mqcFile: mqc_report, mqcMaxSize: params.max_multiqc_email_size.toBytes() ]
    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 (email_address) {
        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/atacseq] Sent summary e-mail to $email_address (sendmail)"
        } catch (all) {
            // Catch failures and try with plaintext
            [ 'mail', '-s', subject, email_address ].execute() << email_txt
            log.info "[nf-core/atacseq] Sent summary e-mail to $email_address (mail)"
        }
    }

    // Write summary e-mail HTML to a file
    def output_d = new File("${params.outdir}/pipeline_info/")
    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 }

    c_green = params.monochrome_logs ? '' : "\033[0;32m";
    c_purple = params.monochrome_logs ? '' : "\033[0;35m";
    c_red = params.monochrome_logs ? '' : "\033[0;31m";
    c_reset = params.monochrome_logs ? '' : "\033[0m";

    if (workflow.stats.ignoredCount > 0 && workflow.success) {
        log.info "-${c_purple}Warning, pipeline completed, but with errored process(es) ${c_reset}-"
        log.info "-${c_red}Number of ignored errored process(es) : ${workflow.stats.ignoredCount} ${c_reset}-"
        log.info "-${c_green}Number of successfully ran process(es) : ${workflow.stats.succeedCount} ${c_reset}-"
    }
    if (workflow.success) {
        log.info "-${c_purple}[nf-core/atacseq]${c_green} Pipeline completed successfully${c_reset}-"
    } else {
        checkHostname()
        log.info "-${c_purple}[nf-core/atacseq]${c_red} Pipeline completed with errors${c_reset}-"
    }

}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
/* --                                                                     -- */
/* --                       NF-CORE HEADER                                -- */
/* --                                                                     -- */
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

def nfcoreHeader() {
    // Log colors ANSI codes
    c_black = params.monochrome_logs ? '' : "\033[0;30m";
    c_blue = params.monochrome_logs ? '' : "\033[0;34m";
    c_cyan = params.monochrome_logs ? '' : "\033[0;36m";
    c_dim = params.monochrome_logs ? '' : "\033[2m";
    c_green = params.monochrome_logs ? '' : "\033[0;32m";
    c_purple = params.monochrome_logs ? '' : "\033[0;35m";
    c_reset = params.monochrome_logs ? '' : "\033[0m";
    c_white = params.monochrome_logs ? '' : "\033[0;37m";
    c_yellow = params.monochrome_logs ? '' : "\033[0;33m";

    return """    -${c_dim}--------------------------------------------------${c_reset}-
                                            ${c_green},--.${c_black}/${c_green},-.${c_reset}
    ${c_blue}        ___     __   __   __   ___     ${c_green}/,-._.--~\'${c_reset}
    ${c_blue}  |\\ | |__  __ /  ` /  \\ |__) |__         ${c_yellow}}  {${c_reset}
    ${c_blue}  | \\| |       \\__, \\__/ |  \\ |___     ${c_green}\\`-._,-`-,${c_reset}
                                            ${c_green}`._,._,\'${c_reset}
    ${c_purple}  nf-core/atacseq v${workflow.manifest.version}${c_reset}
    -${c_dim}--------------------------------------------------${c_reset}-
    """.stripIndent()
}

def checkHostname() {
    def c_reset = params.monochrome_logs ? '' : "\033[0m"
    def c_white = params.monochrome_logs ? '' : "\033[0;37m"
    def c_red = params.monochrome_logs ? '' : "\033[1;91m"
    def c_yellow_bold = params.monochrome_logs ? '' : "\033[1;93m"
    if (params.hostnames) {
        def hostname = "hostname".execute().text.trim()
        params.hostnames.each { prof, hnames ->
            hnames.each { hname ->
                if (hostname.contains(hname) && !workflow.profile.contains(prof)) {
                    log.error "====================================================\n" +
                            "  ${c_red}WARNING!${c_reset} You are running with `-profile $workflow.profile`\n" +
                            "  but your machine hostname is ${c_white}'$hostname'${c_reset}\n" +
                            "  ${c_yellow_bold}It's highly recommended that you use `-profile $prof${c_reset}`\n" +
                            "============================================================"
                }
            }
        }
    }
}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
/* --                                                                     -- */
/* --                        END OF PIPELINE                              -- */
/* --                                                                     -- */
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////