pileup can omit or split counts for valid combined and explicit modification selections

[SuhasSrinivasan]

Summary

Several related count-matrix assumptions in modkit pileup can omit valid modification calls, split one biological base/code category into contradictory rows, or place ancillary counts at a different coordinate from the corresponding canonical and modified calls. These behaviors reproduce on accepted, tag-valid inputs in modkit 0.6.4 and current upstream revision 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c.

The failures share the dynamic/combined pileup count-matrix and option-normalization path, so they are reported together and corrected as independently reviewable commits in one proposed PR.

Severity

Severity: High — scientific correctness

Rationale: The affected commands can exit successfully while omitting passing modification classes, reporting an incorrect N_fail, losing selected canonical bases, or splitting counts for one requested category. The errors directly change bedMethyl coverage and modification percentages and are not necessarily evident from the exit status.

User and scientific impact

• Affected result or workflow: bedMethyl output from combined-C, dynamic, high-depth, explicit base/code, duplicate-selection, and non-CpG strand-combined pileup modes.
• Direction of error: omission, category splitting, undercounting of failure/ancillary categories, or duplicate rows.
• Likely exposure: option-specific but routine for multi-modification and high-depth workflows.
• Detectability or workaround: compare equivalent normal/high-depth runs and inspect exact category conservation; otherwise the command commonly exits zero. Avoiding --high-depth, duplicate selections, or same raw codes across canonical bases is only a partial workaround.

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; indexed FASTA
• Fixture tool: samtools 1.23.1 / HTSlib 1.23.1

Steps to reproduce

The following fixtures contain no private data. Save each SAM/FASTA block with the shown filename, then build its BAM/index with:

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

1. Combined-C accounting omits valid modification classes and filtered calls

input.sam:

@HD	VN:1.6	SO:coordinate
@SQ	SN:chrC	LN:1
canonical	0	chrC	1	60	1M	*	0	0	C	?	MM:Z:C+m?,0;	ML:B:C,0	MN:i:1	NM:i:0
methyl	0	chrC	1	60	1M	*	0	0	C	?	MM:Z:C+m?,0;	ML:B:C,255	MN:i:1	NM:i:0
hydroxy	0	chrC	1	60	1M	*	0	0	C	?	MM:Z:C+h?,0;	ML:B:C,255	MN:i:1	NM:i:0
formyl	0	chrC	1	60	1M	*	0	0	C	?	MM:Z:C+f?,0;	ML:B:C,255	MN:i:1	NM:i:0
chebi	0	chrC	1	60	1M	*	0	0	C	?	MM:Z:C+19228?,0;	ML:B:C,255	MN:i:1	NM:i:0
filtered	0	chrC	1	60	1M	*	0	0	C	?	MM:Z:C+m?,0;	ML:B:C,128	MN:i:1	NM:i:0

reference.fa:

>chrC
C

modkit pileup input.bam combined-c.bed \
  --ref reference.fa \
  --modified-bases C \
  --combine-mods \
  --filter-threshold 0.9 \
  --threads 1 \
  --suppress-progress

Observed with 0.6.4:

chrC  0  1  C  3  +  0  1  255,0,0  3  66.67  2  1  0  0  0  0  0

Only m and h enter N_mod; passing f and ChEBI 19228 are omitted, and the threshold-failing call does not increment N_fail.

Expected:

chrC  0  1  C  5  +  0  1  255,0,0  5  80.00  4  1  0  0  1  0  0

Every passing selected C modification contributes to the combined modified count, and the filtered call contributes exactly once to N_fail.

2. Dynamic and high-depth slots can omit or split selected bases

Replace input.sam with:

@HD	VN:1.6	SO:coordinate
@SQ	SN:chr1	LN:3
modified	0	chr1	1	60	3M	*	0	0	ACT	???	MM:Z:A+m?,0;C+m?,0;T+g?,0;	ML:B:C,255,255,255	MN:i:3	NM:i:0
canonical	0	chr1	1	60	3M	*	0	0	ACT	???	MM:Z:A+m?,0;C+m?,0;T+g?,0;	ML:B:C,0,0,0	MN:i:3	NM:i:0

Replace reference.fa with:

>chr1
ACT

Rebuild and index the BAM/FASTA, then run:

modkit pileup input.bam c-normal.bed \
  --ref reference.fa --modified-bases C --combine-mods \
  --no-filtering --interval-size 1 --threads 1 --suppress-progress

modkit pileup input.bam c-high.bed \
  --ref reference.fa --modified-bases C --combine-mods \
  --no-filtering --interval-size 1 --threads 1 --suppress-progress \
  --high-depth --max-depth 10

modkit pileup input.bam explicit.bed \
  --ref reference.fa --modified-bases A:m C:m \
  --no-filtering --interval-size 1 --threads 1 --suppress-progress

Observed with 0.6.4:

• c-normal.bed contains the expected C row, while c-high.bed is empty.
• A:m C:m emits one correct A row but splits C into separate 100% modified and 0% modified rows, each with coverage 1.

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

Expected:

• Normal and high-depth outputs are byte-identical below saturation.
• A and C each have one 50% row with coverage 2, one modified call, and one canonical call.
• Selecting A C T --combine-mods produces all three rows; 0.6.4 omits T on this fixture.

3. Reverse ancillary counts use a different combined-strand anchor

Use this input.sam:

@HD	VN:1.6	SO:coordinate
@SQ	SN:chrR	LN:4
modified	16	chrR	1	60	4M	*	0	0	GATC	????	MM:Z:A+m?,0;	ML:B:C,255	MN:i:4	NM:i:0
canonical	16	chrR	1	60	4M	*	0	0	GATC	????	MM:Z:A+m?,0;	ML:B:C,0	MN:i:4	NM:i:0
filtered	16	chrR	1	60	4M	*	0	0	GATC	????	MM:Z:A+m?,0;	ML:B:C,128	MN:i:4	NM:i:0
no-call	16	chrR	1	60	4M	*	0	0	GATC	????	MM:Z:C+m?,0;	ML:B:C,255	MN:i:4	NM:i:0
mismatch	16	chrR	1	60	4M	*	0	0	GACC	????	MM:Z:C+m?,0;	ML:B:C,255	MN:i:4	NM:i:1
deletion	16	chrR	1	60	2M1D1M	*	0	0	GAC	???	MM:Z:C+m?,0;	ML:B:C,255	MN:i:3	NM:i:1

Use this reference.fa:

>chrR
GATC

Rebuild/index, then run:

modkit pileup input.bam reverse.bed \
  --ref reference.fa \
  --motif GATC 1 \
  --combine-strands \
  --modified-bases A:m \
  --filter-threshold 0.9 \
  --interval-size 4 \
  --threads 1 \
  --suppress-progress

Observed with 0.6.4:

chrR  1  2  m  2  .  1  2  255,0,0  2  50.00  1  1  0  0  0  0  0

Expected: deletion, failure, difference, and no-call must use the same combined anchor as canonical/modified calls, so the final four category columns are each 1:

chrR  1  2  m  2  .  1  2  255,0,0  2  50.00  1  1  0  1  1  1  1

4. Duplicate explicit selections are not idempotent

Using the two-read fixture from section 2:

modkit pileup input.bam unique.bed \
  --ref reference.fa --modified-bases C:m \
  --no-filtering --threads 1 --suppress-progress --use-dynamic

modkit pileup input.bam duplicate.bed \
  --ref reference.fa --modified-bases C:m C:m \
  --no-filtering --threads 1 --suppress-progress --use-dynamic

Observed with 0.6.4: the unique command emits one row; the duplicate command emits two rows for the same key. The same duplication occurs in high-depth mode.

Expected: exact duplicate (canonical base, modification code) selections are idempotent. The two files are byte-identical, while A:m and C:m remain distinct selections.

5. A valid CpG combined-C command panics in debug builds

From a source checkout at the affected revision:

cargo run -p modkit -- pileup \
  tests/resources/bc_anchored_10_reads.sorted.bam /tmp/cpg-combined.bed \
  --cpg --modified-bases C --combine-mods --no-filtering \
  --ref tests/resources/CGI_ladder_3.6kb_ref.fa --threads 1

Observed: the debug build panics because overlapping internal CG,0 and C,0 mask bits violate a single-hit assertion. Release mode already applies deterministic first-hit routing and remains CpG-only.

Expected: debug and release builds both succeed and emit the same CpG site set.

Root-cause evidence

• The combined-C matrix admits only m and h at pileup_processor.rs:1476-1481, while its filtered path assigns the counter to itself at pileup_processor.rs:1503-1504.
• Dynamic combined slots use the canonical-base enum value as a compact-vector offset, and explicit slots search by raw code without the canonical base at pileup_processor.rs:1885-1894.
• Reverse anchoring is applied to canonical/modified/depth paths but not consistently to delete/filter/difference/no-call paths in pileup_processor.rs:1819-1969.
• Modified-base selections are sorted but not deduplicated before preset and writer construction at subcommand.rs:581-612.
• The debug-only single-hit assertion is at pileup_processor.rs:260-263, although the following code already defines first-hit behavior.

Each behavior has a parent-failing regression and an exact fixed-output oracle in the proposed patch.

Proposed fix scope

• Count every selected modified cytosine in combined-C mode and increment the filtered counter with saturation.
• Use explicit canonical-base-to-compact-slot mappings and key explicit slots by (canonical base, modification code).
• Route every reverse ancillary category and depth lookup through the same combined anchor.
• Deduplicate exact selections stably before preset and writer construction, with a warning; do not merge the same raw code across different canonical bases.
• Remove only the invalid debug single-hit assertion while preserving existing first-hit motif selection.

Non-goals

• No sampling, threshold-estimation, ambiguous-base, motif-capacity, or general worker-error policy changes.
• No performance refactor or count-matrix representation redesign.
• No change to supported modification-code syntax or bedMethyl schema.

Acceptance criteria

• The minimal reproducers fail on the affected parent and pass after the fix.
• Combined-C N_mod, N_canonical, N_fail, and valid coverage exactly conserve the six input observations.
• Normal and high-depth outputs are byte-identical below saturation for C-only, A/C/T, and explicit same-code selections across 1/2/3/8 threads.
• Reverse canonical, modified, deletion, failure, difference, no-call, and max-depth accounting share one combined anchor.
• Exact duplicate selections are idempotent in optimized, dynamic, and high-depth modes; same-code/different-base selections remain distinct.
• Debug and release CpG combined-C commands emit the same CpG-only site set.
• Existing successful pileup output remains unchanged outside the stated scope.
• Focused tests and the full workspace test suite pass.

Reproduction artifacts

Artifact	Size	SHA-256	Notes
Combined-C SAM	469 B	c84c706df586b76a174041c30a07c69cc11c0744e98391a7ef9211946d7e9868	Six one-base observations from section 1
Explicit-slot SAM	233 B	e26fc7c93037cd9b64522cd832ba911a49dbc8accf3e619411779a268c7b2724	Two ACT reads from sections 2 and 4
Reverse-anchor SAM	516 B	b8a66cd0a56bd778303f3065daebc2d00966bdc9ea1f0690b007e7cd84d8b9ef	Six reverse observations from section 3

Related work

• Bug in modkit pileup with --combine-mods v0.6.4 #644 reports a --combine-mods panic with A/T selections and appears to overlap the compact-slot defect. This report adds minimal fixtures, exact scientific count oracles, high-depth parity, same-code/different-base behavior, reverse ancillary anchoring, and duplicate-selection semantics.
• Discrepancy in calls when running individual modified base pileup and combined modbase pileup #606 discusses combined-versus-separate percentage calculation. Its maintainer-verified example conserves counts and is not evidence for these matrix defects.
• Generic pileup can carry a modification call into the next record #647 / Clear generic pileup state between records #648 correct a separate generic-pileup cross-record state leak in the same subsystem.
• Proposed PR: Correct combined and dynamic pileup count slots #656
• Deferred policy or upstream dependency: none.