dmr pair imputes zero-coverage replicates as canonical in balanced statistics

[SuhasSrinivasan]

dmr pair imputes zero-coverage replicates as canonical in balanced single-site statistics

Summary

Single-site modkit dmr pair supports multiple samples per condition by repeating -a or -b. Under the default --min-valid-coverage 0, a bedMethyl row with valid coverage zero is retained as a sample observation. The balanced-replicate calculation divides its zero modified count by zero coverage, converts the resulting non-finite fraction to zero modified pseudo-counts, but still assigns that replicate the target balanced coverage. The no-observation replicate is therefore treated as entirely canonical.

The same row is also counted in pct_a_samples or pct_b_samples, even though the documentation defines those fields as the percentage of samples used in the statistical test. The ordinary pooled counts and effect can remain correct while the advertised balanced effect, balanced MAP-based p-value, and sample-representation percentage are silently wrong. If the zero-coverage row carries a modification code absent from every positive-coverage observation, its zero-valued category can additionally leak into pooled serialization and alter the categorical LLR score.

In the minimal reproduction below, condition A has one zero-coverage row and one 10/10 modified replicate; condition B has one 0/10 modified replicate. modkit reports the correct pooled effect 1, but a balanced effect of 0.5 and pct_a_samples = 100. Treating the zero-coverage row as a missing observation gives balanced effect 1 and pct_a_samples = 50.

Severity

Severity: Medium — silent scientific bias in replicate-balanced single-site results

Rationale: The command exits successfully with finite values, but no-observation replicates are assigned canonical pseudo-observations. This can change both the magnitude and direction of the balanced comparison and overstate sample representation. Exposure requires repeated single-site inputs and a retained zero-valid-coverage row. In the reproduced same-code fixture, pooled totals, fractions, effect, and MAP p-value remain numerically correct because the zero row adds no counts; however, a modification code present only in a zero-coverage row can still leak into pooled count or percentage serialization and can perturb the categorical LLR score. --min-valid-coverage 1 is an effective workaround.

User and scientific impact

• Affected workflow: single-site modkit dmr pair with repeated -a or -b inputs.
• Affected fields: balanced_map_pvalue, balanced_effect_size, pct_a_samples, and/or pct_b_samples.
• Secondary pooled-output exposure: if a modification code occurs only in a zero-coverage row, a zero-valued entry can remain in a_counts/b_counts or the corresponding modification-percentage field, and the extra category can alter score, instead of matching an absent observation.
• Direction of error: a zero-coverage replicate is pulled toward 0% modification. Depending on which condition contains it, this can attenuate or exaggerate the balanced effect.
• Likely exposure: zero-coverage bedMethyl rows are an edge case, but they are accepted by the default filter. Exposure accumulates across sparse multi-replicate tables and does not require a command failure or a non-finite output.
• Detectability or workaround: compare default output with --min-valid-coverage 1, or remove zero-coverage rows before indexing. The ordinary pooled effect does not reveal the balanced-field error.

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: three one-row bgzip/tabix-indexed bedMethyl files, one empty indexed control, and a two-base indexed FASTA.
• Fixture tools: bgzip/tabix 1.23.1 and samtools 1.23.1.

Steps to reproduce

Minimal input

At chr1:0, condition A has one replicate with no valid observations and one replicate with 10/10 modified observations. Condition B has one replicate with 0/10 modified observations.

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

# modified.bed
chr1	0	1	C	10	+	0	1	255,0,0	10	100.00	10	0	0	0	0	0	0

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

# ref.fa
>chr1
CC

Commands

work_dir="$(mktemp -d)"

printf '>chr1\nCC\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\t10\t+\t0\t1\t255,0,0\t10\t100.00\t10\t0\t0\t0\t0\t0\t0\n' \
  | bgzip -c > "$work_dir/modified.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/canonical.bed.gz"

tabix -p bed "$work_dir/zero.bed.gz"
tabix -p bed "$work_dir/modified.bed.gz"
tabix -p bed "$work_dir/canonical.bed.gz"

# An indexed sample with no record at this position is the missing-observation
# control.
bgzip -c /dev/null > "$work_dir/missing.bed.gz"
tabix -p bed "$work_dir/missing.bed.gz"

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

modkit dmr pair \
  -a "$work_dir/missing.bed.gz" \
  -a "$work_dir/modified.bed.gz" \
  -b "$work_dir/canonical.bed.gz" \
  -o "$work_dir/missing-row.bed" \
  --ref "$work_dir/ref.fa" \
  --header \
  --base C \
  --max-coverages 10 10 \
  --threads 1 \
  --io-threads 1 \
  --suppress-progress \
  --force

modkit dmr pair \
  -a "$work_dir/zero.bed.gz" \
  -a "$work_dir/modified.bed.gz" \
  -b "$work_dir/canonical.bed.gz" \
  -o "$work_dir/filtered-row.bed" \
  --ref "$work_dir/ref.fa" \
  --header \
  --base C \
  --min-valid-coverage 1 \
  --max-coverages 10 10 \
  --threads 1 \
  --io-threads 1 \
  --suppress-progress \
  --force

cut -f13-20 "$work_dir/zero-row.bed"
cut -f13-20 "$work_dir/missing-row.bed"
cmp "$work_dir/missing-row.bed" "$work_dir/filtered-row.bed"

Control or independent oracle

The zero-coverage replicate has no defined modification fraction and contributes no observation to an equal-replicate average. The available fractions are therefore:

condition A: [10 / 10] = [1]
condition B: [0 / 10] = [0]
balanced effect = mean(A) - mean(B) = 1 - 0 = 1
pct_a_samples = 1 observed / 2 configured = 50
pct_b_samples = 1 observed / 1 configured = 100

Because the available balanced counts are identical to the pooled counts in this fixture, the balanced MAP-based p-value must also equal the ordinary MAP-based p-value. An explicitly absent record and a record filtered by --min-valid-coverage 1 are independent controls for the same no-observation state.

Availability is determined by total valid coverage, not by the number of modified observations. A positive-coverage canonical sample such as 0/10 is therefore an observed sample: it contributes ten canonical observations and remains in both balancing and the represented-sample numerator. In contrast, 0/0, an absent record, and a record removed by the configured coverage filter are three representations of no observation.

Observed behavior

modkit 0.6.4 exits successfully and reports one successful site. Relevant fields from the zero-row run are:

a_pct_modified  b_pct_modified  map_pvalue                  effect_size  balanced_map_pvalue       balanced_effect_size  pct_a_samples  pct_b_samples
1               0               0.0000006230948043897833    1            0.020888505198542493      0.5                   100            100

Replacing the zero row with an absent observation produces:

a_pct_modified  b_pct_modified  map_pvalue                  effect_size  balanced_map_pvalue       balanced_effect_size  pct_a_samples  pct_b_samples
1               0               0.0000006230948043897833    1            0.0000006230948043897833  1                     50             100

The missing-row control and the --min-valid-coverage 1 output are byte-identical. The erroneous zero-row result is also byte-identical between compute/I/O thread settings 1/1 and 4/2.

Expected behavior

• A sample contributes to pooled or balanced single-site calculations only when its aggregated total valid coverage is positive and it satisfies the configured coverage filter; equivalently, coverage must be at least max(1, min_valid_coverage).
• A zero-coverage row is treated as no observation for pooled aggregation, balanced aggregation, and sample-representation percentages, consistently with an absent or positively filtered row.
• A positive-coverage canonical row remains an observed sample even when its modified count is zero.
• The target balancing coverage is calculated over observed positive-coverage samples only.
• If either condition has no positive-coverage observation, the site emits no result row and follows a recoverable failed-site path without producing NaN or panicking. A retained zero-only condition may report a calculation failure while an absent or filtered condition may report DmrMissing; identical diagnostic categories are not required.
• Numerical pooled totals, fractions, effect, and MAP p-value remain unchanged when a zero row contributed no counts. Zero-only modification-code entries should nevertheless disappear so pooled serialization and categorical LLR scoring match the absent-row state.
• For equal-size matched groups, sample slots remain tied to command-line identity. If any configured pair member is absent, zero coverage, or filtered, both whole replicate_* fields use the existing - sentinel; available samples may still produce the pooled and balanced group row.
• The output schema and behavior for all-positive-coverage inputs remain unchanged.

Root-cause evidence

• Repeated -a and -b arguments are public Vec<PathBuf> inputs at subcommands.rs:91-103.
• --min-valid-coverage defaults to zero at subcommands.rs:269-274, and records are retained when valid_coverage >= min_coverage at tabix.rs:163-178.
• AggregatedCounts::iter_mod_fractions divides each modified count by the sample total, yielding 0 / 0 for this row.
• collapse_counts computes target coverage across both retained records. Here it uses (0 + 10) / 2 = 5, converts the zero row to zero modified counts, and still gives it total coverage 5. Since canonical count is total minus modified count, the missing replicate becomes five canonical pseudo-observations. Combining it with the 5/5 modified replicate yields 5/10 modified and balanced effect 0.5.
• SingleSiteDmrScore::new_multi computes representation from vector length and passes the unfiltered vectors to balancing, so the zero row produces 2 / 2 = 100% representation.
• The documentation says these percentages describe samples used in the statistical test and recognizes that some samples have no modification fraction at a position at intro_dmr.md:222-241.

Proposed fix scope

• Distinguish positive-coverage observations from retained zero-coverage rows before both pooled and balanced group aggregation.
• Compute pooled counts, balanced target coverage, and balanced modification fractions only from positive-coverage observations in each condition.
• Use the number of positive-coverage observations, rather than raw vector length, for pct_a_samples and pct_b_samples.
• Preserve configured sample identity and command-order slots while filtering so zero-coverage removal cannot compact matched samples and create cross-pairs. Retain the matched-replicate all-or-none policy: incomplete pairs suppress both whole replicate fields without suppressing a valid group row.
• When either condition has no positive-coverage observation, emit no result row and use an existing recoverable per-site failure path; do not require zero-only and absent conditions to share the same diagnostic category.
• Add an end-to-end parent-red regression proving that the zero-row, absent-row, and --min-valid-coverage 1 controls have identical successful output fields for the fixture above.

Non-goals

• No change to the pooled estimator formula, MAP statistic, coverage capping, or output schema; the only pooled-input change is exclusion of observations with zero total coverage.
• No change to region-level dmr pair or dmr multi behavior.
• No compact partial replicate vector or new inline missing-value grammar. Matched-pair identity and the all-or-none - policy are specified by the separate matched-replicate report and are preserved when this repair filters group observations.
• No general change to zero-coverage handling in other bedMethyl consumers.

Acceptance criteria

• A zero-coverage replicate cannot contribute canonical or modified pseudo-counts to balanced statistics.
• A positive-coverage canonical replicate remains present in pooled and balanced statistics and in the represented-sample numerator.
• For the fixture above, balanced effect is exactly 1, balanced MAP-based p-value equals the ordinary MAP-based p-value, pct_a_samples is 50, and pct_b_samples is 100.
• For the fixture above, the complete successful zero-row output matches the absent-row and --min-valid-coverage 1 controls.
• A zero-coverage row whose modification code is otherwise absent does not leave a zero-valued code entry in pooled count or percentage serialization and does not change the categorical LLR score.
• Reversing the zero-coverage condition does not introduce a directional imputation bias.
• A condition with no positive-coverage observation emits no result row and follows a recoverable failed-site path; diagnostic category and progress text need not be byte-identical between retained-zero, absent, and filtered cases.
• In matched groups, any absent, zero-coverage, or coverage-filtered pair member makes both entire replicate_map_pvalues and replicate_effect_sizes fields exactly -; no remaining samples are compacted into new pairs, and a valid group row is retained.
• Existing all-positive-coverage multi-sample output remains byte-identical.
• Results are byte-identical across supported compute and I/O thread counts.
• Focused DMR tests and the full workspace suite pass.

Reproduction artifacts

Artifact	Size	SHA-256	Notes
ref.fa	9 bytes	dc95419df161b61d6cb661be877d6fdc3df1fabf3e97bb1bcc3a277d0666d44a	Two-base chr1 reference
ref.fa.fai	13 bytes	6db28c377fc1aa8b697be3bc00fd617943e98090313ac8051ecba98568a7944a	FASTA index
zero.bed.gz	86 bytes	2bba46a459d2063872709db9b5ef357ecca00ea40d0705e7be3b76106c037514	Zero-valid-coverage C row
zero.bed.gz.tbi	109 bytes	c97f9d47a106d4fcedf718a08a9cb42aa36297ab59d69e5f0b03e347c2042205	Tabix index
modified.bed.gz	91 bytes	d0cda187d5c18ce4041008863ed4cfe8884950a7d45d3b3d4bcbef3b3e788c34	10/10 modified C row
modified.bed.gz.tbi	109 bytes	f57b938f9b5433ccd02b21bf3735405177f5af6446ceec0a3f386b11671bc6bd	Tabix index
canonical.bed.gz	90 bytes	6a2a2e515d35f1d130e260165ac9fb9b21b9c892330ddce1e7994aa50f14fb5d	0/10 modified C row
canonical.bed.gz.tbi	109 bytes	894ac2ee3cb92e5e0d79d709acac8ff19e2a3c2c9d5a857d128e424b41db3a6d	Tabix index
missing.bed.gz	28 bytes	d079906378251d29409f7f0f691113fd3b0049b5927b713e2704ca44ac743d48	Empty bgzip missing-observation control
missing.bed.gz.tbi	75 bytes	1a5845ff0c599ef381032688bd8be8b755436ce4721e2ef28fef4690a700705d	Empty Tabix index
zero-row.bed	441 bytes	240e4cd7f63197b566914549db781aee32c5a2a86ca25344d838ddc7fedfbdc0	Erroneous successful output
missing-row.bed	442 bytes	9bbac869ad0de081739f7c0bbaded7a14be72319ae835d44d9ddbff1c0a1200c	Independent missing-observation oracle

Related work

• Duplicate searches used dmr zero coverage balanced, balanced effect zero coverage, pct_a_samples, and zero-coverage replicate. No matching upstream report was found.
• Issue #361 discusses general use and interpretation of balanced DMR outputs but does not report no-observation imputation or representation inflation.
• A separate zero-total-coverage defect concerns a condition with no observations and can produce non-finite output or a worker abort. This report remains separate because the pooled totals here are positive, the command emits finite successful output, and the defect is in replicate balancing and representation accounting.
• Missing matched-pair identity and ordering affect the optional replicate_* fields through a different code path and remain a distinct finding. The two repairs must nevertheless compose: this group-balancing correction preserves configured slots and the separately selected all-or-none matched-pair policy rather than compacting the positive observations into new pairs.