`dmr pair --segment` omits the final successfully scored site from each HMM chunk

[SuhasSrinivasan]

dmr pair --segment omits the final successfully scored site from each HMM chunk

Summary

modkit dmr pair --segment successfully writes every single-site score but decodes only N-1 HMM states for N successful scores. The final genomic position in each HMM chunk is therefore absent from the segment output, and the final decoded state can be assigned to the preceding position.

On the public chromosome 20 test inputs, modkit 0.6.4 exits successfully after reporting 17,271 successfully processed sites, while the sole segment reports only 17,270 sites and ends at the penultimate coordinate.

Severity

Severity: High — scientific correctness and output completeness

Rationale: The command reports success while omitting a successfully scored site from every HMM chunk. Segment boundaries, site counts, aggregated modification counts, effect sizes, and potentially the final state assignment can therefore be wrong without an error or warning. A user can detect the omission only by reconciling the site and segment artifacts.

User and scientific impact

• Affected result or workflow: HMM segmentation requested with modkit dmr pair --segment.
• Direction of error: one omitted genomic site per nonempty HMM chunk, a shortened final interval, undercounted segment aggregates, and possible state misassignment near the chunk end.
• Likely exposure: routine whenever segmentation is enabled; the number of omitted sites grows with the number of chunks created by chromosome or maximum-gap boundaries.
• Detectability or workaround: compare the number of successful site rows with the sum of segment num_sites and compare the last site and segment ends. Running without --segment preserves the single-site rows but does not provide segmented DMRs.

Affected versions and environment

• Released version: modkit 0.6.4.
• Development revision inspected: 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c.
• Operating system and architecture used for reproduction: macOS on Apple Silicon.
• Inputs: bgzip/tabix-indexed bedMethyl files and an indexed FASTA reference from the public test resources.

Steps to reproduce

From the modkit repository:

work_dir="$(mktemp -d)"

modkit dmr pair \
  -a tests/resources/lung_00733-m_adjacent-normal_5mc-5hmc_chr20_cpg_pileup.bed.gz \
  -b tests/resources/lung_00733-m_primary-tumour_5mc-5hmc_chr20_cpg_pileup.bed.gz \
  -o "$work_dir/sites.bed" \
  --segment "$work_dir/segments.bed" \
  --ref tests/resources/GRCh38_chr20.fa \
  --header \
  --base C \
  --max-coverages 100 100 \
  --threads 1 \
  --io-threads 1 \
  --suppress-progress \
  --force

awk 'NR > 1 { n += 1 } END { print "site rows:", n }' "$work_dir/sites.bed"
awk 'NR > 1 { n += $6 } END { print "segment num_sites total:", n }' "$work_dir/segments.bed"
tail -n 1 "$work_dir/sites.bed"
tail -n 1 "$work_dir/segments.bed"
shasum -a 256 "$work_dir/sites.bed" "$work_dir/segments.bed"

Control or independent oracle

There is one emission score for every successfully scored genomic position. The HMM therefore must return exactly one ordered state per score, and converting those states to half-open BED regions must cover every input position exactly once.

For this one-chunk input, the control invariants are:

number of decoded states = number of site rows = 17,271
sum of segment num_sites = 17,271
last segment end = last site end = 10,804,378

An independent exhaustive oracle can verify tiny HMMs by enumerating the hidden start state and every emitted state, scoring each possible path without using the dynamic-programming matrix or its back-pointers.

Observed behavior

modkit 0.6.4 exits with status 0 and reports:

finished, processed 17271 sites successfully, 175 failed
site rows: 17271
segment num_sites total: 17270

The last successful site is:

chr20  10804377  10804378  ...

The only segment ends at the preceding successful coordinate:

chr20  9701049  10804085  same  907.0660975771025  17270  ...

The artifacts are deterministic for this invocation:

f9c0e0566b5a78c0bd868f11b53af23058e0c0506f6d78764be707ff8c7df982  sites.bed
a387a5ae5f7c5721f269e0a9869702e7af54b8df66a8280f876e90617d8adf57  segments.bed

Expected behavior

• Decode exactly one correctly ordered HMM state for every successful score.
• Cover every input genomic position exactly once in the resulting half-open segments, including singleton chunks.
• Make the sum of segment num_sites equal the number of successfully segmented site rows.
• End the final segment at one base past the final represented site.
• Preserve deterministic site and segment bytes across supported compute and I/O thread counts.

Root-cause evidence

• HmmModel::viterbi_path explicitly asserts the incorrect scores.len() - 1 cardinality.
• HmmModel::viterbi_decode obtains the predecessor of the final state but does not append that predecessor before querying the preceding pointer cell. It later removes one path entry, yielding N-1 ordered outputs for N emissions.
• viterbi_forward constructs one hidden start cell plus one dynamic-programming cell and pointer per emission, confirming that the decode must return N emitted states.
• path_to_region_labels repeats the positions.len() - 1 assumption and zips the shortened path with the positions, silently leaving the last position uncovered.
• Segment counts and aggregate statistics are computed from those shortened region bounds at single_site.rs:1107-1151.

The parent behavior is reproduced both by the released binary and by focused parent-red tests: singleton inputs return no emitted state/region, and multi-site inputs return one fewer state than scores.

Proposed fix scope

• Backtrace from the final dynamic-programming state through exactly the N-1 emitted-state transitions, retaining the final state, so the returned path contains N states in genomic order.
• Require path.len() == positions.len() when constructing segment regions.
• Add an independent exhaustive tiny-path oracle that checks both state order and cardinality.
• Add singleton, two-position, multi-position, and state-transition region tests that require every input position to belong to exactly one half-open region.
• Add a real CLI regression that reconciles site rows, segment num_sites, and the final interval end under at least --threads 1 --io-threads 1 and --threads 4 --io-threads 2.

Non-goals

• No change to emission scores, transition probabilities, distance projection, HMM priors, or state labels.
• No output-schema or interval-convention change.
• No performance claim beyond retaining linear-time backtrace and region integration.

Acceptance criteria

• For every nonempty score vector of length N, decoding returns exactly N states in the order selected by an independent exhaustive oracle.
• Singleton segmentation emits one one-base half-open region.
• Every position in multi-site inputs is covered exactly once by a region with the decoded state.
• The public chromosome 20 reproduction emits 17,271 site rows, a segment-site total of 17,271, and final segment end 10,804,378.
• Site and segment outputs are byte-identical between the 1/1 and 4/2 compute/I/O thread configurations.
• Existing DMR error propagation and successful non-segmented output remain unchanged.
• Focused HMM/DMR tests and the full workspace suite pass.

Reproduction artifacts

Artifact	Size	SHA-256	Notes
tests/resources/lung_00733-m_adjacent-normal_5mc-5hmc_chr20_cpg_pileup.bed.gz	253,043 bytes	95b5d56cb595f08a1cc7973c1e7124ecfaf90860b315afaa37d810964edd2c78	Public condition A input
tests/resources/lung_00733-m_primary-tumour_5mc-5hmc_chr20_cpg_pileup.bed.gz	251,296 bytes	3b55a523bf334d68de93c58736147197464c1db96d9888d42a59b2bc355defac	Public condition B input
tests/resources/GRCh38_chr20.fa	65,249,832 bytes	5806b504166c43762e11c809f720b28b6d50462104fea22e7490deffbe6b59cf	Public reference

Related work

• Issue #448 asks broadly about single-site positions that appear absent from segmentation output, and the maintainer noted that such missing positions would sound like a bug. This report isolates a narrower deterministic implementation defect: exactly the final successfully scored site of each nonempty HMM chunk is omitted, with a source-level cause and independent path oracle.
• A completed tests-first fix is prepared. Its PR will be linked only after the prerequisite branch chain is rebased onto current upstream and the exact recomposed head passes a fresh full-workspace gate.
• The separate zero-coverage DMR scoring defect has different triggers, root causes, severity, and test fixtures; it should be reviewed independently.