Fixed nf-core#164 - Introduced 4/5bp shift as it is common for ATAC-s…

…eq data. Fixed nf-core#168 - Always write out genome fa and fai so IGV session file can be opened.

conf/modules.config

@@ -453,6 +453,37 @@ process {
         ]
     }
 
+    withName: '.*:MERGED_LIBRARY_BAM_SHIFT_READS:DEEPTOOLS_ALIGNMENTSIEVE' {
+        ext.args   = '--ATACshift'
+        ext.prefix = { "${meta.id}.mLb.clN.shifted" }
+        publishDir = [
+            path: { "${params.outdir}/${params.aligner}/merged_library/shifted_reads" },
+            mode: params.publish_dir_mode,
+            pattern: '*.bam',
+            enabled: params.save_align_intermeds
+        ]
+    }
+
+    withName: '.*:MERGED_LIBRARY_BAM_SHIFT_READS:SAMTOOLS_SORT' {
+        ext.prefix = { "${meta.id}.mLb.clN.shifted.sorted" }
+        publishDir = [
+            path: { "${params.outdir}/${params.aligner}/merged_library/shifted_reads" },
+            mode: params.publish_dir_mode,
+            pattern: '*.bam',
+            enabled: params.shift_reads
+        ]
+    }
+
+    withName: '.*:MERGED_LIBRARY_BAM_SHIFT_READS:SAMTOOLS_INDEX' {
+        ext.prefix = { "${meta.id}.mLb.clN.shifted.sorted" }
+        publishDir = [
+            path: { "${params.outdir}/${params.aligner}/merged_library/shifted_reads" },
+            mode: params.publish_dir_mode,
+            pattern: '*.bai',
+            enabled: params.shift_reads
+        ]
+    }
+
     withName: '.*:MERGED_LIBRARY_BAM_TO_BIGWIG:BEDTOOLS_GENOMECOV' {
         ext.args   = {
             [
@@ -782,6 +813,37 @@ if (!params.skip_merge_replicates) {
             ]
         }
 
+        withName: '.*:MERGED_REPLICATE_BAM_SHIFT_READS:DEEPTOOLS_ALIGNMENTSIEVE' {
+            ext.args   = '--ATACshift'
+            ext.prefix = { "${meta.id}.mRp.clN.shifted" }
+            publishDir = [
+                path: { "${params.outdir}/${params.aligner}/merged_replicate" },
+                mode: params.publish_dir_mode,
+                pattern: '*.bam',
+                enabled: params.save_align_intermeds
+            ]
+        }
+
+        withName: '.*:MERGED_REPLICATE_BAM_SHIFT_READS:SAMTOOLS_SORT' {
+            ext.prefix = { "${meta.id}.mRp.clN.shifted.sorted" }
+            publishDir = [
+                path: { "${params.outdir}/${params.aligner}/merged_replicate/shifted_reads" },
+                mode: params.publish_dir_mode,
+                pattern: '*.bam',
+                enabled: params.shift_reads
+            ]
+        }
+
+        withName: '.*:MERGED_REPLICATE_BAM_SHIFT_READS:SAMTOOLS_INDEX' {
+            ext.prefix = { "${meta.id}.mRp.clN.shifted.sorted" }
+            publishDir = [
+                path: { "${params.outdir}/${params.aligner}/merged_replicate/shifted_reads" },
+                mode: params.publish_dir_mode,
+                pattern: '*.bai',
+                enabled: params.shift_reads
+            ]
+        }
+
         withName: '.*:MERGED_REPLICATE_BAM_TO_BIGWIG:BEDTOOLS_GENOMECOV' {
             ext.args   = {
                 [
@@ -956,8 +1018,7 @@ if (!params.skip_igv) {
                 [
                     path: { "${params.outdir}/genome" },
                     mode: params.publish_dir_mode,
-                    pattern: '*.{fa,fasta}',
-                    enabled: params.save_reference
+                    pattern: '*.{fa,fasta,fai}'
                 ]
             ]
         }

modules/local/deeptools_alignmentsieve.nf

@@ -0,0 +1,36 @@
+process DEEPTOOLS_ALIGNMENTSIEVE {
+    tag "$meta.id"
+    label 'process_medium'
+
+    conda "bioconda::deeptools=3.5.1"
+    container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ?
+        'https://depot.galaxyproject.org/singularity/deeptools:3.5.1--py_0' :
+        'biocontainers/deeptools:3.5.1--py_0' }"
+
+    input:
+    tuple val(meta), path(bam), path(bai)
+
+    output:
+    tuple val(meta), path("*.bam"), emit: bam
+    path  "versions.yml", emit: versions
+
+    when:
+    task.ext.when == null || task.ext.when
+
+    script:
+    def args = task.ext.args ?: ''
+    def prefix = task.ext.prefix ?: "${meta.id}"
+
+    """
+    alignmentSieve \\
+        $args \\
+        -b $bam \\
+        -o ${prefix}.bam \\
+        --numberOfProcessors $task.cpus
+
+    cat <<-END_VERSIONS > versions.yml
+    "${task.process}":
+        deeptools: \$(alignmentSieve --version | sed -e "s/alignmentSieve //g")
+    END_VERSIONS
+    """
+}

modules/local/igv.nf

@@ -7,6 +7,7 @@ process IGV {
 
     input:
     path fasta
+    path fai
     path ("${bigwig_library_publish_dir}/*")
     path ("${peak_library_publish_dir}/*")
     path ("${consensus_library_publish_dir}/*")
@@ -25,6 +26,7 @@ process IGV {
     path "*files.txt"  , emit: txt
     path "*.xml"       , emit: xml
     path fasta         , emit: fasta
+    path fai           , emit: fai
     path "versions.yml", emit: versions
 
     when:

nextflow.config

@@ -39,6 +39,7 @@ params {
     skip_merge_replicates      = false
     save_align_intermeds       = false
     save_unaligned             = false
+    shift_reads                = true
 
     // Options: Peaks
     narrow_peak                = false

nextflow_schema.json

@@ -289,6 +289,13 @@
                     "hidden": true,
                     "description": "BAMTools JSON file with custom filters for single-end data.",
                     "fa_icon": "fas fa-cog"
+                },
+                "shift_reads": {
+                    "type": "boolean",
+                    "fa_icon": "fas fa-chart-area",
+                    "default": true,
+                    "help_text": "Shift aligned reads as commonly done for ATACseq, +4bp for reads on the + strand, -5 bp for reads on the - strand. This can only be applied if all samples are paired-end.",
+                    "description": "Shift aligned reads (+4bp and -5bp)."
                 }
             }
         },

subworkflows/local/bam_shift_reads.nf

@@ -0,0 +1,40 @@
+include { SAMTOOLS_SORT            } from '../../modules/nf-core/samtools/sort/main'
+include { SAMTOOLS_INDEX           } from '../../modules/nf-core/samtools/index/main'
+include { DEEPTOOLS_ALIGNMENTSIEVE } from '../../modules/local/deeptools_alignmentsieve'
+
+workflow BAM_SHIFT_READS {
+    take:
+    ch_bam_bai                   // channel: [ val(meta), [ bam ], [bai] ]
+
+    main:
+    ch_versions = Channel.empty()
+
+    //
+    // Shift reads
+    //
+    DEEPTOOLS_ALIGNMENTSIEVE (
+        ch_bam_bai
+    )
+    ch_versions = ch_versions.mix(DEEPTOOLS_ALIGNMENTSIEVE.out.versions)
+
+    //
+    // Sort reads
+    //
+    SAMTOOLS_SORT (
+        DEEPTOOLS_ALIGNMENTSIEVE.out.bam
+    )
+    ch_versions = ch_versions.mix(SAMTOOLS_SORT.out.versions)
+
+    //
+    // Index reads
+    //
+    SAMTOOLS_INDEX (
+        SAMTOOLS_SORT.out.bam
+    )
+    ch_versions = ch_versions.mix(SAMTOOLS_INDEX.out.versions)
+
+    emit:
+    bam = SAMTOOLS_SORT.out.bam                     // channel: [ val(meta), [ bam ] ]
+    bai = SAMTOOLS_INDEX.out.bai                    // channel: [ val(meta), [ bai ] ]
+    versions = ch_versions                          // channel: [ versions.yml ]
+}

workflows/atacseq.nf

@@ -68,9 +68,11 @@ include { MULTIQC } from '../modules/local/multiqc'
 //
 // SUBWORKFLOW: Consisting of a mix of local and nf-core/modules
 //
-include { INPUT_CHECK    } from '../subworkflows/local/input_check'
-include { PREPARE_GENOME } from '../subworkflows/local/prepare_genome'
-include { ALIGN_STAR     } from '../subworkflows/local/align_star'
+include { INPUT_CHECK     } from '../subworkflows/local/input_check'
+include { PREPARE_GENOME  } from '../subworkflows/local/prepare_genome'
+include { ALIGN_STAR      } from '../subworkflows/local/align_star'
+include { BAM_SHIFT_READS as MERGED_LIBRARY_BAM_SHIFT_READS } from '../subworkflows/local/bam_shift_reads'
+include { BAM_SHIFT_READS as MERGED_REPLICATE_BAM_SHIFT_READS } from '../subworkflows/local/bam_shift_reads'
 include { BIGWIG_PLOT_DEEPTOOLS as MERGED_LIBRARY_BIGWIG_PLOT_DEEPTOOLS       } from '../subworkflows/local/bigwig_plot_deeptools'
 include { BAM_FILTER_BAMTOOLS as MERGED_LIBRARY_FILTER_BAM                    } from '../subworkflows/local/bam_filter_bamtools'
 include { BAM_BEDGRAPH_BIGWIG_BEDTOOLS_UCSC as MERGED_LIBRARY_BAM_TO_BIGWIG   } from '../subworkflows/local/bam_bedgraph_bigwig_bedtools_ucsc'
@@ -96,6 +98,7 @@ include { PRESEQ_LCEXTRAP as MERGED_LIBRARY_PRESEQ_LCEXTRAP
 include { DEEPTOOLS_PLOTFINGERPRINT as MERGED_LIBRARY_DEEPTOOLS_PLOTFINGERPRINT         } from '../modules/nf-core/deeptools/plotfingerprint/main'
 include { ATAQV_ATAQV as MERGED_LIBRARY_ATAQV_ATAQV                                     } from '../modules/nf-core/ataqv/ataqv/main'
 include { ATAQV_MKARV as MERGED_LIBRARY_ATAQV_MKARV                                     } from '../modules/nf-core/ataqv/mkarv/main'
+include { SAMTOOLS_INDEX } from '../modules/nf-core/samtools/index/main'
 
 include { PICARD_MERGESAMFILES as PICARD_MERGESAMFILES_LIBRARY   } from '../modules/nf-core/picard/mergesamfiles/main'
 include { PICARD_MERGESAMFILES as PICARD_MERGESAMFILES_REPLICATE } from '../modules/nf-core/picard/mergesamfiles/main'
@@ -141,6 +144,24 @@ workflow ATACSEQ {
     )
     ch_versions = ch_versions.mix(INPUT_CHECK.out.versions)
 
+    //
+    // Check if reads are all paired-end if 'shift_reads' parameter is set
+    //
+    if (params.shift_reads) {
+        INPUT_CHECK
+            .out
+            .reads
+            .filter { meta, reads -> meta.single_end }
+            .collect()
+            .map {
+                it ->
+                    def count = it.size()
+                    if (count > 0) {
+                        exit 1, 'The parameter --shift_reads can only be applied if all samples are paired-end.'
+                    }
+            }
+    }
+
     //
     // SUBWORKFLOW: Read QC and trim adapters
     //
@@ -242,8 +263,8 @@ workflow ATACSEQ {
             [],
             []
         )
-        ch_genome_bam        = FASTQ_ALIGN_CHROMAP.out.bam
         ch_genome_bam_index  = FASTQ_ALIGN_CHROMAP.out.bai
+        ch_genome_bam        = FASTQ_ALIGN_CHROMAP.out.bam
         ch_samtools_stats    = FASTQ_ALIGN_CHROMAP.out.stats
         ch_samtools_flagstat = FASTQ_ALIGN_CHROMAP.out.flagstat
         ch_samtools_idxstats = FASTQ_ALIGN_CHROMAP.out.idxstats
@@ -342,11 +363,27 @@ workflow ATACSEQ {
         ch_versions = ch_versions.mix(MERGED_LIBRARY_PICARD_COLLECTMULTIPLEMETRICS.out.versions.first())
     }
 
+    //
+    // SUBWORKFLOW: Shift paired-end reads
+    //
+    ch_merged_library_filter_bam = MERGED_LIBRARY_FILTER_BAM.out.bam
+    ch_merged_library_filter_bai = MERGED_LIBRARY_FILTER_BAM.out.bai
+
+    if (params.shift_reads && params.aligner != 'chromap' ) {
+        MERGED_LIBRARY_BAM_SHIFT_READS (
+            ch_merged_library_filter_bam.join(ch_merged_library_filter_bai, by: [0]),
+        )
+        ch_versions = ch_versions.mix(MERGED_LIBRARY_BAM_SHIFT_READS.out.versions)
+
+        ch_merged_library_filter_bam = MERGED_LIBRARY_BAM_SHIFT_READS.out.bam
+        ch_merged_library_filter_bai = MERGED_LIBRARY_BAM_SHIFT_READS.out.bai
+    }
+
     //
     // SUBWORKFLOW: Normalised bigWig coverage tracks
     //
     MERGED_LIBRARY_BAM_TO_BIGWIG (
-        MERGED_LIBRARY_FILTER_BAM.out.bam.join(MERGED_LIBRARY_FILTER_BAM.out.flagstat, by: [0]),
+        ch_merged_library_filter_bam.join(MERGED_LIBRARY_FILTER_BAM.out.flagstat, by: [0]),
         PREPARE_GENOME.out.chrom_sizes
     )
     ch_versions = ch_versions.mix(MERGED_LIBRARY_BAM_TO_BIGWIG.out.versions)
@@ -366,10 +403,8 @@ workflow ATACSEQ {
     }
 
     // Create channels: [ meta, [bam], [bai] ]
-    MERGED_LIBRARY_FILTER_BAM
-        .out
-        .bam
-        .join(MERGED_LIBRARY_FILTER_BAM.out.bai, by: [0])
+    ch_merged_library_filter_bam
+        .join(ch_merged_library_filter_bai, by: [0])
         .set { ch_bam_bai }
 
     //
@@ -523,24 +558,37 @@ workflow ATACSEQ {
         ch_markduplicates_replicate_metrics  = MERGED_REPLICATE_MARKDUPLICATES_PICARD.out.metrics
         ch_versions = ch_versions.mix(MERGED_REPLICATE_MARKDUPLICATES_PICARD.out.versions)
 
+        //
+        // SUBWORKFLOW: Shift paired-end reads
+        // Shift again, as ch_merged_library_replicate_bam is generated out of unshifted reads
+        //
+        ch_merged_replicate_markduplicate_bam = MERGED_REPLICATE_MARKDUPLICATES_PICARD.out.bam
+        ch_merged_replicate_markduplicate_bai = MERGED_REPLICATE_MARKDUPLICATES_PICARD.out.bai
+
+        if (params.shift_reads && params.aligner != 'chromap' ) {
+            MERGED_REPLICATE_BAM_SHIFT_READS (
+                ch_merged_replicate_markduplicate_bam.join(ch_merged_replicate_markduplicate_bai, by: [0]),
+            )
+            ch_versions = ch_versions.mix(MERGED_REPLICATE_BAM_SHIFT_READS.out.versions)
+
+            ch_merged_replicate_markduplicate_bam = MERGED_REPLICATE_BAM_SHIFT_READS.out.bam
+            ch_merged_replicate_markduplicate_bai = MERGED_REPLICATE_BAM_SHIFT_READS.out.bai
+        }
+
         // SUBWORKFLOW: Normalised bigWig coverage tracks
         //
         MERGED_REPLICATE_BAM_TO_BIGWIG (
-            MERGED_REPLICATE_MARKDUPLICATES_PICARD.out.bam.join(MERGED_REPLICATE_MARKDUPLICATES_PICARD.out.flagstat, by: [0]),
+            ch_merged_replicate_markduplicate_bam.join(MERGED_REPLICATE_MARKDUPLICATES_PICARD.out.flagstat, by: [0]),
             PREPARE_GENOME.out.chrom_sizes
         )
         ch_ucsc_bedgraphtobigwig_replicate_bigwig = MERGED_REPLICATE_BAM_TO_BIGWIG.out.bigwig
         ch_versions = ch_versions.mix(MERGED_REPLICATE_BAM_TO_BIGWIG.out.versions)
 
         // Create channels: [ meta, bam, ([] for control_bam) ]
-        MERGED_REPLICATE_MARKDUPLICATES_PICARD
-            .out
-            .bam
-            .map {
-                meta, bam ->
-                    [ meta , bam, [] ]
-            }
-            .set { ch_bam_replicate }
+        // Create channels: [ meta, [bam], [bai] ]
+        ch_merged_replicate_markduplicate_bam
+        .join(ch_merged_replicate_markduplicate_bai, by: [0])
+        .set { ch_bam_replicate }
 
         //
         // SUBWORKFLOW: Call peaks with MACS2, annotate with HOMER and perform downstream QC
@@ -593,6 +641,7 @@ workflow ATACSEQ {
     if (!params.skip_igv) {
         IGV (
             PREPARE_GENOME.out.fasta,
+            PREPARE_GENOME.out.fai,
             MERGED_LIBRARY_BAM_TO_BIGWIG.out.bigwig.collect{it[1]}.ifEmpty([]),
             MERGED_LIBRARY_CALL_ANNOTATE_PEAKS.out.peaks.collect{it[1]}.ifEmpty([]),
             ch_macs2_consensus_library_bed.collect{it[1]}.ifEmpty([]),