diff --git a/bin/featurecounts_merge.sh b/bin/featurecounts_merge.sh
new file mode 100755
index 00000000..93839275
--- /dev/null
+++ b/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"
diff --git a/conf/modules.config b/conf/modules.config
index 54d8cc62..009eea2d 100644
--- a/conf/modules.config
+++ b/conf/modules.config
@@ -705,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" },
@@ -946,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" },
diff --git a/modules.json b/modules.json
index 717b50ac..94a6afc8 100644
--- a/modules.json
+++ b/modules.json
@@ -164,8 +164,7 @@
                     "subread/featurecounts": {
                         "branch": "master",
                         "git_sha": "6d46786420b4d7bc88eba026eb389c0c5535d120",
-                        "installed_by": ["modules"],
-                        "patch": "modules/nf-core/subread/featurecounts/subread-featurecounts.diff"
+                        "installed_by": ["modules"]
                     },
                     "trimgalore": {
                         "branch": "master",
diff --git a/modules/local/featurecounts_merge.nf b/modules/local/featurecounts_merge.nf
new file mode 100644
index 00000000..48d5113d
--- /dev/null
+++ b/modules/local/featurecounts_merge.nf
@@ -0,0 +1,33 @@
+process FEATURECOUNTS_MERGE {
+    tag "${meta.id}"
+    label 'process_single'
+
+    conda "conda-forge::sed=4.7"
+    container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ?
+        'https://depot.galaxyproject.org/singularity/ubuntu:20.04' :
+        'nf-core/ubuntu:20.04' }"
+
+    input:
+    tuple val(meta), path('featurecounts/*')
+
+    output:
+    tuple val(meta), path("*.featureCounts.tsv"), emit: counts
+    tuple val("${task.process}"), val('sed'), eval("sed --version | sed '1!d;s/.*GNU sed) //'"), topic: versions
+
+    when:
+    task.ext.when == null || task.ext.when
+
+    script:
+    def prefix = task.ext.prefix ?: "${meta.id}"
+    """
+    featurecounts_merge.sh \\
+        ${prefix}.featureCounts.tsv \\
+        \$(ls featurecounts/*.featureCounts.tsv | sort)
+    """
+
+    stub:
+    def prefix = task.ext.prefix ?: "${meta.id}"
+    """
+    touch ${prefix}.featureCounts.tsv
+    """
+}
diff --git a/modules/local/tests/featurecounts_merge.nf.test b/modules/local/tests/featurecounts_merge.nf.test
new file mode 100644
index 00000000..1f191d04
--- /dev/null
+++ b/modules/local/tests/featurecounts_merge.nf.test
@@ -0,0 +1,102 @@
+nextflow_process {
+
+    name "Test Process FEATURECOUNTS_MERGE"
+    script "../featurecounts_merge.nf"
+    process "FEATURECOUNTS_MERGE"
+
+    tag "modules"
+    tag "modules_local"
+    tag "featurecounts_merge"
+
+    test("merge single-end and paired-end count tables") {
+
+        when {
+            process {
+                """
+                input[0] = [
+                    [ id:'consensus_peaks' ],
+                    [
+                        file("\${projectDir}/modules/local/tests/fixtures/se.featureCounts.tsv", checkIfExists: true),
+                        file("\${projectDir}/modules/local/tests/fixtures/pe.featureCounts.tsv", checkIfExists: true)
+                    ]
+                ]
+                """
+            }
+        }
+
+        then {
+            def merged = path(process.out.counts.get(0).get(1)).readLines()
+            def header = merged.find { it.startsWith('Geneid') }.split('\t')
+            def sampleCols = header[6..-1]
+            assertAll(
+                { assert process.success },
+                // one merged table emitted
+                { assert process.out.counts.size() == 1 },
+                // The module merges its inputs in `ls | sort` filename order, so the
+                // column order is deterministic. Here 'pe.featureCounts.tsv' sorts
+                // before 'se.featureCounts.tsv' (in the pipeline the files are
+                // '...PE.featureCounts.tsv' / '...SE.featureCounts.tsv', so PE-first).
+                // Every SE and PE sample column appears exactly once.
+                { assert sampleCols == ['T0_PE.bam', 'T15_PE.bam', 'T100_SE.bam', 'T150_SE.bam'] },
+                { assert sampleCols.toUnique().size() == sampleCols.size() },
+                // Annotation + row order preserved from the first-sorted (PE) table;
+                // SE counts matched to rows by Geneid despite the different row order.
+                { assert merged.contains('peak_3\tIII\t9\t309\t+\t301\t131\t132\t31\t32') },
+                { assert merged.contains('peak_1\tI\t1\t100\t+\t100\t111\t112\t11\t12') },
+                { assert merged.contains('peak_2\tII\t5\t205\t+\t201\t121\t122\t21\t22') },
+                { assert snapshot(process.out.versions).match("versions") }
+            )
+        }
+    }
+
+    test("single-file pass-through (homogeneous cohort)") {
+
+        when {
+            process {
+                """
+                input[0] = [
+                    [ id:'consensus_peaks' ],
+                    [ file("\${projectDir}/modules/local/tests/fixtures/se.featureCounts.tsv", checkIfExists: true) ]
+                ]
+                """
+            }
+        }
+
+        then {
+            def merged = path(process.out.counts.get(0).get(1)).readLines()
+            def header = merged.find { it.startsWith('Geneid') }.split('\t')
+            assertAll(
+                { assert process.success },
+                { assert (header[6..-1]) == ['T100_SE.bam', 'T150_SE.bam'] },
+                { assert merged.contains('peak_1\tI\t1\t100\t+\t100\t11\t12') },
+                { assert merged.contains('peak_3\tIII\t9\t309\t+\t301\t31\t32') }
+            )
+        }
+    }
+
+    test("stub") {
+
+        options "-stub"
+
+        when {
+            process {
+                """
+                input[0] = [
+                    [ id:'consensus_peaks' ],
+                    [
+                        file("\${projectDir}/modules/local/tests/fixtures/se.featureCounts.tsv", checkIfExists: true),
+                        file("\${projectDir}/modules/local/tests/fixtures/pe.featureCounts.tsv", checkIfExists: true)
+                    ]
+                ]
+                """
+            }
+        }
+
+        then {
+            assertAll(
+                { assert process.success },
+                { assert snapshot(process.out).match() }
+            )
+        }
+    }
+}
diff --git a/modules/local/tests/featurecounts_merge.nf.test.snap b/modules/local/tests/featurecounts_merge.nf.test.snap
new file mode 100644
index 00000000..47a579bf
--- /dev/null
+++ b/modules/local/tests/featurecounts_merge.nf.test.snap
@@ -0,0 +1,44 @@
+{
+    "versions": {
+        "content": null,
+        "timestamp": "2026-07-23T03:28:47.894742",
+        "meta": {
+            "nf-test": "0.9.4",
+            "nextflow": "26.04.4"
+        }
+    },
+    "stub": {
+        "content": [
+            {
+                "0": [
+                    [
+                        {
+                            "id": "consensus_peaks"
+                        },
+                        "consensus_peaks.featureCounts.tsv:md5,d41d8cd98f00b204e9800998ecf8427e"
+                    ]
+                ],
+                "1": [
+                    [
+                        "FEATURECOUNTS_MERGE",
+                        "sed",
+                        ""
+                    ]
+                ],
+                "counts": [
+                    [
+                        {
+                            "id": "consensus_peaks"
+                        },
+                        "consensus_peaks.featureCounts.tsv:md5,d41d8cd98f00b204e9800998ecf8427e"
+                    ]
+                ]
+            }
+        ],
+        "timestamp": "2026-07-23T03:28:54.791169",
+        "meta": {
+            "nf-test": "0.9.4",
+            "nextflow": "26.04.4"
+        }
+    }
+}
\ No newline at end of file
diff --git a/modules/local/tests/fixtures/pe.featureCounts.tsv b/modules/local/tests/fixtures/pe.featureCounts.tsv
new file mode 100644
index 00000000..bbd6a111
--- /dev/null
+++ b/modules/local/tests/fixtures/pe.featureCounts.tsv
@@ -0,0 +1,5 @@
+# Program:featureCounts v2.1.1; Command:"featureCounts" "-F" "SAF" "-p"
+Geneid	Chr	Start	End	Strand	Length	T0_PE.bam	T15_PE.bam
+peak_3	III	9	309	+	301	131	132
+peak_1	I	1	100	+	100	111	112
+peak_2	II	5	205	+	201	121	122
diff --git a/modules/local/tests/fixtures/se.featureCounts.tsv b/modules/local/tests/fixtures/se.featureCounts.tsv
new file mode 100644
index 00000000..da0e2690
--- /dev/null
+++ b/modules/local/tests/fixtures/se.featureCounts.tsv
@@ -0,0 +1,5 @@
+# Program:featureCounts v2.1.1; Command:"featureCounts" "-F" "SAF"
+Geneid	Chr	Start	End	Strand	Length	T100_SE.bam	T150_SE.bam
+peak_1	I	1	100	+	100	11	12
+peak_2	II	5	205	+	201	21	22
+peak_3	III	9	309	+	301	31	32
diff --git a/modules/nf-core/subread/featurecounts/environment.yml b/modules/nf-core/subread/featurecounts/environment.yml
index 4b433c9f..23133559 100644
--- a/modules/nf-core/subread/featurecounts/environment.yml
+++ b/modules/nf-core/subread/featurecounts/environment.yml
@@ -4,8 +4,4 @@ channels:
   - conda-forge
   - bioconda
 dependencies:
-  # PINNED: subread held at 2.0.1. 2.1.1 strictly rejects mixed single-end /
-  # paired-end BAMs in one featureCounts call (both -p and no -p error), which
-  # breaks consensus-peak quantification for mixed cohorts. 2.0.1 auto-detects
-  # per BAM. Remove the pin once mixed-library counting is handled explicitly.
-  - bioconda::subread=2.0.1
+  - bioconda::subread=2.1.1
diff --git a/modules/nf-core/subread/featurecounts/main.nf b/modules/nf-core/subread/featurecounts/main.nf
index a8c5facf..7e99581b 100644
--- a/modules/nf-core/subread/featurecounts/main.nf
+++ b/modules/nf-core/subread/featurecounts/main.nf
@@ -4,8 +4,8 @@ process SUBREAD_FEATURECOUNTS {
 
     conda "${moduleDir}/environment.yml"
     container "${workflow.containerEngine in ['singularity', 'apptainer'] && !task.ext.singularity_pull_docker_container
-        ? 'https://depot.galaxyproject.org/singularity/subread:2.0.1--hed695b0_0'
-        : 'quay.io/biocontainers/subread:2.0.1--hed695b0_0'}"
+        ? 'https://depot.galaxyproject.org/singularity/subread:2.1.1--h577a1d6_0'
+        : 'quay.io/biocontainers/subread:2.1.1--h577a1d6_0'}"
 
     input:
     tuple val(meta), path(bams), path(annotation)
diff --git a/modules/nf-core/subread/featurecounts/subread-featurecounts.diff b/modules/nf-core/subread/featurecounts/subread-featurecounts.diff
deleted file mode 100644
index a920bc96..00000000
--- a/modules/nf-core/subread/featurecounts/subread-featurecounts.diff
+++ /dev/null
@@ -1,35 +0,0 @@
-Changes in component 'nf-core/subread/featurecounts'
-Changes in 'subread/featurecounts/environment.yml':
---- modules/nf-core/subread/featurecounts/environment.yml
-+++ modules/nf-core/subread/featurecounts/environment.yml
-@@ -4,4 +4,8 @@
-   - conda-forge
-   - bioconda
- dependencies:
--  - bioconda::subread=2.1.1
-+  # PINNED: subread held at 2.0.1. 2.1.1 strictly rejects mixed single-end /
-+  # paired-end BAMs in one featureCounts call (both -p and no -p error), which
-+  # breaks consensus-peak quantification for mixed cohorts. 2.0.1 auto-detects
-+  # per BAM. Remove the pin once mixed-library counting is handled explicitly.
-+  - bioconda::subread=2.0.1
-
-'modules/nf-core/subread/featurecounts/meta.yml' is unchanged
-Changes in 'subread/featurecounts/main.nf':
---- modules/nf-core/subread/featurecounts/main.nf
-+++ modules/nf-core/subread/featurecounts/main.nf
-@@ -4,8 +4,8 @@
- 
-     conda "${moduleDir}/environment.yml"
-     container "${workflow.containerEngine in ['singularity', 'apptainer'] && !task.ext.singularity_pull_docker_container
--        ? 'https://depot.galaxyproject.org/singularity/subread:2.1.1--h577a1d6_0'
--        : 'quay.io/biocontainers/subread:2.1.1--h577a1d6_0'}"
-+        ? 'https://depot.galaxyproject.org/singularity/subread:2.0.1--hed695b0_0'
-+        : 'quay.io/biocontainers/subread:2.0.1--hed695b0_0'}"
- 
-     input:
-     tuple val(meta), path(bams), path(annotation)
-
-'modules/nf-core/subread/featurecounts/tests/main.nf.test.snap' is unchanged
-'modules/nf-core/subread/featurecounts/tests/nextflow.config' is unchanged
-'modules/nf-core/subread/featurecounts/tests/main.nf.test' is unchanged
-************************************************************
diff --git a/subworkflows/local/bed_consensus_quantify_qc_bedtools_featurecounts_deseq2.nf b/subworkflows/local/bed_consensus_quantify_qc_bedtools_featurecounts_deseq2.nf
index dcf506c5..5dbae1f9 100644
--- a/subworkflows/local/bed_consensus_quantify_qc_bedtools_featurecounts_deseq2.nf
+++ b/subworkflows/local/bed_consensus_quantify_qc_bedtools_featurecounts_deseq2.nf
@@ -6,6 +6,7 @@ include { HOMER_ANNOTATEPEAKS    } from '../../modules/nf-core/homer/annotatepea
 include { SUBREAD_FEATURECOUNTS  } from '../../modules/nf-core/subread/featurecounts/main'
 
 include { MACS3_CONSENSUS        } from '../../modules/local/macs3_consensus'
+include { FEATURECOUNTS_MERGE    } from '../../modules/local/featurecounts_merge'
 include { DESEQ2_QC              } from '../../modules/local/deseq2_qc'
 
 workflow BED_CONSENSUS_QUANTIFY_QC_BEDTOOLS_FEATURECOUNTS_DESEQ2 {
@@ -55,30 +56,71 @@ workflow BED_CONSENSUS_QUANTIFY_QC_BEDTOOLS_FEATURECOUNTS_DESEQ2 {
         ch_homer_annotatepeaks = HOMER_ANNOTATEPEAKS.out.txt
     }
 
-    // Create channels: [ meta, [ bams ], saf ]
-    // The bam list comes from an unordered channel collect, so its order (and
-    // therefore the featureCounts column order and output file) is otherwise
-    // non-deterministic across runs/hosts. Sort by filename so the consensus
-    // count matrix is reproducible and its snapshot is stable.
+    //
+    // Quantify peaks across samples with featureCounts.
+    //
+    // featureCounts (subread >= 2.1.0) applies paired-end mode (-p) to a whole
+    // invocation and aborts when that invocation mixes single-end and paired-end
+    // BAMs. The consensus BAMs can span both library types, so split them by
+    // endedness, count each homogeneous batch with the correct pairing flag
+    // (derived from meta.single_end inside SUBREAD_FEATURECOUNTS), then merge the
+    // per-batch matrices back into one consensus table for DESeq2 and MultiQC.
+    // The join with ch_peaks keeps only samples that contributed peaks; combining
+    // with MACS3_CONSENSUS.out.saf also gates counting on a consensus existing
+    // (>= 2 samples), matching the previous behaviour.
+    //
+    ch_consensus_saf = MACS3_CONSENSUS.out.saf.map { _meta, saf -> saf }
+
+    // The merged-library caller joins in a control-BAM column
+    // ([ meta, bam, control ] -> [ meta, bam, control, peak ]) while the
+    // merged-replicate caller does not ([ meta, bams ] -> [ meta, bams, peak ]),
+    // so the joined tuple arity differs between the two instantiations of this
+    // subworkflow. Index positionally (meta = item[0], bam = item[1]) to stay
+    // tolerant of both shapes, as the pre-split implementation did.
     ch_bams
         .join(ch_peaks)
-        .collect { item -> item[1] }
-        .filter { item -> item.size() > 1 }
-        .map { item -> [ item ] }
-        .concat(MACS3_CONSENSUS.out.saf)
-        .collect()
-        .filter { item -> item.size() == 3 }
-        .map {
-            bam, meta, saf ->
-                [ meta, bam.toSorted { it.name }, saf ]
+        .branch { item ->
+            single_end: item[0].single_end
+            paired_end: !item[0].single_end
         }
-        .set { ch_bam_saf }
+        .set { ch_consensus_bams }
+
+    // Each batch is assembled from an unordered channel collect, so sort by
+    // filename: the BAM order sets the featureCounts column order, and an
+    // unsorted list makes the count matrix (and its snapshot md5) vary between
+    // runs and hosts.
+    ch_se_batch = ch_consensus_bams.single_end
+        .map { item -> item[1] }
+        .collect()
+        .filter { bams -> bams }
+        .map { bams -> [ [ id: 'consensus_peaks', single_end: true ], bams.toSorted { it.name } ] }
+
+    ch_pe_batch = ch_consensus_bams.paired_end
+        .map { item -> item[1] }
+        .collect()
+        .filter { bams -> bams }
+        .map { bams -> [ [ id: 'consensus_peaks', single_end: false ], bams.toSorted { it.name } ] }
+
+    ch_featurecounts_input = ch_se_batch
+        .mix(ch_pe_batch)
+        .combine(ch_consensus_saf)
+
+    SUBREAD_FEATURECOUNTS (
+        ch_featurecounts_input
+    )
 
     //
-    // Quantify peaks across samples with featureCounts
+    // Merge the per-library-type count matrices into a single consensus matrix
     //
-    SUBREAD_FEATURECOUNTS (
-        ch_bam_saf
+    // Sorted for the same reason: the merge script's column order follows the
+    // order of the per-batch matrices it is handed.
+    ch_merged_counts = SUBREAD_FEATURECOUNTS.out.counts
+        .map { _meta, counts -> counts }
+        .collect()
+        .map { counts -> [ [ id: 'consensus_peaks' ], counts.toSorted { it.name } ] }
+
+    FEATURECOUNTS_MERGE (
+        ch_merged_counts
     )
 
     //
@@ -95,7 +137,7 @@ workflow BED_CONSENSUS_QUANTIFY_QC_BEDTOOLS_FEATURECOUNTS_DESEQ2 {
     ch_deseq2_qc_size_factors  = channel.empty()
     if (!skip_deseq2_qc) {
         DESEQ2_QC (
-            SUBREAD_FEATURECOUNTS.out.counts,
+            FEATURECOUNTS_MERGE.out.counts,
             ch_deseq2_pca_header_multiqc,
             ch_deseq2_clustering_header_multiqc
         )
@@ -119,8 +161,8 @@ workflow BED_CONSENSUS_QUANTIFY_QC_BEDTOOLS_FEATURECOUNTS_DESEQ2 {
 
     homer_annotatepeaks     = ch_homer_annotatepeaks            // channel: [ txt ]
 
-    featurecounts_txt       = SUBREAD_FEATURECOUNTS.out.counts  // channel: [ txt ]
-    featurecounts_summary   = SUBREAD_FEATURECOUNTS.out.summary // channel: [ txt ]
+    featurecounts_txt       = FEATURECOUNTS_MERGE.out.counts    // channel: [ val(meta), txt ]
+    featurecounts_summary   = SUBREAD_FEATURECOUNTS.out.summary // channel: [ val(meta), txt ] (one per library type)
 
     deseq2_qc_pdf           = ch_deseq2_qc_pdf                  // channel: [ pdf ]
     deseq2_qc_rdata         = ch_deseq2_qc_rdata                // channel: [ rdata ]
diff --git a/tests/bin/test_featurecounts_merge.sh b/tests/bin/test_featurecounts_merge.sh
new file mode 100755
index 00000000..24987364
--- /dev/null
+++ b/tests/bin/test_featurecounts_merge.sh
@@ -0,0 +1,110 @@
+#!/usr/bin/env bash
+#
+# Unit tests for bin/featurecounts_merge.sh — pure bash, no Docker, milliseconds.
+#
+# The consensus-peak quantification runs featureCounts once per library type
+# (single-end / paired-end) and merges the resulting count tables. These tests
+# pin the merge logic directly: exact column-bind of a mixed SE+PE pair, and the
+# single-file pass-through used when a cohort is all-SE or all-PE.
+#
+# Run: tests/bin/test_featurecounts_merge.sh
+set -euo pipefail
+
+here="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
+merge="$here/../../bin/featurecounts_merge.sh"
+tmp="$(mktemp -d)"
+trap 'rm -rf "$tmp"' EXIT
+
+fail=0
+pass=0
+
+# assert_eq NAME EXPECTED_FILE ACTUAL_FILE
+assert_eq() {
+    local name=$1 exp=$2 act=$3
+    if diff -u "$exp" "$act" >/dev/null; then
+        echo "ok   - $name"
+        pass=$((pass + 1))
+    else
+        echo "FAIL - $name"
+        echo "----- diff (expected vs actual) -----"
+        diff -u "$exp" "$act" || true
+        echo "-------------------------------------"
+        fail=$((fail + 1))
+    fi
+}
+
+# --- fixtures --------------------------------------------------------------
+# Both tables share an identical annotation block (same SAF); they differ only
+# in their sample count columns. The PE table intentionally lists its data rows
+# in a DIFFERENT order to prove the merge keys on Geneid, not on row position.
+
+printf '%s\n' \
+'# Program:featureCounts v2.1.1; Command:"featureCounts" "-F" "SAF"' \
+$'Geneid\tChr\tStart\tEnd\tStrand\tLength\tT100_SE.bam\tT150_SE.bam' \
+$'peak_1\tI\t1\t100\t+\t100\t11\t12' \
+$'peak_2\tII\t5\t205\t+\t201\t21\t22' \
+$'peak_3\tIII\t9\t309\t+\t301\t31\t32' \
+> "$tmp/se.featureCounts.tsv"
+
+printf '%s\n' \
+'# Program:featureCounts v2.1.1; Command:"featureCounts" "-F" "SAF" "-p"' \
+$'Geneid\tChr\tStart\tEnd\tStrand\tLength\tT0_PE.bam\tT15_PE.bam' \
+$'peak_3\tIII\t9\t309\t+\t301\t131\t132' \
+$'peak_1\tI\t1\t100\t+\t100\t111\t112' \
+$'peak_2\tII\t5\t205\t+\t201\t121\t122' \
+> "$tmp/pe.featureCounts.tsv"
+
+# === Case 1: mixed SE + PE merge (SE table first) ==========================
+# Annotation and row order come from the first (SE) file; sample columns are
+# appended SE-then-PE; PE counts are matched to rows by Geneid.
+printf '%s\n' \
+'# Program:featureCounts (merged single-end and paired-end libraries)' \
+$'Geneid\tChr\tStart\tEnd\tStrand\tLength\tT100_SE.bam\tT150_SE.bam\tT0_PE.bam\tT15_PE.bam' \
+$'peak_1\tI\t1\t100\t+\t100\t11\t12\t111\t112' \
+$'peak_2\tII\t5\t205\t+\t201\t21\t22\t121\t122' \
+$'peak_3\tIII\t9\t309\t+\t301\t31\t32\t131\t132' \
+> "$tmp/expected_mixed.tsv"
+
+"$merge" "$tmp/out_mixed.tsv" "$tmp/se.featureCounts.tsv" "$tmp/pe.featureCounts.tsv"
+assert_eq "mixed SE+PE merge column-binds on Geneid" "$tmp/expected_mixed.tsv" "$tmp/out_mixed.tsv"
+
+# === Case 2: all-SE cohort -> single-file pass-through =====================
+# Empty PE branch means featureCounts runs once; the merge must pass the single
+# table through unchanged (except the normalised comment line), order preserved.
+printf '%s\n' \
+'# Program:featureCounts (merged single-end and paired-end libraries)' \
+$'Geneid\tChr\tStart\tEnd\tStrand\tLength\tT100_SE.bam\tT150_SE.bam' \
+$'peak_1\tI\t1\t100\t+\t100\t11\t12' \
+$'peak_2\tII\t5\t205\t+\t201\t21\t22' \
+$'peak_3\tIII\t9\t309\t+\t301\t31\t32' \
+> "$tmp/expected_se_only.tsv"
+
+"$merge" "$tmp/out_se_only.tsv" "$tmp/se.featureCounts.tsv"
+assert_eq "all-SE single-file pass-through" "$tmp/expected_se_only.tsv" "$tmp/out_se_only.tsv"
+
+# === Case 3: all-PE cohort -> single-file pass-through =====================
+# Same pass-through, preserving the PE file's own (unsorted) row order.
+printf '%s\n' \
+'# Program:featureCounts (merged single-end and paired-end libraries)' \
+$'Geneid\tChr\tStart\tEnd\tStrand\tLength\tT0_PE.bam\tT15_PE.bam' \
+$'peak_3\tIII\t9\t309\t+\t301\t131\t132' \
+$'peak_1\tI\t1\t100\t+\t100\t111\t112' \
+$'peak_2\tII\t5\t205\t+\t201\t121\t122' \
+> "$tmp/expected_pe_only.tsv"
+
+"$merge" "$tmp/out_pe_only.tsv" "$tmp/pe.featureCounts.tsv"
+assert_eq "all-PE single-file pass-through" "$tmp/expected_pe_only.tsv" "$tmp/out_pe_only.tsv"
+
+# === Case 4: usage error on too few arguments =============================
+if "$merge" "$tmp/out_none.tsv" >/dev/null 2>&1; then
+    echo "FAIL - merge should exit non-zero with no input files"
+    fail=$((fail + 1))
+else
+    echo "ok   - usage error when no input tables given"
+    pass=$((pass + 1))
+fi
+
+# --- summary ---------------------------------------------------------------
+echo
+echo "featurecounts_merge.sh: $pass passed, $fail failed"
+[ "$fail" -eq 0 ]
diff --git a/tests/consensus_all_pe.nf.test b/tests/consensus_all_pe.nf.test
new file mode 100644
index 00000000..779f3b09
--- /dev/null
+++ b/tests/consensus_all_pe.nf.test
@@ -0,0 +1,40 @@
+nextflow_pipeline {
+
+    name "Test pipeline - homogeneous all-paired-end cohort"
+    script "../main.nf"
+    tag "pipeline"
+    tag "consensus_endedness"
+
+    // Failure-B edge case: an all-paired-end cohort. The single-end branch of the
+    // consensus split is empty, so featureCounts runs exactly once (paired mode)
+    // and FEATURECOUNTS_MERGE takes the single-file pass-through path. The merged
+    // matrix must still be produced and carry only paired-end sample columns.
+    test("all paired-end cohort - single featureCounts batch, pass-through merge") {
+
+        when {
+            params {
+                input  = "${projectDir}/tests/csv/all_pe_samplesheet.csv"
+                outdir = "$outputDir"
+            }
+        }
+
+        then {
+            def mergedFc = []
+            new File("${outputDir}").eachFileRecurse { f ->
+                if (f.name ==~ /consensus_peaks\.m(Lb|Rp)\.clN\.featureCounts\.tsv/) {
+                    mergedFc << f
+                }
+            }
+            assertAll(
+                { assert workflow.success },
+                { assert mergedFc.size() > 0 : "no merged consensus featureCounts matrix was produced" },
+                { mergedFc.each { f ->
+                    def header = f.readLines().find { it.startsWith('Geneid') }.split('\t') as List
+                    def cols = header[6..-1]
+                    assert cols.toUnique().size() == cols.size() : "duplicate sample columns in ${f.name}: ${cols}"
+                    assert cols.every { it.contains('_PE') } : "unexpected non-PE column in ${f.name}: ${cols}"
+                } }
+            )
+        }
+    }
+}
diff --git a/tests/consensus_all_se.nf.test b/tests/consensus_all_se.nf.test
new file mode 100644
index 00000000..bdec3183
--- /dev/null
+++ b/tests/consensus_all_se.nf.test
@@ -0,0 +1,40 @@
+nextflow_pipeline {
+
+    name "Test pipeline - homogeneous all-single-end cohort"
+    script "../main.nf"
+    tag "pipeline"
+    tag "consensus_endedness"
+
+    // Failure-B edge case: an all-single-end cohort. The paired-end branch of the
+    // consensus split is empty, so featureCounts runs exactly once (single-end
+    // mode) and FEATURECOUNTS_MERGE takes the single-file pass-through path. The
+    // merged matrix must still be produced and carry only single-end sample columns.
+    test("all single-end cohort - single featureCounts batch, pass-through merge") {
+
+        when {
+            params {
+                input  = "${projectDir}/tests/csv/all_se_samplesheet.csv"
+                outdir = "$outputDir"
+            }
+        }
+
+        then {
+            def mergedFc = []
+            new File("${outputDir}").eachFileRecurse { f ->
+                if (f.name ==~ /consensus_peaks\.m(Lb|Rp)\.clN\.featureCounts\.tsv/) {
+                    mergedFc << f
+                }
+            }
+            assertAll(
+                { assert workflow.success },
+                { assert mergedFc.size() > 0 : "no merged consensus featureCounts matrix was produced" },
+                { mergedFc.each { f ->
+                    def header = f.readLines().find { it.startsWith('Geneid') }.split('\t') as List
+                    def cols = header[6..-1]
+                    assert cols.toUnique().size() == cols.size() : "duplicate sample columns in ${f.name}: ${cols}"
+                    assert cols.every { it.contains('_SE') } : "unexpected non-SE column in ${f.name}: ${cols}"
+                } }
+            )
+        }
+    }
+}
diff --git a/tests/csv/all_pe_samplesheet.csv b/tests/csv/all_pe_samplesheet.csv
new file mode 100644
index 00000000..6f76672a
--- /dev/null
+++ b/tests/csv/all_pe_samplesheet.csv
@@ -0,0 +1,5 @@
+sample,fastq_1,fastq_2,replicate
+OSMOTIC_STRESS_T0_PE,https://raw.githubusercontent.com/nf-core/test-datasets/atacseq/testdata/SRR1822153_1.fastq.gz,https://raw.githubusercontent.com/nf-core/test-datasets/atacseq/testdata/SRR1822153_2.fastq.gz,1
+OSMOTIC_STRESS_T0_PE,https://raw.githubusercontent.com/nf-core/test-datasets/atacseq/testdata/SRR1822154_1.fastq.gz,https://raw.githubusercontent.com/nf-core/test-datasets/atacseq/testdata/SRR1822154_2.fastq.gz,2
+OSMOTIC_STRESS_T15_PE,https://raw.githubusercontent.com/nf-core/test-datasets/atacseq/testdata/SRR1822157_1.fastq.gz,https://raw.githubusercontent.com/nf-core/test-datasets/atacseq/testdata/SRR1822157_2.fastq.gz,1
+OSMOTIC_STRESS_T15_PE,https://raw.githubusercontent.com/nf-core/test-datasets/atacseq/testdata/SRR1822158_1.fastq.gz,https://raw.githubusercontent.com/nf-core/test-datasets/atacseq/testdata/SRR1822158_2.fastq.gz,1
diff --git a/tests/csv/all_se_samplesheet.csv b/tests/csv/all_se_samplesheet.csv
new file mode 100644
index 00000000..2c436e37
--- /dev/null
+++ b/tests/csv/all_se_samplesheet.csv
@@ -0,0 +1,5 @@
+sample,fastq_1,fastq_2,replicate
+OSMOTIC_STRESS_T100_SE,https://raw.githubusercontent.com/nf-core/test-datasets/atacseq/testdata/SRR1822153_1.fastq.gz,,1
+OSMOTIC_STRESS_T100_SE,https://raw.githubusercontent.com/nf-core/test-datasets/atacseq/testdata/SRR1822154_1.fastq.gz,,2
+OSMOTIC_STRESS_T150_SE,https://raw.githubusercontent.com/nf-core/test-datasets/atacseq/testdata/SRR1822157_1.fastq.gz,,1
+OSMOTIC_STRESS_T150_SE,https://raw.githubusercontent.com/nf-core/test-datasets/atacseq/testdata/SRR1822158_1.fastq.gz,,1
diff --git a/tests/default.nf.test b/tests/default.nf.test
index 15421c16..31d36934 100644
--- a/tests/default.nf.test
+++ b/tests/default.nf.test
@@ -28,7 +28,29 @@ nextflow_pipeline {
                     stable_name,
                     // All files with stable contents
                     stable_path
-                ).match() }
+                ).match() },
+                // Failure-B regression: the -profile test cohort mixes single-end
+                // (OSMOTIC_STRESS_*_SE) and paired-end (OSMOTIC_STRESS_*_PE) libraries.
+                // subread >= 2.1 aborts if one featureCounts call mixes endedness, so
+                // the consensus quantification splits the cohort by library type,
+                // counts each batch, and merges the matrices. Assert the merged
+                // consensus matrix carries every SE and PE sample column exactly once.
+                {
+                    def mergedFc = []
+                    new File("${outputDir}").eachFileRecurse { f ->
+                        if (f.name ==~ /consensus_peaks\.m(Lb|Rp)\.clN\.featureCounts\.tsv/) {
+                            mergedFc << f
+                        }
+                    }
+                    assert mergedFc.size() > 0 : "no merged consensus featureCounts matrix was produced"
+                    mergedFc.each { f ->
+                        def header = f.readLines().find { it.startsWith('Geneid') }.split('\t') as List
+                        def cols = header[6..-1]
+                        assert cols.toUnique().size() == cols.size() : "duplicate sample columns in ${f.name}: ${cols}"
+                        assert cols.any { it.contains('_SE') } : "no single-end sample column in ${f.name}: ${cols}"
+                        assert cols.any { it.contains('_PE') } : "no paired-end sample column in ${f.name}: ${cols}"
+                    }
+                }
             )
         }
     }