Preserve pileup calls around ambiguous query bases

[SuhasSrinivasan]

Fixes #668.

Summary

• Prevents valid BAM IUPAC query symbols from panicking optimized or generic pileup paths.
• Skips only an unclassifiable aligned observation while preserving pending state and valid downstream modification calls from the same record.
• Keeps ambiguity symbols from consuming A/C/G/T MM deltas and complements BAM IUPAC ambiguity pairs for reverse reads.
• Rejects a record containing the distinct SAM/BAM reference-equality symbol = before tallying and reports it through one aggregate diagnostic; resolving = from the reference is not implemented here.
• Adds exact forward/reverse adapter, optimized/generic CLI, packed-sequence, and release direct-RNA hot-loop regressions.

Severity

Severity: High — scientific completeness and reliability

Rationale: Accepted ambiguity can silently remove valid downstream modified calls and coverage while the command exits successfully, or can panic after writing partial output. This can bias pileup counts toward reads without ambiguous observations.

Root cause

The optimized reverse scanner treated every non-ACGT query symbol as unreachable and its complement helper panicked outside A/C/G/T. In both pileup workers, several aligned-position conversions either unwrapped an ambiguous query base or abandoned the entire record. One overrun path also advanced the modification iterator before converting the current query base but saved that pending state only after the fallible conversion, so skipping the observation lost downstream state.

Implementation

• Ignore non-ACGT symbols when calculating optimized reverse-read A/C/G/T frequencies used by MM deltas.
• Explicitly complement the BAM IUPAC ambiguity pairs, including R/Y, M/K, B/V, D/H, S, W, and N.
• Centralize aligned query-base conversion in an inline, allocation-free helper returning Option<DnaBase>.
• In optimized and generic overrun paths, save the already-advanced pending modification state before attempting to classify the current query observation.
• Continue the aligned-position loop when only that observation is ambiguous; do not increment the failed-record count.
• Scan the packed BAM sequence for code zero (=) before either worker parses modifications or changes a tally. Reject that complete record, increment the existing failed-record total and a dedicated approximate counter, and emit one final explanation of the reference-resolution requirement.

Preserved behavior

• A/C/G/T MM deltas, ML probabilities, thresholding, strand semantics, output schema, and worker scheduling remain unchanged.
• Ambiguous observations are excluded from every bedMethyl category, including N_diff; no A/C/G/T identity is invented.
• The tested canonical-only direct-RNA optimized release control is byte-identical.
• No allocation, clone, lock, sort, or synchronization is added per aligned position. The = guard adds one allocation-free packed-byte scan per processed record.

Non-goals

• No support claim for MM groups whose fundamental base is N or U, negative-strand MM groups, or malformed MM/ML recovery.
• No reference-backed resolution of SEQ =; affected records are rejected before tallying and pileup continues. Deciding whether unresolved IUPAC observations should increment N_diff is separate policy work.
• No generic record-error policy, output rollback/atomicity, threshold, schema, scheduler, or broader performance redesign.
• No requirement that existing optimized and generic motif labels be byte-identical where their selected reference-position scopes already differ.
• This branch materially overlaps PR #597: both ignore non-ACGT symbols in reverse-frequency accounting and make several aligned query-base conversions position-local. This branch extends that work by covering q-before-next-mod handling in both workers, generic overrun/pending-state preservation, BAM IUPAC ambiguity complements, explicit = rejection, and exact optimized/generic scientific and performance oracles. The maintainer should choose or combine the implementations rather than merge both overlapping diffs independently.

Behavior before and after

Case	Before	After	Expected oracle
Leading N before a valid downstream 5mC call	Both workers can abandon the record and emit no row while exiting zero	Skip the ambiguous observation and retain the call	Exact downstream row at chr1:1-2, modified count/valid coverage 1
Soft-clipped earlier call, then aligned ambiguity before a later call	Conversion unwrap can panic	Pending state is saved before the ambiguity is skipped	Same exact downstream 5mC row in both workers
Reverse read containing BAM ambiguity	Reverse frequency/complement path can panic	IUPAC is complemented without consuming A/C/G/T deltas	Both paired forward/reverse rows emitted; 102-byte SHA-256 5a315a8f4557d043736415c4f9743619f39f8f7d31974b3ed620708859df1784
Record containing SEQ = plus a valid earlier call	Treating = as ordinary ambiguity can silently include only part of the record	Reject the complete record before any tally; retain an otherwise identical control record and log the reason once in aggregate	Optimized and generic output exactly one modified call with valid coverage 1, not 2; both report one reference-equality rejection
Ordinary direct-RNA input without ambiguous query bases in the sampled scan	Existing output	Byte-identical	60,959 rows, 9,989,028 bytes, SHA-256 a85202b69e79e708656e71858e4243cb1cc51e22a44bb9f65650a23dacb407fa

Testing

Test environment

• Revision tested: bb20328778b2a5e3f0f815b36af2ae09c1df7ae7
• Tree tested and worktree state: 7cc1afd1313f8af725dd34635a4620c704a774fe; clean after testing
• Toolchain: rustc/cargo 1.90.0
• Platform: macOS 26.6, arm64 Apple Silicon
• Reference binary: installed modkit 0.6.4 and exact upstream parent 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c
• External tools: samtools 1.23.1
• Dependency resolution identity: ignored Cargo.lock SHA-256 f9d3389f8fc63c9ad653dad3fbf43a1c26cc54008f0a39e06e71ff6455a1e9b5; rust-htslib 0.46.0, hts-sys 2.2.0

Test layer	Exact command, fixture, or matrix	Result and evidence
Core: parent-red regression	Installed 0.6.4 and exact parent; leading ambiguity, softclip-overrun ambiguity, and paired reverse fixtures	Leading case silently loses the downstream call; softclip/reverse cases exit 101; partial paired output SHA-256 794047e609d4b3993f27bff54b1a6249c93d1d472668b0f766a91fd761eb4f6d
Core: focused adapter regression	cargo test -p mod_kit pileup::base_mods_adapter::base_mods_adapter_tests::test_ambiguous_query_bases_preserve_mm_scanning -- --exact --test-threads=1	1 passed, 0 failed; forward/reverse R/N/Y positions and ML-derived qualities are exact
Adapter module	cargo test -p mod_kit base_mods_adapter_tests --lib -- --test-threads=1	9 passed, 0 failed
Packed-sequence guard	cargo test --offline --locked -p mod_kit pileup::pileup_processor::reference_equal_base_tests::odd_length_padding_is_not_a_reference_equal_base --lib -- --exact --test-threads=1	1 passed, 0 failed; a real = is detected and an odd-length BAM record's unused low padding nibble is not
Core: = no-partial-tally regression	cargo test --offline --locked -p modkit --test test_pileup test_pileup_rejects_reference_equal_query_before_tallying -- --exact --test-threads=1	1 passed, 0 failed; optimized and generic paths retain only the valid control's exact 1/1 call and log both the explicit rejection and failed-record total
Core: affected pileup integration	cargo test -p modkit --test test_pileup -- --test-threads=1	19 passed, 10 ignored, 0 failed; both IUPAC tests and the = rejection test exercise optimized and generic workers
Core: applicable full workspace gate	TMPDIR=/private/tmp/modkit-668-supervisor-full-bb20328 cargo test --offline --locked --workspace --all-targets --no-fail-fast -- --test-threads=1	184 passed, 14 ignored, 0 failed on the exact final head
Installed-versus-fixed comparison	Public inline SAM/BAM fixtures in issue #668	Four optimized/generic downstream outputs are each the exact 50-byte row, SHA-256 41458d8efc760a6849bbc6e95f2322be5d02c4a602ec4e24d24325d942866319; paired reverse output matches the behavior table
= branch guard	Matched two-record optimized/generic fixture on the parent behavior and final branch	The parent includes both valid q0 calls as 2/2 without a reference-equality diagnostic; the final branch deliberately rejects the unresolved = record before tallying and emits exact 1/1 from the valid control. This guards the branch's = boundary rather than reproducing issue #668's original IUPAC defect
Optimized/generic semantic parity	Leading and softclip-overrun fixtures	Both workers emit the same scientific position, modification identity, count, and coverage without failed-record inflation
Exact release artifact	cargo build --offline --locked --release -p modkit at the tested head	23,833,584-byte binary, SHA-256 fb9dff56b0b5ec233e5f65f5fed6eb4e029dd408ca4a314c92ca293307f0a4e1; the earlier pre-= artifact did not match and its timing was discarded
Real-data performance and byte guard	Reconstructed command: modkit pileup rna_top50_tx.bam OUT --reference rna_top50_tx.padded.fa --modified-bases C:m A:a A:17596 T:17802 --no-filtering --threads 8 --io-threads 4 --suppress-progress; 898,198,413-byte/1,632,831-read direct-RNA BAM; one warm-up then U1/F1/U2/F2/U3/F3 release runs	Upstream wall/user/sys medians 6.87/23.73/0.32 s and means 6.8467/23.7133/0.3200 s; final medians 6.75/23.15/0.32 s and means 6.7567/23.1467/0.3200 s. The small apparent improvement is treated as run/build noise, not a speedup claim; no material regression is measurable. All six timed outputs and both warm-ups are byte-identical with the hash/count above
Generic unaffected control	One release generic run per exact upstream/final binary with the same command minus --modified-bases	Both outputs byte-identical: 131,900 rows, 20,887,346 bytes, SHA-256 2c98dc0193764a410a8d84b1227733ec101f47db5541ec7c313e7d07e1a447c8; raw generic/optimized files are intentionally not compared because their selected reference-position scopes differ
Hot-loop source/layout audit	get_query_base, packed = preflight, adapter frequency/complement, and both overrun state paths	O(1) borrowed aligned-base access; one allocation-free sequential scan over approximately half the sequence bytes per fetched record exposure; no lock/atomic/synchronization change. Existing or newly advanced pending state is preserved across each ambiguity skip, with newly advanced state saved before fallible conversion in both overrun branches. A read crossing K chunks is scanned K times, matching existing per-exposure processing; the default 1 Mb interval makes K=1 for every active contig in this control
Core: formatting/diff checks	Direct stable rustfmt --check on all five changed Rust files; git diff --check upstream/master...HEAD	Passed; only upstream nightly-setting warnings were emitted

Tests not performed

• The adapter test samples representative R/N/Y symbols rather than enumerating every IUPAC ambiguity symbol, although the source explicitly covers every supported ambiguity complement pair.
• The sampled 100,000-primary-record real-data scan contained no non-ACGT query bases; it is therefore a normal-hot-loop regression guard, not a real-data correctness exposure test.
• The literal earlier benchmark command was not retained; the command reported above was reconstructed and accepted only after reproducing the prior exact row/byte/SHA oracle.
• The ordinary-IUPAC integration does not parse N_diff directly; exclusion from N_diff is established by both workers skipping before every tally method.
• No output rollback, malformed-tag, negative-strand-MM, or GPU test was run because those paths are outside scope.

Scientific validation

• Population/eligibility invariant: one ambiguous aligned observation is skipped without removing the rest of an otherwise usable record.
• Downstream-call preservation invariant: each IUPAC fixture retains exactly one downstream modified call and its unit valid coverage; ambiguity itself contributes to no bedMethyl category, including N_diff.
• Whole-record rejection invariant: a record containing unresolved = contributes no partial call, while the matched valid control contributes exactly one call and unit coverage in both workers.
• Coordinate/strand/interval invariant: forward/reverse MM deltas count only their fundamental A/C/G/T bases and retain exact query/reference positions.
• Determinism invariant: fixed optimized/generic synthetic outputs agree scientifically; normal-path release artifacts are byte-identical across all A/B runs.
• Independent oracle or specialist review: MM delta positions and ML values are directly enumerable; CPU and supervisor reviews traced every state-save branch and verified the hot-loop cost.

Output and compatibility

• User-visible change: valid IUPAC symbols no longer panic or discard downstream pileup calls.
• Expected output differences: affected IUPAC reads recover scientifically usable downstream calls and coverage; the ambiguous position itself remains unclassified. A record containing = is rejected before tallying and produces an aggregate error log because the base identity requires reference resolution.
• Byte-identical controls: normal direct-RNA optimized and generic release controls are byte-identical to the parent.
• CLI/API/schema compatibility: no option, public API, column, or schema change.
• Partial-output or failure semantics: ordinary ambiguity is position-local; unresolved = is whole-record-fatal but not process-fatal, and the command continues after logging aggregate rejection/failure counts. This PR does not make output transactional for unrelated fatal errors.

Reviewer guide

1. Review reverse frequency/complement handling in base_mods_adapter.rs.
2. Review the packed-nibble = preflight, including its odd-length padding rule, and confirm it precedes all tally mutation in both workers.
3. Review the inline get_query_base helper and confirm every ordinary-IUPAC None branch skips only the current position.
4. Pay particular attention to pending-state assignments before break 'overran in both workers.
5. Rerun the adapter, two IUPAC, packed-padding, and = no-partial-tally regressions.
6. Compare the bounded additions with PR Allow N #597 and retain the direct-RNA performance evidence requested there.

Checklist

• The issue contains reproducible observed and expected behavior.
• The change is limited to the linked issue's approved scope.
• The regressions are demonstrably red on the exact parent/released binary.
• All tests actually performed are listed above with their results.
• Unrun or inapplicable tests are disclosed.
• Scientific counts/statistics and output compatibility are explicitly checked.
• Formatting and diff-hygiene checks pass, or unrelated findings are documented.
• No generated data, private sample identifiers, or unrelated changes are included.

[SuhasSrinivasan]

Closing this overlapping PR to keep #597 as the single active review path for ordinary IUPAC query bases. The two remaining downstream-state cases and the applicable performance evidence are documented in #597 (comment); this branch remains preserved for its exact tests. The separate SEQ = whole-record policy is intentionally deferred and is not being proposed as part of #597.