`entropy --regions` omits singleton summaries and computes medians from genomic rather than value order

[SuhasSrinivasan]

entropy --regions omits singleton summaries and computes medians from genomic rather than value order

Summary

Regional entropy summaries have two independent median defects in the same local statistics function. A region with exactly one successful entropy window is omitted because the shared percentile helper rejects singleton input, while a region with multiple windows can report the wrong median because genomic-order values are passed to a helper that requires sorted values.

The command can exit zero with an empty region file or a plausible but incorrect median.

Severity

Severity: High — scientific/statistical accuracy

Rationale: The unsorted-input case silently changes a reported regional statistic and is not detectable from the summary row alone. The singleton case is narrower and would be Medium by itself, but it silently omits a valid requested summary. Both fixes are small, entropy-local branches in the same function and should be reviewed together.

User and scientific impact

• Affected result or workflow: modkit entropy --regions regional BED summaries.
• Direction of error: omitted singleton region rows and incorrect odd/even medians.
• Likely exposure: every region with exactly one successful window; any multi-window region whose genomic entropy values are not already sorted by value.
• Detectability or workaround: retain and independently sort the per-window bedGraph values, then recompute the regional median. The summary output does not indicate that sorting was skipped.

Affected versions and environment

• Released version: modkit 0.6.4.
• Development revision: 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c.
• Verification environment: macOS 26.6, arm64, Apple M2 Ultra.
• Input formats: coordinate-sorted BAM with BAI, FASTA with FAI, and BED regions.
• Related tool version: samtools 1.23.1.

Steps to reproduce

Singleton region

Save the following as anchor.fa, anchor.sam, and singleton.regions.bed.

>chr1
AACGAA

@HD	VN:1.6	SO:coordinate
@SQ	SN:chr1	LN:6
read-0	0	chr1	1	60	6M	*	0	0	AACGAA	??????	MM:Z:C+m?,0;	ML:B:C,255	MN:i:6	NM:i:0

chr1	2	3	singleton

Use --base C so this reproducer does not depend on flanking motif lookup. Use an explicit threshold to isolate regional statistics from automatic-threshold sampling.

samtools faidx anchor.fa
samtools view -b -o anchor.bam anchor.sam
samtools index anchor.bam

modkit entropy \
  --in-bam anchor.bam \
  --out-bed singleton-output \
  --ref anchor.fa \
  --base C \
  --num-positions 1 \
  --window-size 1 \
  --min-coverage 1 \
  --max-filtered-positions 0 \
  --filter-threshold 0 \
  --threads 1 \
  --io-threads 1 \
  --regions singleton.regions.bed \
  --prefix singleton \
  --suppress-progress

wc -c singleton-output/*

Multi-window statistical oracle

The affected function can be exercised with ordinary regional entropy values in genomic order. Two minimal examples are:

odd genomic order:  [0.75, 0.25, 0.50]
even genomic order: [0.75, 0.00, 0.25, 0.50]

The correct medians after value sorting are:

median([0.25, 0.50, 0.75])       = 0.50
median([0.00, 0.25, 0.50, 0.75]) = 0.375

Control or independent oracle

A singleton sample has mean, median, minimum, and maximum all equal to its only value. For the fixture above the entropy is 0, the accepted-read count is 1, and the regional statistics must be:

mean_entropy=0
median_entropy=0
min_entropy=0
max_entropy=0
mean_num_reads=1
min_num_reads=1
max_num_reads=1
successful_window_count=1
failed_window_count=0

For multiple values, a median is an order statistic. Genomic coordinate order has no statistical meaning for selecting the middle value. With the existing linear-interpolation convention, the two exact medians above are 0.50 and 0.375.

Observed behavior

The singleton command exits zero and writes one successful window, but the region file is empty:

finished, 1 regions processed successfully, 1 windows failed
not enough datapoints, got 1

singleton_windows.bedgraph: 15 bytes
singleton_regions.bed:        0 bytes

The empty region file has SHA-256 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855.

For the multi-window examples, the affected interpolation path uses the values at their current positions rather than their sorted positions:

• [0.75, 0.25, 0.50] reports 0.25 instead of 0.50.
• [0.75, 0.00, 0.25, 0.50] interpolates 0.00 and 0.25 and reports 0.125 instead of 0.375.

Mean, minimum, maximum, read-count statistics, and success/failure counts are not implicated by the ordering defect.

Expected behavior

• A region with one successful entropy window emits one regional summary row.
• Singleton mean/median/min/max equal the only entropy measurement; read-count statistics and success/failure counts remain exact.
• A multi-window regional median is invariant under permutation of the same entropy measurements.
• Odd and even medians use the existing linear-interpolation convention after sorting by numeric value.
• Computing a regional median does not mutate the caller's genomic-order measurements.

Root-cause evidence

• DescriptiveStats::new sends every nonempty measurement slice directly to the shared percentile helper at entropy/mod.rs:1365-1391.
• The shared helper explicitly rejects fewer than two data points and indexes the supplied slice without sorting it at thresholds.rs:16-35.
• Other callers sort threshold probabilities before invoking the helper, but entropy regional values arrive in genomic window order. The parent-regression examples therefore deterministically produce the wrong 0.25 and 0.125 medians above.

Proposed fix scope

• In entropy's DescriptiveStats::new, return the single measurement directly as the singleton median.
• For two or more measurements, copy the slice, sort the copy with f32::total_cmp, and pass that sorted copy to the unchanged interpolation helper.
• Continue computing mean/min/max and read-count statistics from the original input.

Non-goals

• Do not change the shared percentile helper's established two-observation contract for automatic-threshold callers.
• Do not change entropy-window coordinates, region ownership, filtering, thresholds, entropy normalization, or final floating-point reduction.
• Do not redefine the regional aggregation unit or weight windows by coverage.

Acceptance criteria

• A singleton measurement such as entropy 0.25 with seven reads yields mean/median/min/max 0.25, all read-count statistics 7, success count 1, and the supplied failure count unchanged.
• The CLI singleton fixture emits a nonempty region row with the exact statistics above.
• [0.75, 0.25, 0.50] reports median 0.50.
• [0.75, 0.00, 0.25, 0.50] reports median 0.375.
• Permutations of the same odd/even measurement multisets return the same median.
• The input measurement order is unchanged after summary construction.
• Mean/min/max, read-count statistics, successful/failed counts, shared threshold percentiles, and existing multi-window outputs remain unchanged except for previously incorrect medians.
• Focused entropy tests and the full workspace suite pass.

Reproduction artifacts

Artifact	Size	SHA-256	Notes
anchor.sam	122 bytes	a6024e8e4849097b0864a525e80b3df8df84a4c82ad232f4f7a7cf3772742f93	Inline one-read fixture
anchor.bam	262 bytes	c7036578a3538d3683f919b33fd9a8040940052f507a10f23278817f009f443c	samtools 1.23.1 conversion
anchor.bam.bai	96 bytes	5fac9d50548cbfa173e5cfada3df0907b562da155f28acab0e8bc00c938d7eea	BAM index
anchor.fa	13 bytes	f765963cba9bf7f33907ba35f2a569ec9cec427f00c4751d125c2b162e7d2f74	Inline reference
anchor.fa.fai	13 bytes	f8fa22465200794bd9facbf11c8daf763ff02341dc1d2b3d9fc8ac358dba29ed	FASTA index
singleton.regions.bed	19 bytes	04d53c9bec5bdcbe34301c49832eae377370a382d9f6e46518c9822f13937b0b	One requested region

Related work

• Proposed PR: one small regional-statistics PR containing separate singleton and sorted-median commits.
• Coordinate/window ownership should remain a separate issue and PR even though its corrected interval can be composed in the final integration test.
• Final-bit entropy reduction reproducibility is a separate numerical-ordering issue; it is not required for these exact statistical oracles.

[SuhasSrinivasan]

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