dmr pair loses sample identity in matched replicate statistics

[SuhasSrinivasan]

dmr pair loses sample identity and cross-pairs matched replicate statistics

Summary

Single-site modkit dmr pair documents that repeated -a and -b inputs with equal group sizes are paired by command-line position: a1 with b1, a2 with b2, and so on. The input reader initially assigns each sample an integer identity, but the single-site organizer discards those identities while collecting per-position counts into separate vectors for conditions A and B. The replicate calculation then zips those independently ordered vectors.

As a result, the optional replicate_map_pvalues and replicate_effect_sizes fields can describe different cross-sample comparisons from those requested. This occurs with complete positive-coverage inputs: for three configured pairs whose expected effects are [1, -0.5, -0.5], modkit 0.6.4 reports [0.5, -1, 0.5]. Running the three intended comparisons separately produces the documented command-order oracle.

Missing observations expose the same identity loss more directly. If only a2 and b1 have records at a position, there is no complete configured pair, but modkit zips those remaining records and emits a replicate result for the unrequested a2-versus-b1 comparison.

Severity

Severity: Medium — silent scientific misassignment in matched-replicate output

Rationale: The command exits successfully and emits finite, plausible per-replicate statistics under the wrong sample pairing. The ordinary pooled and balanced group-level fields remain correct, limiting the impact to the optional matched-replicate columns. However, complete three-replicate inputs reproduce the defect without malformed or edge-case data, and users may rely on these fields to assess replicate concordance or identify sample-specific effects.

User and scientific impact

• Affected workflow: single-site modkit dmr pair with equal numbers of repeated -a and -b inputs.
• Affected fields: replicate_map_pvalues and replicate_effect_sizes.
• Direction of error: statistics are assigned to unrequested cross-sample comparisons. Signs and magnitudes can both change.
• Likely exposure: routine for at least some three-or-more-replicate layouts; missing observations can expose cross-pairing even with two configured pairs.
• Detectability or workaround: run each intended a_i/b_i pair in a separate command and join the per-pair results externally, or ignore the replicate fields. Filtering missing records alone is not sufficient because sample identity has already been discarded.

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

Steps to reproduce

Minimal input

The three bedMethyl files contain 100%, 50%, and 0% modified observations at the same positive-strand C.

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

# half.bed
chr1	0	1	C	10	+	0	1	255,0,0	10	50.00	5	5	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\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\t50.00\t5\t5\t0\t0\t0\t0\t0\n' \
  | bgzip -c > "$work_dir/half.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/modified.bed.gz"
tabix -p bed "$work_dir/half.bed.gz"
tabix -p bed "$work_dir/canonical.bed.gz"

# Configured pairs and expected effects:
#   modified  - canonical =  1
#   canonical - half      = -0.5
#   half      - modified  = -0.5
modkit dmr pair \
  -a "$work_dir/modified.bed.gz" \
  -a "$work_dir/canonical.bed.gz" \
  -a "$work_dir/half.bed.gz" \
  -b "$work_dir/canonical.bed.gz" \
  -b "$work_dir/half.bed.gz" \
  -b "$work_dir/modified.bed.gz" \
  -o "$work_dir/three-pairs.bed" \
  --ref "$work_dir/ref.fa" \
  --header \
  --base C \
  --max-coverages 10 10 \
  --threads 1 \
  --io-threads 1 \
  --suppress-progress \
  --force

# Independent command-order controls.
modkit dmr pair -a "$work_dir/modified.bed.gz"  -b "$work_dir/canonical.bed.gz" -o "$work_dir/pair1.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/canonical.bed.gz" -b "$work_dir/half.bed.gz"      -o "$work_dir/pair2.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/half.bed.gz"      -b "$work_dir/modified.bed.gz"  -o "$work_dir/pair3.bed" --ref "$work_dir/ref.fa" --header --base C --max-coverages 10 10 --threads 1 --io-threads 1 --suppress-progress --force

cut -f21-22 "$work_dir/three-pairs.bed"
for fp in "$work_dir"/pair{1,2,3}.bed; do
  tail -n 1 "$fp" | cut -f15-16
done

# Crossed-missing extension: neither configured pair has both observations.
bgzip -c /dev/null > "$work_dir/missing.bed.gz"
tabix -p bed "$work_dir/missing.bed.gz"

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

# Aligned control: the first configured pair has both observations.
modkit dmr pair \
  -a "$work_dir/modified.bed.gz" \
  -a "$work_dir/missing.bed.gz" \
  -b "$work_dir/canonical.bed.gz" \
  -b "$work_dir/missing.bed.gz" \
  -o "$work_dir/aligned-pair-control.bed" \
  --ref "$work_dir/ref.fa" --header --base C \
  --max-coverages 10 10 --threads 1 --io-threads 1 \
  --suppress-progress --force

cmp "$work_dir/cross-pair-missing.bed" "$work_dir/aligned-pair-control.bed"

Control or independent oracle

The documentation defines pairing by the relative order of repeated arguments. The three pairwise effects can also be calculated directly from the input fractions:

pair 1: 1.0 - 0.0 =  1.0
pair 2: 0.0 - 0.5 = -0.5
pair 3: 0.5 - 1.0 = -0.5

Running each intended pair separately is independent of the multi-sample organization and produces the expected MAP/effect sequence:

0.0000006230948043897833   1
0.020888505198542493      -0.5
0.01995939262165838       -0.5

Observed behavior

modkit 0.6.4 exits successfully but the combined command reports:

replicate_map_pvalues:
0.01995939262165838,0.0000006230948043897833,0.020888505198542493

replicate_effect_sizes:
0.5,-1,0.5

Those values correspond to cross-pairing the same sample fractions rather than the configured comparisons. The output is byte-identical at compute/I/O thread settings 1/1 and 4/2.

The missing-observation extension is also deterministic. With configured pairs (a1=missing, b1=canonical) and (a2=modified, b2=missing), no configured pair is complete, but modkit emits the a2-versus-b1 result. Its current output is byte-identical to an aligned incomplete-pair control where pair 1 is complete and pair 2 is missing. That equality demonstrates the cross-pairing mechanism; under the proposed all-or-none policy, neither incomplete site should emit replicate vectors, and both replicate fields should be -.

The two current outputs are byte-identical, SHA-256 d5dbd66cce750a62b28480322e7e7938abd3156aaaa2480003cb1105542a55c3.

The current fabricated cross-pair fields are:

replicate_map_pvalues = 0.0000006230948043897833
replicate_effect_sizes = 1

Expected behavior

• Preserve each sample's configured identity through per-position aggregation.
• Pair the ith -a input only with the ith -b input and serialize complete-pair statistics in command-line order.
• A configured pair is complete only when both members have positive aggregated coverage after all input filters. An absent, zero-coverage, or min_valid_coverage-filtered member makes that pair unavailable.
• Emit replicate statistics only when every configured pair is complete. If any configured pair is unavailable, serialize both entire fields, replicate_map_pvalues and replicate_effect_sizes, as the existing - sentinel.
• Do not emit a compact list of complete pairs, introduce an inline missing marker, or substitute another available sample.
• Preserve the group-level row whenever pooled and balanced calculations remain valid; those calculations use all independently available positive-coverage observations and do not require complete pairs.

Root-cause evidence

• The input data structure is initially keyed by sample ID, and sample IDs are assigned from argument position at tabix.rs:18-30 and tabix.rs:177-191.
• Parallel input collection retains those IDs in an FxHashMap at tabix.rs:54-100.
• organize_bedmethy_lines consumes sample.into_values(), discarding each sample ID, and pushes the resulting counts into a per-position Vec<AggregatedCounts>. Conditions A and B are organized independently.
• SingleSiteDmrScore::new_multi considers equal-length vectors matched and calculates replicate statistics with counts_a.iter().zip(counts_b). Vector order is therefore being treated as pair identity after that identity has been discarded.
• The public documentation explicitly says these statistics are calculated according to command-line order at intro_dmr.md:241-265.

Proposed fix scope

• Retain the integer sample ID together with each AggregatedCounts value through per-position organization.
• Construct each configured pair from its original relative A/B argument index, rather than vector or hash-table iteration order.
• Serialize replicate statistics in command-line order.
• Retain fixed configured sample slots at each position. If every A/B slot pair contains positive-coverage observations, calculate all replicate statistics in configured order; otherwise leave both replicate vectors empty so the existing serializer emits - for both whole fields.
• Add complete three-pair, aligned-incomplete, crossed-missing, zero-coverage-member, and coverage-filtered parent-red regressions. Use independent single-pair commands as the complete-pair numerical oracle.

Non-goals

• No change to MAP estimation, balanced weighting, pooled aggregation, coverage capping, or output field definitions.
• No change to unequal-group behavior, where replicate-pair fields are not emitted.
• No change to region-level dmr pair or dmr multi.
• Canonical imputation of zero-coverage rows in balanced group statistics is a separate defect. For this identity repair, however, a zero-coverage pair member is unavailable and therefore triggers the all-or-none - policy.

Acceptance criteria

• The complete three-pair fixture emits effects exactly [1, -0.5, -0.5] and MAP-based p-values exactly matching the three independent commands, in argument order.
• Changing hash-map insertion, traversal, worker, compute-thread, or I/O-thread order cannot change sample pairing.
• If any configured pair is incomplete, including the aligned-missing and crossed-missing fixtures, both replicate fields are exactly -; no partial compact vector is emitted.
• A site containing only a2 and b1 never emits an a2-versus-b1 statistic.
• A zero-coverage or coverage-threshold-filtered pair member is treated as unavailable for replicate completeness.
• The group-level row is retained when each condition still has positive observations. Its pooled, balanced, and representation fields are computed from the independently available positive-coverage observations and are not suppressed merely because a configured pair is incomplete.
• Existing pooled and balanced output fields remain byte-identical for complete positive-coverage inputs.
• 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
modified.bed.gz	91 bytes	d0cda187d5c18ce4041008863ed4cfe8884950a7d45d3b3d4bcbef3b3e788c34	10/10 modified C row
modified.bed.gz.tbi	109 bytes	f57b938f9b5433ccd02b21bf3735405177f5af6446ceec0a3f386b11671bc6bd	Tabix index
half.bed.gz	92 bytes	7118f965db6d273c4d02ba8f00d0634af4f91672593ebf0dafa7e52d193de979	5/10 modified C row
half.bed.gz.tbi	109 bytes	47f8313737fe6cdc0178b117685280a960698552e7843cde53bde591995230dd	Tabix index
canonical.bed.gz	90 bytes	6a2a2e515d35f1d130e260165ac9fb9b21b9c892330ddce1e7994aa50f14fb5d	0/10 modified C row
canonical.bed.gz.tbi	109 bytes	894ac2ee3cb92e5e0d79d709acac8ff19e2a3c2c9d5a857d128e424b41db3a6d	Tabix index
three-pairs.bed	495 bytes	1431de11b030a61d79e7a737b4fa41396833c4be6731bc2944b3dd890bf36e61	Erroneously cross-paired complete output
pair1.bed	339 bytes	2f8ec43cfc80eb0bbb1ded99b0eaf9934ebafc86c32ae4f45e30813fdd59c345	Independent pair 1 control
pair2.bed	342 bytes	0f56474b26eb45da75f6575d5ab15378929fccff9ed2c7f58c6a611f3bf3d3d8	Independent pair 2 control
pair3.bed	343 bytes	64fd87eb880d92d9816468a7424bb118c4916b54e16c5552893f338b014f821a	Independent pair 3 control
missing.bed.gz	28 bytes	d079906378251d29409f7f0f691113fd3b0049b5927b713e2704ca44ac743d48	Empty indexed missing-observation extension
missing.bed.gz.tbi	75 bytes	1a5845ff0c599ef381032688bd8be8b755436ce4721e2ef28fef4690a700705d	Empty Tabix index
cross-pair-missing.bed	513 bytes	d5dbd66cce750a62b28480322e7e7938abd3156aaaa2480003cb1105542a55c3	Fabricated cross-pair replicate result
aligned-pair-control.bed	513 bytes	d5dbd66cce750a62b28480322e7e7938abd3156aaaa2480003cb1105542a55c3	Current aligned incomplete-pair control; byte-identical to the fabricated cross-pair result, although both replicate fields should be - under the proposed policy

Related work

• Duplicate searches used dmr replicate pairing order, replicate_map_pvalues, matched samples DMR, per-replicate p-values, and dmr sample order. No matching upstream report was found.
• Issue #361 discusses general use and interpretation of balanced DMR outputs but does not report loss of command-order pair identity.
• Issue #279 uses repeated DMR inputs but concerns filtered bedMethyl validity rather than replicate pairing.