From 3ac37f35a8b551d19208515f8e096b0ab15e3fe6 Mon Sep 17 00:00:00 2001 From: Chris Fields Date: Mon, 22 Jun 2026 13:49:53 -0500 Subject: [PATCH] feat(dada): add --failed-uniques diagnostic (issue #60) Emit the unique sequences that failed to denoise (final p-value < omega_c, `map == null`) as a tidy long-format TSV (`sequencesamplereads`, one row per failed unique per sample), so users can answer "what failed to denoise and how many reads did it cost?" without hand-joining the dada `map` to the derep uniques. Motivated by benjjneb/dada2#1899. Available on `dada`, `dada-pseudo`, and `dada-pooled`, with semantics that follow the subcommand: `dada`/`dada-pseudo` report per-sample failures (`map == null`, decided in pseudo's round-2 pass); `dada-pooled` reports global failures (`result.map == null` on the merged unique table), expanded to one row per sample the failed merged unique appears in. The pooled per-sample-`map`-null ambiguity noted in the issue is unreachable: a unique present in a sample always carries >=1 read into its cluster, so a per-sample `map == null` always means a genuine denoising failure. - new module src/failed_uniques.rs (Row + write_tsv, deterministic sort) - denoise_and_serialize returns (json, Vec) gated by collect_failed - integration test: dada_failed_uniques_matches_map_nulls - docs/diagnostics.md: new --failed-uniques section Co-Authored-By: Claude Opus 4.8 --- docs/diagnostics.md | 63 +++++++++++++++++++++++ src/cli.rs | 25 +++++++++ src/failed_uniques.rs | 50 ++++++++++++++++++ src/main.rs | 117 ++++++++++++++++++++++++++++++++++++++++-- tests/dada_pseudo.rs | 46 +++++++++++++++++ 5 files changed, 296 insertions(+), 5 deletions(-) create mode 100644 src/failed_uniques.rs diff --git a/docs/diagnostics.md b/docs/diagnostics.md index 0eaec5c..abd4371 100644 --- a/docs/diagnostics.md +++ b/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 ` 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** diff --git a/src/cli.rs b/src/cli.rs index d68b218..69cf31c 100644 --- a/src/cli.rs +++ b/src/cli.rs @@ -344,6 +344,14 @@ pub enum Commands { #[arg(long)] output_dir: Option, + /// 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: + /// `sequencesamplereads`, 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, + /// 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, + /// Write a TSV of uniques that failed to denoise to this file. Tidy long + /// format with a header: `sequencesamplereads`. 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, + /// 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, + /// Write a TSV of uniques that failed to denoise (per-sample `map == null`) + /// to this file. Tidy long format with a header: + /// `sequencesamplereads`, 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, + /// Output compact (minified) JSON instead of pretty-printed #[arg(long)] compact: bool, diff --git a/src/failed_uniques.rs b/src/failed_uniques.rs new file mode 100644 index 0000000..07b298f --- /dev/null +++ b/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 `sequencesamplereads`, 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) -> io::Result { + 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()) +} diff --git a/src/main.rs b/src/main.rs index f87c68e..ebfe1ac 100644 --- a/src/main.rs +++ b/src/main.rs @@ -17,6 +17,7 @@ mod derep; mod error; mod error_models; mod evaluate; +mod failed_uniques; mod filter; mod filter_trim; mod kdist_calibrate; @@ -453,6 +454,7 @@ fn main() -> io::Result<()> { trace_min_abund, output, output_dir, + failed_uniques: failed_uniques_path, compact, verbose, } => { @@ -560,6 +562,9 @@ fn main() -> io::Result<()> { input.len() ); } + let collect_failed = failed_uniques_path.is_some(); + let failed_rows: std::sync::Mutex> = + std::sync::Mutex::new(Vec::new()); for_each_sample_concurrent(input.len(), jobs, threads, |i, sub_pool| { let path = &input[i]; let (derep, json_sample) = @@ -598,7 +603,7 @@ fn main() -> io::Result<()> { } } let sample = json_sample.unwrap_or_else(|| fastq_stem(path)); - let json = denoise_and_serialize( + let (json, failed) = denoise_and_serialize( "dada", &sample, &file_basename(path), @@ -607,8 +612,12 @@ fn main() -> io::Result<()> { &resolved.run, sub_pool, compact, + collect_failed, verbose, )?; + if collect_failed { + failed_rows.lock().unwrap().extend(failed); + } let out_path = output_dir.join(format!("{sample}.json")); std::fs::write(&out_path, &json)?; if verbose { @@ -616,6 +625,16 @@ fn main() -> io::Result<()> { } Ok(()) })?; + if let Some(ref fu_path) = failed_uniques_path { + let rows = failed_rows.into_inner().unwrap(); + let n = failed_uniques::write_tsv(fu_path, rows)?; + if verbose { + eprintln!( + "[dada] wrote {n} failed-unique row(s) to {}", + fu_path.display() + ); + } + } return Ok(()); } @@ -809,6 +828,29 @@ fn main() -> io::Result<()> { let sample = sample_name .or(json_sample) .unwrap_or_else(|| fastq_stem(input)); + + // ---- Optional failed-to-denoise unique diagnostic (issue #60) ---- + if let Some(ref fu_path) = failed_uniques_path { + let rows: Vec = result + .map + .iter() + .enumerate() + .filter(|(_, m)| m.is_none()) + .map(|(i, _)| failed_uniques::Row { + sequence: raw_inputs[i].seq.clone(), + sample: sample.clone(), + reads: raw_inputs[i].abundance, + }) + .collect(); + let n = failed_uniques::write_tsv(fu_path, rows)?; + if verbose { + eprintln!( + "[dada] wrote {n} failed-unique row(s) to {}", + fu_path.display() + ); + } + } + let total_reads: u32 = result.clusters.iter().map(|c| c.reads).sum(); let asvs: Vec = result @@ -920,6 +962,7 @@ fn main() -> io::Result<()> { kdist_cutoff, kmer_size, no_kmer_screen, + failed_uniques: failed_uniques_path, compact, verbose, } => { @@ -1182,12 +1225,25 @@ fn main() -> io::Result<()> { map: Vec>, } + // Failed-to-denoise uniques (issue #60). Pooled denoising runs once on + // the merged unique table, so "failed" is a global property + // (`result.map[mu] == None`); rows are collected per sample the failed + // merged unique appears in, carrying that sample's read count. + let collect_failed = failed_uniques_path.is_some(); + let mut failed_rows: Vec = Vec::new(); + for (s, sample_name) in sample_names.iter().enumerate() { // Sum per-cluster reads for this sample by walking its local uniques. let mut cluster_reads: Vec = vec![0u32; result.clusters.len()]; for (lu, &mu) in local_to_merged[s].iter().enumerate() { if let Some(c) = result.map[mu] { cluster_reads[c] += sample_unique_counts[s][lu]; + } else if collect_failed { + failed_rows.push(failed_uniques::Row { + sequence: raw_inputs[mu].seq.clone(), + sample: sample_name.clone(), + reads: sample_unique_counts[s][lu], + }); } } @@ -1240,6 +1296,16 @@ fn main() -> io::Result<()> { ); } } + + if let Some(ref fu_path) = failed_uniques_path { + let n = failed_uniques::write_tsv(fu_path, failed_rows)?; + if verbose { + eprintln!( + "[dada-pooled] wrote {n} failed-unique row(s) to {}", + fu_path.display() + ); + } + } if verbose { let t_output = t_output.elapsed(); let total = (t_derep + t_merge + t_dada + t_output) @@ -1294,6 +1360,7 @@ fn main() -> io::Result<()> { kdist_cutoff, kmer_size, no_kmer_screen, + failed_uniques: failed_uniques_path, compact, verbose, } => { @@ -1529,6 +1596,8 @@ fn main() -> io::Result<()> { if low_memory { ", streaming" } else { "" }, ); } + let collect_failed = failed_uniques_path.is_some(); + let failed_rows: Mutex> = Mutex::new(Vec::new()); if !low_memory { // Cached: mark priors serially (cheap; scoped so the &mut borrow // is released), then denoise concurrently with shared read access. @@ -1547,7 +1616,7 @@ fn main() -> io::Result<()> { let sample_raws = sample_raws_opt.as_ref().unwrap(); for_each_sample_concurrent(n_samples, jobs, threads, |s, sub_pool| { let sample_name = &sample_names[s]; - let json = denoise_and_serialize( + let (json, failed) = denoise_and_serialize( "dada-pseudo", sample_name, &file_basename(&input[s]), @@ -1556,8 +1625,12 @@ fn main() -> io::Result<()> { &resolved.run, sub_pool, compact, + collect_failed, verbose, )?; + if collect_failed { + failed_rows.lock().unwrap().extend(failed); + } let out_path = output_dir.join(format!("{sample_name}.json")); std::fs::write(&out_path, &json)?; if verbose { @@ -1579,7 +1652,7 @@ fn main() -> io::Result<()> { raws.len(), ); } - let json = denoise_and_serialize( + let (json, failed) = denoise_and_serialize( "dada-pseudo", sample_name, &file_basename(&input[s]), @@ -1588,8 +1661,12 @@ fn main() -> io::Result<()> { &resolved.run, sub_pool, compact, + collect_failed, verbose, )?; + if collect_failed { + failed_rows.lock().unwrap().extend(failed); + } let out_path = output_dir.join(format!("{sample_name}.json")); std::fs::write(&out_path, &json)?; if verbose { @@ -1598,6 +1675,16 @@ fn main() -> io::Result<()> { Ok(()) })?; } + if let Some(ref fu_path) = failed_uniques_path { + let rows = failed_rows.into_inner().unwrap(); + let n = failed_uniques::write_tsv(fu_path, rows)?; + if verbose { + eprintln!( + "[dada-pseudo] wrote {n} failed-unique row(s) to {}", + fu_path.display() + ); + } + } } Commands::MergePairs { @@ -4013,6 +4100,9 @@ fn to_json(value: &T, compact: bool) -> io::Result { } #[allow(clippy::too_many_arguments)] +/// Denoise one sample and serialize its dada JSON. When `collect_failed` is set, +/// also returns the uniques that failed to denoise (`map == null`) as +/// [`failed_uniques::Row`]s tagged with `sample`; otherwise the row vec is empty. fn denoise_and_serialize( tag: &'static str, sample: &str, @@ -4022,8 +4112,9 @@ fn denoise_and_serialize( run_params: &DadaRunParams, pool: &rayon::ThreadPool, compact: bool, + collect_failed: bool, verbose: bool, -) -> io::Result { +) -> io::Result<(String, Vec)> { #[derive(Serialize)] struct DadaOutput { sample: String, @@ -4061,6 +4152,22 @@ fn denoise_and_serialize( let mut run_params = *run_params; run_params.n_prior = raw_inputs.iter().filter(|r| r.prior).count(); + let failed = if collect_failed { + result + .map + .iter() + .enumerate() + .filter(|(_, m)| m.is_none()) + .map(|(i, _)| failed_uniques::Row { + sequence: raw_inputs[i].seq.clone(), + sample: sample.to_string(), + reads: raw_inputs[i].abundance, + }) + .collect() + } else { + Vec::new() + }; + let out = DadaOutput { sample: sample.to_string(), input_file: input_file.to_string(), @@ -4075,7 +4182,7 @@ fn denoise_and_serialize( map: result.map, }; - to_json(&Tagged::new(tag, out), compact) + Ok((to_json(&Tagged::new(tag, out), compact)?, failed)) } /// Build a [`derep::Derep`] for `dada` / `dada-pooled` from either a FASTQ file diff --git a/tests/dada_pseudo.rs b/tests/dada_pseudo.rs index 2a33402..c289593 100644 --- a/tests/dada_pseudo.rs +++ b/tests/dada_pseudo.rs @@ -772,3 +772,49 @@ fn dada_from_fastq_matches_dada_from_derep_json() { "dada output from derep JSON differs from dada output from FASTQ", ); } + +/// `--failed-uniques` (issue #60) must emit a header plus exactly one row per +/// `map == null` unique (a unique that failed to denoise), with that unique's +/// sequence and in-sample abundance. The single-sample dada `map` is the clean +/// reference signal, so the TSV row count must equal the JSON null count and +/// every TSV sequence must be a `map == null` input unique. +#[test] +fn dada_failed_uniques_matches_map_nulls() { + let dir = scratch("failed_uniques"); + let err = shared_err_model(); + let s1 = fixture("sam1F.fastq.gz"); + + let out_json = dir.join("d.json"); + let fu_tsv = dir.join("failed.tsv"); + run(&[ + "dada", + s1.to_str().unwrap(), + "--error-model", + err.to_str().unwrap(), + "-o", + out_json.to_str().unwrap(), + "--failed-uniques", + fu_tsv.to_str().unwrap(), + ]); + + let v: serde_json::Value = serde_json::from_slice(&std::fs::read(&out_json).unwrap()).unwrap(); + let map = v["map"].as_array().unwrap(); + let null_count = map.iter().filter(|m| m.is_null()).count(); + assert!(null_count > 0, "fixture should have some failed uniques"); + + let tsv = std::fs::read_to_string(&fu_tsv).unwrap(); + let mut lines = tsv.lines(); + assert_eq!( + lines.next().unwrap(), + "sequence\tsample\treads", + "TSV must start with the header row", + ); + let rows: Vec<&str> = lines.collect(); + assert_eq!(rows.len(), null_count, "one TSV row per map==null unique",); + for row in rows { + let cols: Vec<&str> = row.split('\t').collect(); + assert_eq!(cols.len(), 3, "row must have sequence/sample/reads"); + assert_eq!(cols[1], "sam1F", "sample column"); + assert!(cols[2].parse::().unwrap() >= 1, "reads >= 1"); + } +}