`entropy` motif rows depend on window boundaries and overlapping motifs can collapse incompatible strand pairs

[SuhasSrinivasan]

entropy motif rows depend on window boundaries and overlapping motifs can collapse incompatible strand pairs

Summary

Accepted motif-based modkit entropy runs can omit sites or emit invalid BED intervals when only --window-size changes. The same scanner loses anchors whose motif context lies outside a requested region or a 10,000-base scan chunk, and multiple selected motifs can duplicate one biological anchor or silently collapse incompatible positive/negative partner relationships.

These behaviors change the set and coordinates of entropy windows without any change to the reads, reference, motif, or filtering threshold.

Severity

Severity: High — scientific correctness and output integrity

Rationale: Affected commands can exit zero with plausible but incomplete entropy rows, and some rows extend past the reference end. Window size is an execution/search bound, not a biological site-selection parameter for one-position windows. A single successful output does not indicate that a boundary motif was omitted or counted twice.

User and scientific impact

• Affected result or workflow: motif-selected entropy windows, region summaries, combined-strand entropy, and downstream comparisons based on their coordinates or counts.
• Direction of error: site omission, duplicate biological anchors, invalid/over-wide half-open intervals, region-versus-whole-contig disagreement, and order-dependent collapse of combined-strand partners.
• Likely exposure: routine for small --window-size values, exact/short regions, motifs near a 10,000-base scan seam or contig end, and overlapping motif selections.
• Detectability or workaround: rerunning several window sizes can expose some omissions, but it is not a complete workaround for region, chunk-seam, overlap, or conflicting-partner cases.

Affected versions and environment

• Released version: modkit 0.6.4.
• Development revision: 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c.
• Verification environment: macOS 26.6, arm64, Apple M2 Ultra.
• Input formats: coordinate-sorted BAM with BAI and reference FASTA with FAI.
• Related tool version used to materialize the fixture: samtools 1.23.1.

Steps to reproduce

Window-size and half-open-coordinate fixture

Save the following files as motif.fa and motif.sam.

>chr1
CGTACG

@HD	VN:1.6	SO:coordinate
@SQ	SN:chr1	LN:6
read-0	0	chr1	1	60	6M	*	0	0	CGTACG	??????	MM:Z:C+m?,0,0;	ML:B:C,255,255	MN:i:6	NM:i:0
read-1	16	chr1	1	60	6M	*	0	0	CGTACG	??????	MM:Z:C+m?,0,0;	ML:B:C,255,255	MN:i:6	NM:i:0

Materialize and index the fixture, then vary only --window-size. The explicit threshold isolates motif geometry from automatic-threshold behavior.

samtools faidx motif.fa
samtools view -b -o motif.bam motif.sam
samtools index motif.bam

for window in 1 2 3 6 100; do
  modkit entropy \
    --in-bam motif.bam \
    --out-bed "entropy-${window}.bed" \
    --ref motif.fa \
    --motif CG 0 \
    --num-positions 1 \
    --window-size "$window" \
    --min-coverage 1 \
    --max-filtered-positions 0 \
    --filter-threshold 0 \
    --threads 1 \
    --io-threads 1 \
    --suppress-progress
done

wc -l entropy-*.bed
shasum -a 256 entropy-*.bed

Anchor-only region fixture

The reference and read are AACGAA; the requested BED interval owns the C anchor at coordinate 2 but not the following G context.

>chr1
AACGAA

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

chr1	2	3	anchor-only

samtools faidx anchor.fa
samtools view -b -o anchor.bam anchor.sam
samtools index anchor.bam

modkit entropy \
  --in-bam anchor.bam \
  --out-bed anchor-output \
  --ref anchor.fa \
  --motif CG 0 \
  --num-positions 1 \
  --window-size 1 \
  --min-coverage 1 \
  --filter-threshold 0 \
  --threads 1 \
  --io-threads 1 \
  --regions anchor.regions.bed \
  --prefix anchor \
  --suppress-progress

Conflicting combined-strand motifs

Save the following four-base reference/read as conflict.fa and conflict.sam, then select both CG 0 and CGCG 0:

>chr1
CGCG

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

samtools faidx conflict.fa
samtools view -b -o conflict.bam conflict.sam
samtools index conflict.bam

modkit entropy \
  --in-bam conflict.bam \
  --out-bed conflict.bed \
  --ref conflict.fa \
  --motif CG 0 \
  --motif CGCG 0 \
  --combine-strands \
  --num-positions 1 \
  --window-size 4 \
  --min-coverage 1 \
  --filter-threshold 0 \
  --threads 1 \
  --io-threads 1 \
  --suppress-progress

Here the positive C anchor at coordinate 0 has partner 1 under CG 0 and partner 3 under CGCG 0. There is no single bijective combined-strand interpretation.

Control or independent oracle

CGTACG has exactly two CG occurrences. With one requested position, each forward and reverse anchor owns one one-base half-open output interval. Every window size must therefore produce exactly these four rows:

chr1	0	1	0	+	1
chr1	1	2	0	-	1
chr1	4	5	0	+	1
chr1	5	6	0	-	1

The anchor-only BED region owns coordinate 2. Motif discovery may inspect clipped right context outside the region, but output ownership remains [2,3), so the window file must contain chr1 2 3 0 + 1 and the region must report one successful window.

For overlapping motifs, identical (strand, canonical base, anchor) hits represent one biological position and must be unioned. In combined mode, exact identical pairs may deduplicate, but both positive-to-negative and negative-to-positive relationships must be one-to-one. A conflict must return nonzero before output creation rather than select a partner by motif order or hash insertion order.

Observed behavior

All five whole-contig commands exit zero, but row membership and bytes depend on --window-size:

Window size	Rows	SHA-256	Relevant output
1	2	bc1ed7dfa1fa75c647860a59930d7302396f418deb50461273e14f52364b352d	only the two reverse rows
2	4	5802882211f8cf7bfbf5f73afae797b81381493c440aefbf7cb7b4d1bbf859cb	all anchors, but every row has width 2
3	3	4ae81aa0e8144826be8efd97026df300a5b4ee97b6c489f509cff3abbb200c29	one forward anchor omitted
6	4	5802882211f8cf7bfbf5f73afae797b81381493c440aefbf7cb7b4d1bbf859cb	all anchors, width 2
100	4	5802882211f8cf7bfbf5f73afae797b81381493c440aefbf7cb7b4d1bbf859cb	all anchors, width 2

The emitted rows include chr1 5 7 even though the six-base contig has exclusive end 6.

The anchor-only region exits nonzero with region anchor-only has no valid positions, skipping and creates two empty output files. Whole-contig and region discovery therefore disagree about an anchor owned by the region.

The conflicting-motif command exits zero and emits two rows instead of rejecting the incompatible pair relation:

chr1	0	2	0	+	1
chr1	2	4	0	+	1

The whole-contig behavior was deterministic across repeated one-thread runs. Focused parent regressions additionally reproduce a missed motif spanning the 10,000-base scan seam and an incorrect global offset for motifs first found in a later scan chunk.

Expected behavior

• Entropy motif anchors and output bytes are invariant to --window-size when the requested biological windows are otherwise identical.
• Every one-position row is exactly [p,p+1), including a reverse anchor at the final reference base.
• Motif matching receives clipped context on both sides, translates local hits to global coordinates, and emits each anchor only from its half-open owner interval.
• Whole-contig, nonzero-start region, and anchor-only region discovery agree on all owned anchors.
• Motifs at and across the 10,000-base scan seam are neither omitted nor duplicated.
• Overlapping selected motifs form a deterministic union of biological anchors.
• Combined-strand pairs are bijective. Conflicting mappings fail before creating or modifying output.

Root-cause evidence

• Multi-position intervals are constructed as min..max, while later output code adds another base, mixing inclusive and exclusive endpoint conventions at entropy/mod.rs:73-95 and entropy/mod.rs:521-539.
• Window discovery searches left context but stops the reference slice at curr_position + window_size; right-context motifs can disappear as the window size changes at entropy/mod.rs:1058-1135.
• Initial motif search scans disjoint 10,000-base chunks, does not overlap them for motif context, and returns a chunk-local coordinate without adding the chunk origin at entropy/mod.rs:1187-1201.
• Hits from every selected motif are concatenated without a biological-anchor union, and combined pairs are collected into an FxHashMap, allowing duplicate weighting or overwrite/cardinality collapse at entropy/mod.rs:939-980 and entropy/mod.rs:1090-1133.

Proposed fix scope

• Represent every entropy window and BAM fetch interval once as a global half-open range [min_anchor,max_anchor+1); remove downstream endpoint increments.
• Scan disjoint owner stripes with longest-motif-minus-one clipped context on both sides, translate hits globally, and retain only anchors owned by each stripe/region.
• Produce monotonic strand-specific hit streams with bounded lookahead rather than materializing and cloning all hits for a dense contig.
• Sort and deduplicate identical biological anchors deterministically.
• Validate both directions of the combined-partner bijection over every requested reference owner before writer construction.

Non-goals

• No change to filtering or automatic-threshold semantics.
• No entropy-state representation or final floating-point reduction change.
• No change to the separate pileup motif scanner.
• No claim of transactional rollback after runtime write/flush failure.

Acceptance criteria

• The six-base fixture emits the exact four-row oracle for window sizes 1, 2, 3, 6, and 100.
• One- and multi-position intervals are valid global half-open ranges; a final-base reverse anchor emits exactly [3,4) on a four-base reference.
• Whole-contig, nonzero-start, and anchor-only region runs produce the same owned-anchor map.
• Motifs spanning and following a 10,000-base seam are found at their exact global coordinates once.
• Identical and overlapping motifs deduplicate by (strand, base, anchor); exact combined pairs deduplicate once.
• Conflicting positive-to-negative or negative-to-positive mappings return a deterministic nonzero error before output creation, including a conflict on a later contig.
• Exact hit order and output are stable across scan chunk sizes and thread counts.
• A dense multi-base scan has an explicit live-hit high-water bound determined by requested lookahead and fixed stripe/context settings, not contig length or motif density.
• Existing entropy controls, focused regressions, and the full workspace suite pass.

Reproduction artifacts

Artifact	Size	SHA-256	Notes
motif.sam	215 bytes	7bff3f176b98eb21516c56aabfbe4a94f486f89d91522736d6f5931ea83adb7c	Inline two-read six-base fixture
motif.bam	273 bytes	de4ae3ca28bf992607017820b43b055076fe1afc9e17ba5bf485cf3be8a5925e	samtools 1.23.1 conversion
motif.bam.bai	96 bytes	141d0e5aa09cdf5f5a97ad328ea75b05881c9a4be37d706b193a70d58b501b15	BAM index
motif.fa	13 bytes	13c23b1cbfeeecb091fe1b52b3e6a715fa6224265b0fae90d487d85a155a78bc	Inline reference
motif.fa.fai	13 bytes	f8fa22465200794bd9facbf11c8daf763ff02341dc1d2b3d9fc8ac358dba29ed	FASTA index
anchor.sam	122 bytes	a6024e8e4849097b0864a525e80b3df8df84a4c82ad232f4f7a7cf3772742f93	Anchor-only region fixture
anchor.regions.bed	21 bytes	18fe0d1321077378e6c3d03fd8fbb115079c966daa5c8fed08532a7205b87555	[2,3) owner interval
conflict.sam	124 bytes	c3f4fb6dd716a9895217f6e41364f9a2f603a46782958fd1e94c89a20576a20b	Conflicting combined-partner fixture

Related work

• Related issue: the analogous pileup motif-ownership path is separate because entropy uses its own in-memory reference scanner.
• Proposed PR: to be linked after the focused geometry/ownership commits are rebuilt on current upstream and fresh gates pass.
• Deferred work: zero-valued resource validation, entropy state cardinality, final-bit floating reduction, automatic filtering, and broad runtime lifecycle/output transactionality remain separate review units.

[SuhasSrinivasan]

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