Evaluate alignment paradigms that solve the global problem differently (pair-HMM, WFA)

[cjfields]

Background

A brief lit review confirms that Needleman-Wunsch / Gotoh global DP is the dominant pairwise-alignment paradigm, and that the research frontier is mostly acceleration of the same recurrence (WFA, KSW2, Parasail/striped SW, Edlib bit-parallel, GASAL2/GPU) rather than new alignment models.

The alignment paradigm in dada2-rs is deliberately fixed to global / ends-free NW because the error-model contract requires a positional alignment: al2subs compresses it into a Sub (per-position map + substitution list with both nucleotides and qualities), which the p-value model in pval.rs consumes to classify each difference as biological divergence vs base-calling error. Distance-only or local-only methods break this pipeline.

Rule of thumb for what's worth evaluating: a paradigm should address the global alignment problem from a genuinely different perspective — not change the problem.

Candidates that clear that bar

1. Pair-HMMs — fuse alignment and the error model into a single likelihood layer. DADA2 currently bolts pval.rs onto an NW alignment as two separate steps; a pair-HMM models error-vs-signal natively. Richer, but higher cost/complexity — assess whether it's justified for amplicon data.

2. WFA (wavefront alignment) — changes the complexity regime to O(n·s) (sequence length × edit distance), exploiting the same "sequences are very similar" assumption DADA2 already relies on (k-mer pre-screen + banded DP). Strongest candidate by the "different perspective, same problem" criterion. Caveat: we need a variant that retains the positional traceback, not score/edit-distance only, to feed al2subs.

Scope

• Literature/feasibility evaluation, not a commitment to implement.
• For any candidate: confirm it can produce the positional substitution record pval.rs needs, and benchmark vs the current anti-diagonal auto-vectorized banded NW (align_vectorized_with_buf).
• Out of scope: switching away from a global paradigm (local/Smith-Waterman would wrongly trim divergent-but-real read ends).

[cjfields]

Scoped note: pair-HMM as a learn-errors deliverable (not just a faster aligner)

Separating two distinct ideas under the pair-HMM heading, since they carry different risk/reward:

• Pair-HMM as a faster aligner (WFA's lane): lower conceptual risk; must preserve the positional Sub → pval.rs contract.
• Pair-HMM as the error model itself (the more interesting avenue): higher reward, higher risk. It changes what learn-errors produces. Today that's a transition × quality-binned error matrix consumed pointwise by pval.rs, with alignment (NW) and error modeling as separate steps. A pair-HMM fuses them — indels and substitutions in one generative model — which is more principled than DADA2's comparatively crude gap handling. Payoff is the possibility of a model that is efficient in the dada stage and potentially shared across sequencing runs.

The crux: cross-run transferability

The win only materializes if the model travels across runs. Run-to-run variation, library-prep and chemistry-batch differences mean a single frozen model won't generalize. So viability hinges on structure, not just retraining speed:

• A hierarchical / partially-pooled parameterization that separates technology-intrinsic error structure (relatively stable) from run-specific nuisance (cluster density, chemistry, prep): shared backbone + small per-run adaptation. The model needs to be retrainable / flexible to this variation.
• Must keep quality-score conditioning — DADA2's discriminating power comes from binning errors by quality. A quality-agnostic pair-HMM would regress; conditioning emissions on quality multiplies parameter count, feeding back into overfitting/retrainability.
• Mechanically already expressible: learned params are serialized in JSON alongside the dada/learn-errors params, so "train on run A, apply to run B" is at least representable in the current design.

Suggested framing

Evaluate the pair-HMM specifically as a learn-errors deliverable, with viability gated on a cross-run-transferable (hierarchical/adaptable) parameterization — not on raw retraining speed. This is a separate, larger avenue than the WFA acceleration track and should be assessed on its own.

[cjfields]

For future note: any approach outlined here would be very exploratory at this point. It's worth noting here in case this is something that we want to delve into at a later point, and may actually become a completely separate project since one (using paired HMM or similar approaches) would deviate significantly from the original DADA2 implementation.

[cjfields]

WFA track: dependency = COMBINE-lab/wfa2lib-rs

Decision: build the WFA acceleration track on COMBINE-lab/wfa2lib-rs — a from-scratch pure-Rust WFA implementation (no C WFA2-lib FFI; no cc/bindgen, so no cross-compile complication).

Availability

Not on crates.io (sparse index checked: wfa2lib-rs / wfa2 / wfa all absent). Consumable only as a git dependency for now:

wfa2lib-rs = { git = "https://github.com/COMBINE-lab/wfa2lib-rs", rev = "<pin-a-sha>" }

Caveats — all resolved

1. License: BSD-3-Clause, © 2026 COMBINE lab. Clean. Only gap: not declared in its Cargo.toml.
2. Library, not just a binary: src/lib.rs exposes WavefrontAligner; the CLI is isolated in src/bin/. Same shape as the original WFA2 C++ library.
3. edition 2024 / MSRV rustc 1.91: we're on rustc 1.95 and already edition 2024. Nuance — our Cargo.toml declares rust-version = "1.87"; depending on wfa2lib-rs raises our effective floor to 1.91, so we'll decide whether to bump the advertised MSRV. Not a build blocker.

API fit with the Sub → pval.rs contract

WavefrontAligner + WavefrontPenalties::new_affine() (match/mismatch/gap), align_end2end(pattern, text), CIGAR in the cigar field. CIGAR → positional [al0, al1] pair → al2subs satisfies the positional-substitution contract pval.rs requires. Supports global/semi-global/extension, score-only or full traceback.

Plan

• Contribute upstream rather than maintain a divergent fork (we have a line to COMBINE-lab). Working on a fork + local checkout in the meantime.
• Target PRs back to COMBINE-lab: declare license/description/repository in Cargo.toml; feature-gate clap/tracing/mimalloc to the bin so library consumers don't inherit the CLI dep tree; optionally nudge toward a crates.io release to drop the git pin later.
• In dada2-rs: depend via pinned git rev on the fork branch initially.

[cjfields]

feat/wfa-backend — first integration result: WFA ends-free is not yet a drop-in

Wired wfa2lib-rs in (path dep, default-features = false to avoid its mimalloc-alloc global allocator; bumped our rust-version 1.87 → 1.91 for its MSRV) and wrote a CIGAR→[al0,al1] adapter (align_wfa_endsfree_with_buf) so it can be compared head-to-head against align_endsfree.

Penalty mapping confirmed working (Eizenga adjustment handles the non-zero match): AffinePenalties { match_: -5, mismatch: 4, gap_opening: 0, gap_extension: 8 }. Two integration gotchas: the scope defaults to ComputeScore (empty CIGAR) — must set alignment_scope = ComputeAlignment; and op semantics are M/X = advance both, I = gap in pattern (al0), D = gap in text (al1).

Finding: divergence at gap/homopolymer boundaries and ends

Four hand-picked realistic cases (identical, 1 sub, 1 gap, diff-length) pass exact parity. But a 10k randomized sweep shows WFA-via-Eizenga ends-free does not reproduce DADA2's score-maximizing optimum — it returns alignments 1–8 points worse, concentrated at gap/homopolymer boundaries and sequence ends. Failure rate scales with divergence: 0% at 0 edits, ~1% at 1 edit, ~10% at 8 edits.

Concrete (1 edit, extra base in a leading AAA run):

EF : -AACTAGAC...   (free leading end-gap, score 275)
     AAACTAGAC...
WFA: AA-CTAGAC...   (penalized interior gap, score 267)
     AAACTAGAC...

Also the reverse: cases where EF's interior-gap solution beats WFA's mismatch + free-end-gap. So it's specifically the free end-gap crediting in the fork's ends-free path, not the adapter (which is exact on end2end-equivalent cases).

Implication & next step

WFA cannot be a drop-in for align_endsfree as-is — it would violate the concordance guardrail. Next step is to investigate the fork's align_endsfree begin/end-free handling (the begin-free skip isn't being credited at cost 0 in these cases). Likely a genuine upstream issue worth raising with COMBINE-lab. The branch is parked at a clean checkpoint: adapter + exact-parity unit tests green, plus an ignored sweep_wfa_parity (bucketed by edit count) and a wfa_diagnose test that prints counterexamples.

[cjfields]

Does affine gap scoring mitigate the divergence? No — isolation experiment

Added wfa_global_isolation to test whether the gap model is implicated. It compares WFA global (align_end2end) against the scalar global align_standard, alongside the ends-free comparison, on 10k random pairs:

Comparison	Disagree
WFA global vs align_standard	9 / 10000 (0.09%)
WFA ends-free vs align_endsfree	550 / 10000 (5.5%)

So ~98% of the divergence is in ends-free free-end-gap crediting, independent of the interior gap model. Conclusions:

• Affine gap scoring would not mitigate this — the gap model is not the cause. Affine also makes single-base gaps costlier (open + extend ≥ linear), which nudges WFA further toward the "mismatch + free-end-gap" resolutions that already lose to align_endsfree's interior-gap solutions. And switching WFA to affine diverges from DADA2's linear-gap reference by construction, so it can't improve parity-vs-DADA2.
• The fix belongs in the fork's align_endsfree begin/end-free handling, exactly where wfa_diagnose points.
• Residual 9/10000 in the global case is minor tie-breaking/edge behavior worth a later look, but it's noise next to the ends-free effect.

Branch feat/wfa-backend: commits e169a6e (backend + adapter) and e53167e (isolation experiment).

[cjfields]

Upstream packaging PRs opened against COMBINE-lab/wfa2lib-rs

Two of the integration caveats are now upstream PRs (branches on the HPCBio fork, targeting COMBINE-lab/wfa2lib-rs:main):

• COMBINE-lab/wfa2lib-rs#2 — feature-gate CLI deps: clap/tracing/tracing-subscriber made optional behind a cli feature so library consumers can use default-features = false and skip the CLI/logging tree. (branch feat/gate-cli-deps)
• COMBINE-lab/wfa2lib-rs#3 — add a public set_alignment_scope setter (default scope is ComputeScore, which yields an empty CIGAR; the setter makes ComputeAlignment discoverable). (branch feat/set-alignment-scope)

A third change — declaring license = "BSD-3-Clause" in Cargo.toml — is being submitted separately by @cjfields.

Still outstanding

• Cosmetic: switch the dada2-rs adapter from the direct alignment_scope = field write to set_alignment_scope(...) once Benchmarking checklists #3 lands.
• The real drop-in blocker is unchanged: the ends-free free-end-gap crediting fix in the fork's align_endsfree (the 5.5% divergence isolated above). Reproduce via the wfa_diagnose test.

[cjfields]

Considered and rejected: the libwfa crate

For the record (so it doesn't resurface): libwfa (chfi/rs-wfa) is the one WFA Rust crate actually published on crates.io, but it's a non-starter for us:

• FFI bindings to the original C WFA, not pure Rust — vendors ~3,600 lines of C across 25 files and requires cc/bindgen at build time. That reintroduces exactly the cross-compilation friction (musl/aarch64, C toolchain in CI) that pure-Rust wfa2lib-rs lets us avoid.
• Unmaintained: last release 0.1.2 in Oct 2020.
• Pre-dates WFA2's Eizenga non-zero-match support — the original WFA can't express a positive match score, so DADA2's scoring would not even map onto it. That generalization is the whole reason WFA is feasible here.

Being "on crates.io" doesn't outweigh being worse on every axis that mattered. The packaging gap is better closed by nudging wfa2lib-rs toward a crates.io release. We stick with wfa2lib-rs.

[cjfields]

End-to-end real-data A/B: WFA gives an identical ASV table on Illumina

Wired an opt-in dispatch (DADA2RS_ALIGN_BACKEND=wfa, read once via LazyLock; replaces only the vectorized/ends-free path, leaves the gapless and homopolymer paths untouched; zero-cost when unset — commit 257ce28) and ran a full forward-reads A/B.

Setup: data/dada2/sam1F + sam2F → filter-and-trim (trunc-len 240, maxEE 2) → learn-errors (default backend, held constant) → dada twice (default vs WFA backend).

Result: the ASV table is identical between backends —

	sam1	sam2
ASV sequences equal	✅	✅
abundance diffs	0	0
churn (compare_asvs.py)	0	0

The only differences anywhere in the output are the diagnostic birth_fold / birth_e / birth_pval stats, which wobble sub-threshold. That wobble is useful: it confirms the WFA flag genuinely engaged and the alignments did differ — they just never changed which sequences become ASVs or their counts.

Takeaway: WFA's 5.5% per-alignment ends-free divergence (on synthetic random pairs) does not propagate to ASVs on screened, band-limited real data. The boundary suboptimality is academic for this Illumina case.

Caveats: forward reads only, 2 Illumina samples, error model held constant. PacBio is untested and would be a tougher case (longer reads, more indels, homopolymers) and only a partial swap (the homopolymer path stays scalar). So this lowers the urgency of the ends-free fix for Illumina, but we should run a PacBio A/B before generalizing.

[cjfields]

PacBio HiFi A/B: also an identical ASV table (full swap)

Ran the same A/B on PacBio HiFi (data/dada2/samPB, 500 CCS reads, ~1472 bp full-length 16S; filter min-len 1000 / max-len 1600 / maxN 0 / maxEE 2; band 32; homo-gap-p left unset).

Important: for HiFi we do not use homopolymer settings, so homo_gap_p defaults to gap_p and use_homo is false — meaning the WFA backend handles every alignment. This is a full swap, correcting my earlier "partial swap" caveat for PacBio.

Result:

	samPB (PacBio HiFi)
ASVs (default / wfa)	11 / 11
sequences equal	✅
abundance diffs	0
churn	0

Again the only differences are the diagnostic birth_e / birth_fold stats (sub-threshold).

Conclusion across both platforms

WFA's 5.5% synthetic per-alignment ends-free divergence does not propagate to the ASV table on real data — confirmed on Illumina short reads (full pipeline) and PacBio HiFi full-length 16S (full swap).

So the ends-free free-end-gap fix in the fork is low priority for production output — worth doing as correctness hygiene, but not a blocker to using WFA as a backend. WFA is ASV-equivalent today on both platforms.

[cjfields]

CI: WFA backend now gated by the concordance guardrail

The concordance workflow now runs a second time with the WFA backend switched on, gated against the same static R DADA2 reference.

• run_illumina.sh / run_pacbio.sh honor an optional ALIGN_BACKEND env var (nw|wfa2), threaded as --align-backend through every alignment-using subcommand: learn-errors, dada, and remove-bimera-denovo. (The chimera step matters — it dispatches on its own BimeraAlignParams.backend, not the env-var fallback, so it would otherwise stay on NW.) Unset = default nw, so existing behavior is unchanged.
• concordance.yml adds "Illumina (WFA)" and "PacBio (WFA)" steps that run the full pipeline with ALIGN_BACKEND=wfa2 and compare to the same reference with the same thresholds.

Verified locally — the WFA Illumina run clears the gate identically to NW:

recall    1.000  (>= 0.85) PASS
precision 1.000  (>= 0.80) PASS
count corr 1.000 (>= 0.95) PASS

So WFA is now regression-guarded at the ASV level in CI, on the same footing as the default NW path. (PacBio WFA step is guarded on the not-yet-committed PacBio fixture, like the existing NW PacBio step.)

[cjfields]

Opened PR #50 (feat/wfa-backend) with the experimental WFA backend: opt-in --align-backend nw|wfa2 across dada/learn-errors/chimera, error-model provenance, and an Illumina WFA concordance gate (green). ASV-equivalent to NW on Illumina + PacBio HiFi; PacBio WFA gate deferred on perf (cached aligner is ~4.3× NW on long reads). Git-pinned to the public HPCBio wfa2lib-rs fork.