`dmr pair` emits NaN results or aborts when input or coverage capping creates a zero DMR denominator

[SuhasSrinivasan]

dmr pair emits NaN results or aborts when input or coverage capping creates a zero DMR denominator

Summary

modkit dmr pair retains bedMethyl records whose valid coverage is zero under its default filtering policy, then computes modification fractions by dividing by that coverage. If one condition has a zero-coverage record and the other has positive canonical-only coverage at the same site, the default prior produces a successful output row containing NaN scientific values. The accepted boundary prior --prior 0.5 0.5 sends the same undefined site into a singular beta-density calculation, panics inside a Rayon worker, and aborts the process with status 134.

The same zero-denominator path is reachable even when both scored input records have positive coverage. If either effective maximum coverage is zero—through an accepted --max-coverages value or an automatically estimated 95th percentile of zero—coverage resizing converts the positive count to (0, 0) and again computes 0 / 0. A repair limited to rejecting raw zero coverage in Counts::new is therefore incomplete because resizing happens afterward.

Severity

Severity: Medium — scientific output validity and reliability on accepted edge-case input

Rationale: With the default prior, the command exits successfully and labels an undefined comparison as a successfully processed site while emitting NaN scientific values. With a boundary prior, the same zero-denominator state can abort the process. Exposure includes retained raw zero-coverage rows, explicit maximum-coverage vectors containing zero, and sparse automatic-cap inputs whose 95th percentile floors to zero. These are different routes to the same undefined DMR denominator, so severity remains Medium and this extension does not add another finding. --min-valid-coverage 1 handles the original raw-zero fixture but is not a universal automatic-cap workaround because too few positive observations can remain; explicitly positive maximum coverages are the deterministic control.

User and scientific impact

• Affected result or workflow: single-site modkit dmr pair scoring when either condition has raw total coverage zero or either final effective maximum coverage is zero.
• Direction of error: false-success output containing NaN scientific fields under the default prior, or worker panic and process abort under a boundary prior.
• Likely exposure: uncommon but accepted by the current CLI because --min-valid-coverage defaults to zero, explicit maximum-coverages accept zero, and automatic estimation retains zero coverage values.
• Detectability or workaround: the boundary-prior panic is visible, but default-prior results are reported as successful and require explicit NaN or resolved-cap checks to detect. Removing zero-coverage rows or passing --min-valid-coverage 1 avoids the original raw-zero path, while explicit positive maximum coverages avoid the cap path.

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: raw-zero and positive one-row bgzip/tabix-indexed bedMethyl files, a 21-site sparse automatic-cap fixture, and indexed FASTA references.
• Fixture tools used: bgzip, tabix, and samtools faidx from the modkit conda environment.

Steps to reproduce

Minimal input

The two bedMethyl records describe the same positive-strand C at chr1:0. The first has valid coverage zero; the second has one canonical observation and no modified observations.

# zero.bed
chr1	0	1	C	0	+	0	1	255,0,0	0	0.00	0	0	0	0	0	0	0

# canonical.bed
chr1	0	1	C	1	+	0	1	255,0,0	1	0.00	0	1	0	0	0	0	0

# ref.fa
>chr1
C

Commands

work_dir="$(mktemp -d)"

printf '>chr1\nC\n' > "$work_dir/ref.fa"
samtools faidx "$work_dir/ref.fa"

printf 'chr1\t0\t1\tC\t0\t+\t0\t1\t255,0,0\t0\t0.00\t0\t0\t0\t0\t0\t0\t0\n' \
  | bgzip -c > "$work_dir/zero.bed.gz"
printf 'chr1\t0\t1\tC\t1\t+\t0\t1\t255,0,0\t1\t0.00\t0\t1\t0\t0\t0\t0\t0\n' \
  | bgzip -c > "$work_dir/canonical.bed.gz"
tabix -p bed "$work_dir/zero.bed.gz"
tabix -p bed "$work_dir/canonical.bed.gz"

modkit dmr pair \
  -a "$work_dir/zero.bed.gz" \
  -b "$work_dir/canonical.bed.gz" \
  -o "$work_dir/default-prior.bed" \
  --ref "$work_dir/ref.fa" \
  --base C \
  --max-coverages 1 1 \
  --threads 1 \
  --io-threads 1 \
  --suppress-progress \
  --force

set +e
RUST_BACKTRACE=1 modkit dmr pair \
  -a "$work_dir/zero.bed.gz" \
  -b "$work_dir/canonical.bed.gz" \
  -o "$work_dir/boundary-prior.bed" \
  --ref "$work_dir/ref.fa" \
  --base C \
  --prior 0.5 0.5 \
  --max-coverages 1 1 \
  --threads 1 \
  --io-threads 1 \
  --suppress-progress \
  --force
boundary_exit_code=$?
set -e

cat "$work_dir/default-prior.bed"
printf 'boundary exit status: %s\n' "$boundary_exit_code"
wc -c "$work_dir/boundary-prior.bed"

Positive-input coverage-cap extension

The same estimator receives a zero denominator after coverage resizing when the source rows are both positive. Create 10/10 modified and 0/10 modified inputs:

printf 'chr1\t0\t1\tC\t10\t+\t0\t1\t255,0,0\t10\t100.00\t10\t0\t0\t0\t0\t0\t0\n' \
  | bgzip -c > "$work_dir/modified10.bed.gz"
printf 'chr1\t0\t1\tC\t10\t+\t0\t1\t255,0,0\t10\t0.00\t0\t10\t0\t0\t0\t0\t0\n' \
  | bgzip -c > "$work_dir/canonical10.bed.gz"
tabix -p bed "$work_dir/modified10.bed.gz"
tabix -p bed "$work_dir/canonical10.bed.gz"

for caps in "0 0" "0 10" "10 0" "10 10"; do
  read -r a_cap b_cap <<< "$caps"
  modkit dmr pair \
    -a "$work_dir/modified10.bed.gz" \
    -b "$work_dir/canonical10.bed.gz" \
    -o "$work_dir/caps-${a_cap}-${b_cap}.bed" \
    --ref "$work_dir/ref.fa" \
    --header \
    --base C \
    --max-coverages "$a_cap" "$b_cap" \
    --threads 1 \
    --io-threads 1 \
    --suppress-progress \
    --force
done

For the automatic-cap route, create twenty zero-coverage sites followed by one positive canonical observation. Use the same file for both conditions to isolate cap estimation from effect-size direction:

printf '>chr1\nCCCCCCCCCCCCCCCCCCCCC\n' > "$work_dir/ref21.fa"
samtools faidx "$work_dir/ref21.fa"

for pos in {0..19}; do
  end=$((pos + 1))
  printf 'chr1\t%s\t%s\tC\t0\t+\t%s\t%s\t255,0,0\t0\t0.00\t0\t0\t0\t0\t0\t0\t0\n' \
    "$pos" "$end" "$pos" "$end"
done > "$work_dir/automatic-zero-cap.bed"
printf 'chr1\t20\t21\tC\t1\t+\t20\t21\t255,0,0\t1\t0.00\t0\t1\t0\t0\t0\t0\t0\n' >> "$work_dir/automatic-zero-cap.bed"
bgzip -c "$work_dir/automatic-zero-cap.bed" > "$work_dir/automatic-zero-cap.bed.gz"
tabix -p bed "$work_dir/automatic-zero-cap.bed.gz"

modkit dmr pair \
  -a "$work_dir/automatic-zero-cap.bed.gz" \
  -b "$work_dir/automatic-zero-cap.bed.gz" \
  -o "$work_dir/automatic-zero-cap-output.bed" \
  --ref "$work_dir/ref21.fa" \
  --header \
  --base C \
  -N 21 \
  --interval-size 4 \
  --threads 1 \
  --io-threads 1 \
  --suppress-progress \
  --force

Control or independent oracle

With a shared prior Beta(0.5, 0.5), posterior alpha parameters are:

alpha_A = 0.5 + n_mod_A = 0.5
alpha_B = 0.5 + n_mod_B = 0.5
alpha_A + alpha_B = 1

The closed-form density at zero evaluates ln_beta(alpha_A + alpha_B - 1, ...), so this boundary would call ln_beta(0, ...). The function requires strictly positive arguments and must not be called at equality.

When both conditions have positive canonical-only coverage, their empirical effect is exactly zero and the estimator returns before evaluating this density. The zero-coverage condition instead computes 0 / 0 = NaN; comparisons with NaN are false. With the default Beta(0.55, 0.55) prior, the later calculations propagate NaN into a nominally successful result. With Beta(0.5, 0.5), the null-density call also reaches the singular pair-sum boundary. Both the default --delta 0.05 and --delta 1 reproduce that abort.

A positive input count is resized when its coverage exceeds the effective cap. At a cap of zero, the resized modified count and total are both zero, so the resized fraction is 0 / 0; no MAP estimand exists. Zero is therefore outside the valid effective-cap domain. It must be rejected, not converted to one or treated as an ordinary observation filter.

With A=10/10 modified and B=0/10 modified, explicit cap vectors [0,0], [0,10], and [10,0] each produce a nominally successful row whose MAP p-value and effect size are NaN. [10,10] produces MAP p-value 0.0000006230948043897833 and effect 1.

Both conditions use the same 21-row fixture containing twenty 0/0 records followed by one positive-coverage canonical 0/1 record. The fixture isolates automatic maximum-coverage estimation: its 95th percentile is zero in both conditions. It is separate from the explicit A=10/10 versus B=0/10 numerical control. On the intermediate carrier that already rejects raw-zero source counts, the automatic-cap command reports twenty recoverable raw-zero failures plus one successful NaN row for the positive-coverage site. With explicit caps [10,10], that carrier retains the twenty failures and emits one finite canonical-versus-canonical result with MAP p-value 1 and effect 0.

Observed behavior

With the default prior, modkit 0.6.4 exits with status 0 and reports one successful site:

finished, processed 1 sites successfully, 0 failed
chr1  0  1  .  0  +  C:0  0  C:0  1  C:NaN  C:0.00  NaN  0  NaN  NaN  NaN  NaN  NaN

With --prior 0.5 0.5, modkit 0.6.4 prints the following panic and then aborts:

thread '<unnamed>' panicked at .../statrs-0.16.1/src/function/beta.rs:21:27:
called `Result::unwrap()` on an `Err` value: ArgMustBePositive("a")
...
statrs::function::beta::ln_beta
mod_kit::dmr::beta_diff::PMapEstimator::calc_beta_diff
mod_kit::dmr::beta_diff::PMapEstimator::predict
...
Rayon: detected unexpected panic; aborting
exit status: 134
0 .../boundary-prior.bed

The boundary-prior abort was reproduced with the default delta and with --delta 1.

The exact source repair that rejects raw zero coverage still accepts every explicit zero-cap vector under the default prior, exits zero, reports one successful site, and emits NaN for the MAP p-value and effect. The three outputs are byte-identical, SHA-256 db1ba76e9a02d00017940d087e2f9cff8f24df3743252829123e4f1a0571cb28. Its [10,10] control is finite and exact. On the same-file automatic fixture, that repair rejects the twenty source-zero sites but still resizes the remaining positive site to zero coverage and emits it as a successful NaN result. The explicit-cap signatures are byte-identical across compute/I/O thread settings 1/1 and 4/2.

Expected behavior

A condition whose raw pooled coverage is zero should continue through the recoverable per-site calculation-error path, producing no row for that site. For the minimal fixture, the command should complete normally, report one failed site, and emit no data row under both the default and boundary priors.

An effective maximum-coverage component of zero is instead an invalid global estimator configuration. Whether supplied explicitly or derived automatically, it should cause a normal nonzero analysis-construction error before scoring begins. No scientific row may be emitted. Existing behavior remains byte-identical for representable, strictly positive effective caps. Replicate-scaled caps that overflow usize now return a normal construction error instead of wrapping or panicking.

Valid positive-coverage inputs using --prior 0.5 0.5 must remain accepted: the CLI prior domain is alpha > 0, beta > 0, and alpha + beta >= 1, while the closed-form density at zero independently requires both paired posterior sums to be strictly greater than one.

Root-cause evidence

• --min-valid-coverage defaults to zero at subcommands.rs:269-274, and the tabix reader retains records with valid_coverage >= min_coverage at tabix.rs:163-178.
• --max-coverages accepts zero-valued components at subcommands.rs:307-310. Automatic estimation collects every retained valid_coverage, including zero, and can floor its 95th percentile to zero at single_site.rs:858-977.
• PMapEstimator::new stores either explicit or automatically resolved zero caps without validating their denominator domain.
• Counts::new permits (n_mod, coverage) = (0, 0) and computes 0 / 0, producing NaN.
• Counts::new is invoked before coverage capping. Therefore, even after adding a raw-zero guard there, Counts::resize can still construct (0, 0) directly and divide by zero when the effective cap is zero.
• The estimator's early return and later null-density evaluation are at beta_diff.rs:226-265. NaN.abs() <= rope is false, and later probability calculations either propagate NaN or reach the null-density boundary.
• The boundary guard checks < 1 rather than <= 1 before passing alpha_1 + alpha_2 - 1 and beta_1 + beta_2 - 1 to ln_beta at beta_diff.rs:168-189. At equality, statrs 0.16.1 unwraps an invalid nonpositive argument and panics.
• The CLI accepts a positive prior whose parameters sum to one at single_site.rs:76-87, although its error message states a different boundary.

Proposed fix scope

• Reject zero total coverage in the DMR count-to-fraction boundary before performing 0 / 0. Return the condition through the established per-site calculation-error path so no scientific NaN row is emitted.
• Make PMapEstimator::new reject either final effective maximum coverage equal to zero and return a propagated construction error.
• Perform that check after replicate scaling and before the estimator is returned or used. Use checked multiplication while resolving replicated coverages, so every constructed estimator guarantees strictly positive, non-overflowed effective caps.
• Retain raw zero-total rejection in Counts::new; the source-count and resolved-cap guards protect different boundaries.
• Also reject posterior alpha-pair or beta-pair sums equal to one before evaluating the zero-difference closed form; the finite-domain guard must use <= 1.
• Keep the user-specified prior contract distinct: each parameter must be positive and their sum may equal one.
• Align CLI help, error text, and scoring documentation with those two domains.
• Add parent-red end-to-end regressions using the exact zero-coverage/canonical-only fixture above: the default prior currently produces a successful NaN row, and --prior 0.5 0.5 currently aborts. Both must become the same normal failed-site result.
• Add explicit [0,0], [0,10], [10,0], and positive [10,10] unit/CLI controls plus the automatic twenty-zero-and-one-positive fixture.
• Retain a positive-coverage --prior 0.5 0.5 CLI control and an independently calculated interior-density control.

Non-goals

• No change to the MAP estimand, posterior update, effect-size adjustment, quadrature, or results for representable, strictly positive effective coverage caps.
• Do not clamp a zero effective cap to one.
• Do not silently change automatic-cap sampling to discard zero values or redefine its percentile/minimum-datapoint policy in this repair.
• Replicate-scaled caps that overflow usize are outside the preserved positive-cap domain and return a construction error.
• No general policy change applies to other bedMethyl consumers or records below a user-selected positive minimum. Raw-zero DMR observations are rejected at the count-to-fraction boundary, while zero resolved coverage caps are rejected at estimator construction.
• No change to HMM segmentation or its site-completeness repair.

Acceptance criteria

• The default-prior public CLI regression reports zero successful sites and one failed site, emits no data row, and contains no NaN scientific output.
• The boundary-prior version of the same regression no longer panics or aborts under the default delta or --delta 1; it follows the same failed-site policy.
• Zero total coverage is rejected before any modification fraction or downstream DMR statistic is computed.
• Explicit cap vectors [0,0], [0,10], and [10,0] return a normal nonzero constructor error and emit no scientific data row.
• With [10,10], A=10/10 versus B=0/10 retains MAP p-value 0.0000006230948043897833, effect 1, and existing output bytes.
• The twenty-zero-plus-one-positive automatic fixture reports the derived zero caps and exits cleanly with a constructor error; it cannot report the positive site as a successful NaN result.
• The same sparse same-file fixture with explicit [10,10] retains twenty recoverable raw-zero failures and one finite canonical-versus-canonical result with MAP p-value 1 and effect 0.
• The original raw-zero fixture remains a recoverable per-site failure rather than becoming a command-level configuration error.
• No scientific output contains NaN, and neither failure path panics or aborts the process.
• Equality at either posterior alpha-pair or beta-pair boundary returns an ordinary error before calling ln_beta with a nonpositive argument.
• Positive-coverage inputs with shared prior 0.5 0.5 continue to succeed.
• Priors with positive parameters and sum below one are rejected; positive parameters with sum equal to or above one retain the documented behavior.
• An interior closed-form value is checked against an independent exact oracle; for posterior parameters (2, 3) and (4, 5), the log-density ratio is ln(16/11).
• Existing DMR output on valid positive-coverage fixtures remains byte-identical.
• Focused tests and the full workspace suite pass.

Reproduction artifacts

Artifact	Size	SHA-256	Notes
ref.fa	8 bytes	bdaa0b566f868a175805d61a3e04ad59ab51951543111703233aae0841f0f923	One-base chr1 reference
ref.fa.fai	13 bytes	930ebf9ac3bd56acfde9b26962ede12cc26cf7584454838addbdaa8ed7120810	samtools faidx index
zero.bed.gz	86 bytes	2bba46a459d2063872709db9b5ef357ecca00ea40d0705e7be3b76106c037514	One zero-valid-coverage C record, bgzip 1.23.1
zero.bed.gz.tbi	109 bytes	c97f9d47a106d4fcedf718a08a9cb42aa36297ab59d69e5f0b03e347c2042205	Tabix BED index
canonical.bed.gz	87 bytes	133b02b5baffafaac6566abd5ea9a325b180dbd32d062c1708a18f219fb2bf24	One canonical-only C record, bgzip 1.23.1
canonical.bed.gz.tbi	109 bytes	bc3b84bd63132fac713edee64984b160c4e0026ac1ef21cfed5d04acd06ffbd0	Tabix BED index
default-prior.bed	66 bytes	f06354a4f868e126e9010d008a1f2f848ca844f6866486da8f7cdf8a900aab6b	Successful installed output containing NaN fields
boundary-prior.bed after abort	0 bytes	e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855	Empty output created before worker abort

Related work

• Proposed PR: to be linked after publication. The initial Counts::new and posterior-boundary carrier is incomplete until the effective-cap constructor guard and its focused regressions are integrated, independently reviewed, and fully regated.
• The HMM final-site omission is a separate High-severity correctness issue with unrelated code, inputs, and acceptance criteria.