Seeded indexed fractional sampling changes with worker and interval settings

[SuhasSrinivasan]

Summary

Seeded indexed fractional sampling is not reproducible across worker scheduling and interval geometry. Modern probability collectors use deterministic but mutable RNG streams local to each worker, so the same alignment can receive a different decision when it is assigned to another worker or refetched in another interval. The legacy mapped indexed path also fails to apply the user seed consistently, and sparse retained intervals can use the wrong predecessor boundary.

As a result, changing only --threads, --interval-size, or --sampling-interval-size can change the sampled population, probability histograms, automatic thresholds, command output, or even whether threshold estimation succeeds.

Severity

Severity: High — scientific reproducibility

Rationale: A fixed input, sampling fraction, and explicit seed do not identify a fixed sampled population. Downstream calls and summary statistics can change with internal scheduling rather than biological data or a user-visible statistical choice.

User and scientific impact

• Affected workflows include indexed fractional sampling in sample-probs, Summary, standard pileup threshold estimation, call-mods, and duplex/legacy automatic-threshold consumers.
• The direction is not predictable: different runs can include different alignments, change per-base histograms and thresholds, and thereby change pass/fail classification.
• Long alignments fetched in more than one interval and sparse BED/motif selections are especially sensitive to ownership errors.
• Repeating a command with one thread is not a general workaround because interval geometry, repeated fetch ownership, and the legacy path's seed handling can still matter.

Affected versions and environment

• Released version: modkit 0.6.4.
• Development revision: 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c.
• Reproduced on macOS arm64 with indexed BAM and the repository's ten-read fixture.
• Reproduced with workers 1 and 8, intervals 20 and 1000, seed 7, and fraction 0.5.

Steps to reproduce

From the repository's modkit package directory:

tmp=$(mktemp -d /tmp/modkit-fraction.XXXXXX)
input=../tests/resources/bc_anchored_10_reads.sorted.bam

modkit sample-probs "$input" --sample-frac 0.5 --seed 7 \
  --threads 1 --interval-size 20 --hist --out-dir "$tmp" \
  --prefix w1_i20 --force --suppress-progress

modkit sample-probs "$input" --sample-frac 0.5 --seed 7 \
  --threads 8 --interval-size 20 --hist --out-dir "$tmp" \
  --prefix w8_i20 --force --suppress-progress

modkit sample-probs "$input" --sample-frac 0.5 --seed 7 \
  --threads 8 --interval-size 1000 --hist --out-dir "$tmp" \
  --prefix w8_i1000 --force --suppress-progress

for name in w1_i20 w8_i20 w8_i1000; do
  printf '%s rows=' "$name"
  awk 'END { print NR - 1 }' "$tmp/${name}_probabilities.tsv"
  shasum -a 256 "$tmp/${name}_thresholds.tsv" \
    "$tmp/${name}_probabilities.tsv"
done

The legacy path can be checked with the same fixture:

modkit call-mods "$input" "$tmp/call-w1.sam" --output-sam \
  --sampling-frac 0.5 --seed 7 --threads 1 \
  --sampling-interval-size 20 --filter-percentile 0.25 \
  --suppress-progress

modkit call-mods "$input" "$tmp/call-w8.sam" --output-sam \
  --sampling-frac 0.5 --seed 7 --threads 8 \
  --sampling-interval-size 1000 --filter-percentile 0.25 \
  --suppress-progress

awk '!/^@/' "$tmp/call-w1.sam" | shasum -a 256
awk '!/^@/' "$tmp/call-w8.sam" | shasum -a 256

Control or independent oracle

For a Bernoulli fraction sampler, one alignment identity plus one master seed must produce one immutable decision. Worker number, completion order, interval width, and repeated interval fetches are not statistical inputs and must not affect membership. Fraction 0 must include none and fraction 1 must include all.

Observed behavior

One installed-0.6.4 run produced:

Case	Probability rows	Threshold SHA-256	Probability SHA-256
workers 1, interval 20	50	d1295a2a354a5166a54f53a3edbafde34b5bfbf4981c3ea744a9c153a548461a	31b57ff26efa889ccac91da8c3ede89990fe937f327e34ad5918293410b56ac3
workers 8, interval 20	66	81ec4043bda3b0c7a84522379b56cafd4ddb9192f0e55ad96344661187b7395d	8b7bae3034f8cc61a0d900d9bd7ef7d145b313ca58eeb1b972f3390f34c724a8
workers 8, interval 1000	60	c39577b69c3def5925fa3f39d53a78e44b5ca93564a1dd3e8298e55e9609e336	264dc2a57500ead8154058860f29bba9984dd04554f94f907756388506465f76

The two legacy commands both wrote ten records, but their header-normalized SAM bodies had different SHA-256 values (29cbf9af... versus 78cf8d33...) because their sampled thresholds differed.

On the repository's one-site include-BED control, fraction 0.5, seed 7, and sampling interval 20 failed with one datapoint, while interval 1000 succeeded with threshold 0.5839844. A separate spanning/gap fixture shows that retained sparse intervals can otherwise omit gap-starting reads or process spanning reads more than once.

Expected behavior

For indexed fractional sampling, the same file, explicit seed, fraction, and scientific selectors should produce the same alignment-level membership across worker counts, interval sizes, repeated processes, and interval refetches. Different seeds should normally diverge. Sparse selectors must retain every eligible alignment exactly once.

Root-cause evidence

The modern implementation initializes mutable RNG state per worker. Scheduling therefore determines which random stream is applied to a record. The legacy mapped indexed implementation derives interval schedules without preserving the explicit Bernoulli seed contract. Its ownership cutoff is also tied to physical interval traversal rather than the previous retained interval when sparse selectors remove intervening chunks.

These are upstream defects, not Rayon defects: parallel reduction combines already different samples correctly.

Proposed fix scope

• Resolve one master seed for the complete indexed fractional job.
• Use a versioned stateless score over stable alignment identity fields and compare it with an exact integer fraction threshold.
• Give every worker the same sampler and include QNAME in selective CRAM decoding.
• Keep fraction 0 and 1 exact and reject non-finite/out-of-range fractions before starting workers.
• Filter and annotate retained intervals serially before parallel execution so each retained interval receives the previous retained end.
• Update seed help to distinguish stable indexed alignment decisions from seeded encounter-order streaming on unindexed input.

Non-goals

• Do not choose or change indexed fixed-count sampling semantics; that is the separate policy question in issue Decide indexed --num-reads semantics: reference-span-balanced approximate sampling vs exact record-global sampling #654.
• Do not promise identical membership after reference-header reordering or arbitrary BAM/CRAM re-encoding.
• Byte-identical duplicate alignments may share a decision because no portable stable occurrence ordinal exists.
• Do not change unmapped-record eligibility, same-QNAME aggregation, or coordinate-selector policy in this repair.

Acceptance criteria

• Per-alignment inclusion decisions agree across workers 1/2/3/8, intervals 20/100/1000, repeats, and interval traversal order for one seed. Exact threshold/probability output agrees on the frozen unique-QNAME fixture and whenever the post-aggregation probability population is identical.
• sample-probs, Summary, standard pileup, call-mods, and duplex threshold consumers are covered.
• Long records fetched in multiple intervals receive one decision and contribute complete calls once.
• Sparse BED/motif regressions cover both a spanning record and a record beginning in a gap between retained intervals.
• Fraction 0 includes none; fraction 1 matches the all-record control; distinct seeds diverge on a nondegenerate fixture.
• Existing unindexed encounter-order and fixed-count behavior remains unchanged.
• End-to-end runtime and peak RSS are measured on a representative indexed direct-RNA modBAM.

Reproduction artifacts

Artifact	SHA-256	Notes
tests/resources/bc_anchored_10_reads.sorted.bam	4441acbe1ad59caf6ab5d56bf18097d8999c08351395219c919cb863652d85f0	Ten mapped reads with modification probabilities
tests/resources/bc_anchored_10_reads.sorted.bam.bai	e2a74713867df4a01785d614e066fd310fea08293fc0b86488ec336d89fe0c35	BAM index
tests/resources/include-pos-1-site.bed	repository fixture	Sparse-selector control

Related work

• Issue Decide indexed --num-reads semantics: reference-span-balanced approximate sampling vs exact record-global sampling #654 concerns fixed-count estimands and is intentionally not addressed here.
• Issue Probability sampling applies --edge-filter inconsistently across execution paths #667 concerns edge-filter population consistency, not sampling membership.
• Issue Probability sampling can hang or return partial histograms after reader and worker failures #677 concerns fatal probability-worker/reader errors, not deterministic selection.
• A separate explicit-canonical reduction defect can change thresholds even after membership is stabilized.
• Distinct alignments sharing one QNAME are governed by a separate legacy aggregation-key/order contract. This repair stabilizes their per-alignment sampling decisions but does not redefine downstream raw-QNAME deduplication.