`pileup` can panic or lose downstream calls on valid ambiguous query bases

[SuhasSrinivasan]

Summary

Valid IUPAC ambiguity symbols in BAM SEQ can make modkit pileup panic or abandon the rest of an otherwise usable record. The optimized reverse-read scanner assumes every query symbol is A/C/G/T, while both optimized and generic aligned-position paths either unwrap conversion of an ambiguous symbol or skip the entire record when only that one observation is unusable.

As a result, modkit 0.6.4 can exit 101 after writing a partial bedMethyl file, or exit 0 after silently losing a valid downstream modification call and its coverage. These behaviors also reproduce at current upstream revision 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c.

Severity

Severity: High — scientific completeness and reliability

Rationale: Accepted query ambiguity can selectively remove downstream modified calls from affected reads without invalidating their MM/ML groups. Some paths fail visibly, but others count the record as failed, write no row, and exit successfully. The failure can therefore bias counts and coverage toward reads without ambiguous base observations.

User and scientific impact

• Affected result or workflow: optimized and generic pileup, including reverse-aligned reads and records with ambiguity before a later valid modification call.
• Direction of error: loss of valid downstream modified calls and coverage; in panic paths, a nonzero exit with a partial output file.
• Likely exposure: data-dependent. It affects BAMs whose query sequences contain N or another BAM IUPAC symbol, including workflows that synthesize or transform modBAM records.
• Detectability or workaround: panics are visible, but record-local loss can appear only as a failed-record log and lower counts. Pre-filtering every ambiguous read avoids the crash but discards the same scientifically usable calls.

Affected versions and environment

• Released version: modkit 0.6.4
• Development revision: 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c
• Operating system and architecture: macOS 26.6, arm64
• Input format and index: coordinate-sorted BAM with BAI
• Related tools: samtools 1.23.1 for the synthetic fixtures

Steps to reproduce

These fixtures are synthetic and contain no private data.

Create reference.fa:

>chr1
CCTGTACTTC

Create leading-ambiguity.sam. The N at query position 0 precedes a valid modified C at query/reference position 1; the trailing N exercises the no-pending-call path.

@HD	VN:1.6	SO:coordinate
@SQ	SN:chr1	LN:10
leading-ambiguity	0	chr1	1	60	10M	*	0	0	NCTGTACTTN	IIIIIIIIII	MM:Z:C+m?,0;	ML:B:C,255	MN:i:10

Create softclip-ambiguity.sam. The first modified C is soft-clipped, then an aligned N occurs before the second valid modified C.

@HD	VN:1.6	SO:coordinate
@SQ	SN:chr1	LN:10
softclip-ambiguity	0	chr1	1	60	1S9M	*	0	0	CNCTGTACTT	IIIIIIIIII	MM:Z:C+m?,0,0;	ML:B:C,255,255	MN:i:10

Build and index both BAMs:

samtools faidx reference.fa
samtools view -b -o leading-ambiguity.bam leading-ambiguity.sam
samtools index leading-ambiguity.bam
samtools view -b -o softclip-ambiguity.bam softclip-ambiguity.sam
samtools index softclip-ambiguity.bam

Run both the optimized (--modified-bases C:m) and generic paths:

modkit pileup leading-ambiguity.bam leading-optimized.bed \
  --no-filtering --ref reference.fa --threads 1 --modified-bases C:m
modkit pileup leading-ambiguity.bam leading-generic.bed \
  --no-filtering --ref reference.fa --threads 1

modkit pileup softclip-ambiguity.bam softclip-optimized.bed \
  --no-filtering --ref reference.fa --threads 1 --modified-bases C:m
modkit pileup softclip-ambiguity.bam softclip-generic.bed \
  --no-filtering --ref reference.fa --threads 1

The optimized reverse-scanner panic can be reproduced independently with this paired fixture.

Create paired-reference.fa:

>chrF3
CAN
>chrR3
NTG

Create paired-reverse.sam:

@HD	VN:1.6	SO:coordinate
@SQ	SN:chrF3	LN:3
@SQ	SN:chrR3	LN:3
forward	0	chrF3	1	60	3M	*	0	0	CAN	III	MM:Z:C+m?,0;	ML:B:C,255	MN:i:3
reverse	16	chrR3	1	60	3M	*	0	0	NTG	III	MM:Z:C+m?,0;	ML:B:C,255	MN:i:3

samtools faidx paired-reference.fa
samtools view -b -o paired-reverse.bam paired-reverse.sam
samtools index paired-reverse.bam

modkit pileup paired-reverse.bam paired-optimized.bed \
  --no-filtering --ref paired-reference.fa --threads 1 \
  --modified-bases C:m

Control or independent oracle

MM deltas count occurrences of the group's stated fundamental base in the original read orientation. An N, R, Y, or other ambiguity symbol is not a C and must not consume a C+m delta. An ambiguous aligned observation cannot be classified as a specific A/C/G/T mismatch or canonical call, but it does not invalidate later MM/ML observations in the same read.

For each of the first two fixtures, both pileup workers must therefore emit exactly this downstream 5mC row:

chr1	1	2	m	1	+	1	2	255,0,0	1	100.00	1	0	0	0	0	0	0

Observed behavior

On modkit 0.6.4:

• Both leading-ambiguity commands exit 0, report approximately one failed record, and write an empty output file. The valid downstream 5mC call is silently lost.
• Both softclip-ambiguity commands exit 101 after unwrapping InvalidDnaBase in pileup_processor.rs; each leaves an empty output file.
• The paired optimized command writes only the forward row, then exits 101 from base_mods_adapter.rs after the reverse record reaches an unreachable!() branch. The 51-byte partial output has SHA-256 794047e609d4b3993f27bff54b1a6249c93d1d472668b0f766a91fd761eb4f6d.

With the bounded repair, optimized and generic execution both emit the exact expected row for each downstream-call fixture. All four files are 50 bytes with SHA-256 41458d8efc760a6849bbc6e95f2322be5d02c4a602ec4e24d24325d942866319. The paired optimized command exits 0 and emits both rows; its 102-byte output has SHA-256 5a315a8f4557d043736415c4f9743619f39f8f7d31974b3ed620708859df1784.

Expected behavior

• Valid BAM IUPAC query symbols do not panic either pileup worker.
• Ambiguity symbols do not consume A/C/G/T MM-group deltas.
• Only the ambiguous aligned observation is skipped; pending state and every valid downstream modification remain available to the rest of the record.
• Ambiguous observations are not misclassified as a canonical call, mismatch, deletion, or No Call for an invented A/C/G/T identity.
• Valid ambiguity does not increment the failed-record count.
• Optimized and generic workers agree on scientific positions, modification identities, counts, and valid coverage.

Root-cause evidence

• Reverse-read A/C/G/T frequency accounting treats every other BAM query symbol as unreachable at base_mods_adapter.rs:68-87.
• The optimized adapter's complement helper panics for every non-ACGT symbol at base_mods_adapter.rs:466-474.
• The optimized worker either abandons the whole record on a failed query-base conversion or unwraps the same conversion while advancing pending modification state at pileup_processor.rs:320-420.
• The generic worker has an analogous pending-state unwrap path at pileup_processor.rs:830-900, plus record-level skips at neighboring conversion sites.

Proposed fix scope

• Ignore non-ACGT query symbols when calculating optimized reverse-read A/C/G/T frequencies and MM-delta offsets.
• Complement the complete BAM IUPAC alphabet explicitly so reverse scanning remains orientation-correct.
• Centralize fallible A/C/G/T conversion for aligned query observations in the two pileup workers.
• Skip only the ambiguous observation, saving any already-advanced pending modification state before continuing to the next aligned position.
• Add exact forward/reverse adapter and optimized/generic end-to-end regressions.

Non-goals

• No support claim for MM groups whose fundamental base is N or U; this change concerns ordinary A/C/G/T MM groups on reads containing ambiguous query symbols.
• No new support for negative-strand MM groups.
• No malformed MM/ML recovery, generic record-error policy, or output rollback/atomicity change.
• No classification of an ambiguous observation as canonical, mismatch, deletion, or No Call.
• No requirement that optimized and generic bedMethyl column-4 motif labels be byte-identical where their existing documented labels differ.
• No threading, scheduling, threshold, schema, or performance redesign.

Acceptance criteria

• The two downstream-call fixtures emit exactly one 5mC row with modified count and valid coverage equal to one in both workers.
• The paired forward/reverse fixture emits both rows and exits 0 in optimized mode.
• A direct adapter regression proves exact forward and reverse positions and ML-derived qualities when ambiguity symbols occur between canonical bases.
• All supported BAM IUPAC complement pairs are covered without consuming unrelated MM deltas.
• Existing canonical-only pileup output remains byte-identical.
• Focused adapter, pileup integration, and full workspace tests pass.
• A representative direct-RNA output/performance guard shows no material hot-loop regression.

Reproduction artifacts

Artifact	Size	SHA-256	Notes
reference.fa	17 B	f3e4974762aaceac59f73162525a2c93a27baf1051ef31663fe2530dad06cd72	Ten-base downstream-call reference
leading-ambiguity.sam	137 B	6bf5f06412eeab65fd7a92c2990ef09587f30a7e4963075a77573e4f4101eee2	Ambiguity before and after one valid call
leading-ambiguity.bam	283 B	80dbb07cf3f5e412a84550ca9f33333aff0de25819ea5ba02c57b1b03f074187	Indexed fixture data
softclip-ambiguity.sam	145 B	5156b5d8c74e5a8c5ac894e2cc7e851e5cbe5a7a35843343e8c3824eeb64b04c	Unaligned call, ambiguity, downstream call
softclip-ambiguity.bam	290 B	4d17850fdb57cbb03c2411728634881d4fc2c2dbde007571f737266179c7136c	Indexed fixture data
paired-reference.fa	22 B	9bc64067eb65501ab4db210fc61868bc2785dd27de2ec0fd000b6bd07dda553f	Forward/reverse adapter reference
paired-reverse.sam	200 B	53be83b534c050a525197b0b148d047fbe97ddf825927a30f61abd1d50add0c7	Optimized reverse-scanner fixture
paired-reverse.bam	296 B	4979d3dbbfd695fcee45d7f53383e67c2e636d509e2b52ac3dd0da1175f22e14	Indexed paired fixture data

Related work

• PR #597, Allow N, by @killercup, correctly identifies and addresses the initial reverse-frequency crash and several query-base conversions. As of its current head a6115966351539b1f30d917eaf450f1bd4658e43, it remains open and materially overlaps this finding. The exact downstream-call fixtures above additionally exercise optimized and generic q < next modification handling and generic pending-state advancement, which still abandon the record or panic in that patch. This report is intended to make the remaining behavior and exact scientific oracle explicit, not to discount that contribution.
• The maintainer noted on Allow N #597 that this adapter is a hot pileup loop and requested performance evidence. The bounded repair has an existing representative direct-RNA byte/performance guard and should retain that evidence in any completion PR.
• Fresh dependency resolution selects incompatible hts-sys bindings #665 / Pin hts-sys to bindings compatible with rust-htslib 0.46 #666 address the separate fresh-build dependency incompatibility and do not change the ambiguity semantics.