Generic pileup can carry a modification call into the next record

[SuhasSrinivasan]

Summary

When the generic pileup worker processes an MM group with no listed modification calls immediately after a record that does contain a modification call, it can reuse the preceding record's cached modification state. This attributes the earlier call to the second read.

Severity

Severity: High — scientific correctness

Rationale: The command accepts both records and exits successfully, but inflates modified counts and reduces canonical counts. The error is silent in ordinary bedMethyl output and can directly change estimated modification fractions.

User and scientific impact

• Affected workflow: Generic fallback modkit pileup processing.
• Direction of error: Modified counts can be inflated while canonical counts are reduced.
• Exposure: Data-dependent; reached when an empty explicit/implicit MM group follows a record with a call at the same processed position.
• Detectability: There is no diagnostic. Detection requires an independent per-read accounting oracle.

Affected versions and environment

• Reproduced with modkit 0.6.4.
• Present at assessed source revision 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c.
• Reproduced on macOS/Apple silicon with an indexed BAM generated by samtools.

Steps to reproduce

Minimal input

Save as two-records.sam:

@HD	VN:1.6	SO:coordinate
@SQ	SN:chr1	LN:5
modified	0	chr1	1	60	5M	*	0	0	AAAAA	IIIII	MM:Z:A+a.,4;	ML:B:C,255
canonical	0	chr1	1	60	5M	*	0	0	AAAAA	IIIII	MM:Z:A+a.;	ML:B:C

Commands

samtools view -b -o two-records.bam two-records.sam
samtools index two-records.bam
modkit pileup two-records.bam pileup.bed \
  --no-filtering \
  --threads 1 \
  --io-threads 1 \
  --interval-size 100
status=$?
awk '$2 == 4' pileup.bed
printf 'exit=%s\n' "$status"

Control or independent oracle

At query/reference position 4, the first read contributes one explicit modified-A call. The second read's empty implicit-mode group contributes one canonical A. Therefore:

coverage = 2
modified = 1
canonical = 1
percent_modified = 50.00

Observed behavior

The command exits 0, creates pileup.bed, and reports 100% modified with two modified calls and no canonical calls:

chr1	4	5	a	2	+	4	5	255,0,0	2	100.00	2	0	0	0	0	0	0

Expected behavior

The command should exit 0 and create the same output except that the final position must contain one modified and one canonical call:

chr1	4	5	a	2	+	4	5	255,0,0	2	50.00	1	1	0	0	0	0	0

Root-cause evidence

• update_mods_iter2 populates mod_pos, canonical_base, pos_base_mod_call, and mod_strand when the iterator yields a call.
• When the iterator returns None, those fields are left unchanged and can retain values from the preceding record.
• The focused regression on the parent revision reports n_modified = 2; clearing the exhausted state produces the exact 1/1 oracle.

Proposed fix scope

Clear every per-record cached modification field when the iterator is exhausted, and add a worker-level two-record regression for the populated-to-empty transition.

Non-goals

• Does not redesign MM/ML parsing or generic-pileup data structures.
• Does not change optimized/high-depth pileup behavior.
• Does not introduce a CPU-performance refactor.

Acceptance criteria

• The focused regression is red on the assessed parent and green on the proposed fix.
• Position 4 has coverage 2, modified 1, canonical 1, and 50.00% modified.
• Every other count category remains zero and unrelated output remains unchanged.
• Focused, library, and workspace tests pass.

Related work

• Proposed fix: Clear generic pileup state between records #648