`dmr pair` can reorder contigs and split segments at internal batch boundaries

[SuhasSrinivasan]

dmr pair can reorder contigs and split segments at internal batch boundaries

Summary

On accepted multi-contig bedMethyl inputs, modkit dmr pair can emit contig blocks in a different order from the command's already-established lexical contig schedule. When --segment is enabled, a contig can leave and later re-enter the stream, which makes the stateful HMM finish early and split one biologically continuous segment into multiple rows.

With explicit fixed coverage caps, the output depends on the internal --interval-size/--batch-size geometry even when the scored sites, statistical parameters, and scientific inputs are identical. The behavior is deterministic for a given geometry; this report does not claim worker-thread nondeterminism.

Severity

Severity: High — scientific correctness and reproducibility

Rationale: With --segment, accepted input can produce different segment boundaries, row counts, num_sites, aggregate counts, scores, effect statistics, and state continuity solely because an internal batching option changed. The command exits successfully and gives no indication that a contig was artificially closed and reopened. Without segmentation, the scientific values for individual sites are retained but contig row order can still change; that site-only impact is Medium reproducibility/order rather than High.

User and scientific impact

• Affected result or workflow: multi-contig modkit dmr pair, especially --segment output.
• Direction of error: contig reordering/re-entry, premature HMM flushes, and artificial segment splitting; no site loss from this defect itself.
• Likely exposure: data-dependent but ordinary. It occurs when one DmrBatchOfPositions spans a contig transition and its hash-map iteration order differs from the lexical scheduler order.
• Detectability or workaround: inspect site/segment contig blocks for re-entry and compare outputs across interval sizes. Choosing an interval size that happens not to combine contigs can avoid the trigger, but this is not a reliable data-independent workaround.

Affected versions and environment

• Released version: reproduced with modkit 0.6.4.
• Development revision inspected: 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c.
• Deterministic regression/contract revision: d80b35c7b59cb269f9b03519035bf0b94d7fed75, directly atop issue dmr pair --segment omits the final successfully scored site from each HMM chunk #686 repair ba63c29527fde8c25e6f916ebac1fd62e9adb028. It extracts the real handoff into a no-op helper and adds a deliberately scrambled cc, aa, bb unit oracle, so the parent-red result no longer depends on an accidental hash iteration order. Its production behavior remains equivalent to ba63c29; dmr pair --segment omits the final successfully scored site from each HMM chunk #686 restores final-site completeness but does not repair this defect.
• Exact CLI artifact carrier: superseded tests/fixtures commit 0a966eeb214fc206025c93ea48b72a877f735e59, also directly atop ba63c29. The geometry-dependent outputs and hashes below were captured from this exact build and are retained as evidence rather than silently relabeled as artifacts from the rewritten regression commit.
• Operating system and architecture: macOS 26.6, arm64.
• Relevant input format and index: synthetic bgzip/tabix-indexed bedMethyl inputs plus FASTA/FAI.
• Related tools: samtools, bgzip, and tabix 1.23.1.

Steps to reproduce

Minimal input

The fixture contains 6 positive sites on aa, 15 on bb, and 6 on cc. Condition A is fully modified and condition B is fully canonical at every site. Each reference contig contains only C, so every row is unambiguous.

Commands

work_dir="$(mktemp -d)"
cd "$work_dir"

: > a.bed
: > b.bed
for spec in aa:6 bb:15 cc:6; do
    chrom="${spec%%:*}"
    count="${spec##*:}"
    awk -v chrom="$chrom" -v count="$count" 'BEGIN {
        OFS="\t"
        for (i = 0; i < count; i++) {
            print chrom, i, i + 1, "C", 10, "+", i, i + 1, \
                  "255,0,0", 10, "100.00", 10, 0, 0, 0, 0, 0, 0
        }
    }' >> a.bed
    awk -v chrom="$chrom" -v count="$count" 'BEGIN {
        OFS="\t"
        for (i = 0; i < count; i++) {
            print chrom, i, i + 1, "C", 10, "+", i, i + 1, \
                  "255,0,0", 10, "0.00", 0, 10, 0, 0, 0, 0, 0
        }
    }' >> b.bed
done

printf '>aa\nCCCCCC\n>bb\nCCCCCCCCCCCCCCC\n>cc\nCCCCCC\n' > ref.fa
samtools faidx ref.fa
bgzip -c a.bed > a.bed.gz
bgzip -c b.bed > b.bed.gz
tabix -p bed a.bed.gz
tabix -p bed b.bed.gz

modkit dmr pair \
  -a a.bed.gz -b b.bed.gz \
  -o sites-i10.bed --segment segments-i10.bed \
  --ref ref.fa --header --base C --max-coverages 10 10 \
  --batch-size 1 --interval-size 10 \
  --threads 1 --io-threads 1 --suppress-progress --force

modkit dmr pair \
  -a a.bed.gz -b b.bed.gz \
  -o sites-i3.bed --segment segments-i3.bed \
  --ref ref.fa --header --base C --max-coverages 10 10 \
  --batch-size 1 --interval-size 3 \
  --threads 1 --io-threads 1 --suppress-progress --force

awk 'NR > 1 && $1 != previous { print $1; previous = $1 }' sites-i10.bed
awk 'NR > 1 && $1 != previous { print $1; previous = $1 }' sites-i3.bed
cut -f1-3,6 segments-i10.bed
cut -f1-3,6 segments-i3.bed
shasum -a 256 sites-i10.bed segments-i10.bed sites-i3.bed segments-i3.bed

Repeat both commands with --threads 4 --io-threads 2 to distinguish batch geometry from thread scheduling.

Control or independent oracle

The scheduler constructs the input contig queue in Rust lexical order: aa, bb, cc. Within each contig, the positions are already ordered and identical between the two runs. Therefore:

site contig blocks = aa, bb, cc, with no contig re-entry
site counts = aa:6, bb:15, cc:6, total:27
segment rows = aa:[0,6), bb:[0,15), cc:[0,6)
segment num_sites = aa:6, bb:15, cc:6, total:27
interval-size 10 bytes = interval-size 3 bytes
threads 1/1 bytes = threads 4/2 bytes

The final-site cardinality in this oracle is supplied by the independent issue #686 repair. This issue concerns only contig delivery order and the artificial HMM reset.

Observed behavior

The released 0.6.4 binary and the earlier freshly built exact 0a966ee CLI artifact carrier both emit these site contig blocks at interval size 10:

bb
aa
bb
cc

At interval size 3 they emit the expected single-pass blocks:

aa
bb
cc

On that exact artifact carrier, all 27 sites are conserved by #686 but interval size 10 produces:

bb  0   10  10
aa  0    6   6
bb 10   15   5
cc  0    6   6

Interval size 3 instead produces one continuous bb segment:

aa  0    6   6
bb  0   15  15
cc  0    6   6

Exact 0a966ee CLI artifact hashes are:

interval 10 sites:    f039943e7c5c73187be1375fe8720082215dba6481617abcdd61dbd85d0a64f8
interval 10 segments: 4a87999334d6007d394065ad181e1245e7939b1d1b9d8c5ecd5786ac01ba0432
interval 3 sites:     4714c16b976cc7c4df346bf7cf33846274eebca556ec152bbfc22021d053b235
interval 3 segments:  835bacf28ab382f6e002b066f604c6eceb283e6ca1643b470daad4c9871ca766

For each geometry, 1/1 and 4/2 compute/I/O-thread outputs are byte-identical. The geometry-dependent difference is therefore deterministic and not evidence of thread scheduling nondeterminism.

Released 0.6.4 shows the same contig re-entry and artificial flushes, plus the separate final-site omission reported in #686. Its interval-10 site artifact is byte-identical to the exact 0a966ee site artifact, confirming that #686 does not affect the site-order evidence.

Because the source consumes an unordered map, a different binary build of the same unfixed production source happened to emit the three keys in lexical order for this small fixture. That variability is part of the reproducibility defect; the hashes above are attributed only to the cited fresh exact build rather than claimed as universal for every build. The installed eight-contig control independently emits bb, aa, dd, cc, ee, ff, gg, hh instead of lexical aa through hh, reducing reliance on a favorable three-key iteration order.

Expected behavior

• Preserve the lexical contig order already established by SingleSiteBatches when delivering each completed batch.
• Emit each contig in one contiguous block without leaving and later re-entering it.
• Keep site and segment bytes invariant to internal interval/batch geometry when the selected sites and scoring parameters are identical.
• Never finish or reset the HMM merely because a hash map reordered contigs inside one batch.
• Preserve byte identity across supported compute and I/O thread counts.

Root-cause evidence

• SingleSiteBatches::new constructs a lexically sorted contig queue.
• SingleSiteBatches::get_next_batch retains current_batch when update() advances to the next contig, so one valid batch can contain adjacent contigs.
• DmrBatchOfPositions and the organized condition maps use FxHashMap for their outer contig maps.
• process_batch_of_positions consumes the organized map with into_iter() and collects it directly into a vector without restoring scheduler order.
• The receiver sends that vector first to the segmenter and then to the site writer. HmmDmrSegmenter::add finishes the current HMM chunk whenever the next vector entry has a different contig.

Rayon's indexed collect::<Vec<_>>() preserves the order of the batch vector, and the single producer sends completed results sequentially. The unstable order originates inside the per-batch contig map, not in worker completion order.

At exact deterministic regression/contract parent d80b35c, the scrambled-key unit oracle is red with observed order cc, aa, bb rather than lexical aa, bb, cc. The focused integration is also red at the first cross-geometry comparison (site output changed for interval 3, one thread); the same-geometry 1/1-versus-4/2 comparison has already passed at that point. The unit oracle makes the tests-first failure independent of allocator, build, hash, or scheduling behavior.

Proposed fix scope

Restore the scheduler's established lexical order at the final per-batch handoff before the vector can reach either the site writer or the stateful HMM segmenter. The expected implementation is one in-place sort of the unique contig keys after collection.

For C contigs in a batch, this adds O(C log C) string comparisons and O(log C) stack space; C is normally one or a few, and the existing score vectors are moved by handle rather than cloned. No site arithmetic, within-contig position order, floating-point reduction, worker schedule, or I/O operation changes.

Non-goals

• No natural-chromosome or FASTA/header-order policy. This report preserves the Rust lexical order already selected by the scheduler.
• No change to HMM probabilities, transitions, state labels, or maximum-gap behavior.
• No change to the independent final-site decoder defect in dmr pair --segment omits the final successfully scored site from each HMM chunk #686; the exact 15-site segment oracle depends on that correction.
• No change to automatic maximum-coverage sampling. Its sampling budget can have a separate interval/super-batch overshoot policy; this regression supplies explicit caps so only output delivery order changes.
• No thread-pool, channel, batching, or hash-map redesign.
• No claim that the current reproducer varies across thread counts.

Acceptance criteria

• The minimal interval-10 reproducer is red before the ordering repair and green afterward.
• A deterministic unit oracle supplies cc, aa, bb at the real post-collection ordering seam and requires lexical aa, bb, cc; it is red at the regression/contract commit and green after the repair.
• With explicit maximum coverages, interval sizes 10 and 3 produce byte-identical site and segment artifacts from the same selected sites and scoring parameters.
• Compute/I/O configurations 1/1 and 4/2 produce byte-identical artifacts for both interval geometries.
• Site output contains exactly 27 unique positions in lexical aa, bb, cc blocks with no contig re-entry and per-contig counts 6, 15, and 6.
• Segment output contains exactly one row per contig: [0,6), [0,15), and [0,6), with num_sites 6, 15, and 6.
• Existing single-contig DMR output and the issue dmr pair --segment omits the final successfully scored site from each HMM chunk #686 final-site regression remain unchanged.
• Focused DMR tests and the applicable full workspace suite pass.

Reproduction artifacts

Artifact	Size	SHA-256	Notes
dmr_contig_order_a.bed	1,399 bytes	8e19f4d5256f95c9c4e4b3ace52343932cd7d1d20188abc48f28d4ae66b4cf7a	Human-readable fully modified condition A
dmr_contig_order_a.bed.gz	237 bytes	d8afab2481cd3ecdf833f742070d1043f53d0586179585d923d5b0f890f791a3	BGZF condition A
dmr_contig_order_a.bed.gz.tbi	142 bytes	56fb19cf633800b132c85e2d97e9b4dbb324bf695ff2b39d9e4bebfd1e6ef2a3	Tabix index for condition A
dmr_contig_order_b.bed	1,345 bytes	a42024e37e2b4b19c9fb6b208f90a12b78e771b6b88c5b2a677d4acc9ae273f9	Human-readable fully canonical condition B
dmr_contig_order_b.bed.gz	234 bytes	9422da7f32af55929c2c9c98a86c4cc4a8ce876bddec313a9315c3002a1bfb61	BGZF condition B
dmr_contig_order_b.bed.gz.tbi	142 bytes	9f599173c84ae808deb9a9089461052954eb871ecec537c5e16499e5ad85c7f1	Tabix index for condition B
dmr_contig_order.fa	42 bytes	d64e0cbf2234c99951a85805e327ab8c15cca286daf9768b636844a25e9a7ae2	Three-contig reference
dmr_contig_order.fa.fai	38 bytes	f4db8a10d8e0f05577d8981c8fd3c1caf995faa813bb4c84b298713c03f27f82	FASTA index

The human-readable BED and FASTA files completely define the synthetic fixture; the compressed and indexed files are included only so the Rust integration test does not require external tools at test runtime.

Related work

• Issue #686 covers the independent N-1 final-site decoder defect. Its repair is used only to make the segment-count oracle exact; it does not change or repair contig order.
• Issue #448 is broader discussion of missing positions in segmentation and does not cover multi-contig batch order or geometry-dependent HMM resets.
• Proposed PR: to be linked only after the focused implementation passes parent-red, fix-green, thread/geometry, and full-workspace gates and after the dmr pair --segment omits the final successfully scored site from each HMM chunk #686 dependency is available.

[SuhasSrinivasan]

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