Optimized `pileup` can silently lose calls from independent same-base MM groups

[SuhasSrinivasan]

Optimized pileup can silently lose calls from independent same-base MM groups

Summary

The optimized pileup adapter does not advance independent MM groups for the same canonical base independently. While locating one group's explicit event, another group's delta can reach zero; the adapter then consumes both groups at the first position, assigns the later group's probability there, and silently loses its valid downstream call while the command exits successfully.

Severity

Severity: High — scientific correctness

Rationale: A valid accepted MM/ML tag can produce the wrong modification identity/quality at one genomic position and omit a true call at another without an error. This directly changes pileup rows and biological call counts. Exposure requires separate MM groups for the same canonical base with interleaved delta progression.

User and scientific impact

• Affected result or workflow: optimized pileup paths selected by explicit modified-base/all-context processing.
• Direction of error: downstream call omission and upstream modification-code/probability misassignment.
• Likely exposure: data dependent; requires independent same-base MM groups rather than one combined multi-code group.
• Detectability or workaround: force the generic pileup processor with equivalent multi-base motif selections and compare nonzero calls. There is no diagnostic in the affected optimized run.

Affected versions and environment

• Affected release: modkit 0.6.4.
• Development revision: 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c.
• Operating system and architecture used for verification: macOS on Apple Silicon.
• Input format: valid coordinate-sorted/indexed modBAM with separate same-base MM groups.

Steps to reproduce

Create a two-base reference and one valid alignment:

printf '>chr1\nCC\n' > ref.fa
samtools faidx ref.fa

printf '@HD\tVN:1.6\tSO:coordinate\n@SQ\tSN:chr1\tLN:2\nindependent-groups\t0\tchr1\t1\t60\t2M\t*\t0\t0\tCC\t??\tMM:Z:C+m?,0;C+h?,1;\tML:B:C,200,250\tMN:i:2\tNM:i:0\n' > input.sam
samtools view -b -o input.bam input.sam
samtools index input.bam

Run the optimized path and print only rows with a nonzero modified count:

modkit pileup input.bam optimized.bed \
  --ref ref.fa \
  --modified-bases C:m C:h \
  --no-filtering \
  --threads 1 \
  --io-threads 1 \
  --suppress-progress

awk -F '\t' '$12 > 0 {print $2, $4, $12}' optimized.bed

Observed behavior

0 h 1

The q0 event should be m with raw ML 200, but it is reported as h using ML 250. The valid q1 h event is absent.

Control or independent oracle

Force the generic path with two equivalent motif anchors:

modkit pileup input.bam generic.bed \
  --ref ref.fa \
  --motif CC 0 \
  --motif CC 1 \
  --no-filtering \
  --threads 1 \
  --io-threads 1 \
  --suppress-progress

awk -F '\t' '$12 > 0 {print $2, $4, $12}' generic.bed

The generic path retains both events (motif suffixes in column 4 are expected):

0 m,CC,0 1
1 h,CC,1 1

The direct adapter oracle is q0 m/200 and q1 h/250.

Expected behavior

Optimized and generic pileup retain the same two biological events: q0/genomic position 0 is m with raw probability 200, and q1/position 1 is h with raw probability 250. Each independent MM group advances only when its own explicit event is returned.

Root-cause evidence

• modkit-core/src/pileup/base_mods_adapter.rs: scanning decrements every matching canonical-base group's pending delta. Once any group is ready, aggregation and move_forward treat every group now at zero as explicit at that returned position.
• On C+m?,0;C+h?,1;, scanning q0 makes the second group's delta reach zero even though its explicit event belongs to q1. The parent consumes both at q0.
• The focused parent regression emits only q0 as h/250; the corrected adapter emits q0 m/200 and q1 h/250.

Proposed fix scope

Record which code slots were already ready at the returned position, use that readiness set for probability aggregation, and advance only those groups. Keep the repair isolated from cardinality/storage changes.

Non-goals

• No ambiguous/IUPAC-base handling; that is tracked separately in issue pileup can panic or lose downstream calls on valid ambiguous query bases #668 and upstream PR Allow N #597.
• No support for more than 16 modification-code slots or more than eight motifs.
• No generic-worker redesign, malformed MM/ML/MN policy change, negative-strand MM policy, threshold/count-schema change, threading change, or performance refactor.

Acceptance criteria

• The minimal fixture emits exactly the q0 m/200 and q1 h/250 events in the optimized adapter.
• Forced-generic and optimized scientific events agree after motif-label normalization.
• One-group, combined-code, implicit-call, reverse, and ordinary same-base controls remain unchanged.
• Focused adapter/CLI tests and the applicable full workspace suite pass.

Reproduction artifacts

Artifact	Size	SHA-256	Notes
input.sam	137 bytes	12929266ea4c1a2d278ca1aee2f40c6c15fd047dbedee41b83da85f18537e8f3	Inline SAM above.
input.bam	278 bytes	6ba4ccc9992767a6e5655b0f91c5edae76383c957ac799d0745d70663a02cfe5	samtools 1.23.1 conversion.
input.bam.bai	96 bytes	826d99d8857c284b8ca88b391b52174e1c40052d35c89a0dd7ca9e3d32f898c9	samtools 1.23.1 index.
ref.fa	9 bytes	dc95419df161b61d6cb661be877d6fdc3df1fabf3e97bb1bcc3a277d0666d44a	Inline FASTA above.

Related work

• Related issues/PRs: issue pileup can panic or lose downstream calls on valid ambiguous query bases #668 and PR Allow N #597 cover ambiguity handling in the same adapter but not independent MM-group progression; issue pileup can omit or split counts for valid combined and explicit modification selections #655/PR Correct pileup record state and count accounting #656 cover separate pileup slot/count mapping.
• Proposed PR: to be linked after dependency-resolution PR Pin hts-sys to bindings compatible with rust-htslib 0.46 #666 and the overlapping adapter work are integrated, rebased, and freshly regated.
• Deferred policy question or upstream dependency: none for this scientific behavior.

[SuhasSrinivasan]

The complete reproducer, fix, regression oracle, and test evidence are now consolidated in ready-for-review PR #709. Closing this standalone issue so the maintainer can review the fix directly in the PR, per the requested PR-only workflow.