Automatic thresholding can silently accept samples with no probability observations

[SuhasSrinivasan]

Summary

Automatic thresholding silently accepts a sample that contains no usable modification-probability observations. The legacy path converts an empty probability map into a caller whose missing per-base thresholds default to zero. Modern QualHist consumers can likewise report successful empty results without distinguishing parsed records from actual probability observations.

In the legacy pileup-hemi reproduction, this makes an empty automatic sample scientifically indistinguishable from an explicit --no-filtering request. Other modern consumers can instead return a successful empty report without distinguishing an undefined estimate from an intentional explicit/no-filter mode. Fraction zero, empty inputs, and selector/filter combinations that sample no eligible calls can therefore disable filtering or falsely report successful automatic analysis.

Severity

Severity: High — scientific correctness and reporting

Rationale: An automatic threshold request can silently become threshold zero and emit results equivalent to no filtering. Users receive successful output without any estimated population, while explicit thresholds and true no-filter modes have materially different intent.

User and scientific impact

• Affected legacy automatic-threshold consumers include duplex/entropy-style paths using the shared probability map.
• Affected modern consumers include Summary, standard pileup, and extract calls paths using QualHist.
• Empty samples can arise from fraction 0, empty modification tags, mapped/selector restrictions, or a valid sample containing records but no eligible probabilities.
• The failure is easy to miss because the command may exit zero and produce plausible-looking no-filter output.

Affected versions and environment

• Released version: modkit 0.6.4.
• Development revision: 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c.
• Reproduced on macOS arm64 with checked-in duplex and empty-tag BAM fixtures.

Steps to reproduce

Legacy empty sample:

tmp=$(mktemp -d /tmp/modkit-empty-threshold.XXXXXX)
input=../tests/resources/duplex_modcalls_sort.bam
ref=../tests/resources/GRCh38_chr20.fa

modkit pileup-hemi "$input" --out-bed "$tmp/fraction-zero.bed" \
  --ref "$ref" --motif CG 0 \
  --region chr20:22,613,835-22,640,468 \
  --sampling-frac 0 --suppress-progress

modkit pileup-hemi "$input" --out-bed "$tmp/no-filter.bed" \
  --ref "$ref" --motif CG 0 \
  --region chr20:22,613,835-22,640,468 \
  --no-filtering --suppress-progress

wc -l "$tmp/fraction-zero.bed" "$tmp/no-filter.bed"
shasum -a 256 "$tmp/fraction-zero.bed" "$tmp/no-filter.bed"

Modern empty-observation sample:

modkit summary ../tests/resources/empty-tags.sorted.bam \
  --num-reads 10 --threads 1 --io-threads 1 --tsv \
  --suppress-progress

Control or independent oracle

Automatic thresholding requires at least one canonical or separated-modification probability observation. Record count alone is not evidence of an estimable distribution. Explicit --filter-threshold 0 and --no-filtering are the controls that may intentionally accept an empty population.

Observed behavior

Installed modkit 0.6.4 exits zero for both legacy commands. Each writes 341 rows, and both outputs have the same SHA-256:

336479a1351d5594dd46295d65e6e5894620acdeb7208ca9d5973b2a16fe2b0a

The automatic Summary command also exits zero and reports:

mod_bases
total_reads_used	0

No automatic threshold was estimable in either case.

Expected behavior

Automatic mode must return nonzero and emit no scientific rows when the sampled population contains no usable modification probabilities. Paths that already resolve thresholds before writer creation must also preserve absent outputs and existing sentinels. Explicit thresholds and true no-filter modes must retain their documented empty-output behavior, and a sample containing only real modified-state probability histograms remains valid.

Root-cause evidence

The legacy threshold helper accepts an empty HashMap and constructs a caller whose absent per-base entries fall back to zero. The modern path does not consistently check whether canonical or separated-modification histograms contain observations before building thresholds or Summary results; an ok_records count can be nonzero even when the probability population is empty.

Proposed fix scope

• Add one shared legacy guard before percentile/caller construction that rejects an empty probability map.
• Add an explicit modern has_probability_observations predicate over canonical and separated-modification histograms.
• Apply the guard only to automatic modes.
• Preserve pre-output failure behavior on commands that already resolve automatic thresholds before writer creation; do not absorb unrelated writer-lifecycle changes.
• Ensure entropy's documented no-filter mode truly bypasses sampling before the legacy guard is enabled there.

Non-goals

• Fraction 0 remains a valid numeric fraction; it is the combination with automatic threshold estimation that cannot produce a threshold.
• Do not reject explicit zero thresholds, --no-filtering, or valid all-modified samples.
• Do not redefine malformed-record skip policy or selector eligibility.
• Do not make a general transactional-output claim for later runtime failures.

Acceptance criteria

• Fraction-zero and genuinely empty/filter-empty automatic samples fail with cannot calculate automatic thresholds because no modification probabilities were sampled.
• The automatic failure emits no scientific rows. On paths whose setup is already bounded before writer creation, it neither creates a new output nor mutates an existing sentinel.
• Summary, pileup, extract calls, duplex, and entropy automatic paths are covered across legacy and modern helpers.
• Explicit threshold zero and no-filter controls succeed with their current empty outputs.
• A sample containing only valid modified-state observations still computes a threshold.
• BAM/CRAM and worker-count controls preserve the same outcome.

Reproduction artifacts

Artifact	SHA-256	Notes
tests/resources/duplex_modcalls_sort.bam	21c2286ff7b2972a9c94ce24003bd3e20ba7c4fb38f416898a48eb5c2c789f94	Legacy fraction-zero reproduction
tests/resources/GRCh38_chr20.fa	5806b504166c43762e11c809f720b28b6d50462104fea22e7490deffbe6b59cf	Duplex reference
tests/resources/empty-tags.sorted.bam	29ab421d7f924e681737904eaadd10bd9a392f4587a4dfb96ae175c39012d845	Modern zero-observation reproduction
tests/resources/empty-tags.sorted.bam.bai	3696c4b4f0fe549db35b87e81e974967a727edfbf8f89595ed03faa29d9c1dd9	BAM index

Related work

• Entropy's no-filter bypass and pre-output setup ordering are prerequisites for applying the legacy guard safely to that command.
• Deterministic fractional sampling ensures a stable nonempty population but does not define the empty-population contract.
• Issue Probability sampling applies --edge-filter inconsistently across execution paths #667 concerns which edge-filtered calls enter a nonempty threshold population, not how an empty population is handled.
• Standard pileup currently opens its writer before automatic threshold resolution; moving that writer is a separate setup/output-lifecycle repair and is not claimed here.