Automatic thresholds can use an arbitrary explicit-canonical probability cap

[SuhasSrinivasan]

Summary

Automatic threshold calculation does not reliably reduce explicit-canonical probabilities across every sampled read. Within each Rayon fold, the per-base value is overwritten for every HashMap read, so the fold contributes the last arbitrary iteration value rather than its maximum. The final threshold cap can therefore vary across processes or thread layouts and be lower than the true global explicit-canonical maximum.

Severity

Severity: High — scientific correctness and reproducibility

Rationale: When the requested percentile lies above the explicit-canonical cap, an arbitrary lower cap makes filtering more permissive. The same sampled calls can receive different thresholds and pass/fail classifications solely because HashMap iteration and Rayon folding changed.

User and scientific impact

• Affected workflows are automatic-threshold consumers using the legacy ReadIdsToBaseModProbs aggregation path.
• A spuriously low cap can retain lower-confidence modification calls that should have failed the intended threshold.
• The defect is schedule/process dependent and can be hidden when the requested percentile is already below the arbitrary cap.
• It is independent of which reads were sampled; deterministic sampling alone does not repair it.

Affected versions and environment

• Released version: modkit 0.6.4.
• Development revision: 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c.
• Reproduced on macOS arm64 with the repository's ten-read indexed BAM.

Steps to reproduce

First obtain an independent maximum from unfiltered extracted calls:

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

modkit extract calls "$bam" "$tmp/calls.tsv" --no-filtering \
  --threads 1 --io-threads 1 --suppress-progress --force

awk -F '\t' '
  NR == 1 { for (i = 1; i <= NF; i++) c[$i] = i; next }
  $c["canonical_base"] == "C" && $c["call_code"] == "-" &&
    $c["inferred"] == "false" {
      n++
      if ($c["call_prob"] + 0 > max) max = $c["call_prob"] + 0
    }
  END { print "rows=" n, "global_max=" max }
' "$tmp/calls.tsv"

Then force the percentile above that cap and repeat fresh processes:

for run in 1 2 3; do
  modkit call-mods "$bam" "$tmp/run-$run.sam" --output-sam \
    --sampling-frac 1 --seed 7 --threads 1 \
    --sampling-interval-size 20 --filter-percentile 0.99 \
    --suppress-progress --log-filepath "$tmp/run-$run.log"
  grep 'estimated pass threshold' "$tmp/run-$run.log"
done

Control or independent oracle

For each primary base, calculate max(call_prob) across every explicit, non-inferred canonical call in the unchanged sampled population. The final cap is mathematically min(percentile_threshold, global_explicit_canonical_max). Maximum is associative, commutative, and idempotent, so partitioning and iteration order cannot change it.

Observed behavior

The fixture contains ten explicit C-canonical rows and has a true global maximum of 0.8671875. Fresh installed-0.6.4 processes returned caps including:

0.58984375
0.77734375
0.78515625

The uncorrected percentile before capping was 0.9977344. A corrected build consistently caps it at 0.8671875.

An independent 99-read unit fixture with a true maximum of 0.99 produced a parent value of 0.71000004 before the same one-line reduction repair.

Expected behavior

Every fold must retain the maximum explicit-canonical probability seen so far for each base, and the cross-fold reduction must take the maximum again. The result must be the global maximum for the sampled population regardless of read order, process hash seed, or Rayon thread count.

Root-cause evidence

The fold currently performs an unconditional canonical_probs.insert(base, p). self.inner is a HashMap, so the retained value is the last arbitrary read visited by that fold. The later cross-fold maximum cannot recover larger values already overwritten inside a fold.

Proposed fix scope

Replace the fold-local overwrite with entry-wise maximum reduction:

canonical_probs
    .entry(base)
    .and_modify(|current| {
        if p > *current {
            *current = p;
        }
    })
    .or_insert(p);

Retain the existing cross-fold maximum, sampled calls, percentile algorithm, fallback threshold, and per-base separation.

Non-goals

• Do not change the sampling population or sampling algorithm.
• Do not change percentile interpolation or the scientific meaning of the explicit-canonical cap.
• Do not merge different canonical bases or modification codes.
• Do not change raw-QNAME aggregation or claim downstream repeatability when that separate aggregation-key/order contract produces a different post-aggregation population.

Acceptance criteria

• A multi-read fixture with a known per-base maximum returns that maximum across fresh processes and Rayon thread counts.
• The repository fixture caps 0.9977344 at exactly 0.8671875 in every run.
• A multi-base control proves each base is reduced independently.
• Threshold and downstream output are repeatable when the aggregated ReadIdsToBaseModProbs population is unchanged; the executable fixture uses unique QNAMEs to isolate this reduction invariant.
• Existing percentile/fallback tests remain byte-identical where the cap is inactive.

Reproduction artifacts

Artifact	SHA-256	Notes
tests/resources/bc_anchored_10_reads.sorted.bam	4441acbe1ad59caf6ab5d56bf18097d8999c08351395219c919cb863652d85f0	Public CLI reproduction; global explicit C-canonical maximum 0.8671875
tests/resources/bc_anchored_10_reads.sorted.bam.bai	e2a74713867df4a01785d614e066fd310fea08293fc0b86488ec336d89fe0c35	BAM index

Related work

• Deterministic fractional sampling is a prerequisite for stable membership but is a separate root cause.
• Fixed-count sampling policy in issue Decide indexed --num-reads semantics: reference-span-balanced approximate sampling vs exact record-global sampling #654 does not change this reduction invariant.
• Edge-filter population consistency in issue Probability sampling applies --edge-filter inconsistently across execution paths #667 is separate from reducing one already selected population.