`pileup` motif output depends on chunking/preload mode, and invalid `--combine-strands` requests can truncate output

[SuhasSrinivasan]

pileup motif output depends on chunking/preload mode, and invalid --combine-strands requests can truncate output

Summary

Accepted motif-filtered modkit pileup commands can emit different bedMethyl rows when only --interval-size or --preload-references changes. The same ownership defect can silently omit a valid non-CpG combined-strand motif pair while the command exits successfully.

Validation of --combine-strands is also split across processor-specific and post-writer paths. A missing motif or non-palindromic motif can create or truncate the requested output before returning an error, while a palindromic motif whose reverse anchor precedes its forward anchor is accepted by the generic processor but rejected by the optimized processor.

Dense overlapping palindromic motifs expose a related implementation constraint: exact pair closure must not repeatedly refetch and rescan the complete growing owner prefix. A bounded regression protects the correctness repair from turning a long tandem repeat into quadratic repeated scan work.

Severity

Severity: High — scientific correctness and output integrity

Rationale: Affected commands can exit zero while omitting motif sites or accepting an invalid strand-combination geometry, and the result can change solely with execution chunking or FASTA loading mode. These are silent changes to reported genomic coordinates and modification counts. Late validation is separately destructive: even when the command exits nonzero, it can replace an existing result with an empty file.

User and scientific impact

• Affected result or workflow: motif-filtered and non-CpG strand-combined bedMethyl output from modkit pileup.
• Direction of error: silent site omission, interval-dependent site membership, preload-dependent site membership, inconsistent processor acceptance, and destructive truncation on invalid input.
• Likely exposure: option-specific but routine when --motif is used with small intervals, --preload-references, --modified-bases, or --combine-strands.
• Detectability or workaround: comparing large-interval, small-interval, preloaded, faidx-backed, generic, and optimized runs exposes the problem, but a single successful output has no completeness marker. Using a large interval is not a general correctness guarantee. Write invalid-request output to a fresh temporary path to avoid losing an existing file.

Affected versions and environment

• Released version: modkit 0.6.4.
• Development revision: 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c.
• Operating system and architecture used for verification: macOS 26.6 (Darwin 25.6.0), arm64, Apple M2 Ultra.
• Relevant input format and index: coordinate-sorted BAM with BAI; FASTA with FAI.
• Related tool versions: Rust 1.90.0; samtools/HTSlib is needed only to materialize the inline six-base BAM fixture.

Steps to reproduce

Minimal input

The interval/preload failure can be reproduced with this six-base synthetic fixture.

motif-boundaries.sam:

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

motif-boundaries.fa:

>chr1
ACGTCG

The combined-pair and invalid-input reproductions use only public files already tracked in the modkit repository:

tests/resources/CG_5mC_20230207_1700_6A_PAG66026_3c0abf27_oligo_741_adapters_modcalls_0th_sort_10_reads.bam
tests/resources/CG_5mC_20230207_1700_6A_PAG66026_3c0abf27_oligo_741_adapters_modcalls_0th_sort_10_reads.bam.bai
tests/resources/CGI_ladder_3.6kb_ref.fa
tests/resources/CGI_ladder_3.6kb_ref.fa.fai

Run the commands below from the repository root with modkit 0.6.4 on PATH.

Commands

Set up a fresh output directory, materialize the synthetic fixture, and point to the tracked public fixture:

issue_dir="$(mktemp -d)"
public_bam="tests/resources/CG_5mC_20230207_1700_6A_PAG66026_3c0abf27_oligo_741_adapters_modcalls_0th_sort_10_reads.bam"
public_reference="tests/resources/CGI_ladder_3.6kb_ref.fa"
boundary_bam="$issue_dir/motif-boundaries.bam"
boundary_reference="$issue_dir/motif-boundaries.fa"

# Save the two inline blocks above as these files, then index them.
samtools view -b -o "$boundary_bam" "$issue_dir/motif-boundaries.sam"
samtools index "$boundary_bam"
samtools faidx "$boundary_reference"
modkit --version

First vary only interval size and FASTA loading mode for one CG 0 motif:

for interval in 1 2 3 6 100; do
  modkit pileup "$boundary_bam" \
    "$issue_dir/boundary-i${interval}-faidx.bed" \
    --ref "$boundary_reference" --motif CG 0 --modified-bases C:m \
    --no-filtering --interval-size "$interval" \
    --threads 1 --suppress-progress

  modkit pileup "$boundary_bam" \
    "$issue_dir/boundary-i${interval}-preload.bed" \
    --ref "$boundary_reference" --motif CG 0 --modified-bases C:m \
    --no-filtering --interval-size "$interval" \
    --threads 1 --suppress-progress --preload-references
done

wc -c "$issue_dir"/boundary-*.bed
shasum -a 256 "$issue_dir"/boundary-*.bed

Then compare optimized and generic CGCG 0 --combine-strands processing while varying only interval size:

modkit pileup "$public_bam" "$issue_dir/cgcg0-optimized-i1.bed" \
  --ref "$public_reference" --motif CGCG 0 --combine-strands \
  --modified-bases C:m --no-filtering \
  --region oligo_741_adapters:22-62 \
  --interval-size 1 --threads 1 --suppress-progress

modkit pileup "$public_bam" "$issue_dir/cgcg0-optimized-i40.bed" \
  --ref "$public_reference" --motif CGCG 0 --combine-strands \
  --modified-bases C:m --no-filtering \
  --region oligo_741_adapters:22-62 \
  --interval-size 40 --threads 1 --suppress-progress

modkit pileup "$public_bam" "$issue_dir/cgcg0-generic-i1.bed" \
  --ref "$public_reference" --motif CGCG 0 --combine-strands \
  --no-filtering --region oligo_741_adapters:22-62 \
  --interval-size 1 --threads 1 --suppress-progress

modkit pileup "$public_bam" "$issue_dir/cgcg0-generic-i40.bed" \
  --ref "$public_reference" --motif CGCG 0 --combine-strands \
  --no-filtering --region oligo_741_adapters:22-62 \
  --interval-size 40 --threads 1 --suppress-progress

wc -c "$issue_dir"/cgcg0-*.bed
shasum -a 256 "$issue_dir"/cgcg0-*.bed

Finally exercise invalid combine-strand configurations against absent and existing output paths. set +e is intentional so every failing command is observed:

set +e

modkit pileup "$public_bam" "$issue_dir/no-motif-absent.bed" \
  --ref "$public_reference" --combine-strands --no-filtering \
  --region oligo_741_adapters:22-62 \
  --interval-size 40 --threads 1 --suppress-progress
echo "no motif, absent output: $?"

modkit pileup "$public_bam" "$issue_dir/cat0-absent.bed" \
  --ref "$public_reference" --motif CAT 0 --combine-strands --no-filtering \
  --region oligo_741_adapters:22-62 \
  --interval-size 40 --threads 1 --suppress-progress
echo "CAT 0, absent output: $?"

printf 'existing output must remain unchanged\n' > "$issue_dir/no-motif-existing.bed"
printf 'existing output must remain unchanged\n' > "$issue_dir/cat0-existing.bed"
shasum -a 256 "$issue_dir"/*-existing.bed

modkit pileup "$public_bam" "$issue_dir/no-motif-existing.bed" \
  --ref "$public_reference" --combine-strands --no-filtering \
  --region oligo_741_adapters:22-62 \
  --interval-size 40 --threads 1 --suppress-progress
echo "no motif, existing output: $?"

modkit pileup "$public_bam" "$issue_dir/cat0-existing.bed" \
  --ref "$public_reference" --motif CAT 0 --combine-strands --no-filtering \
  --region oligo_741_adapters:22-62 \
  --interval-size 40 --threads 1 --suppress-progress
echo "CAT 0, existing output: $?"

modkit pileup "$public_bam" "$issue_dir/cgcg2-generic.bed" \
  --ref "$public_reference" --motif CGCG 2 --combine-strands --no-filtering \
  --region oligo_741_adapters:22-62 \
  --interval-size 40 --threads 1 --suppress-progress
echo "CGCG 2, generic: $?"

modkit pileup "$public_bam" "$issue_dir/cgcg3-generic.bed" \
  --ref "$public_reference" --motif CGCG 3 --combine-strands --no-filtering \
  --region oligo_741_adapters:22-62 \
  --interval-size 40 --threads 1 --suppress-progress
echo "CGCG 3, generic: $?"

modkit pileup "$public_bam" "$issue_dir/cgcg2-optimized.bed" \
  --ref "$public_reference" --motif CGCG 2 --combine-strands \
  --modified-bases C:m --no-filtering \
  --region oligo_741_adapters:22-62 \
  --interval-size 40 --threads 1 --suppress-progress
echo "CGCG 2, optimized: $?"

wc -c "$issue_dir"/no-motif-*.bed "$issue_dir"/cat0-*.bed \
  "$issue_dir"/cgcg2-generic.bed "$issue_dir"/cgcg3-generic.bed
shasum -a 256 "$issue_dir"/no-motif-*.bed "$issue_dir"/cat0-*.bed \
  "$issue_dir"/cgcg2-generic.bed "$issue_dir"/cgcg3-generic.bed

Control or independent oracle

--interval-size controls work partitioning, and --preload-references controls FASTA access. Neither option changes the motif, reads, thresholds, or biological counting rule. The six-base fixture contains two CG occurrences with one positive and one negative anchor each, so every run must emit the same four exact rows. The fixed output is 200 bytes with SHA-256 a6641668973c23e4d2b4a96fa13dd89320466681c13d5afd6770b191ad703735 for all ten interval/preload combinations.

chr1	1	2	m	1	+	1	2	255,0,0	1	100.00	1	0	0	0	0	0	0
chr1	2	3	m	1	-	2	3	255,0,0	1	100.00	1	0	0	0	0	0	0
chr1	4	5	m	1	+	4	5	255,0,0	1	100.00	1	0	0	0	0	0	0
chr1	5	6	m	1	-	5	6	255,0,0	1	100.00	1	0	0	0	0	0	0

The generic CGCG 0 processor independently emits the same one-row result at interval sizes 1 and 40. The palindromic motif has a forward anchor at offset 0 and its paired reverse anchor at offset 3; both must be discovered and owned atomically by the positive anchor:

oligo_741_adapters	38	39	m	11	.	38	39	255,0,0	11	100.00	11	0	0	1	0	0	0

That row is 71 bytes with SHA-256 588ced97ec269b522ecd8ec09d96a1ea96e8e48f6a9f73f3c496f77756237d57.

For a reverse-complement palindromic motif of length L and forward focus offset f, the paired reverse offset is L - f - 1. CGCG 0 therefore has offsets 0 and 3, and GATC 1 has offsets 1 and 2; both have a representable forward-owned pair. CGCG 2 and CGCG 3 have reverse offsets 1 and 0, respectively, before their forward anchors. Those unsupported geometries must be rejected consistently. CAT 0 is not reverse-complement palindromic, and combining without a motif is undefined.

As an early-validation control, the already-supported multiple-motif rejection (--motif CG 0 --motif GATC 1 --combine-strands) exits 1 without creating an absent output and preserves the 38-byte sentinel above unchanged.

Observed behavior

Every accepted interval/preload command exits 0, but output membership changes on the six-base fixture:

Motif command	FASTA mode	Interval	Rows	Bytes	SHA-256
CG 0	faidx	1	2	100	71120e71b3ef66531fedbf719819e4bd730bb8f2f1a4985d18486c97d816e03f
CG 0	preload	1	0	0	e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855
CG 0	faidx	2	3	150	2bc6ad31df0c92e50008bf90c2760e48b2663c6f98eb32e57702e1b3d25aba1c
CG 0	preload	2	2	100	d7085daf54a206ee522a066e5a942161fef47be84f1013705fa4c43218a1283b
CG 0	faidx/preload	3	4	200	a6641668973c23e4d2b4a96fa13dd89320466681c13d5afd6770b191ad703735
CG 0	faidx/preload	6	4	200	a6641668973c23e4d2b4a96fa13dd89320466681c13d5afd6770b191ad703735
CG 0	faidx/preload	100	4	200	a6641668973c23e4d2b4a96fa13dd89320466681c13d5afd6770b191ad703735

For CGCG 0 --combine-strands, both generic runs emit the 71-byte oracle. The optimized interval-40 run also emits it, but the optimized interval-1 run exits 0 with an empty file (SHA-256 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855).

Invalid-input handling is also processor- and path-dependent:

• With no motif, the generic command exits 1 with cannot combine strands without a motif, but an absent output is created as a zero-byte file and the 38-byte existing sentinel is truncated to zero bytes.
• With CAT 0, the generic command exits 1 with cannot combine strands with a motif that is not a palindrome, but likewise creates or truncates the output. The optimized command instead exits 0 and writes an empty file.
• Generic CGCG 2 exits 0 and writes two rows (135 bytes; SHA-256 175e4baaac5736bbd1f16e3631bad77aa80230832ca30a0e40c5b1ed23816885). Generic CGCG 3 exits 0 with an empty output. The optimized path rejects both with invalid palindromic motif before creating an absent output.
• The pre-run sentinel is 38 bytes with SHA-256 c890f44a8b775c9c416cfebdf74147ff67bbc453c48093fdd58aa1d247bf3911; after the late no-motif or CAT0 error, its checksum is the empty-file SHA-256 above.

The behavior was deterministic across repeated single-threaded runs.

Expected behavior

• Motif site membership, coordinates, counts, and output bytes are invariant to --interval-size and --preload-references for the same biological request.
• Motif matching uses enough clipped flanking reference context to discover every anchor owned by a half-open requested interval, without duplicating anchors between intervals.
• A combined palindromic pair is owned by its positive anchor. Boundary extension includes its paired reverse anchor atomically, including chained overlapping motifs; a selected region or stranded position filter cannot admit an orphan reverse anchor whose positive owner is excluded.
• Generic and optimized processing agree on motif-anchor membership and ordinary canonical/modified counts for supported non-CpG strand-combined motifs, including CGCG 0 and GATC 1. Exceptional reverse ancillary categories are covered separately by issue pileup can omit or split counts for valid combined and explicit modification selections #655 and PR Correct combined and dynamic pileup count slots #656.
• --combine-strands requires exactly one reverse-complement palindromic motif whose reverse anchor does not precede its forward anchor. Missing, multiple, non-palindromic, and negative-displacement configurations return a causal nonzero error before writer construction.
• An invalid request does not create an absent output path and does not modify an existing output file.

Root-cause evidence

• The preloaded FASTA path uses a half-open substring while the faidx path passes the same end to an inclusive fetch API at fasta.rs:72-94. This gives the two access modes different sequence slices.
• Motifs are searched only inside each requested sequence slice, without left/right context or global half-open anchor ownership, at fasta.rs:177-218. Motifs crossing a chunk boundary can therefore disappear or be attributed inconsistently.
• Combined-strand extension examines mask bits at a chunk-end index rather than discovering complete positive/reverse anchor pairs and closing the owner interval over them at fasta.rs:221-283.
• When an overlapping motif chain reaches beyond the fixed lookahead buffer, that extension loop refetches and rescans a growing prefix. A direct fixed-point rewrite can make this happen once per newly reached pair; the repair therefore needs a bounded, forward-only closure scan rather than repeated full-prefix work.
• The optimized non-CpG DynamicAllContext preset is not marked as strand-combined in the processor selection match at subcommand.rs:1025-1042, so generic and optimized execution can use different motif ownership semantics.
• Negative motif displacement is checked only inside the optimized modified-base preset path at subcommand.rs:636-647; the generic path can accept the same unsupported motif.
• File and phased writers are constructed at subcommand.rs:1056-1184 before the later palindrome lookup validation at subcommand.rs:1185-1191 and no-motif feeder validation at interval_chunks.rs:200-211. File::create has already truncated the target when those checks fail.
• Parent-revision regressions reproduce every ownership mismatch. Separate absent/existing-output tests for no motif and CAT 0 fail on the affected parent because the invalid generic request creates or truncates output; CGCG 2/CGCG 3 regressions demonstrate processor inconsistency. The supported CGCG 0 and GATC 1 controls establish that validation is about representable pair ownership, not motif length.

Proposed fix scope

• Give both FASTA backends one checked half-open coordinate contract.
• Fetch clipped flanking context around an owner interval, translate hits to global coordinates, and retain only anchors owned by that interval.
• For strand combination, discover complete palindromic pairs from their positive anchor, extend interval ownership through the paired reverse anchor to a fixed point, and apply region/stranded filtering to the positive owner.
• Compute that fixed point with disjoint, geometrically growing lookahead windows and build the mask only through the exact closure, bounding repeated scan work on dense overlapping motifs.
• Preserve strand-combination semantics for the optimized non-CpG preset so generic and optimized processors consume the same pair mask.
• Normalize and validate combine-strand motifs immediately after parsing: require exactly one palindrome and reject reverse-before-forward anchor displacement before constructing any writer.
• Retain existing downstream palindrome and no-motif checks as defense in depth.
• Document the exactly-one, palindrome, and anchor-order contract in CLI help and the pileup guide.

Non-goals

• No support for multiple motifs in one strand-combined invocation.
• No new ownership convention for reverse-before-forward displacement; those configurations are rejected rather than silently re-anchored.
• No independent change to optimized reverse-strand deletion, filtered, difference, no-call, or max-depth accounting. That Medium defect is addressed by issue #655 and PR #656; this motif work must be composed after that correction and rerun its asymmetric GATC 1 matrix.
• No change to threshold estimation, modification classification, bedMethyl schema, or output ordering.
• No entropy, pileup-hemi, open-chromatin/GPU, or broader owner-span/processing-span performance redesign. A contig-spanning transitive motif chain can still form one large work item; decoupling ownership from processing is separate work.

Acceptance criteria

• The six-base CG 0 fixture is byte-identical across faidx/preloaded reference access and interval sizes 1, 2, 3, 6, and 100.
• Generic and optimized CGCG 0 --combine-strands emit the exact 71-byte ordinary modified-call oracle across interval sizes 1, 2, 3, and 40.
• A region beginning after the positive CGCG 0 owner does not emit an orphan row for its reverse anchor.
• GATC 1 preserves the expected positive/reverse anchor pair across interval boundaries and both FASTA modes.
• A deterministic 512-base CG tandem repeat plus two terminal bases returns all 255 CGCG pairs (510 anchors) in exact coordinate/strand order while requiring no more than 12 sequence fetches and three reference-lengths of fetched bases in either FASTA mode.
• Missing motif, multiple motifs, CAT 0, CGCG 2, and CGCG 3 requests return nonzero consistently in generic and optimized paths before writer construction.
• Every invalid configuration is tested with both an absent output path and a pre-existing sentinel; the absent path remains absent and the sentinel remains byte-identical.
• Existing downstream motif/feeder validation remains in place, supported CpG and non-CpG controls remain unchanged, and CLI/user-guide text states the actual contract.
• Focused motif/FASTA and pileup integration tests pass; cargo test --workspace --all-targets passes with 193 tests passed, 14 ignored, and zero failed.

Reproduction artifacts

Artifact	Size	SHA-256	Notes
motif-boundaries.sam	217 bytes	32413dc27f75a5fdeb0a622fa7909c6044ee2a2fe4cd5cee7ca65853ca800483	Inline two-read, six-base boundary fixture
motif-boundaries.bam	291 bytes	7cb938315a35726ea4ebf7e0a9d381f63bb9143874a8562798543cc1513d3a42	BAM materialized from the inline SAM
motif-boundaries.bam.bai	96 bytes	372169b941872d66ef041e0c1c287462d3d6a7628215878891c73f5b05951647	Synthetic BAM index
motif-boundaries.fa	13 bytes	6e5423d16aec8429b887c7279a5b5945912e3ffaa4546fa580eaf412dc10e3b8	Inline six-base FASTA
motif-boundaries.fa.fai	13 bytes	f8fa22465200794bd9facbf11c8daf763ff02341dc1d2b3d9fc8ac358dba29ed	Synthetic FASTA index
tests/resources/CG_5mC_20230207_1700_6A_PAG66026_3c0abf27_oligo_741_adapters_modcalls_0th_sort_10_reads.bam	5,496 bytes	f32816ce26c0e165cce243d6e614f18c673939d0165b1cd152547eb2ffdff765	Public coordinate-sorted ten-read fixture
tests/resources/CG_5mC_20230207_1700_6A_PAG66026_3c0abf27_oligo_741_adapters_modcalls_0th_sort_10_reads.bam.bai	192 bytes	f662d3220eb99a42e9b9ec1d84f3c6d94924312b2c0a52cb67d04535cb81ca25	BAM index
tests/resources/CGI_ladder_3.6kb_ref.fa	9,606 bytes	7c02f1e6dc3a3dae186ae47696c822efafe849b8c8526bc131237bf7c3a64570	Public reference FASTA
tests/resources/CGI_ladder_3.6kb_ref.fa.fai	1,037 bytes	e2ce9558bb7a0a1898a210e342ba7033e64460802b4602a31dfce2a01347a20e	FASTA index
boundary-i1-preload.bed	0 bytes	e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855	Affected small-interval/preloaded result
Fixed boundary output	200 bytes	a6641668973c23e4d2b4a96fa13dd89320466681c13d5afd6770b191ad703735	Four-row oracle for every interval/preload combination
cgcg0-optimized-i1.bed	0 bytes	e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855	Affected combined-pair omission
cgcg0-generic-i1.bed	71 bytes	588ced97ec269b522ecd8ec09d96a1ea96e8e48f6a9f73f3c496f77756237d57	Exact combined-pair oracle
cgcg2-generic.bed	135 bytes	175e4baaac5736bbd1f16e3631bad77aa80230832ca30a0e40c5b1ed23816885	Two rows from an unsupported negative-displacement motif
Existing-output sentinel	38 bytes	c890f44a8b775c9c416cfebdf74147ff67bbc453c48093fdd58aa1d247bf3911	Becomes the empty-file hash after late no-motif or CAT0 validation

Related work

• Related issue and prerequisite: issue #655 / PR #656 fixes optimized non-CpG combined-strand reverse ancillary counts. The eventual motif PR must be rebuilt after Correct combined and dynamic pileup count slots #656 and rerun its asymmetric filtered/deletion/difference/no-call/max-depth matrix.
• Proposed PR: to be linked after publication.
• Deferred policy questions and follow-up: supporting multiple combined motifs or reverse-before-forward pair ownership is intentionally deferred. Separating motif owner span from processing span for a pathological contig-spanning transitive chain is also outside this correction. PR Correct combined and dynamic pileup count slots #656 is a hard composition dependency for PR publication.