validate can omit later truth sets and eligible observations

[SuhasSrinivasan]

Summary

modkit validate can silently omit accepted observations in two related ways: when one BAM is paired with more than one truth BED, only the first BED sees the BAM records; and within a truth set, eligible sites without a truth-overlapping emitted modification call can be absent from the No Call, mismatch, and deletion totals.

Adjacent fallback-classification cases also treat CIGAR reference skips as deletions, discard a whole record when one aligned query base is ambiguous, allow edge-filtered sites to re-enter as fallback observations, and fail when a valid input contains no confident calls. These behaviors reproduce with modkit 0.6.4 and current upstream revision 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c.

Severity

Severity: High — scientific correctness and reproducibility

Rationale: The command commonly exits successfully while omitting eligible truth observations, changing contingency-table denominators and accuracy statistics. Reusing one BAM with multiple truth BEDs can make the result depend on truth-set order. The omissions are not evident from the exit status and can be mistaken for properties of the data.

User and scientific impact

• Affected result or workflow: modkit validate full contingency tables, No Call/mismatch/deletion totals, raw and filtered accuracy, and multi-truth-set comparisons.
• Direction of error: silent undercounting, misclassification, order dependence, or failure on a valid all-No-Call population.
• Likely exposure: routine when one BAM is paired with multiple truth BEDs; data-dependent for sparse explicit MM groups, CIGAR N, ambiguous aligned bases, or --edge-filter.
• Detectability or workaround: run every BAM/BED pair separately and independently audit every truth-overlapping alignment. There is no complete output-only check for the within-record omissions.

Affected versions and environment

• Released version: modkit 0.6.4
• Development revision: 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c
• Operating system and architecture: macOS 26.6, arm64
• Input format: coordinate-sorted BAM with valid MM/ML/MN tags and BED6 truth annotations
• Related tools: samtools 1.23.1 / HTSlib 1.23.1

Steps to reproduce

The fixtures below are synthetic and contain no private data.

1. One BAM paired with two truth BEDs

Create split.sam:

@HD	VN:1.6	SO:coordinate
@SQ	SN:chr1	LN:2
both-sites	0	chr1	1	60	2M	*	0	0	CC	??	MM:Z:C+m?,0,0;	ML:B:C,255,255	MN:i:2	NM:i:0

Create truth-a.bed and truth-b.bed:

# truth-a.bed
chr1	0	1	m	.	+

# truth-b.bed
chr1	1	2	-	.	+

Build the BAM and run both input orders:

samtools view -b -o split.bam split.sam

modkit validate \
  --bam-and-bed split.bam truth-a.bed \
  --bam-and-bed split.bam truth-b.bed \
  --canonical-base C --filter-threshold 0 --threads 1 \
  --suppress-progress --out-filepath split-ab.tsv

modkit validate \
  --bam-and-bed split.bam truth-b.bed \
  --bam-and-bed split.bam truth-a.bed \
  --canonical-base C --filter-threshold 0 --threads 1 \
  --suppress-progress --out-filepath split-ba.tsv

Observed with modkit 0.6.4:

• The log reports Processed 1 mapping records for the first BED and Processed 0 mapping records for the second.
• split-ab.tsv is byte-identical to validating truth-a.bed alone: SHA-256 fd0cb870360160f7c8d53b6d24c1fe065e4d810f561ce195e99495422faa3754.
• split-ba.tsv is byte-identical to validating truth-b.bed alone: SHA-256 e0a57a68351324ec066a417229d212a06a09f46fc6e0315405d7ce66bf843f6.

The first truth BED therefore determines which ground-truth row survives.

Expected: both truth sets are evaluated independently against the complete BAM. Both input orders produce the same 326-byte report, containing one called m observation for ground-truth m and one called m observation for ground-truth canonical C.

2. Eligible unseeded sites are omitted from No Call totals

Create unseeded.sam:

@HD	VN:1.6	SO:coordinate
@SQ	SN:chr1	LN:2
called-1	0	chr1	1	60	2M	*	0	0	CC	??	MM:Z:C+m?,1;	ML:B:C,255	MN:i:2	NM:i:0
called-2	0	chr1	1	60	2M	*	0	0	CC	??	MM:Z:C+m?,1;	ML:B:C,255	MN:i:2	NM:i:0
called-3	0	chr1	1	60	2M	*	0	0	CC	??	MM:Z:C+m?,1;	ML:B:C,255	MN:i:2	NM:i:0
uncalled-1	0	chr1	1	60	2M	*	0	0	CC	??	MM:Z:C+m?,0;	ML:B:C,255	MN:i:2	NM:i:0
uncalled-2	0	chr1	1	60	2M	*	0	0	CC	??	MM:Z:C+m?,0;	ML:B:C,255	MN:i:2	NM:i:0
uncalled-3	0	chr1	1	60	2M	*	0	0	CC	??	MM:Z:C+m?,0;	ML:B:C,255	MN:i:2	NM:i:0
uncalled-4	0	chr1	1	60	2M	*	0	0	CC	??	MM:Z:C+m?,0;	ML:B:C,255	MN:i:2	NM:i:0
uncalled-5	0	chr1	1	60	2M	*	0	0	CC	??	MM:Z:C+m?,0;	ML:B:C,255	MN:i:2	NM:i:0
uncalled-6	0	chr1	1	60	2M	*	0	0	CC	??	MM:Z:C+m?,0;	ML:B:C,255	MN:i:2	NM:i:0

Create unseeded.bed:

chr1	1	2	m	.	+

samtools view -b -o unseeded.bam unseeded.sam

modkit validate \
  --bam-and-bed unseeded.bam unseeded.bed \
  --canonical-base C --filter-threshold 0 --threads 1 \
  --suppress-progress --out-filepath unseeded.tsv

head -n 1 unseeded.tsv

Observed with modkit 0.6.4:

full_contingency_table: [["ground_truth_label","m"],["m",3]]

Expected:

full_contingency_table: [["ground_truth_label","m","No Call"],["m",3,6]]

Additional exact classification cases

Accepted case	Observed behavior	Expected behavior
1M3N1M across five truth positions	Three CIGAR N positions are counted as deletions	Reference-skip positions contribute no molecular observation; matched 1M3D1M still contributes three deletions
One aligned record containing modified C, canonical C, N, IUPAC R, and A mismatch	Parsing R discards the record	Exact categories are modified=1, canonical=1, No Call=2, mismatch-A=1
A truth-aligned fallback site removed by normal or inverted --edge-filter	It re-enters as No Call or mismatch	It contributes no observation; a deletion without query position retains the existing explicit policy
Three records with an eligible empty explicit descriptor C+m?;	The complete table is printed, then the command fails with No minimum value found	Success with No Call=3 and undefined metrics represented as NA
Filtering removes every confident call	Human cells and filtered accuracy can contain NaN	Undefined percentages are NA; defined values remain numeric

Each case has a focused parent-failing regression and exact category oracle in the proposed patch.

Expected behavior

• Every BAM/BED pair is evaluated against a fresh complete BAM traversal, independent of pair order.
• Every truth-overlapping site eligible under its MM descriptor and edge filter is classified exactly once as canonical, modified, No Call, mismatch, deletion, or intentionally omitted reference skip.
• Ambiguous aligned bases affect only their own sites.
• Valid zero-confident-call populations succeed with a complete machine-readable schema and NA for undefined metrics.

Root-cause evidence

• The BAM reader is created once per path outside the truth-BED loop, then consumed repeatedly at subcommand.rs:967-994.
• Fallback traversal is seeded only from emitted truth-overlapping calls at subcommand.rs:309-374, so empty, sparse, and opposite-strand descriptors can leave no eligible strand to traverse.
• Missing query positions are all classified as deletion and ambiguous bases propagate a record-level parse error at subcommand.rs:379-405.
• Ground-truth balancing requires a nonempty called population at subcommand.rs:624-637.

Proposed fix scope

• Open a fresh BAM reader for each BAM/BED pair.
• Retain private normalized MM-descriptor strand eligibility, including empty and N groups, before positional expansion.
• Classify fallback sites with reference-skip, alignment orientation, ambiguous-base, and edge-filter semantics applied site-locally.
• Render zero-denominator statistics as unavailable and allow an empty confident-call population to complete successfully.

Non-goals

• No change to ground-truth balancing, threshold selection, read-filter eligibility, or supported MM/ML syntax.
• No change to deletion-without-query-position policy under edge filtering.
• No parallel scheduling, I/O architecture, or performance redesign.
• Strict truth-BED parsing and conflicting-label rejection are tracked separately.

Acceptance criteria

• The two minimal reproducers fail on the affected parent and pass after the fix.
• Two truth BEDs using one BAM equal the union of separate invocations in both input orders.
• Exact m/No Call/mismatch/deletion totals match the stated oracles for explicit, implicit, empty-descriptor, both-strand, reference-skip, ambiguous-base, and edge-filter fixtures.
• All-No-Call and all-filtered reports retain the complete seven-line machine schema with NA only where the denominator is undefined.
• Existing successful validate output remains unchanged outside the stated scope.
• Focused tests and the full workspace test suite pass.

Reproduction artifacts

Artifact	Size	SHA-256	Notes
split.sam	124 B	c59dbb79a197258e693c3306fe16d73284cf8f89821c079eeda0de8e0e9cd967	One two-site explicit modification record
truth-a.bed	15 B	6f0ea56683eef12e8b03b3ad02572921aca474b58732645489b61063acd301b8	Modified truth at position 0
truth-b.bed	15 B	94ca61bfba6a303e510d6e658ce9d5c0fc3c029697a382c002470c4e898c6cc0	Canonical truth at position 1
unseeded.sam	720 B	4ca69eb4a652c6a91e306219d6fc6a478cfd15d4eccf01e90d8632d4e42d21c3	Three called and six eligible unseeded records
unseeded.bed	15 B	5696c154d255eba8bcc6d65080fd2d9f448e6340bf1b6978e094aab5a415620b	One modified truth site

Related work

• Related issue: validate silently accepts invalid or conflicting truth BED rows #658 covers strict truth-BED parsing and conflicting-label rejection as a distinct input-validation review surface.
• Proposed PR: Correct validate truth parsing and observation accounting #659
• Deferred policy or upstream dependency: none.

[SuhasSrinivasan]

The confirmed validate fix is now included in the consolidated, ready-for-review PR #659. Closing this standalone issue in line with the maintainer's request to keep the identified bug and its code change together in the PR.