nebiolabs · bwlang · Sep 8, 2026 · Sep 8, 2026 · Sep 8, 2026 · Sep 8, 2026
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
@@ -8,7 +8,8 @@ on:
     branches:
         - master
   schedule:
-    - cron: '0 0 1 * *'
+    - cron: '0 6 * * *'
+  workflow_dispatch:
 
 jobs:
   test:
@@ -17,21 +18,52 @@ jobs:
 
     steps:
       - name: Checkout
-        uses: actions/checkout@v4
+        uses: actions/checkout@v7
 
       - name: Setup Micromamba
-        uses: mamba-org/setup-micromamba@v1
+        uses: mamba-org/setup-micromamba@v3
         with:
           micromamba-version: '1.5.10-0' # any version from https://github.com/mamba-org/micromamba-releases
           post-cleanup: 'all'
 
-      - name: Setup Nextflow latest-edge
-        uses: nf-core/setup-nextflow@v2
+      - name: Setup Nextflow (minimum supported version)
+        uses: nf-core/setup-nextflow@v3
+        with:
+          version: '25.04.0'
+
+      - name: Install nf-test
+        uses: nf-core/setup-nf-test@v2
+        with:
+          version: 0.9.5
+
+      - name: Run Tests
+        run: nf-test test --verbose
+
+  test-latest-nextflow:
+    # Runs on a schedule (and manually), not on every push/PR: this catches a future
+    # Nextflow release breaking us, without blocking normal PRs on an upstream change
+    # unrelated to them.
+    if: github.event_name == 'schedule' || github.event_name == 'workflow_dispatch'
+    runs-on: ubuntu-latest
+    timeout-minutes: 60
+
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v7
+
+      - name: Setup Micromamba
+        uses: mamba-org/setup-micromamba@v3
+        with:
+          micromamba-version: '1.5.10-0' # any version from https://github.com/mamba-org/micromamba-releases
+          post-cleanup: 'all'
+
+      - name: Setup Nextflow (latest)
+        uses: nf-core/setup-nextflow@v3
         with:
           version: latest
 
       - name: Install nf-test
-        uses: nf-core/setup-nf-test@v1
+        uses: nf-core/setup-nf-test@v2
         with:
           version: 0.9.5
 

diff --git a/BENCHMARK_RESULTS.md b/BENCHMARK_RESULTS.md
@@ -0,0 +1,40 @@
+# Benchmark results: bamadap + fgumi zipper + picodup vs the current pipeline
+
+Performance changes in branch: `bamslice_zip`
+
+**Bottom line**: ~2.4x faster end-to-end on a clean dedicated node (16m/6m45s), with indistinguishable output. Correctness verified in depth (§2);
+
+## 1. What changed, and the measured impact
+
+| Change | Impact |
+| --- | --- |
+| `fastp`+`bwa`+Picard-dedup pipeline → `bamadap` + `fgumi zipper` + `picodup`, fused into single-pipe `trimAndAlign` + streamed `mergeAndPicodup` | End-to-end: **16m → 6m45s (~2.4x)** |
+| Dedup: Picard `MarkDuplicates` → `picodup`, streamed straight off `samtools merge` (no intermediate merged BAM) | **2m02s → 27s (~4.6x)**, peak RSS 8.5GB → <1GB |
+| Fixed a real `bwameth.py` bug (naive interleave-detection failed on `bamadap`'s `read/1`,`read/2` mate names, silently mis-converting bisulfite reads); fixed via `bamadap --no-mate-suffix` rather than a downstream patch | Removes a 3x slowdown *and* a silent correctness bug in the align stage |
+| `convert_methylkit_to_bed`: single-threaded `gawk` → `mawk` under `parallel --pipepart` | **7m53s → 35s (~14x)** — biggest single bottleneck in either pipeline; benefits `master` too |
+| `tasmanian` 1.x (single-threaded, capped at 2M reads, errored on some real reads) → `tasmanian-mismatch` 2.x (parallel, full library) | **3m42s → 17.5s (~13x), on more data** |
+| `fastqc` → `falco` | 57.8s → 30.8s (~2x) |
+| Fixed `gc_bias` (silently broken: `CollectGcBiasMetrics` requires an `Rscript` on PATH even when no chart is produced; `picard-slim` excludes R) | No-op `Rscript` shim; restores a previously-silent-failing QC step |
+| Right-sized cpu reservations for QC steps that don't read `task.cpus` (new `single_threaded_qc` label, fixed at 2 cpus instead of scaling with `--max_cpus`) | Frees queue slots on shared executors; no effect on single-task time |
+| Kept `bwameth --read-group` (populated from the uBAM's first `@RG`) instead of dropping it as originally planned; `zipper` still fixes the per-read `RG:Z:` tag | Avoids `picodup` mislabeling every metric row "Unknown Library" |
+| Fixed `bamadap`'s EL8 build (glibc-2.34-only symbols, `target-cpu=native` SIGILL'd on other CPUs) | Portability now solved at deploy time, not build time: Capistrano (`capistrano-rust-buildcache`) compiles a `-C target-cpu=native` binary per distinct machine type on first run and caches it, rather than shipping one generic binary |
+
+Chunk-size sweep on the real SGE cluster (18/37/55/90MB): smaller chunks buy wall-clock speed at
+the cost of more aggregate CPU-hours (7m14s/7.1 CPU-h at 18MB vs. 10m18s/4.4 CPU-h at 90MB); no
+change made to the 37MB default — it's a reasonable middle ground, not a correctness question.
+
+## 2. Correctness (final `.md.bam`, same test uBAM)
+
+| | Baseline | New | Δ |
+| --- | --- | --- | --- |
+| Total reads | 10,017,473 | 10,017,757 | +0.003% |
+| Duplicates | 1,362,823 | 1,362,959 | +0.01% |
+| PERCENT_DUPLICATION | 0.136415 | 0.136426 | matches to 4th decimal |
+| CpG methylation (Pearson r) | — | **0.999998** | 5 sites (0.0001%) differ by >1pp |
+
+Residual differences trace to expected causes (`bamadap` vs `fastp` trimming, corrected
+per-read `RG` moving optical-duplicate grouping, each dedup tool's own library-size extrapolation
+formula) — no read loss, no regression. `ngs-aggregate_results` compatibility for
+`tasmanian-mismatch` 2.x's new output format was verified against the real deployed parser
+(PR #932, merged and deployed 2026-09-07) and an end-to-end real (non-stubbed) aggregation run
+succeeded.
diff --git a/README.md b/README.md
@@ -106,7 +106,7 @@ Modify the conf/references.config file to specify your genome files
 ## Reference Genomes
 Pre-built reference genomes with methylation spike-in controls + CpG Island Annotation files:
 - **T2T CHM13**: [T2T_chm13v2.0+  meth_controls](https://neb-em-seq-sra.s3.us-east-1.amazonaws.com/T2T_chm13v2.0%2Bmeth_controls.fa.gz) [CpG Islands](https://neb-em-seq-sra.s3.us-east-1.amazonaws.com/human_t2t_cpg_islands.gtf)
-- **GRCh38**: [GRCh39 + meth controls](https://neb-em-seq-sra.s3.us-east-1.amazonaws.com/grch38_core%2Bmeth_controls.fa.gz) [CpG Islands](https://neb-em-seq-sra.s3.us-east-1.amazonaws.com/human_grch38_cpg_islands.gtf)
+- **GRCh38**: [GRCh38 + meth controls](https://neb-em-seq-sra.s3.us-east-1.amazonaws.com/grch38_core%2Bmeth_controls.fa.gz) [CpG Islands](https://neb-em-seq-sra.s3.us-east-1.amazonaws.com/human_grch38_cpg_islands.gtf)
 - **GRCm39**: [GRCm39 + meth_controls](https://neb-em-seq-sra.s3.us-east-1.amazonaws.com/grcm39%2Bmeth_controls.fa.gz) [CpG Islands](https://neb-em-seq-sra.s3.us-east-1.amazonaws.com/grcm39_cpg_islands.gtf)
 - Create your own reference by appending the [control sequences](assets/methylation_controls.fa) to your preferred genome fasta (e.g. `cat genome.fa methylation_controls.fa > genome+methylation_controls.fa`)
    | Sequence | Methylation State                                | Purpose |
@@ -157,6 +157,6 @@ nf-test test --updateSnapshot
 ```
 
 ## Upgrade
-As of July, 2026, Nextflow v24 or older is not supported anymore. 
-If using Nextlfow version 24 is your only option, 
-please uncomment `nextflow.preview.topic = true` as the top line in *main.nf* 
+This pipeline requires Nextflow >=25.04, where topic channels (used throughout for
+version reporting) are a stable feature rather than a preview one. `main.nf` checks
+this at startup and fails immediately with a clear message on older versions. 
diff --git a/conf/base.config b/conf/base.config
@@ -23,6 +23,14 @@ process {
         cpus   = { params.max_cpus ?: 8 }
     }
 
+    withLabel: dedup {
+        // Deliberately not capped by params.max_cpus: this label exists so the dedup stage
+        // can be tuned independently of the trim/align concurrency (see BENCHMARK_PLAN.md phase 6).
+        // Still capped to whatever the machine actually has, so it doesn't request more cpus
+        // than any executor can ever grant (e.g. GitHub Actions' 4-core runners).
+        cpus   = { Math.min(16, Runtime.runtime.availableProcessors()) }
+    }
+
     withLabel: low_cpu {
         cpus   = 2
         memory = { params.max_memory ?: 16.GB }
@@ -38,6 +46,14 @@ process {
         memory = { params.max_memory ?: 8.GB }
     }
 
+    withLabel: single_threaded_qc {
+        // For fastqc (falco), gc_bias, insert_size_metrics, picard_metrics, idx_stats:
+        // none of them read task.cpus, so scaling with params.max_cpus like medium_cpu
+        // would only reserve idle cpu slots, starving other queued tasks on the local
+        // executor for no speed benefit.
+        cpus   = 2
+    }
+
     executor = 'local' // 'sge'. 'slurm', etc... 
 
 

diff --git a/conf/references.config b/conf/references.config
@@ -1,10 +1,11 @@
 params {
     genomes = [
-        'T2T_chm13v2.0+bs_controls': [
-            bwa_index: "/mnt/galaxy/galaxyworks/tool-data/T2T_chm13v2.0+bs_controls/bwameth_index/T2T_chm13v2.0+bs_controls/T2T_chm13v2.0+bs_controls.fa",
-            genome_fa: "/mnt/galaxy/galaxyworks/tool-data/T2T_chm13v2.0+bs_controls/bwameth_index/T2T_chm13v2.0+bs_controls/T2T_chm13v2.0+bs_controls.fa",
-            genome_fai: "/mnt/galaxy/galaxyworks/tool-data/T2T_chm13v2.0+bs_controls/bwameth_index/T2T_chm13v2.0+bs_controls/T2T_chm13v2.0+bs_controls.fa.fai",
-            target_bed: "/bioinfo/ref/t2t_chm13_v2+meth_controls+m13+phix/folded_windows/T2T_chm13v2.0+meth_controls_hairpins.s18w200nogt.bed"
+        'T2T_chm13v2.0+meth_controls': [
+            bwa_index:  "/bioinfo/ref/t2t_chm13v2+meth_controls/bwameth_index/T2T_chm13v2.0+bs_controls.fa",
+            genome_fa:  "/bioinfo/ref/t2t_chm13v2+meth_controls/bwameth_index/T2T_chm13v2.0+bs_controls.fa",
+            genome_fai: "/bioinfo/ref/t2t_chm13v2+meth_controls/bwameth_index/T2T_chm13v2.0+bs_controls.fa.fai",
+            genome_dict:"/bioinfo/ref/t2t_chm13v2+meth_controls/bwameth_index/T2T_chm13v2.0+bs_controls.fa.dict",
+            target_bed: "/bioinfo/ref/t2t_chm13v2+meth_controls/folded_windows/T2T_chm13v2.0+meth_controls_hairpins.s18w200nogt.bed"
         ],
         'grcm39+meth_controls': [
             bwa_index: "/mnt/galaxy/data/genome/grcm39+meth_controls/bwameth_index/grcm39+meth_controls/grcm39+meth_controls.fa",
@@ -40,6 +41,7 @@ params {
             bwa_index: "${projectDir}/tests/fixtures/reference_files/reference.fa",
             genome_fa: "${projectDir}/tests/fixtures/reference_files/reference.fa",
             genome_fai: "${projectDir}/tests/fixtures/reference_files/reference.fa.fai",
+            genome_dict: "${projectDir}/tests/fixtures/reference_files/reference.dict",
             target_bed: "${projectDir}/tests/fixtures/target_bed/emseq_test_regions.bed"
         ]
     ]

diff --git a/conf/test.config b/conf/test.config
@@ -24,6 +24,12 @@ params {
     max_time   = '6.h'
 
     outputDir = "test_output"
+
+    // Not yet on bioconda. nf-test's pipeline tests run trimAndAlign/mergeAndPicodup under
+    // -stub-run (see their stub: blocks), so CI never actually invokes these paths -- they
+    // only matter for a manual, non-stub `nextflow run main.nf -profile test`.
+    bamadap_bin = '/appdev/langhorst/src/bamadap/target-el8/release/bamadap'
+    picodup_bin = '/appdev/langhorst/src/picodup/target-el8/release/picodup'
     }
 
     trace.overwrite = true

diff --git a/main.nf b/main.nf
@@ -1,12 +1,14 @@
-// nextflow.preview.topic = true
+// Minimum version: topic channels are stable here, preview (and unsupported) before it.
+if (!nextflow.version.matches('>=25.04')) {
+    error "This pipeline requires Nextflow >=25.04. Found: ${nextflow.version}."
+}
 
 include { registerEmailNotifications }                        from './lib/notifications.nf'
 include { createVersionsFile }                                from './lib/versions.nf'
 include { format_ngs_agg_opts }                               from './modules/aggregate_results'
-include { fastp }                                             from './modules/fastp'
 include { mergeFastpJson }                                    from './modules/merge_fastp_json'
-include { alignReads }                                        from './modules/align_reads'
-include { mergeAndMarkDuplicates }                            from './modules/merge_and_mark_duplicates'
+include { trimAndAlign }                                      from './modules/trim_and_align'
+include { mergeAndPicodup }                                   from './modules/merge_and_picodup'
 include { methylDackel_mbias }                                from './modules/methyldackel_mbias'
 include { methylDackel_extract }                              from './modules/methyldackel_extract'
 include { extract_cytosine_report }                           from './modules/extract_cytosine_report'
@@ -54,6 +56,7 @@ workflow {
          def reference_list = params.genomes[params.genome]
          genome_fa = channel.value(reference_list.genome_fa)
          genome_fai = channel.value(reference_list.genome_fai)
+         genome_dict = channel.value(reference_list.genome_dict)
 
         // Validate reference indices exist before running
         def ref = reference_list.bwa_index
@@ -63,6 +66,9 @@ workflow {
         if (!file("${reference_list.genome_fa}.fai").exists()) {
             exit 1, "Fasta index (.fai) not found for ${reference_list.genome_fa}. Run: samtools faidx ${reference_list.genome_fa}"
         }
+        if (!reference_list.genome_dict || !file(reference_list.genome_dict).exists()) {
+            exit 1, "Sequence dictionary not found for ${reference_list.genome_fa}. Run: samtools dict ${reference_list.genome_fa} -o <ref>.dict"
+        }
 
         passed_bams = bams.filter { _library, bam -> checkFileSize(bam) }
         failed_bams = bams.filter { _library, bam -> !checkFileSize(bam) }
@@ -99,23 +105,18 @@ workflow {
             }
         }
 
-        fastp( passed_bams.combine(bam_chunks, by:0), adapter_fasta )
-        mergeFastpJson( fastp.out.fastp_json.groupTuple() )
+        trimAndAlign(
+            passed_bams.combine(bam_chunks, by:0),
+            adapter_fasta,
+            reference_list.bwa_index,
+            genome_fa,
+            genome_fai,
+            genome_dict
+        )
+        mergeFastpJson( trimAndAlign.out.fastp_json.groupTuple() )
 
-        fastq_chunks = fastp.out.trimmed_fastq.map { library, chunk_name, fq_files ->
-            def fq_list = fq_files instanceof List ? fq_files : [fq_files]
-            if (params.single_end) {
-                tuple(library, chunk_name, fq_list[0])
-            }
-            else {
-                def r1 = fq_list.find { fq -> fq.name.contains('.1.trimmed.fastq') }
-                def r2 = fq_list.find { fq -> fq.name.contains('.2.trimmed.fastq') }
-                tuple(library, chunk_name, [r1, r2])
-            }
-        }
-        alignReads( passed_bams.combine(fastq_chunks, by:0), reference_list.bwa_index )
-        mergeAndMarkDuplicates( alignReads.out.bam_files.groupTuple() )
-        md_bams = mergeAndMarkDuplicates.out.md_bams
+        mergeAndPicodup( trimAndAlign.out.bam_files.groupTuple() )
+        md_bams = mergeAndPicodup.out.md_bams
 
         ///////// Methylation Calling //////////
         methylDackel_extract( md_bams, genome_fa, genome_fai )
@@ -157,18 +158,18 @@ workflow {
         fastqc( md_bams )
         picard_metrics( md_bams, genome_fa, genome_fai )
         tasmanian( md_bams, genome_fa, genome_fai )
-        combine_nonconverted_counts( alignReads.out.nonconverted_counts.groupTuple() )
+        combine_nonconverted_counts( trimAndAlign.out.nonconverted_counts.groupTuple() )
         find_switchback_reads( md_bams, genome_fa )
 
         //////// Collect files for internal summaries //////////
         agg_opts = [
-        ['--bam', mergeAndMarkDuplicates.out.md_bams.map{ row -> tuple(row[0], row[1]) }],
-        ['--bai', mergeAndMarkDuplicates.out.md_bams.map{ row -> tuple(row[0], row[2]) }],
+        ['--bam', mergeAndPicodup.out.md_bams.map{ row -> tuple(row[0], row[1]) }],
+        ['--bai', mergeAndPicodup.out.md_bams.map{ row -> tuple(row[0], row[2]) }],
         ['--metadata_bam_file', bams],
         ['--fastp', mergeFastpJson.out.merged_json],
         ['--aln', picard_metrics.out.for_agg ],
         ['--gc', gc_bias.out.for_agg ],
-        ['--dup', mergeAndMarkDuplicates.out.log],
+        ['--dup', mergeAndPicodup.out.log],
         ['--idx_stats', idx_stats.out.for_agg],
         ['--flagstat', flagstats.out.for_agg],
         ['--fastqc', fastqc.out.for_agg],

diff --git a/modules/align_reads.nf b/modules/align_reads.nf