Single-pass partitioning experiment (#323)

dev/concordance/run_illumina.sh

@@ -138,6 +138,13 @@ VERBOSE="${VERBOSE:-}"
 verbose_arg=()
 [ -n "$VERBOSE" ] && verbose_arg=(--verbose)
 METRICS="${METRICS:-}"
+# metrics_for <step>: sets metrics_arg to --metrics-json $OUT/<step>.metrics.json
+# when METRICS is set, else to nothing.
+metrics_for() {
+  metrics_arg=()
+  [ -n "$METRICS" ] && metrics_arg=(--metrics-json "$OUT/$1.metrics.json")
+  return 0
+}
 
 mkdir -p "$OUT"/{filtered,dada_fwd,dada_rev,control_persample}
 
@@ -301,14 +308,16 @@ if [ "$POOL" = "pseudo" ] || [ "$POOL" = "pseudo-fixed-err" ]; then
   else
     echo "==> dada-pseudo (fwd, rev; one error model across both rounds -- our default)"
   fi
+  metrics_for dada-pseudo.fwd
   run_step dada-pseudo.fwd "$BIN" dada-pseudo "${filtFs[@]}" --error-model "$OUT/errF.json" \
       -o "$OUT/dada_fwd" --priors-out "$OUT/priors_fwd.fasta" \
       ${reest_arg[@]+"${reest_arg[@]}"} \
-      --threads "$THREADS" ${backend_arg[@]+"${backend_arg[@]}"} ${screen_arg[@]+"${screen_arg[@]}"} ${end_arg[@]+"${end_arg[@]}"} ${verbose_arg[@]+"${verbose_arg[@]}"}
+      --threads "$THREADS" ${backend_arg[@]+"${backend_arg[@]}"} ${screen_arg[@]+"${screen_arg[@]}"} ${end_arg[@]+"${end_arg[@]}"} ${verbose_arg[@]+"${verbose_arg[@]}"} ${metrics_arg[@]+"${metrics_arg[@]}"}
+  metrics_for dada-pseudo.rev
   run_step dada-pseudo.rev "$BIN" dada-pseudo "${filtRs[@]}" --error-model "$OUT/errR.json" \
       -o "$OUT/dada_rev" --priors-out "$OUT/priors_rev.fasta" \
       ${reest_arg[@]+"${reest_arg[@]}"} \
-      --threads "$THREADS" ${backend_arg[@]+"${backend_arg[@]}"} ${screen_arg[@]+"${screen_arg[@]}"} ${end_arg[@]+"${end_arg[@]}"} ${verbose_arg[@]+"${verbose_arg[@]}"}
+      --threads "$THREADS" ${backend_arg[@]+"${backend_arg[@]}"} ${screen_arg[@]+"${screen_arg[@]}"} ${end_arg[@]+"${end_arg[@]}"} ${verbose_arg[@]+"${verbose_arg[@]}"} ${metrics_arg[@]+"${metrics_arg[@]}"}
 
   # Positive control: did the priors actually CHANGE anything?
   #
@@ -355,19 +364,23 @@ elif [ "$POOL" = "true" ]; then
   # MEMORY: pooling holds every sample's uniques at once. Budget accordingly --
   # docs/findings has pooled peaks in the tens of GB for diverse soil runs.
   echo "==> dada-pooled (fwd, rev; full pooling)"
+  metrics_for dada-pooled.fwd
   run_step dada-pooled.fwd "$BIN" dada-pooled "${filtFs[@]}" --error-model "$OUT/errF.json" \
       -o "$OUT/dada_fwd" --threads "$THREADS" \
-      ${backend_arg[@]+"${backend_arg[@]}"} ${screen_arg[@]+"${screen_arg[@]}"} ${end_arg[@]+"${end_arg[@]}"} ${verbose_arg[@]+"${verbose_arg[@]}"}
+      ${backend_arg[@]+"${backend_arg[@]}"} ${screen_arg[@]+"${screen_arg[@]}"} ${end_arg[@]+"${end_arg[@]}"} ${verbose_arg[@]+"${verbose_arg[@]}"} ${metrics_arg[@]+"${metrics_arg[@]}"}
+  metrics_for dada-pooled.rev
   run_step dada-pooled.rev "$BIN" dada-pooled "${filtRs[@]}" --error-model "$OUT/errR.json" \
       -o "$OUT/dada_rev" --threads "$THREADS" \
-      ${backend_arg[@]+"${backend_arg[@]}"} ${screen_arg[@]+"${screen_arg[@]}"} ${end_arg[@]+"${end_arg[@]}"} ${verbose_arg[@]+"${verbose_arg[@]}"}
+      ${backend_arg[@]+"${backend_arg[@]}"} ${screen_arg[@]+"${screen_arg[@]}"} ${end_arg[@]+"${end_arg[@]}"} ${verbose_arg[@]+"${verbose_arg[@]}"} ${metrics_arg[@]+"${metrics_arg[@]}"}
 else
   # Per-sample denoising (pool=FALSE analog) — matches R dada() default.
   echo "==> dada (fwd, rev; per-sample)"
+  metrics_for dada.fwd
   run_step dada.fwd "$BIN" dada "${filtFs[@]}" --error-model "$OUT/errF.json" \
-      --output-dir "$OUT/dada_fwd" --threads "$THREADS" ${backend_arg[@]+"${backend_arg[@]}"} ${screen_arg[@]+"${screen_arg[@]}"} ${end_arg[@]+"${end_arg[@]}"} ${verbose_arg[@]+"${verbose_arg[@]}"}
+      --output-dir "$OUT/dada_fwd" --threads "$THREADS" ${backend_arg[@]+"${backend_arg[@]}"} ${screen_arg[@]+"${screen_arg[@]}"} ${end_arg[@]+"${end_arg[@]}"} ${verbose_arg[@]+"${verbose_arg[@]}"} ${metrics_arg[@]+"${metrics_arg[@]}"}
+  metrics_for dada.rev
   run_step dada.rev "$BIN" dada "${filtRs[@]}" --error-model "$OUT/errR.json" \
-      --output-dir "$OUT/dada_rev" --threads "$THREADS" ${backend_arg[@]+"${backend_arg[@]}"} ${screen_arg[@]+"${screen_arg[@]}"} ${end_arg[@]+"${end_arg[@]}"} ${verbose_arg[@]+"${verbose_arg[@]}"}
+      --output-dir "$OUT/dada_rev" --threads "$THREADS" ${backend_arg[@]+"${backend_arg[@]}"} ${screen_arg[@]+"${screen_arg[@]}"} ${end_arg[@]+"${end_arg[@]}"} ${verbose_arg[@]+"${verbose_arg[@]}"} ${metrics_arg[@]+"${metrics_arg[@]}"}
 fi
 
 echo "==> merge-pairs"

dev/unoise/README.md

@@ -0,0 +1,116 @@
+# Single-pass and UNOISE arm tools (#323)
+
+Development scripts for the experimental `--partition single-pass` and
+`--birth-rule unoise` (issue #323). They are not installed and not part of the
+CLI. The flags are documented in `docs/commands/dada.md`.
+
+| Script | Purpose |
+| --- | --- |
+| `run_arms.sh` | Run the arm grid on one dataset: learn once, then denoise, merge, remove chimeras and optionally score each arm against a truth set |
+| `vsearch_to_seqtab.py` | Convert the vsearch arm's ZOTUs and OTU table into the files the other arms produce |
+| `summarize_arms.py` | One row per arm: ASV counts, set identity and Bray-Curtis against a baseline arm, dada time, and `reference-eval` scores |
+| `phase_share.py` | Phase 0: the share of denoising time that the single pass could remove, from `--metrics-json` |
+
+Each script's header documents its inputs, outputs and environment.
+
+## The arms
+
+| Arm | Flags | `_indel` variant |
+| --- | --- | --- |
+| `divisive` | none (`dada` as released) | `--indel-rates` |
+| `sp` | `--partition single-pass` | `--indel-rates` |
+| `sp_reassign` | `--partition single-pass --reassign` | `--indel-rates` |
+| `unoise` | `--partition single-pass --birth-rule unoise` | none: d counts gaps by definition |
+| `vsearch` | UNOISE3 as vsearch runs it (external control) | none |
+
+The `unoise` arm uses UNOISE's rule to pick centres and counts every absorbed
+member's reads (`--omega-c` does not apply), but keeps DADA2's merging, chimera
+removal and table. The `vsearch` arm is the UNOISE3 workflow end to end: merge
+pairs, filter at `--fastq_maxee 1` (`VSEARCH_MAXEE`), dereplicate,
+`--cluster_unoise`, `--uchime3_denovo`, then map every read, filtered or not,
+to the ZOTUs at 97% identity for the table. It uses no error model. Its
+`--minsize` follows the `unoise` arm's `--min-abund`. HiFi quality scores go
+up to 93, so every vsearch step reads FASTQ with `--fastq_qmax 93`.
+
+On the pooled MiSeq SOP forward reads, the two arms differ in three separable
+steps:
+
+| Step | `unoise` | `vsearch` |
+| --- | --- | --- |
+| Centres, `VSEARCH_MAXEE` off | 449 | 449, identical |
+| Centres, `--fastq_maxee 1` | 449 | 427, all among the 449 |
+| After chimera removal | 284 (DADA2 consensus) | 387 (UCHIME3); our 284 are all among UCHIME3's 406 at maxEE off |
+| Reads in the table | 127,286 | 139,114 (97% mapping) |
+
+Compare `sp_indel` with `unoise` to test the learned error model against the
+fixed rule with indels counted in both (`summarize_arms.py --pair
+sp_indel:unoise`). The indel λ counts indel events and the UNOISE rule counts
+gap columns, so the two still count indels differently; see the command page.
+
+One `learn-errors` run serves every arm. Learning with `--indel-rates-out`
+leaves the substitution model byte-identical, so the arms with and without
+indels share it.
+
+`run_arms.sh` checks that each arm's output JSON records the partition, birth
+rule, `--reassign`, indel λ and `--min-abund` it was asked for, and stops if
+not. Arms that silently ran alike would look like a null result.
+
+## `--min-abund`
+
+UNOISE needs an abundance cutoff (vsearch `--cluster_unoise` defaults to
+`--minsize 8`); DADA2 defaults to 1. Run each dataset twice, sharing one
+`LEARN_DIR`:
+
+- `MIN_ABUND` unset: each arm at its own default, `unoise` at
+  `UNOISE_MIN_ABUND` (8). This compares the methods as each would be used.
+- `MIN_ABUND=8`: every arm at 8. This isolates the partition and birth rule
+  from the cutoff.
+
+## Running a dataset
+
+The manifest is tab-separated with no header: `sample  fwd  [rev]`. Inputs must
+already be trimmed and filtered, for example the `filtered/` directory of a
+`dev/concordance/run_illumina.sh` or `run_pacbio.sh` run. Row order selects the
+training set when `--nbases` truncates.
+
+The `vsearch` arm needs vsearch (`VSEARCH`, default `vsearch` on `PATH`); the
+script checks before any work. Drop it from `ARMS` where vsearch is not
+installed. Every DADA2-side arm removes chimeras with `remove-bimera-denovo
+--method consensus`; `CHIMERA_EXTRA` passes it flags, such as the
+`--min-fold-parent-over-abundance 3.5` recommended for PacBio. The vsearch arm
+uses UCHIME3 regardless. For ITS2, `VSEARCH_MERGE_EXTRA="--fastq_allowmergestagger"` may be
+needed when reads overrun each other.
+
+```bash
+BIN=target/release/dada2-rs
+
+# MiSeq V4
+THREADS=24 bash dev/unoise/run_arms.sh $BIN out/miseq miseq.tsv
+THREADS=24 MIN_ABUND=8 LEARN_DIR=out/miseq/learn \
+  bash dev/unoise/run_arms.sh $BIN out/miseq_ma8 miseq.tsv
+
+# NovaSeq, binned quality scores (bins from `dada2-rs summary --report`)
+LEARN_EXTRA="--errfun binned-qual --binned-quals 2,11,25,37" \
+  THREADS=24 bash dev/unoise/run_arms.sh $BIN out/novaseq novaseq.tsv
+
+# PacBio HiFi, MSA-1003 (k=7, the PacBio recommendation)
+LEARN_EXTRA="--errfun pacbio" EXTRA="--band 32 --kmer-size 7" \
+  CHIMERA_EXTRA="--min-fold-parent-over-abundance 3.5" \
+  REFERENCE=msa1003_truth.fa THREADS=24 \
+  bash dev/unoise/run_arms.sh $BIN out/msa1003 msa1003.tsv
+
+# ITS2: no fixed-length truncation in the prep (cutadapt only)
+THREADS=24 bash dev/unoise/run_arms.sh $BIN out/its2 its2.tsv
+
+python3 dev/unoise/summarize_arms.py out/msa1003 --pair sp_indel:unoise --tsv out/msa1003.tsv
+```
+
+Timing: `summarize_arms.py` reports each dada step's wall time (`wall_s`) and
+peak memory (`peak_mb`) from its `--metrics-json`, which needs no external
+timer. The vsearch arm's come from `time_fwd.txt` and cover `--cluster_unoise`
+only. Where `/usr/bin/time` is missing, set `TIME_BIN` to a GNU time; without
+one, bash's `time` keyword records wall and CPU seconds but no memory.
+
+A run stopped part-way resumes: finished arms (those with `<arm>/.done`) are
+skipped unless `FORCE=1`, and learning is skipped when `LEARN_DIR` already holds
+the error models and indel rates.

dev/unoise/phase_share.py

@@ -0,0 +1,60 @@
+#!/usr/bin/env python3
+"""Share of denoising time a single-pass partition could remove (#323).
+
+Reads --metrics-json files from `dada`, `dada-pooled` or `dada-pseudo` and sums
+the per-run phase timers. Single pass keeps every comparison and drops the
+shuffle, the bud search and the p-value refresh, so the sum of those three is
+the most it can save. Shares are of the summed phases: `loop_total` excludes
+the initial compare, so it is not used as the denominator.
+
+A share is only meaningful with the thread count it was measured at: the
+removable phases are serial, so their share grows with threads.
+
+Usage: dev/unoise/phase_share.py <metrics.json> [...]
+"""
+import json
+import sys
+
+PHASES = ["compare", "shuffle", "bud", "p_update", "index_add"]
+REMOVABLE = ["shuffle", "bud", "p_update"]
+
+
+def runs(o):
+    """Yield every object carrying a phases block (one per run_dada call)."""
+    if isinstance(o, dict):
+        if isinstance(o.get("phases"), dict) and "shuffle" in o["phases"]:
+            yield o
+        for v in o.values():
+            yield from runs(v)
+    elif isinstance(o, list):
+        for v in o:
+            yield from runs(v)
+
+
+def main(paths):
+    if not paths:
+        sys.exit(__doc__)
+    for path in paths:
+        with open(path) as fh:
+            rows = list(runs(json.load(fh)))
+        if not rows:
+            print(f"{path}: no phases block found")
+            continue
+        t = {k: sum(r["phases"].get(k) or 0.0 for r in rows) for k in PHASES}
+        total = sum(t.values())
+        rem = sum(t[k] for k in REMOVABLE)
+        threads = rows[0].get("run", {}).get("threads", "?")
+
+        def pct(x):
+            return f"{x / total:.1%}" if total else "n/a"
+
+        print(
+            f"{path}: runs={len(rows)} threads={threads} total={total:.2f}s "
+            f"compare={t['compare']:.2f}s ({pct(t['compare'])}) "
+            f"shuffle={t['shuffle']:.2f}s bud={t['bud']:.2f}s "
+            f"p_update={t['p_update']:.2f}s -> removable {rem:.2f}s ({pct(rem)})"
+        )
+
+
+if __name__ == "__main__":
+    main(sys.argv[1:])