count mixed SE/PE libraries per-batch instead of pinning subread

.nf-core.yml

@@ -1,6 +1,10 @@
 lint:
   files_unchanged:
     - .gitignore
+  merge_markers:
+    # False positive: SAM base-quality strings contain runs of '<' that trip the
+    # merge-marker check. Snapshot shipped by nf-core/modules bowtie2/align.
+    - modules/nf-core/bowtie2/align/tests/main.nf.test.snap
   nextflow_config:
     - config_defaults:
         - params.bamtools_filter_se_config

CHANGELOG.md

@@ -7,6 +7,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Enhancements & fixes
 
+- [[#452]](https://github.com/nf-core/atacseq/pull/452) - Count consensus peaks separately for single-end and paired-end libraries and merge the per-batch matrices, so mixed cohorts work under `subread` 2.1.1 and the 2.0.1 pin added in [[#448]](https://github.com/nf-core/atacseq/pull/448) can be removed.
+- [[#452]](https://github.com/nf-core/atacseq/pull/452) - Publish the per-library-type consensus count matrices (`*.SE.featureCounts.tsv` / `*.PE.featureCounts.tsv`) alongside the merged matrix.
+- [[#448]](https://github.com/nf-core/atacseq/pull/448) - Update nf-core modules and subworkflows to their latest versions and reconcile the resulting call-signature changes.
+- [[#448]](https://github.com/nf-core/atacseq/pull/448) - Add read groups to Chromap alignments, which Picard 3.4.0 `MarkDuplicates` now requires.
+- [[#448]](https://github.com/nf-core/atacseq/pull/448) - Quote `--seq_center` when building read-group arguments, so a sequencing-centre name containing whitespace no longer breaks Chromap (Picard) or Bowtie2 alignment.
+- [[#448]](https://github.com/nf-core/atacseq/pull/448) - Pin `subread` to 2.0.1 so consensus peak counting keeps supporting mixed single-end/paired-end cohorts in a single `featureCounts` invocation.
+- [[#448]](https://github.com/nf-core/atacseq/pull/448) - Sort the consensus `featureCounts` BAM inputs so the count matrix column order is deterministic.
 - [[#446]](https://github.com/nf-core/atacseq/pull/446) - Make pipeline code compliant with strict Nextflow v2 syntax parser, with no behaviour change.
 - [[#407]](https://github.com/nf-core/atacseq/pull/407) to add filtering reads according fragment size to help to focus on NFR, MNR, DNR, TNR
 - [[#164]](https://github.com/nf-core/atacseq/issues/164) and partly [[#91]](https://github.com/nf-core/atacseq/issues/91) with code from [[#301]](https://github.com/nf-core/atacseq/pull/301) to address shifting of reads as an option that is turned off by default.
@@ -22,7 +29,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Updated pipeline template to [nf-core/tools 3.1.1](https://github.com/nf-core/tools/releases/tag/3.1.1)
 - Updated pipeline template to [nf-core/tools 3.4.1](https://github.com/nf-core/tools/releases/tag/3.4.1)
 - [[#427](https://github.com/nf-core/atacseq/issues/427)] - Implements default nf-test at the pipeline level.
-- [[436](https://github.com/nf-core/atacseq/issues/437)] - Fix strick syntax.
+- [[#436]](https://github.com/nf-core/atacseq/pull/436) - Fix strict syntax.
 - [[437](https://github.com/nf-core/atacseq/issues/437)] - Follow up to 436.
 - [[#438](https://github.com/nf-core/atacseq/issues/438)] - Add `checkIfExists` to file inputs in `PREPARE_GENOME` and prevent S3 access errors during index validation.
 - [[PR #443](https://github.com/nf-core/atacseq/pull/443)] - Updated pipeline template to [nf-core/tools 4.0.2](https://github.com/nf-core/tools/releases/tag/4.0.2)

bin/featurecounts_merge.sh

@@ -0,0 +1,60 @@
+#!/usr/bin/env bash
+# Merge multiple featureCounts count tables into a single matrix.
+#
+# The consensus-peak quantification runs featureCounts once per library type
+# (single-end / paired-end), so each input table shares an identical annotation
+# block (Geneid, Chr, Start, End, Strand, Length) computed from the same SAF, and
+# differs only in its per-sample count columns. This column-binds those sample
+# columns back together, keyed on Geneid, and reproduces the featureCounts output
+# layout expected downstream by deseq2_qc.r:
+#
+#   line 1        : a "# Program:featureCounts" comment (skipped via read.delim skip=1)
+#   line 2        : header  Geneid  Chr  Start  End  Strand  Length  <sample cols...>
+#   remaining rows: annotation columns 1-6 followed by one count per sample
+#
+# With a single input table this is an order-preserving pass-through.
+#
+# Usage: featurecounts_merge.sh OUTFILE INPUT1 [INPUT2 ...]
+set -euo pipefail
+
+if [ "$#" -lt 2 ]; then
+    echo "Usage: $(basename "$0") OUTFILE INPUT1 [INPUT2 ...]" >&2
+    exit 1
+fi
+
+out=$1
+shift
+
+awk '
+    BEGIN { FS = OFS = "\t" }
+
+    # Skip the leading "# Program:featureCounts ..." comment of every input file.
+    FNR == 1 { fidx++; next }
+
+    # Header row: keep the 6 annotation columns from the first file only,
+    # then append every input file`s sample columns (7..NF) in argument order.
+    FNR == 2 {
+        if (fidx == 1) { hdr = $1; for (i = 2; i <= 6; i++) hdr = hdr OFS $i }
+        for (i = 7; i <= NF; i++) hdr = hdr OFS $i
+        next
+    }
+
+    # Data rows: index on Geneid (column 1). Preserve the first file`s row order
+    # and annotation; append sample counts from each file for the matching Geneid.
+    {
+        key = $1
+        if (fidx == 1) {
+            order[++n] = key
+            a = $1; for (i = 2; i <= 6; i++) a = a OFS $i; ann[key] = a
+            v = "";  for (i = 7; i <= NF; i++) v = v OFS $i; val[key] = v
+        } else {
+            for (i = 7; i <= NF; i++) val[key] = val[key] OFS $i
+        }
+    }
+
+    END {
+        print "# Program:featureCounts (merged single-end and paired-end libraries)"
+        print hdr
+        for (j = 1; j <= n; j++) print ann[order[j]] val[order[j]]
+    }
+' "$@" > "$out"

conf/modules.config

@@ -88,7 +88,7 @@ process {
         ]
     }
 
-    withName: 'CUSTOM_GETCHROMSIZES' {
+    withName: 'SAMTOOLS_FAIDX' {
         publishDir = [
             path: { "${params.outdir}/genome" },
             mode: params.publish_dir_mode,
@@ -246,7 +246,7 @@ process {
             ext.args   = {
                 [
                     meta.read_group ? "--rg-id ${meta.id} --rg SM:${meta.id - ~/_T\d+$/} --rg PL:ILLUMINA --rg LB:${meta.id} --rg PU:1" : '',
-                    params.seq_center ? "--rg CN:${params.seq_center}" : ''
+                    params.seq_center ? "--rg 'CN:${params.seq_center}'" : ''
                 ].join(' ').trim()
             }
             ext.prefix = { "${meta.id}.Lb" }
@@ -293,6 +293,28 @@ process {
         }
     }
 
+    process {
+        withName: '.*:FASTQ_ALIGN_CHROMAP:PICARD_ADDORREPLACEREADGROUPS' {
+            ext.args   = {
+                [
+                    "--RGID ${meta.id}",
+                    "--RGSM ${meta.id - ~/_T\d+$/}",
+                    "--RGPL ILLUMINA",
+                    "--RGLB ${meta.id}",
+                    "--RGPU 1",
+                    params.seq_center ? "--RGCN '${params.seq_center}'" : ''
+                ].join(' ').trim()
+            }
+            ext.prefix = { "${meta.id}.Lb.rg" }
+            publishDir = [
+                path: { "${params.outdir}/${params.aligner}/library" },
+                mode: params.publish_dir_mode,
+                saveAs: { filename -> filename.equals('versions.yml') ? null : filename },
+                enabled: false
+            ]
+        }
+    }
+
     process {
         withName: 'STAR_ALIGN' {
             ext.args   = [
@@ -683,6 +705,15 @@ process {
 
         withName: '.*:MERGED_LIBRARY_CONSENSUS_PEAKS:SUBREAD_FEATURECOUNTS'  {
             ext.args   = '-F SAF -O --fracOverlap 0.2'
+            ext.prefix = { "consensus_peaks.mLb.clN.${meta.single_end ? 'SE' : 'PE'}" }
+            publishDir = [
+                path: { "${params.outdir}/${params.aligner}/merged_library/macs3/${params.narrow_peak ? 'narrow_peak' : 'broad_peak'}/consensus" },
+                mode: params.publish_dir_mode,
+                saveAs: { filename -> filename.equals('versions.yml') ? null : filename }
+            ]
+        }
+
+        withName: '.*:MERGED_LIBRARY_CONSENSUS_PEAKS:FEATURECOUNTS_MERGE'  {
             ext.prefix = "consensus_peaks.mLb.clN"
             publishDir = [
                 path: { "${params.outdir}/${params.aligner}/merged_library/macs3/${params.narrow_peak ? 'narrow_peak' : 'broad_peak'}/consensus" },
@@ -924,6 +955,15 @@ process {
 
             withName: '.*:MERGED_REPLICATE_CONSENSUS_PEAKS:SUBREAD_FEATURECOUNTS'  {
                 ext.args   = '-F SAF -O --fracOverlap 0.2'
+                ext.prefix = { "consensus_peaks.mRp.clN.${meta.single_end ? 'SE' : 'PE'}" }
+                publishDir = [
+                    path: { "${params.outdir}/${params.aligner}/merged_replicate/macs3/${params.narrow_peak ? 'narrow_peak' : 'broad_peak'}/consensus" },
+                    mode: params.publish_dir_mode,
+                    saveAs: { filename -> filename.equals('versions.yml') ? null : filename }
+                ]
+            }
+
+            withName: '.*:MERGED_REPLICATE_CONSENSUS_PEAKS:FEATURECOUNTS_MERGE'  {
                 ext.prefix = "consensus_peaks.mRp.clN"
                 publishDir = [
                     path: { "${params.outdir}/${params.aligner}/merged_replicate/macs3/${params.narrow_peak ? 'narrow_peak' : 'broad_peak'}/consensus" },

docs/CONTRIBUTING.md

@@ -121,10 +121,10 @@ Please also refer to the [pipeline-specific contribution guidelines](#pipeline-s
 
 - [ ] Define the corresponding [input channel](#channel-naming-schemes) into your new process from the expected previous process channel.
 - [ ] Install a module with nf-core/tools, or write a local module (see [default processes resource requirements](#default-processes-resource-requirements)), and add it to the target `<workflow>.nf`.
-- [ ] Define the output channel if needed. Mix the version output channel into `ch_versions` and relevant files into `ch_multiqc`.
+- [ ] Define the output channel if needed. Emit software versions to the `versions` topic channel (`topic: versions`) from the module's `output:` block — they are collated automatically in `workflows/atacseq.nf`. Add relevant report files to the MultiQC inputs.
 - [ ] Add new or updated parameters to `nextflow.config` with a [default value](#default-parameter-values).
 - [ ] Add new or updated parameters and relevant help text to `nextflow_schema.json` with [nf-core/tools](#default-parameter-values).
-- [ ] Add validation for relevant parameters to the pipeline utilisation section of `utils_nfcore_\_pipeline/main.nf` subworkflow.
+- [ ] Add validation for relevant parameters to the pipeline utilisation section of `utils_nfcore_atacseq_pipeline/main.nf` subworkflow.
 - [ ] Perform local tests to validate that the new code works as expected.
   - [ ] If applicable, add a new test in the `tests` directory.
 - [ ] Update `usage.md`, `output.md`, and `citation.md` as appropriate.

docs/output.md

@@ -209,7 +209,9 @@ Various QC plots per sample including number of peaks, fold-change distribution,
 - `<ALIGNER>/merged_library/macs3/<PEAK_TYPE>/consensus/`
   - `*.bed`: Consensus peak-set across all samples in BED format.
   - `*.saf`: Consensus peak-set across all samples in SAF format. Required by featureCounts for read quantification.
-  - `*.featureCounts.txt`: Read counts across all samples relative to consensus peak-set.
+  - `*.featureCounts.tsv`: Read counts across all samples relative to consensus peak-set.
+  - `*.featureCounts.tsv.summary`: featureCounts assignment summary, used by MultiQC.
+  - `*.SE.featureCounts.tsv`, `*.PE.featureCounts.tsv` (and their `.summary` files): Per-library-type counts, before they are merged into the matrix above. featureCounts applies paired-end mode to a whole invocation, so single-end and paired-end libraries are counted in separate batches and then combined.
   - `*.annotatePeaks.txt`: HOMER peak-to-gene annotation file for consensus peaks.
   - `*.boolean.annotatePeaks.txt`: Spreadsheet representation of consensus peak-set across samples **with** gene annotation columns. The columns from individual peak files are included in this file along with the ability to filter peaks based on their presence or absence in multiple replicates/conditions.
   - `*.boolean.txt`: Spreadsheet representation of consensus peak-set across samples **without** gene annotation columns. Same as file above but without annotation columns.

main.nf

@@ -51,8 +51,6 @@ include { PIPELINE_COMPLETION     } from './subworkflows/local/utils_nfcore_atac
 workflow NFCORE_ATACSEQ {
 
     main:
-    ch_versions = channel.empty()
-
     // SUBWORKFLOW: Prepare genome files
     PREPARE_GENOME (
         params.genome,
@@ -73,7 +71,6 @@ workflow NFCORE_ATACSEQ {
         params.macs_gsize,
         params.read_length
     )
-    ch_versions = ch_versions.mix(PREPARE_GENOME.out.versions)
 
     //
     // WORKFLOW: Run nf-core/atacseq workflow
@@ -103,7 +100,6 @@ workflow NFCORE_ATACSEQ {
 
     emit:
     multiqc_report = ATACSEQ.out.multiqc_report // channel: /path/to/multiqc_report.html
-    versions       = ch_versions                // channel: [version1, version2, ...]
 }
 
 /*