Optimize the WFA backend on long (PacBio) reads

[cjfields]

Follow-up to #49 / #50. The experimental WFA backend (--align-backend wfa2, off by default) is ASV-equivalent to the default NW aligner on both Illumina and PacBio HiFi, but it is slower than NW on long reads, which is the opposite of what WFA should deliver. This issue tracks digging into why and closing the gap.

Measurements so far (PacBio HiFi, ~1472 bp full-length 16S, learn-errors, band 32)

Variant	Time	vs NW
NW (default)	~29 s	1×
WFA, fresh aligner per call	>6 min (killed, unfinished)	>13×
WFA, cached per-thread aligner ([#50])	~124 s	~4.3×

Illumina (short reads) shows no meaningful perf difference. The pathology is specific to long reads.

Why this is surprising / worth digging into

WFA's complexity is O(n·s) (length × alignment score), so on similar sequences it should beat NW's O(n·band). The pure-Rust port (COMBINE-lab/wfa2lib-rs) documents parity with — and in some cases (e.g. gap-linear) a slight edge over — the original C++ WFA2-lib. So a >4× slowdown on long reads points at how we're driving it, not the algorithm itself.

Hypotheses to investigate

1. Full ends-free free-gap span. align_wfa_endsfree_with_buf calls set_alignment_free_ends(plen, plen, tlen, tlen) — i.e. the entire sequence length is free at every end. On ~1500 bp reads this may make WFA explore a large free-end space each alignment. A bounded free span (e.g. ~band, matching the amplicon overhang we actually expect) may be much faster — and could also help the known ends-free crediting divergence.
2. CIGAR (ComputeAlignment) vs score-only wavefront retention. ComputeAlignment keeps non-modular wavefronts (all scores alive) for traceback; on long reads that is more memory/time. Check whether the retention/reap_cigar path is the cost.
3. BiWFA (align_biwfa). O(s) memory, designed for long sequences with full traceback. Worth A/B-ing against the plain affine aligner on the PacBio fixture.
4. Heuristics off by default. wfa2lib supports adaptive/Z-drop pruning (set_heuristic). Could speed long reads, but must be validated for ASV-equivalence before use.
5. Penalty magnitude / Eizenga scaling. Our Eizenga-adjusted penalties roughly double mismatch/gap internally (mismatch=18, gap_ext=21 for the default scheme), inflating the score s that drives WFA's work. Worth checking whether a smaller equivalent scheme reduces wavefront growth.

Suggested tasks

• Add a WFA kernel micro-benchmark (mirror the existing bench_align_vectorized in nwalign.rs) on a fixed long-read pair, to isolate per-alignment cost from the pipeline.
• Profile a PacBio learn-errors WFA run (e.g. cargo flamegraph) to see where time goes (free-ends exploration vs wavefront compute vs traceback vs allocation).
• A/B bounded vs full free-ends span; bounded vs BiWFA.
• Re-enable the PacBio WFA concordance gate in concordance.yml once WFA is competitive (it is deliberately deferred in Experimental WFA alignment backend (wfa2lib-rs), opt-in via --align-backend (#49) #50).

Pointers

• Adapter + cached aligner: src/nwalign.rs (align_wfa_endsfree_with_buf, the thread_local cache).
• Divergence reproducer (related, ends-free): wfa_diagnose test.
• Default stays nw, so this is not user-facing-urgent — it gates whether WFA is worth promoting beyond experimental.

[cjfields]

Root cause found + first fix (PR #52): WFA was running unbanded

Kernel benchmarks (bench_wfa_long, bench_wfa_divergence, added in #52) show the WFA kernel is ~60× faster than NW on a near-identical 1.5 kb pair — the algorithm isn't the problem. The problem: WFA is O(n·s), and our adapter ran it unbanded while NW is banded. On divergent-but-screened pairs WFA's cost crosses NW's around ~130 edits / 1450 bp and grows from there (1.97× at 200 edits). PacBio mock-community cross-species pairs sit past that crossover.

This was also a latent semantic bug (WFA globally-optimal vs NW band-limited).

Fix (#52): pass the band to WFA as BandedStatic ±band. PacBio learn-errors 124 s → 68 s (4.3× → 2.4× NW). Concordance still passes both platforms (Illumina 1.000/1.000/1.000, PacBio 1.000/1.000/0.994).

Still open (remaining ~2.4×)

The banded kernel is ≈NW even at 200 edits, so the residual gap is likely per-call overhead not captured by the microbench, or pairs beyond the band width. Next: flamegraph a PacBio learn-errors WFA run to attribute the remaining time (candidates: set_alignment_free_ends full-length span recompute, wavefront init_alignment setup per call, CIGAR→al0/al1 conversion). Then decide whether to re-enable the PacBio WFA concordance gate (currently deferred to avoid adding ~minutes to CI).

[cjfields]

Illumina benchmark: WFA is faster than NW (opposite of PacBio)

Benchmarked on the MiSeq SOP (20 samples, paired V4 16S, ~250 bp), per-sample denoising, 4 threads:

step	NW	WFA
learn-errors (F+R)	35.5 s	33.1 s
dada (F+R)	8.1 s	7.9 s
total	43.6 s	41.1 s (~6% faster)

This confirms the read-length story: WFA's O(n·s) cost + per-call overhead hurt on 1500 bp PacBio (2.4× slower) but are a net win on 250 bp Illumina — short, similar reads are WFA's sweet spot.

Tooling: added --align-backend nw|wfa2 to dev/benchmark/bench_pooled.py (branch perf/wfa-illumina-bench), injected into every aligning step (learn-errors, dada*, remove-bimera-denovo); R always uses NW so --align-backend wfa2 --run-r is an honest R-NW vs dada2-rs-WFA comparison.

ASV equivalence + one divergent call

• Forward reads: byte-identical ASV set (360 = 360, churn 0).
• Reverse reads: 304/305 shared; WFA splits one extra abundance-4 ASV in sample F3D144. Mechanism confirmed: both backends give 64 ASVs in F3D144, but WFA's slightly different ends-free alignment of those 4 reads to their cluster center shifted the substitution counts → the abundance p-value dropped below OMEGA_A (birth_pval=6.95e-42) in WFA but not NW → the partition decision flipped. Same nearest neighbor (Hamming-11) in both.

This single low-abundance divergence traces directly to the open ends-free free-end-gap crediting gap in the fork (the wfa_diagnose finding); fixing that makes WFA alignments byte-identical to NW → identical partitions. Recall stays 1.000; precision 304/305 = 0.997. No ground truth (F3D144 is a real mouse-gut sample, not the Mock).

[cjfields]

Resolved via an edit-budget cap — and the slowdown turns out to be k=5-only

Profiling attributed the residual PacBio slowdown to genuine WFA extend compute (O(n·s)), not per-call overhead. The fix is to cap the edit budget and fall back to NW for pairs that exceed it (set_max_alignment_steps → on MaxStepsReached the CIGAR is empty, so we re-align that pair with the banded NW path — NW-identical for exactly those divergent pairs).

The slowdown is a k=5-only phenomenon

PacBio concordance fixture (6740 filtered reads, band=32, 1 thread), learn-errors:

regime	NW	WFA uncapped	WFA cap=320 (≈40 edits)
k=5 (matches R's KMER_SIZE)	61.7s	141.8s (2.3× slower)	59.1s (beats NW)
k=7 (recommended PacBio)	13.1s	7.3s (1.8× faster)	8.7s

At k=7 the k-mer screen already removes divergent pairs, so WFA wins outright and the cap is unnecessary (slightly hurts). At k=5 the looser screen lets divergent non-error-copy pairs through and WFA pays full O(n·s) on them; the cap bounds that. Banding alone does not — k=5 uncapped was 141s with band=32; the cost cap is the real bound.

Zero ASV cost

Full pipeline at k=5 + cap, vs the static R reference: 93=93 ASVs, recall 1.000, precision 1.000, count-corr 0.994 — identical to the NW baseline. The cap only truncates pairs that were never error-copies; real error-copies sit far under budget given the ~99.9% denoising-identity assumption.

Productized

• --wfa-max-edits on all 6 alignment commands, default 50 edits, baked into the error-model JSON (with the usual consistency warning), mirroring --align-backend.
• The budget is an absolute edit count, not a fraction of read length — error-copies stay a few edits apart regardless of length. Internally edits·|gap_p| converts to the WFA cost unit (a scoring-robust upper bound).
• DADA2RS_WFA_MAX_STEPS env still overrides for ad-hoc sweeps.
• Re-enabled the previously-deferred PacBio WFA concordance gate in CI (k=5, default cap).
• WFA + the cap remain labelled EXPERIMENTAL in --help.

Branch: perf/wfa-illumina-bench.