Fragment length filtering not applied when shifting merged_replicate BAMs (--shift_reads)

[fbnrst]

Description of the bug

When running the pipeline with --shift_reads, the fragment-length filtering options (--minFragmentLength / --maxFragmentLength) are applied to DEEPTOOLS_ALIGNMENTSIEVE for the merged_library BAMs, but not for the merged_replicate BAMs. This causes the merged_replicate shifted BAMs (and everything downstream of them, e.g. MERGED_REPLICATE_CONSENSUS_PEAKS/SUBREAD_FEATURECOUNTS, bigWigs, etc.) to look substantially different from the merged_library BAMs and from the underlying individual libraries.

In conf/modules.config, the MERGED_LIBRARY_BAM_SHIFT_READS:DEEPTOOLS_ALIGNMENTSIEVE process config is:

withName: '.*:MERGED_LIBRARY_BAM_SHIFT_READS:DEEPTOOLS_ALIGNMENTSIEVE' {
    ext.args   = [
        '--ATACshift',
        "--minFragmentLength ${params.minFragmentLength}",
        "--maxFragmentLength ${params.maxFragmentLength}"
    ].join(' ').trim()
    ...
}

while the corresponding MERGED_REPLICATE_BAM_SHIFT_READS:DEEPTOOLS_ALIGNMENTSIEVE only sets:

withName: '.*:MERGED_REPLICATE_BAM_SHIFT_READS:DEEPTOOLS_ALIGNMENTSIEVE' {
    ext.args   = '--ATACshift'
    ...
}

The --minFragmentLength/--maxFragmentLength arguments are missing here, so the merged_replicate shifting step does not filter reads by fragment length the way the merged_library step does.

This was identified and discussed on the nf-core Slack: https://nfcore.slack.com/archives/CE5EL6326/p1749127852180939. In that thread, MERGED_REPLICATE_BAM_SHIFT_READS was confirmed to run (so the shift step itself is not skipped), but the resulting peaks/BAMs for merged_replicate looked visibly different from merged_library and from the individual libraries in IGV, whereas without --shift_reads the merged_replicate peaks matched the individual libraries closely. Merging the raw FASTQs of the replicates upfront (bypassing the merged_replicate BAM-level path) reproduced a peak shape consistent with the individual libraries, pointing at the merged_replicate shifting step (rather than the shifting logic itself) as the source of the discrepancy.

Command used and terminal output

$ nextflow run nf-core/atacseq -r dev --shift_reads ...

[15/b580c0] process > NFCORE_ATACSEQ:ATACSEQ:MERGED_REPLICATE_BAM_SHIFT_READS:DEEPTOOLS_ALIGNMENTSIEVE (4d38_f_Prom1)  [100%] 1 of 1 ✔
[42/bb36f5] process > NFCORE_ATACSEQ:ATACSEQ:MERGED_REPLICATE_BAM_SHIFT_READS:SAMTOOLS_SORT (4d38_f_Prom1)             [100%] 1 of 1 ✔
[40/195747] process > NFCORE_ATACSEQ:ATACSEQ:MERGED_REPLICATE_BAM_SHIFT_READS:SAMTOOLS_INDEX (4d38_f_Prom1)            [100%] 1 of 1 ✔
[ae/b0a507] process > NFCORE_ATACSEQ:ATACSEQ:MERGED_REPLICATE_BAM_SHIFT_READS:SAMTOOLS_FLAGSTAT (4d38_f_Prom1)         [100%] 1 of 1 ✔
[a7/e6dbe9] process > NFCORE_ATACSEQ:ATACSEQ:MERGED_REPLICATE_BAM_TO_BIGWIG:BEDTOOLS_GENOMECOV (4d38_f_Prom1)          [100%] 1 of 1 ✔

Relevant files

Peak shapes for the same region, viewed in IGV, differ substantially between the merged_library (blue) and merged_replicate (red) tracks when --shift_reads is used; without --shift_reads the merged_replicate peak matches the individual libraries much more closely. Screenshots and full discussion are in the linked Slack thread.

System information

• Version of nf-core/atacseq: dev, commit 5960178 (bug still present on dev @ 4177ea9)
• Nextflow version: n/a (not the source of the bug)
• Container engine: n/a
• OS: n/a
• Executor: n/a
• Hardware: n/a