Decide indexed --num-reads semantics: reference-span-balanced approximate sampling vs exact record-global sampling

[SuhasSrinivasan]

Summary

Two indexed fixed-count paths are unstable. Modern probability consumers derive an allowance independently for each processing interval, while legacy BAM/CRAM threshold and extract consumers allocate a schedule using different index information. Changing only internal interval geometry can therefore change selected alignments, exceed or underfill the requested approximate maximum, and alter downstream thresholds and output.

Repairing that instability requires an explicit statistical contract: preserving reference-span-balanced approximate sampling and making it interval-invariant is not equivalent to selecting an exact record-global top N.

Severity

Severity: High — scientific reproducibility and statistical contract

Rationale: Accepted inputs can produce different sampled populations, thresholds, counts, and transformed output when only interval size, reference layout, or BAM-versus-CRAM scheduling changes. Direct-RNA transcript abundance can be highly uneven, so changing from reference-span weighting to occurrence-global weighting also changes the scientific estimand.

User and scientific impact

• Affected result or workflow: Summary, sample-probs, pileup, extract, and legacy automatic-threshold consumers using fixed-count indexed sampling.
• Direction of error: unstable membership, over- or under-selection, changed thresholds/output, and format-dependent exhaustion behavior.
• Likely exposure: fixed-count sampling on indexed BAM/CRAM, especially uneven transcriptomes, multiple contigs, motifs, sparse BEDs, or region-restricted workflows.
• Detectability or workaround: repeat across interval sizes and formats and compare normalized alignment fingerprints/histograms. Fractional sampling is a separate path and is not the fixed-count workaround proposed here.

Affected versions and environment

• Released version: modkit 0.6.4.
• Development revision: 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c.
• Operating system and architecture: macOS 26.6, arm64.
• Input formats and indexes: coordinate-sorted BAM/BAI and equivalent externally referenced CRAM/CRAI.
• Fixture tools: Python 3 and samtools 1.23.1 / HTSlib 1.23.1.
• Reproduced with workers 1 and 8, multiple interval sizes, AllPositions, motif, sparse BED, region, multi-contig, and mixed mapped/unmapped inputs.

Steps to reproduce

Save the following standalone generator as make_fixture.py:

from pathlib import Path

out = Path("/tmp/modkit-fixed-count")
out.mkdir(parents=True, exist_ok=True)
qualities = [144, 160, 176, 192, 208, 224, 240, 248]
with (out / "fixture.sam").open("w") as writer:
    writer.write("@HD\tVN:1.6\tSO:coordinate\n")
    writer.write("@SQ\tSN:chr1\tLN:160\n@SQ\tSN:chr2\tLN:160\n")
    for tid in range(2):
        for index in range(4):
            probability = qualities[tid * 4 + index]
            fields = [
                f"mapped-{tid}-{index}", "0", f"chr{tid + 1}", "6",
                "60", "80M", "*", "0", "0", "C" * 80, "?" * 80,
                "MM:Z:C+m?," + ",".join(["0"] * 80) + ";",
                "ML:B:C," + ",".join([str(probability)] * 80),
                "MN:i:80",
            ]
            writer.write("\t".join(fields) + "\n")

Build and index the coordinate-sorted fixture, then run the same sampling request with only the processing interval changed:

modkit --version
python make_fixture.py
samtools view -b -o /tmp/modkit-fixed-count/fixture.bam \
  /tmp/modkit-fixed-count/fixture.sam
samtools index /tmp/modkit-fixed-count/fixture.bam
cd /tmp/modkit-fixed-count

modkit sample-probs fixture.bam --num-reads 3 --threads 1 \
  --io-threads 1 --interval-size 20 --hist --out-dir . \
  --prefix i20 --force --suppress-progress

modkit sample-probs fixture.bam --num-reads 3 --threads 1 \
  --io-threads 1 --interval-size 160 --hist --out-dir . \
  --prefix i160 --force --suppress-progress

for file in i20_probabilities.tsv i160_probabilities.tsv; do
  echo "$file"
  awk -F '\t' '
    NR == 1 { for (i = 1; i <= NF; i++) column[$i] = i; next }
    $column["code"] == "m" && $column["primary_base"] == "C" &&
      $column["count"] > 0 {
        print $column["range_start"], $column["count"]
      }
  ' "$file"
done

Control or independent oracle

Each alignment has a unique raw ML value repeated at all 80 cytosines. Each nonzero histogram bin therefore identifies one contributing alignment (range_start = ML/256), while count 80 proves that the selected alignment contributed completely. Repeating at workers 1 and 8 separates interval-geometry dependence from worker scheduling.

Observed behavior

On modkit 0.6.4, interval 20 shows only ML 144 and 208:

0.5625 80
0.8125 80

Interval 160 instead shows ML 144, 160, 208, and 224:

0.5625 80
0.625 80
0.8125 80
0.875 80

Both worker counts reproduce the corresponding interval result. Broader frozen regressions reproduce different membership/counts for motif, sparse BED, region, multi-contig, and mixed mapped/unmapped modes.

A separate legacy six-occurrence BAM/CRAM fixture shows that requested N=6 or N=99 can retain only one of six extract occurrences and give automatic thresholding only one datapoint. An all-unmapped mapped-only BAM errors with zero reads found in bam index, while the equivalent CRAM succeeds with empty output.

The modern implementation computes an inflated estimate from requested count, interval width, and total header/reference length, then resets and reapplies that allowance per worker interval. Empty header contigs, region length, focus pruning, and mapped/unmapped fallback can therefore change selection independently of the user-visible count.

Expected behavior

Changing internal interval size or worker count must not grossly alter the sampled population or cause the command-level count to be reapplied independently for every interval. The normative weighting/selection contract must be chosen explicitly:

1. preserve reference-span-balanced approximate sampling and repair its interval dependence;
2. retain that default and add an explicit exact record-global strategy; or
3. replace the default with an exact reference-stratified strategy in a compatibility-breaking release.

Root-cause evidence

The modern fixed-count allowance is derived per interval and state is recreated for each processing call. Parallel reduction only adds completed interval results; it has no immutable global selection plan or deficit redistribution. Equal-probability record-global stable top N would fix exactness but would change the estimand from reference-coordinate span to observed alignment-occurrence abundance.

The legacy BAM scheduler can read raw mapped-alignment counts from BAI/CSI, while CRAI lacks those counts. Proving the same exhausted-sample condition for CRAM requires a record scan and second pass in the exhausted case. Raw index counts also do not define post-filter consumer eligibility.

Equal stable identities create another explicit contract: duplicate alignment occurrences must retain multiplicity and deterministic same-file tie order, while byte-identical membership cannot reasonably be promised after BAM/CRAM conversion or re-encoding.

All listed regressions fail on the affected parent. Separate local repairs for fractional sampling and probability aggregation do not resolve this fixed-count estimand question.

Proposed fix scope

The recommended lowest-risk contract is option 2:

• Preserve existing fractional sampling behavior.
• Preserve the approximate fixed-count mode's chosen reference-span estimand while making its plan interval- and worker-invariant, with a documented approximation tolerance.
• Add an explicit exact same-file sampling strategy.
• Define its selection unit as an alignment occurrence rather than raw QNAME and define consumer-specific eligibility before promising exactness.
• Promise exact min(N, eligible) only for the same indexed input after eligibility is defined.
• Include mapped and unmapped records in one budget only for unconstrained Summary; exclude unmapped records when a coordinate selector is active.
• Preserve established handling of zero-contribution records and CIGAR N/D spans unless separately changed.

For the legacy path, please confirm whether a conservative BAM-only raw-count fast path is acceptable or whether format-independent exhaustion decisions are required. If format agreement is required, both formats need the same bounded record-level predicate; the BAM index can remain an additional sufficient shortcut.

Non-goals

• Do not conflate distinct alignments by raw QNAME in a new occurrence-level strategy.
• Do not promise identical membership across BAM/CRAM conversion or arbitrary re-encoding.
• Do not change unmapped modification-class preservation or coordinate-filter policy.
• Do not hide a change in transcript/reference weighting inside a scheduler refactor.
• Do not claim exact consumer eligibility from BAM/CRAM index counts.

Acceptance criteria

For either modern strategy:

• Normalized scientific membership/fingerprints are invariant across workers 1 and 8 and narrow/wide internal intervals for the same indexed input.
• AllPositions, motif, sparse BED, region, and multiple-contig controls are covered.
• Every selected occurrence contributes every call retained by the active selectors.

For an opt-in exact strategy:

• After eligibility is explicitly defined, exactly min(N, eligible) occurrences are selected; N > eligible selects each eligible occurrence once.
• Same-QNAME distinct occurrences and equal-identity duplicate multiplicity are preserved with deterministic same-file tie-breaking.
• Coordinate-constrained modes exclude unmapped records; unconstrained Summary follows the chosen single-budget policy.
• Zero-contribution focused records and CIGAR deletion/reference-skip controls preserve the selected eligibility contract.

For legacy exhaustion:

• On the frozen six-occurrence fixture, N=6 and N=99 match each format's all-record control across clustered starts, workers, and intervals.
• Extract retains both same-QNAME rows because its all-record control does; thresholding matches its existing all-record QNAME aggregation unless separately changed.
• All-unmapped mapped-only N=0,3,99 behavior is consistent between BAM and CRAM for each consumer contract.

Existing fractional modes remain byte-identical. Planning time, extraction time, plan size, and peak RSS are measured on a representative large indexed modBAM.

Reproduction artifacts

Artifact	Size	SHA-256	Notes
make_fixture.py	Inline	N/A	Deterministic generator embedded above
fixture.sam	Generated	Deterministic from inline script	Eight mapped 80-base alignment occurrences on two contigs
Legacy BAM/CRAM fixture	Tiny synthetic	Maintained with the parent-red regression	Six primary mapped occurrences including two same-QNAME alignments

Related work

• Separate local repairs stabilize seeded fractional sampling, sparse-coordinate metadata ownership, coordinate-selector handling, CRAM reference plumbing, and unmapped modification classification. They are prerequisites but do not choose the fixed-count estimand.
• Proposed PR: none until the maintainer selects the statistical contract.
• Deferred policy question: approximate reference-span weighting versus an explicit exact occurrence-global strategy is the purpose of this issue.