Experiment: single-pass denoising in order of abundance (UNOISE-like) with DADA2's error model, across sequencing technologies

Origin

#303 (wetSpring) prompted this. It is a GPU denoiser that adds every significant sequence as a new center in one pass. We aren't pursuing that design: it skips alignment, re-learns errors inside each sample, and seems to borrows from several methods without publishing results checking agreement with them. Its starting idea, denoising without DADA2's iterative splitting, has a cleaner published form in UNOISE3 (Edgar 2016), and that's what this issue tests.

Background

UNOISE3 makes one pass over the unique sequences, most abundant first. Each sequence becomes a new center unless it is close to an existing, more abundant center and its abundance is low enough relative to that center. "Low enough" is a fixed rule, a maximum abundance ratio β(d) = 1/2^(αd+1), tuned on Illumina data and with no per-run error model (UNOISE2 preprint, which also reports results on ITS).

The author notes that UNOISE "was designed and tested on Illumina reads" and does less well on 454 and Ion Torrent data. The main problem is indel errors from miscalled homopolymer lengths, which "cause two issues: incorrect alignments and inflated abundance of incorrect reads" (USEARCH FAQ). PacBio HiFi data has the same kind of error, albeit at much lower rates.

Deblur has a different limit: it assumes reads of one fixed length. Tedersoo et al. 2022 (doi:10.1111/mec.16460) state that it can't be used for ITS reads because their length varies. That statement is quoted in a QIIME 2 forum thread on fungal ITS workflows.

DADA2 swaps the fixed rule for a learned error model and a Poisson abundance test. It then splits clusters one new center at a time, re-assigning members and recomputing p-values after each. That splitting loop is a large part of the cost, and it can't be run in parallel within a sample.

Question

Does a single pass in order of abundance, with DADA2's learned error model and abundance p-value replacing UNOISE3's fixed ratio, match dada's accuracy across technologies compared to DADA2's iterative approach? And what does it save?

Arms

All arms use the same derep input, learn-errors output, global NW alignment, k-mer screen, and OMEGA_A.

dada (current): splits one new center at a time.
Single pass, error model: process uniques in decreasing order of abundance. For each one, compute lambda against each existing center that passes the screen. Take the best (largest-lambda) center and apply DADA2's p-value test. If p < OMEGA_A, the sequence becomes a new center; otherwise it joins that center. No re-assignment afterwards.
Arm 2 plus one re-assignment pass: after the pass, move each member to its most likely center once. This separates what the error model contributes from what DADA2's re-assignment contributes.
Single pass, UNOISE3 rule: the same loop as arm 2, but with the β(d) ratio rule instead of the error model. This is the control for the error model. Compare it with vsearch --cluster_unoise as a check that our loop reproduces UNOISE3.

Add arms 2–4 as an experimental switch on dada, not a new subcommand.

Data

Each dataset tests a different way the fixed rule could fail:

Illumina MiSeq 16S V4: UNOISE3's home turf; all arms should agree here.
Illumina NovaSeq/NextSeq with binned qualities: a newer chemistry with a different error profile.
PacBio HiFi full-length 16S (MSA-1003 truth set): long reads, homopolymer indels, and a different error profile.
ITS2: reads of varying length with more variation between real sequences.

Measures

Accuracy against a known mock community with reference-eval, wherever a truth set exists.
ASV set identity against dada on the same input: shared, gained, and lost ASVs, with the abundance moved. Don't use ASV counts alone; they hide churn.
Wall time, alignments performed, and how well the work parallelises. A single pass compares each unique with every existing center, roughly the same number of alignments as dada, so any savings will come from dropping the splitting loop, not from fewer alignments.

Predictions, stated so they can be disproved

On MiSeq V4, arms 2–4 come close to dada.
Arm 4 drifts from dada and from the truth on PacBio and binned-quality data, while arm 2 does not. That would confirm the error model, not the single pass, is what makes the method work across technologies.
Arm 2's differences from dada are concentrated in low-abundance sequences that sit between two centers, and arm 3 removes most of them.

If arm 2 holds up, it could be a cheaper mode for very large or pooled runs (see #106), and its comparisons could be batched for #17. GPU work is out of scope here.

Plan

Approach: a switchable backend

run_dada (src/dada.rs) returns the cluster structure B, and everything downstream reads only that: JSON output, --failed-uniques, chimera removal, reference-eval. A single-pass backend that builds the same B gets all of that unchanged. It also works in dada-pooled without extra effort, since that mode calls the same function.

Single pass reuses DADA2's comparison step instead of adding a new one:

Walk the unique sequences once, most abundant first.
When a unique becomes a center, run the existing b_compare_parallel from that center against all the uniques, exactly as dada does after each split. The existing greedy option already skips uniques more abundant than the center.
When the walk reaches unique j, its comparisons to every earlier center are already stored. Deciding its fate is then a lookup: take its best center, run the birth test, and either create a center or assign it to that one.
The shuffle (b_shuffle_converge), the search for the next sequence to split off (b_bud_incremental), and the p-value refresh after each split (b_p_update) all drop out.

The k-mer screen, alignment, lambda and p-value code stay identical to dada, which is what makes the arms comparable.

Switches

Two independent switches, following --align-backend and --screen-backend. Both default to today's behaviour, are written to the output JSON only when not the default, and are marked experimental:

--partition divisive|single-pass: divisive is today's dada and stays the default.
--birth-rule pvalue|unoise: the test that decides whether a unique becomes a new center. pvalue is DADA2's Poisson test with the learned error model. unoise is UNOISE3's fixed ratio rule β(d), with d the mismatch count from the stored alignment. It's only valid with single-pass.
An option for one final reassignment pass (arm 3), implemented as a single shuffle call.

learn-errors stays on divisive and doesn't get the flag, so every arm uses the same error model.

Phases

Estimate the possible speed gain. --metrics-json already times the shuffle, split search and p-value refresh separately (shuffle, bud, p_update). The comparisons are unchanged in single pass, so those three timings are the most it can save. If they're small, this issue is purely about accuracy.
Add the switch with no behaviour change. Add the enums and plumbing, and record the choice in the JSON. Add a check that the default output is byte-identical to today's.
Single pass with the p-value rule (arm 2), with unit tests on small inputs where processing order changes the result.
The final reassignment option (arm 3).
The UNOISE rule (arm 4), checked against vsearch --cluster_unoise on MiSeq V4 before it's used as a control.
Run the four arms on the four datasets on the cluster, scored with reference-eval and ASV set identity. The dev/compare_arms.py harness designed for Re-assess the LOESS surface default at scale, and make the benchmark harness error-model-aware #205 does the same job.
Write up the result in docs/findings/.

Decisions to settle before phase 2

Center abundance in the birth test. The center's own unique abundance doesn't change during the run, so the result depends less on processing order (UNOISE-like). The cluster's running read total is closer to DADA2 but depends on which members joined first. Leaning toward the center's own abundance.
Choosing a unique's best center. Highest expected reads (lambda × center abundance), as the DADA2 shuffle assigns, or highest lambda. Leaning toward expected reads.
DADA2's other birth gates. Keep the MIN_FOLD, MIN_HAMMING and singleton checks on in arm 2, so that only the partitioning differs between arms.
dada-pseudo with priors. Out of scope for the first pass.

Phase 0 result: the most single pass could save

Cluster runs, pooled, with --metrics-json (branch single-pass-323, #324). "Removable" is shuffle + bud + p_update, the three phases single pass drops. Its share is of the summed phases.

Run	Threads	Compare	Removable	Share
MiSeq 362, R1	24	107.1 s	7.2 s	6.3%
MiSeq 362, R1	96	39.3 s	7.3 s	15.7%
MiSeq 362, R2	24	65.9 s	12.2 s	15.5%
MiSeq 362, R2	96	27.3 s	12.5 s	31.2%
PacBio 95, k=7	24	486.5 s	11.4 s	2.3%
PacBio 95, k=7	96	171.6 s	11.5 s	6.3%
PacBio 95, k=5	24	3482.4 s	12.3 s	0.4%
PacBio 95, k=5	96	1013.3 s	9.7 s	0.9%

The removable phases take a fixed 7–12 s per run regardless of thread count, because they run single-threaded. Only compare scales with threads, so the share grows with threads.
The best case is MiSeq R2 at 96 threads: at most 1.45× faster, or 12 s on a 40 s run. PacBio tops out at 6% with k=7.
These are upper bounds. Compare time scales with the number of centers, so a single pass that produces more ASVs would spend more time comparing.

Conclusion: speed doesn't justify single pass. #323 continues as an accuracy question: does a single pass with DADA2's error model match dada across technologies, and how does it compare with UNOISE3's fixed rule? Any work on speed is better aimed at the shuffle and p-value refresh within the current algorithm.

Side result: on pooled PacBio, k=7 compares 6–7× faster than k=5 (3482 → 487 s at 24 threads, 1013 → 172 s at 96). run_pacbio.sh still defaults to PACBIO_KMER=5.

Next: phase 1 (switch and plumbing, no behaviour change), marked highly experimental alongside wfa2 and the indel work.

feat: experimental --partition and --birth-rule switches, no behaviou…

feat: experimental single-pass partition, p-value birth rule (#323)

feat: greedy locks in the single-pass partition (#323)

Phase 2: single pass with the p-value rule, plus greedy locks

On single-pass-323 (#324): 0f65157 (single pass) and 0b0f60c (greedy locks). Uses the choices settled above: expected reads are λ × the centre's own reads, the best centre is the one with the most expected reads, and the MIN_FOLD, MIN_HAMMING, MIN_ABUNDANCE and singleton gates apply. Greedy locks are applied to each newly born centre as dada does, for consistency with dada. Nothing is reassigned after the pass.

First look: MiSeq SOP, forward reads, 20 samples, 8 threads (local)

Pooled	Clusters	Alignments	Aligned pairs	Compare time
Divisive	518	9.36M	2.24M	5.94 s
Single pass, no locks	534	12.08M	2.96M	7.91 s
Single pass, with locks	534	9.65M	2.33M	6.65 s

ASV sets, single pass (with locks) against divisive:

	Shared	Only divisive	Only single pass
Pooled	512	6 ASVs, 133 reads	22 ASVs, 768 reads (median 25, max 216)
Per-sample	356	4 ASVs, 41 reads	3 ASVs, 61 reads

Without locks, single pass made 29% more alignments for 3% more clusters. With them it is 3% above divisive, in line with the cluster count.
The locks don't change the ASV set: the same 534 ASVs, with 52 reads moving between clusters.
Pooled single pass keeps a few more low-abundance variants. That's what you'd expect when the cluster's members can't add to the reads a variant is tested against.
The output is identical at 1 and 8 threads. Default (divisive) output remains byte-identical to main.

This is one small run that only checks single pass works. Accuracy will be judged against truth sets in phase 5.

Next: phase 3 (an optional reassignment pass after single pass) and phase 4 (the UNOISE3 rule as a control).

feat: --reassign, one DADA2 shuffle after the single pass (#323)

Plan update: an indel axis, and the UNOISE rule's definition

Why

UNOISE counts gaps by definition. The UNOISE2 preprint (Edgar 2016, doi:10.1101/081257) says: "Let d be the Levenshtein distance (number of differences including both substitutions and gaps) between M and C". DADA2's defaults count substitutions only, and our indel-aware λ (#291) is opt-in. Comparing UNOISE only against substitution-only arms would mix two differences, the birth rule and whether indels count. So the indel model becomes an axis of its own.

Arm	Indel model off	Indel model on (`--indel-rates`)
1. Divisive `dada`	default	yes
2. Single pass, p-value	`--partition single-pass`	yes
3. Single pass + reassign	`--reassign`	yes
4. Single pass, UNOISE rule	n/a: d counts gaps by definition	d = substitutions + internal gap columns

Comparing arm 2 with indels against arm 4 tests the learned error model against the fixed rule, with indels counted in both. The arms with indels off show what substitution-only DADA2 does.

Single pass already honours --indel-rates with no changes. On the MiSeq SOP (forward reads, pooled), the indel λ adds 2 and removes 1 ASV under divisive, and adds 2 under single pass, which is expected on Illumina. Learning indel rates leaves the error model unchanged (err_in and err_out are identical), so these comparisons share one model.

The UNOISE rule (arm 4)

β(d) = 1/2^(αd+1), with α = 2 (the preprint's default; vsearch --unoise_alpha also defaults to 2.0). A unique is absorbed by a centre if skew = its reads / the centre's reads ≤ β(d); otherwise it becomes a new centre.
d counts substitutions plus internal gap columns, with end gaps excluded. The preprint raises the end-gap ambiguity itself and avoids it by global trimming; excluding end gaps also matches our ends-free alignment and vsearch's --iddef 2 default.
If several centres qualify, the unique joins the one with the smallest d, with ties going to the more abundant centre. The preprint doesn't say which.
The abundance cutoff is left to --min-abund, set the same for every arm. UNOISE2 used γ = 4, and vsearch --cluster_unoise defaults to --minsize 8.

Counting differences to keep in mind:

The indel λ counts indel events toward hamming, while Levenshtein counts gap columns, so a 3-base deletion is 1 event but d = 3.
The indel λ ignores indels within --indel-end-zone (10) bases of the 3′ end, whereas arm 4 excludes end gaps.

Each arm keeps its own definition.

Cross-check: vsearch --cluster_unoise --minsize 8 --unoise_alpha 2 against arm 4 with --min-abund 8 on MiSeq V4. This will be documented on the command page and in the findings write-up.

feat: UNOISE birth rule for the single pass (#323)

dev: phase 5 arm-grid driver and summary under dev/unoise (#323)

feat: UNOISE rule counts every member's reads; vsearch UNOISE3 contro…

dev: CHIMERA_EXTRA for run_arms.sh chimera removal (#323)

dev: run_arms.sh timing without /usr/bin/time; wall and peak columns (#…

Navigation Menu

Repository navigation

Experiment: single-pass denoising in order of abundance (UNOISE-like) with DADA2's error model, across sequencing technologies #323

Origin

Background

Question

Arms

Data

Measures

Predictions, stated so they can be disproved

Plan

Approach: a switchable backend

Switches

Phases

Decisions to settle before phase 2

Phase 0 result: the most single pass could save

Phase 2: single pass with the p-value rule, plus greedy locks

First look: MiSeq SOP, forward reads, 20 samples, 8 threads (local)

Plan update: an indel axis, and the UNOISE rule's definition

Why

The UNOISE rule (arm 4)

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Experiment: single-pass denoising in order of abundance (UNOISE-like) with DADA2's error model, across sequencing technologies #323

Description

Origin

Background

Question

Arms

Data

Measures

Predictions, stated so they can be disproved

Activity

cjfields commented on Oct 31, 2026

Plan

Approach: a switchable backend

Switches

Phases

Decisions to settle before phase 2

cjfields commented on Nov 1, 2026

Phase 0 result: the most single pass could save

cjfields commented on Nov 1, 2026

Phase 2: single pass with the p-value rule, plus greedy locks

First look: MiSeq SOP, forward reads, 20 samples, 8 threads (local)

cjfields commented on Nov 1, 2026

Plan update: an indel axis, and the UNOISE rule's definition

Why

The UNOISE rule (arm 4)

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions