feat(dada): add --failed-uniques diagnostic (issue #60)

docs/diagnostics.md

@@ -4,6 +4,69 @@ This page documents **experimental diagnostic tooling** built into `dada2-rs`, a
 
 ---
 
+## `--failed-uniques`: which sequences failed to denoise
+
+A recurring question (e.g. [benjjneb/dada2#1899](https://github.com/benjjneb/dada2/issues/1899)) is *which* unique sequences fail to denoise — particularly for high-diversity samples (soil, seawater) that shed substantial reads through denoising. DADA2's answer lives in the per-unique `map`: a unique whose final abundance p-value falls below `OMEGA_C`, and that is not corrected to any cluster center, gets a `null` map entry (R's `$map == NA`). Those nulls are the uniques that failed to denoise, and they can be traced back to the reads they cost.
+
+`dada2-rs` exposes this directly. The `--failed-uniques <FILE>` flag — available on `dada`, `dada-pseudo`, and `dada-pooled` — writes a TSV of the dropped uniques and their abundances, so you don't have to hand-join the `map` to the derep uniques.
+
+### Output
+
+A tidy long-format TSV with a header, one row per failed unique per sample it appears in:
+
+```
+sequence	sample	reads
+TACGAAGG…AACA	soil_A	7
+TACGGAGG…AATC	soil_A	3
+TACGGAGG…AAAC	soil_B	2
+```
+
+| Column | Meaning |
+|---|---|
+| `sequence` | The unique that failed to denoise (`map == null`). |
+| `sample` | Sample the unique (and its reads) belongs to. |
+| `reads` | Reads this unique contributed in that sample (i.e. reads lost to the failure). |
+
+Rows are sorted by `sample`, then descending `reads`, then `sequence`, so output is deterministic regardless of how many samples were denoised concurrently. Sum the `reads` column (optionally per sample) to get the total reads lost to failed denoising.
+
+### Semantics follow the subcommand
+
+The notion of "failed" matches how each subcommand denoises:
+
+| Subcommand | "Failed" decided | `reads` is |
+|---|---|---|
+| `dada`, `dada-pseudo` | **per sample** (`map == null`; for pseudo, in the round-2 pass) | the unique's in-sample abundance |
+| `dada-pooled` | **globally** — pooled denoising runs once on the merged unique table, so failure is a property of the merged index (`result.map == null`) | the failed merged unique's read count *in that sample* (one row per sample it appears in) |
+
+!!! note "Per-sample `null` is unambiguous"
+
+    In the pooled per-sample output JSON, a `map` entry could in principle be
+    `null` either because the unique failed denoising or because its cluster has
+    zero reads in that sample. The latter never actually happens for a unique
+    that is present in the sample (a present unique always carries ≥1 read into
+    its cluster), so a per-sample `map == null` always means a genuine denoising
+    failure. The `--failed-uniques` TSV reports only real failures.
+
+### Invocation
+
+```bash
+# single sample
+dada2-rs dada sampleA.fastq.gz --error-model err.json \
+    -o sampleA.dada.json --failed-uniques sampleA.failed.tsv
+
+# multi-sample (per-sample failures, one combined TSV with a sample column)
+dada2-rs dada *.fastq.gz --error-model err.json \
+    --output-dir dada_out/ --failed-uniques failed.tsv
+
+# pooled (global failures, expanded per sample)
+dada2-rs dada-pooled *.fastq.gz --error-model err.json \
+    -o pooled_out/ --failed-uniques failed.tsv
+```
+
+To trace a failed unique's divergence from the nearest surviving ASV (is it a distant artifact, or a real low-abundance variant the abundance test shed?), feed the same `dada` run to [`kdist-calibrate --from-dada`](#3-post-inference-mode-from-dada), whose `failed` class is exactly this population.
+
+---
+
 ## `kdist-calibrate`: k-mer screen, `KDIST_CUTOFF`, `BAND_SIZE`
 
 Currently the k-mer **screen** (`KDIST_CUTOFF`) and the alignment **band size** 

src/cli.rs

@@ -344,6 +344,14 @@ pub enum Commands {
         #[arg(long)]
         output_dir: Option<PathBuf>,
 
+        /// Write a TSV of uniques that failed to denoise (final p-value < omega_c,
+        /// `map == null`) to this file. Tidy long format with a header:
+        /// `sequence<TAB>sample<TAB>reads`, one row per failed unique per sample
+        /// it appears in. Answers "what failed to denoise and how many reads did
+        /// it cost?" without hand-joining the `map` to the derep uniques.
+        #[arg(long)]
+        failed_uniques: Option<PathBuf>,
+
         /// Output compact (minified) JSON instead of pretty-printed
         #[arg(long)]
         compact: bool,
@@ -499,6 +507,15 @@ pub enum Commands {
         #[arg(long)]
         no_kmer_screen: Option<bool>,
 
+        /// Write a TSV of uniques that failed to denoise to this file. Tidy long
+        /// format with a header: `sequence<TAB>sample<TAB>reads`. Because pooled
+        /// denoising runs once on the merged unique table, "failed" is a global
+        /// property (`result.map == null` on the merged index); for each failed
+        /// merged unique a row is emitted per sample it appears in, carrying that
+        /// sample's read count.
+        #[arg(long)]
+        failed_uniques: Option<PathBuf>,
+
         /// Output compact (minified) JSON instead of pretty-printed
         #[arg(long)]
         compact: bool,
@@ -684,6 +701,14 @@ pub enum Commands {
         #[arg(long)]
         no_kmer_screen: Option<bool>,
 
+        /// Write a TSV of uniques that failed to denoise (per-sample `map == null`)
+        /// to this file. Tidy long format with a header:
+        /// `sequence<TAB>sample<TAB>reads`, one row per failed unique per sample.
+        /// Failures are decided in the per-sample round-2 pass, so this is a
+        /// per-sample signal (same semantics as `dada`).
+        #[arg(long)]
+        failed_uniques: Option<PathBuf>,
+
         /// Output compact (minified) JSON instead of pretty-printed
         #[arg(long)]
         compact: bool,

src/failed_uniques.rs

@@ -0,0 +1,50 @@
+//! Failed-to-denoise unique diagnostic (issue #60).
+//!
+//! DADA2 reports which input uniques fell out of denoising via the per-unique
+//! `map`: a unique whose final abundance p-value is below `omega_c` and that was
+//! not corrected to any cluster center gets a `null` (`None`) map entry — Ben
+//! Callahan's guidance on benjjneb/dada2#1899 is that these NAs are how you trace
+//! reads that failed to denoise. This module turns that signal into a direct
+//! artifact: a tidy long-format TSV of the dropped uniques and the reads they
+//! cost, so users don't have to hand-join the dada `map` to the derep uniques.
+//!
+//! Format: a header line `sequence<TAB>sample<TAB>reads`, then one row per failed
+//! unique per sample it appears in. For per-sample modes (`dada`, `dada-pseudo`)
+//! each row's `reads` is the unique's in-sample abundance; for `dada-pooled` the
+//! "failed" decision is global (made once on the merged unique table) and a row
+//! is emitted per sample the failed merged unique appears in.
+
+use std::io::{self, Write};
+use std::path::Path;
+
+/// One failed-to-denoise unique, scoped to a single sample.
+pub struct Row {
+    pub sequence: String,
+    pub sample: String,
+    pub reads: u32,
+}
+
+/// Write `rows` as the tidy long-format TSV (with header) to `path`, creating
+/// parent directories as needed. Rows are sorted by sample, then descending
+/// reads, then sequence, so output is deterministic regardless of the order in
+/// which concurrently denoised samples contributed them.
+pub fn write_tsv(path: &Path, mut rows: Vec<Row>) -> io::Result<usize> {
+    if let Some(parent) = path.parent()
+        && !parent.as_os_str().is_empty()
+    {
+        std::fs::create_dir_all(parent)?;
+    }
+    rows.sort_by(|a, b| {
+        a.sample
+            .cmp(&b.sample)
+            .then(b.reads.cmp(&a.reads))
+            .then(a.sequence.cmp(&b.sequence))
+    });
+    let mut w = io::BufWriter::new(std::fs::File::create(path)?);
+    writeln!(w, "sequence\tsample\treads")?;
+    for r in &rows {
+        writeln!(w, "{}\t{}\t{}", r.sequence, r.sample, r.reads)?;
+    }
+    w.flush()?;
+    Ok(rows.len())
+}