Probability sampling applies --edge-filter inconsistently across execution paths

[SuhasSrinivasan]

Summary

Probability sampling applies --edge-filter inconsistently across execution paths. Without a reference-position selector, both indexed argmax handlers ignore the filter entirely. With --include-bed or another aligned selector, the indexed handlers apply asymmetric start/end trims to raw alignment-orientation query positions, so reverse reads are filtered from the wrong molecular ends. The aligned predicate also subtracts an unchecked end trim from the read length. Separately, the serial histogram path swaps asymmetric ends but always keeps the interior, ignoring EdgeFilter.inverted; extract calls --invert-edge-filter can therefore estimate an automatic threshold from a population disjoint from its emitted end calls.

These behaviors reproduce on accepted indexed modBAM input with modkit 0.6.4 and current upstream revision 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c.

Severity

Severity: High — scientific correctness, threshold populations, and reliability

Rationale: Indexed and serial execution can select different modification-probability populations from the same records and options, and inverted serial thresholding can sample read interiors while output retains read ends. These differences change histograms and estimated percentiles, and can therefore change automatic thresholds and downstream calls solely because an index, reference-position selector, or threshold-collection path is present. Overlong trims also have an unsafe arithmetic path in indexed aligned sampling.

User and scientific impact

• Affected result or workflow: sample-probs, Summary probability statistics, and commands that use indexed or serial probability extraction for automatic threshold estimation, including extract calls --invert-edge-filter.
• Direction of error: unfiltered end calls can be retained, the wrong end of reverse reads can be retained under asymmetric filtering, or thresholds for emitted end calls can be estimated from interior calls.
• Likely exposure: routine whenever --edge-filter is combined with an indexed BAM; reverse-orientation differences require an asymmetric filter to be visible.
• Detectability or workaround: compare against --ignore-index; users can force serial sampling, but that gives up indexed parallelism and is easy to omit in downstream threshold consumers.

Affected versions and environment

• Released version: modkit 0.6.4
• Development revision: 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c
• Operating system and architecture: macOS 26.6, arm64
• Input format and index: coordinate-sorted BAM with BAI
• Related tools: samtools 1.23.1 for the synthetic fixture

Steps to reproduce

These fixtures are synthetic and contain no private data.

Create edge-filter.sam:

@HD	VN:1.6	SO:coordinate
@SQ	SN:chr1	LN:100
forward	0	chr1	1	60	5M	*	0	0	AAAAA	IIIII	MM:Z:A+a?,0,0,0,0,0;	ML:B:C,201,202,203,204,205	MN:i:5
reverse	16	chr1	1	60	5M	*	0	0	TTTTT	IIIII	MM:Z:A+a?,0,0,0,0,0;	ML:B:C,201,202,203,204,205	MN:i:5

Create the indexed BAM and a BED covering every aligned position:

samtools view -b -o edge-filter.bam edge-filter.sam
samtools index edge-filter.bam
printf 'chr1\t0\t5\n' > all-positions.bed

Compare indexed and serial histogram output without a position selector:

mkdir indexed serial

modkit sample-probs edge-filter.bam \
  --hist --out-dir indexed --prefix result --force \
  --no-sampling --edge-filter 1,2 \
  --threads 1 --io-threads 1 --suppress-progress

modkit sample-probs edge-filter.bam \
  --ignore-index \
  --hist --out-dir serial --prefix result --force \
  --no-sampling --edge-filter 1,2 \
  --io-threads 1 --suppress-progress

cat indexed/result_probabilities.tsv
cat serial/result_probabilities.tsv

Then compare the aligned handlers selected by an all-position BED:

mkdir indexed-aligned serial-aligned

modkit sample-probs edge-filter.bam \
  --include-bed all-positions.bed \
  --hist --out-dir indexed-aligned --prefix result --force \
  --no-sampling --edge-filter 1,2 \
  --threads 1 --io-threads 1 --suppress-progress

modkit sample-probs edge-filter.bam \
  --ignore-index --include-bed all-positions.bed \
  --hist --out-dir serial-aligned --prefix result --force \
  --no-sampling --edge-filter 1,2 \
  --io-threads 1 --suppress-progress

cat indexed-aligned/result_probabilities.tsv
cat serial-aligned/result_probabilities.tsv

Finally, reproduce the inverted serial automatic-threshold mismatch. Create invert-edge.sam:

@HD	VN:1.6	SO:coordinate
@SQ	SN:chr1	LN:9
invert-edge	0	chr1	1	60	9M	*	0	0	CCCCCCCCC	IIIIIIIII	MM:Z:C+m?,0,0,0,0,0,0,0,0,0;	ML:B:C,140,150,220,220,220,220,220,230,240	MN:i:9

Create an unindexed BAM and run automatic extraction:

samtools view -b -o invert-edge.bam invert-edge.sam

modkit extract calls invert-edge.bam inverted.tsv \
  --edge-filter 2 --invert-edge-filter \
  --filter-percentile 0.5 \
  --threads 1 --io-threads 1 \
  --suppress-progress --force

Control or independent oracle

--edge-filter 1,2 is defined on the original molecular read: discard the first quality and the last two qualities. Each five-base record must therefore contribute ML qualities 202 and 203, irrespective of alignment orientation. Across the two records, the histogram must contain quality 202 twice and quality 203 twice.

The serial --ignore-index path provides an independent implementation of that documented molecular-orientation rule.

For the nine-base inverted fixture, output correctly retains query positions 0, 1, 7, and 8 with ML qualities 140, 150, 230, and 240. The automatic threshold must be estimated from those same four calls: percentile rank 2 gives quality 150, probability 0.5878906. Position 1 is equal to that threshold and must pass.

Observed behavior

Without --include-bed, indexed mode ignores the filter and retains all five qualities from both records:

code	primary_base	range_start	range_end	count	frac	percentile_rank
a	A	0.78515625	0.7890625	2	0.2	10
a	A	0.7890625	0.79296875	2	0.2	30
a	A	0.79296875	0.796875	2	0.2	50
a	A	0.796875	0.80078125	2	0.2	70
a	A	0.80078125	0.8046875	2	0.2	90

Serial mode retains only the expected four observations:

code	primary_base	range_start	range_end	count	frac	percentile_rank
a	A	0.7890625	0.79296875	2	0.5	25
a	A	0.79296875	0.796875	2	0.5	75

With the all-position BED, indexed mode applies the trim but does not swap the molecular ends for the reverse alignment. It retains quality 202 once, quality 203 twice, and quality 204 once. Serial mode again retains qualities 202 and 203 twice each.

The no-selector indexed and serial probability TSVs had SHA-256 values bee3adc7decd91060d7f72f24704a90100de583334374dfdd0e7572b39b79474 and 09c68f6af2c57329f5c87459f6d2e25af102e90cf51437f68222ffac638b21f5, respectively. The aligned indexed TSV had SHA-256 04667f73abac11e7b6b2e654b2a857ab9ffc5340a610a4ccfbcc0035b6f7b0da; its serial counterpart again had 09c68f6af2c57329f5c87459f6d2e25af102e90cf51437f68222ffac638b21f5.

For the unindexed inverted fixture, modkit reports that it estimated from five probabilities and sets the threshold to quality 220, probability 0.8613281. Those are the five interior qualities at query positions 2–6, not the four end calls written to inverted.tsv. Position 1 is consequently marked fail=true; using the correct end-call population makes it fail=false. With --pass-only, the incorrect threshold writes two calls instead of the expected three.

Expected behavior

• Indexed and --ignore-index execution must retain the same probability multiset for the same records and edge filter.
• Start and end trims must refer to original molecular-read ends on both forward and reverse alignments.
• The five-base 1,2 fixture must retain ML qualities 202 and 203 twice each, with or without the all-position BED.
• Normal and inverted filtering must collect the exact same population used by output, including unindexed input and indexed unmapped fallback.
• The nine-base inverted fixture must set the automatic C threshold to quality 150 (0.5878906) and mark query position 1 as passing.
• An overlong or overlapping filter must retain no calls without panic, integer wrap, or path-dependent output.
• Omitting --edge-filter must preserve existing output.

Root-cause evidence

• Both unaligned indexed argmax handlers accept the edge filter as _edge_filter and discard it at sample_probs/mod.rs:1056-1185.
• The aligned indexed iterator compares raw query positions to edge_filter_start and unchecked read_length - edge_filter_end, without converting reverse-alignment query positions back to molecular orientation, at sample_probs/mod.rs:839-877 and sample_probs/mod.rs:1190-1337.
• The serial histogram implementation extracts only start/end sizes, swaps them by alignment orientation, and always applies the interior conjunction; it never consults EdgeFilter.inverted at sample_probs/mod.rs:394-545.
• EdgeFilter::keep_position provides checked/saturating filter semantics at mod_bam.rs:1668-1705.

Proposed fix scope

• Retain the optional EdgeFilter in the indexed probability extractor instead of reducing or discarding it.
• Use one shared molecular-orientation predicate for all four indexed handler variants and both serial histogram branches.
• Reject out-of-range query positions, map reverse-alignment positions back to original molecular-read coordinates, and delegate the actual trim/inversion decision to EdgeFilter::keep_position.
• Apply the same predicate to indexed unmapped fallback, which reuses serial record collection.
• Preserve record ownership, sampling, selectors, counts, and output schemas.

Non-goals

• No change to fractional or fixed-count sampling, RNG scheduling, or interval ownership; indexed fixed-count policy remains the separate discussion in Decide indexed --num-reads semantics: reference-span-balanced approximate sampling vs exact record-global sampling #654.
• No change to unmapped-record ownership, budget, or reference-position selector semantics; only the calls retained from each already-selected record change.
• No addition of a sample-probs --invert-edge-filter option.
• No change to output-side edge filtering, which already honors inversion.
• No change to --ignore-index threshold-reader selection, explicit thresholds, --no-filtering, pass-only row-count reporting, or log wording.
• No threshold-policy or performance redesign.

Acceptance criteria

• All four indexed extractor variants apply one checked molecular-orientation edge predicate.
• Both serial histogram branches use the same checked predicate and honor normal/inverted filtering.
• Indexed and serial histogram directories are byte-identical for forward/reverse asymmetric-filter fixtures, both with and without a position selector, and mapped plus unmapped fallback retains the correct population once.
• The nine-base inverted extraction fixture produces the quality-150 threshold and expected q0/q1/q7/q8 pass/fail states.
• Overlong and overlapping trims retain zero probabilities without panic or wrapping.
• No-filter indexed output remains byte-identical to the parent behavior.
• Record sampling, ownership, and reported record counts are unchanged.
• Focused, affected-package, and full workspace tests pass.

Reproduction artifacts

Artifact	Size	SHA-256	Notes
edge-filter.sam	237 B	e54f3a14ac617ae450f026bac02231322ccb72cc503be420b29823869943a7aa	Forward/reverse five-base synthetic reads
edge-filter.bam	292 B	ba1a3acd3100480d0bb450edc382610afd69b1fae741e7572929be833f21b5b6	Coordinate-sorted modBAM
edge-filter.bam.bai	96 B	734b1f6c6d1224eb0e1cc494ce1f7399c018cb71cbee9fe3d52535a1d4983358	BAM index
all-positions.bed	9 B	7c0b71d52436cfeb4bad2ee5baa96636a33ea1e135781a93bb3356f794080a72	Covers every aligned fixture position
invert-edge.sam	174 B	86ce8fa6f2fb8cb1cf309b54c13a37b6a346bb06af8ec41316737ccdc3948ec0	Distinct interior/end threshold populations
invert-edge.bam	276 B	d469c4a43ea08c408bf67f25f9399193c8bc682f987bdf379d2ca643b2b8db8f	Unindexed inverted-filter fixture

Related work

• Decide indexed --num-reads semantics: reference-span-balanced approximate sampling vs exact record-global sampling #654 concerns the independent policy for indexed fixed-count sampling and is not changed here.
• Fresh dependency resolution selects incompatible hts-sys bindings #665 / Pin hts-sys to bindings compatible with rust-htslib 0.46 #666 address a fresh-build dependency incompatibility; the proposed code fix remains independent of that manifest-only change.
• The broader extract-threshold population repair for --mapped-only and --ignore <mod> is separate; it does not correct edge-filter inversion.

[SuhasSrinivasan]

The complete code change is in ready-for-review PR #696. Closing this standalone catalog issue in line with the maintainer's request that these submissions be represented by PRs rather than separate AI-generated issues.