`localize` can shift feature-centered coordinates and silently ignore invalid BED rows

[SuhasSrinivasan]

localize can shift feature-centered coordinates and silently ignore invalid BED rows

Summary

modkit localize constructs a feature window one base too far toward lower reference coordinates, discards the original feature anchor, and reports offsets with the opposite sign from the documented reference-axis convention. Its region loader can also discard malformed BED rows and continue successfully when at least one row parses, producing an aggregate over an unintended subset of features.

Severity

Severity: High — scientific correctness and input integrity

Rationale: The window defect changes which genomic observations enter every interior feature profile, while the reconstructed anchor and reversed sign place retained observations in incorrect offset bins. Separately, an accepted mixed-validity region file can produce a successful but incomplete aggregate without identifying the omitted features. These behaviors directly change biological metaprofiles and cannot reliably be corrected from the output alone.

User and scientific impact

• Affected result or workflow: feature-centered modification profiles generated by localize, including base-resolution and strand-filtered analyses.
• Direction of error: one extra lower-coordinate position, one omitted higher-coordinate position, shifted or mirrored offset bins, edge-dependent anchor drift, and silent omission of invalid feature rows.
• Likely exposure: the geometry and sign defects affect ordinary localize use; partial region acceptance is data-dependent.
• Detectability or workaround: users must compare output with the original BED and reference coordinates. Older profiles should be regenerated; simple sign negation is insufficient because the older query window and reconstructed anchor were also shifted.

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; plaintext BED regions; tab-separated genome sizes.
• Related tool versions: htslib bgzip and tabix 1.23.1 for the fixture below.

Steps to reproduce

Minimal input

printf 'chr1\t9\t10\tm\t4\t+\t9\t10\t255,0,0\t4\t25.00\t1\t3\t0\t0\t0\t0\t0\nchr1\t10\t11\tm\t4\t+\t10\t11\t255,0,0\t4\t50.00\t2\t2\t0\t0\t0\t0\t0\nchr1\t11\t12\tm\t4\t+\t11\t12\t255,0,0\t4\t75.00\t3\t1\t0\t0\t0\t0\t0\n' > geometry-counts.bed
bgzip -k -f geometry-counts.bed
tabix -f -p bed geometry-counts.bed.gz

printf 'chr1\t10\t11\n' > geometry-region.bed
printf '# valid row followed by a malformed row and another valid row\nchr1\t9\t10\nnot-a-valid-bed-row\nchr1\t11\t12\n' > mixed-regions.bed
printf 'chr1\t1000\n' > genome-sizes.tsv

Commands

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

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

Control or independent oracle

The midpoint of BED feature chr1:10-11 is reference position 10. An inclusive one-base window is the half-open query [9,12), and the reference-axis offset is bedMethyl start - 10. The exact expected table is therefore:

mod_code	offset	n_valid	n_mod	percent_modified
m	-1	4	1	25
m	0	4	2	50
m	1	4	3	75

Every non-comment region row is part of the requested scientific population. If any row is malformed, the command should reject the file rather than silently change that population.

Observed behavior

Both commands exit zero with modkit 0.6.4. The geometry run emits:

mod_code	offset	n_valid	n_mod	percent_modified
m	-1	4	2	50
m	0	4	1	25

The feature observation at position 10 is assigned offset -1, the lower observation at position 9 is assigned offset 0, and the higher observation at position 11 is omitted. The file is 70 bytes with SHA-256 10b337439dfcd2690d5438dd6520d997b002f531563311fd010c8086239eb1b9.

For mixed-regions.bed, the command reports that two regions were loaded, silently ignores the malformed physical line 3, and writes an aggregate with SHA-256 b0661bc358b949c07c6634d1f7b496c0dab8ecae77ad6b53ee286e0e0d6d3653.

The behavior is deterministic in repeated one-thread runs. Focused regressions also cover both contig edges, --window 0, BED3 through extended BED records, leading/interleaved blanks and comments, strict numeric/score/strand parsing, and output preservation.

Expected behavior

• Preserve the original feature midpoint as the anchor.
• Query [anchor-window, anchor+window+1), clipped only at contig bounds.
• Report bedMethyl start - original anchor; negative means a lower reference coordinate and positive means a higher coordinate, independent of feature strand.
• Validate every non-comment BED row with filepath and physical-line context before opening or truncating output.
• Preserve accepted tab-delimited BED names and opaque columns while rejecting inconsistent schemas, incomplete consumed fields, invalid coordinates, scores, and strands.
• Explain that older affected profiles must be regenerated. Documentation must not claim that negating old offsets alone repairs output, because the old query membership and reconstructed anchor can also differ.

Root-cause evidence

• Relevant code: modkit-core/src/localise/subcommand.rs and modkit-core/src/localise/util.rs at development revision 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c.
• Window mechanism: the query is constructed as [anchor-(window+1), anchor+window), shifting both intended bounds one base lower.
• Anchor mechanism: the original feature is overwritten by the clipped query interval, and the query midpoint is later reused as the biological anchor.
• Offset mechanism: aggregation calculates anchor - position rather than position - anchor.
• Input mechanism: parsing folds successful regions and errors into separate collections but returns success whenever the successful collection is nonempty; output is opened before validation.
• Parent-red/fix-green evidence: the parent emits the two incorrect geometry rows above and accepts the mixed file; the reviewed correction emits all three oracle rows, rejects physical line 3, and preserves an existing 29-byte output sentinel with SHA-256 23cdc9ad2b5b5f71b423f3d602e7b2e0aa4de668a359a71b8a2d3c077c866eb6.

Proposed fix scope

Introduce a localize-specific focus-region value that retains the original anchor separately from its clipped query interval. Parse the region file once with a consistent delimiter/schema, propagate every read/parse/coordinate error with physical-line context, and defer output creation until region validation succeeds. Document and test the reference-axis offset contract explicitly.

Non-goals

• Do not infer feature-oriented reversal for negative-strand BED features; that would require a separate explicit option and compatibility decision.
• Do not implement or redefine the currently separate --batch-size policy.
• Do not include chart scaling/sparse-rendering changes or regional bedMethyl query-error propagation in this correction.
• Continue to skip syntactically valid regions whose contigs are absent from the genome-sizes file or Tabix index, with existing diagnostics.
• Do not redesign bedMethyl aggregation or parallel scheduling for performance.

Acceptance criteria

• Interior, coordinate-zero, and both contig-edge fixtures use the exact clipped half-open window around the original anchor.
• --window 0 queries the anchor itself, and --window N owns exactly N positions on each available side.
• Positive-, negative-, and unstranded features use the same documented reference-axis offsets.
• A mixed valid/invalid region file exits nonzero with filepath, physical line, and reason before creating or modifying output.
• BED3 through extended tab-delimited records, spaces in names, opaque columns, benign blanks/comments, and supported legacy whitespace-delimited records have explicit compatibility coverage.
• Existing valid aggregation is deterministic across thread counts after accounting for the intentional coordinate correction.
• Focused localize tests and the applicable full workspace test suite pass.

Reproduction artifacts

Artifact	Size	SHA-256	Notes
geometry-counts.bed	157 bytes	049100275283a1b1c26bd2a4b6bffb8b7f1e371a856f06957b79d0b3f7b7c6ee	Generated inline above
geometry-counts.bed.gz	133 bytes	955a36248d98224661886caa2111aadf10723cc9800bf9f638571924b1c1475b	bgzip 1.23.1
geometry-counts.bed.gz.tbi	108 bytes	35cd5a6249f79be1d925ccb7a7539af1caec5ac831aa9b806a55e356a6f8a649	tabix 1.23.1
geometry-region.bed	11 bytes	20ab665a6aeab8a0c788ccf4f3181c8d805ebb0fdf2652ccb36e4cda1c002fb8	Single feature with anchor 10
mixed-regions.bed	103 bytes	9256f01cac70c6c98da8e03f108dab74c006fcc3d4803e2b624a111df6c59e3d	Valid, malformed, valid rows
genome-sizes.tsv	10 bytes	b52521d37f5929b00afbf5de9e8dda3acee1f83770eecdfed4c180c4a16b116d	Generated inline above

Related work

• Related issues: minimum-coverage issue localize --min-coverage does not filter bedMethyl records #651 and chart-rendering issue localize charts use fractional values on a percent axis and fail for sparse data #674 touch adjacent localize behavior but have independent contracts.
• Proposed PR: to be linked after prerequisite localize work is ordered, the migration documentation is corrected, and the exact rebased package is regated.
• Deferred policy question or upstream dependency: localize --batch-size is tracked in issue localize --batch-size is accepted but has no effect #653; dependency-resolution PR Pin hts-sys to bindings compatible with rust-htslib 0.46 #666 affects publication timing but not this scientific contract.