`stats` and `localize` can succeed with incomplete results after regional query errors

[SuhasSrinivasan]

stats and localize can succeed with incomplete results after regional query errors

Summary

modkit stats and modkit localize discard errors returned while fetching or decoding individual indexed bedMethyl regions. Both commands can therefore exit successfully and write apparently valid output that omits one or more requested regions.

Severity

Severity: High — scientific correctness and output reliability

Rationale: A malformed or otherwise unreadable regional query can silently remove its observations from counts and percentages while the process returns success. The resulting TSV does not identify the omitted region, so downstream analyses can treat materially incomplete scientific output as complete.

User and scientific impact

• Affected result or workflow: regional modification summaries from stats and feature-centered aggregates from localize.
• Direction of error: omitted region rows, reduced aggregate coverage and modified counts, altered percentages, and false-success output.
• Likely exposure: data-dependent; any Tabix fetch, BGZF decode, UTF-8, or bedMethyl parse failure reached inside a requested region can trigger it.
• Detectability or workaround: debug logging may contain a transient message, but normal exit status and output do not reveal the missing scientific population. Users must independently validate every requested region.

Affected versions and environment

• Reproduced release: modkit 0.6.4.
• Development revision: 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c.
• Operating system and architecture: macOS 26.6 on Apple Silicon.
• Input format: bgzip-compressed, tabix-indexed bedMethyl plus BED regions; localize also uses genome sizes.
• Related tool versions: htslib bgzip and tabix 1.23.1 for the fixture below.

Steps to reproduce

Minimal input

printf 'chr1\t100\t101\tm\t4\t+\t100\t101\t255,0,0\t4\t25.00\t1\t3\t0\t0\t0\t0\t0\nchr1\t200\t201\tm\tnot-a-number\t+\t200\t201\t255,0,0\t6\t50.00\t3\t3\t0\t0\t0\t0\t0\nchr1\t300\t301\tm\t8\t-\t300\t301\t255,0,0\t8\t50.00\t4\t4\t0\t0\t0\t0\t0\nchr1\t400\t401\tm\t10\t-\t400\t401\t255,0,0\t10\t20.00\t2\t8\t0\t0\t0\t0\t0\n' > query-malformed.bed
bgzip -k -f query-malformed.bed
tabix -f -p bed query-malformed.bed.gz

printf 'chr1\t100\t101\nchr1\t200\t201\nchr1\t300\t301\nchr1\t400\t401\nabsent\t0\t1\n' > query-regions.bed
printf 'chr1\t1000\nabsent\t10\n' > genome-sizes.tsv

Commands

modkit stats query-malformed.bed.gz \
  --regions query-regions.bed \
  --threads 1 --io-threads 1 \
  --out-table stats.tsv --force

modkit localize query-malformed.bed.gz \
  --regions query-regions.bed \
  --genome-sizes genome-sizes.tsv \
  --window 1 --min-coverage 1 \
  --threads 1 --io-threads 1 \
  --out-file localize.tsv --force

Repeat with --threads 4 to exercise parallel reduction.

Control or independent oracle

Replace not-a-number with 6, rebuild the compressed file and index, and rerun both commands. The valid input has four requested chr1 observations with total valid coverage 28 and total modified count 10. stats must contain all four region rows; localize must aggregate all four observations. The syntactically valid absent contig intentionally remains an accepted skip control.

Observed behavior

At both one and four threads, modkit 0.6.4 exits zero.

stats silently omits chr1:200-201 and writes only three data rows:

#chrom	start	end	name	strand	count_m	count_valid_m	percent_m
chr1	100	101	.	.	1	4	25
chr1	300	301	.	.	4	8	50
chr1	400	401	.	.	2	10	20

The file is 134 bytes with SHA-256 81e73ccc8f6df0e90758ff0aa911aa27f0108b661aa3cd2990d0c3cf5f441e2e.

localize silently aggregates only the three decodable observations, reporting n_valid=22, n_mod=7, and 31.818182% instead of 28, 10, and 35.714287%. Its 67-byte output has SHA-256 0682a22be0a9cdd110036f0d4b30c16e634938e0faaffb6743783ee6608894b9. The affected release also labels that aggregate offset -1 because of the separate localize geometry defect; the missing counts are independent of that label.

With the valid control, stats emits all four rows byte-identically at one and four threads, SHA-256 6b42943c39a9ff8e58cda0c48ac9158fe666b590e17b78ab51b6de42eb6e0ca5.

Expected behavior

• A fetch/decode/parse failure for a requested region is fatal and produces a nonzero exit.
• The error retains command, region, and causal decode context.
• Failure occurs before creating a new output or truncating an existing output.
• Syntactically valid regions on contigs absent from the input index retain the existing documented skip behavior.
• Valid inputs preserve exact counts, rows, and deterministic thread behavior.

Root-cause evidence

• Relevant code: modkit-core/src/localise/subcommand.rs and modkit-core/src/stats/subcommand.rs at revision 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c.
• localize folds Result<LocalizedModCounts> values but logs and discards every Err before infallible reduction.
• stats uses filter_map to convert each failed Result<MethylationStats> into None before infallible fold/reduce.
• Both output destinations are opened before regional processing is known to have succeeded.
• Parent-red/fix-green evidence: the parent produces the incomplete successful outputs above. The reviewed correction returns failed to localize region chr1:199-202: invalid-bedmethyl-data or failed to calculate stats for region chr1:200-201: invalid-bedmethyl-data, leaves an absent output absent, and preserves a 29-byte existing sentinel with SHA-256 23cdc9ad2b5b5f71b423f3d602e7b2e0aa4de668a359a71b8a2d3c077c866eb6 at both one and four threads.

Proposed fix scope

Keep each regional result fallible through Rayon try_fold/try_reduce, attach command and region context, and return the failure to the CLI. Build the complete table in memory and open the destination only after all requested regional work succeeds.

Non-goals

• Do not change the established skip policy for contigs absent from the Tabix index.
• Do not combine this fix with localize BED parsing or window/offset geometry changes; those have a separate input/coordinate contract.
• Do not claim deterministic selection of one exact error when several regions fail concurrently; every returned error must retain a real causal region and source.
• Do not claim transactional rollback after output writing begins, compressed-writer finalization guarantees, or a broader cross-command failure-policy rewrite.
• Do not redesign Tabix reader reuse or parallel scheduling for performance.

Acceptance criteria

• A malformed requested bedMethyl record makes both commands exit nonzero with region and causal context at one and multiple threads.
• Failed runs do not create absent outputs and do not change forced existing outputs.
• Valid controls retain all four requested chr1 observations and have exact one-thread/multi-thread output equality.
• Missing-index-contig controls retain the existing skip behavior.
• Injected Tabix fetch and BGZF/decode failures follow the same fatal path.
• Focused stats, localize, and bedMethyl parsing tests plus the applicable full workspace suite pass.

Reproduction artifacts

Artifact	Size	SHA-256	Notes
query-malformed.bed	241 bytes	3e12a38e60fe124d8a710cf77991b192acec0d6eb393524b12c842b73ece5add	Generated inline above
query-malformed.bed.gz	162 bytes	ba5190520812e42c49c6f97ba73db389e17017fbc41c646319d0b7c54bd22ea7	bgzip 1.23.1
query-malformed.bed.gz.tbi	108 bytes	38d2d078e5ec5f1919df0f752aa742ded9a5fc8344bc9d1ac7b28aced785f316	tabix 1.23.1
query-regions.bed	63 bytes	63945f31407c114d68326307099d50e29cc270b0bdc1fea7eabca92da7b0be12	Four indexed regions plus one missing-contig control
genome-sizes.tsv	20 bytes	908fa762ee27315ee92ad28a4a59590731aa3d9834e251a4deb92c71f08e7ec8	Used by localize
query-valid.bed.gz	150 bytes	2ee01f500b65bb63e883a66ad79a141b3fac9dbb7ae77f15b61f44e69dd007c3	Control with score field corrected to 6

Related work

• Related issues: the localize input/coordinate contract is reported separately; writer finalization and broader partial-output policy are independent lifecycle work.
• Proposed PR: to be linked after the localize source overlap is ordered and the exact rebased package is regated.
• Deferred policy question or upstream dependency: none for fatal requested-region decode errors; dependency-resolution PR Pin hts-sys to bindings compatible with rust-htslib 0.46 #666 affects publication timing only.