bedmethyl merge can silently skip inputs and combine distinct motif rows

[SuhasSrinivasan]

Summary

modkit bedmethyl merge can silently drop a requested input when its data file or Tabix index cannot be opened, yet continue with exit status 0. In the remaining data, it also reduces the full bedMethyl column-4 name to only the modification code and counts duplicate rows from one file as independent samples.

Together, these behaviors can combine biologically distinct motif rows and allow one surviving input to satisfy --min-samples all after another requested input was skipped. Both reproduce on accepted indexed bedMethyl input with modkit 0.6.4 and current upstream revision 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c.

Severity

Severity: High — scientific correctness, sample membership, and output integrity

Rationale: The command can exit successfully with biologically distinct rows combined, an incorrect sample-support tally, or a requested sample omitted. The output remains syntactically valid, so downstream analysis cannot reliably detect or reconstruct the lost row identity or true contributing-file count.

User and scientific impact

• Affected result or workflow: merging replicate or sample bedMethyl tables, particularly motif-annotated pileups and workflows using --min-samples.
• Direction of error: aggregation across distinct motif identities, false inclusion under a sample-support filter, and silent omission of unreadable inputs.
• Likely exposure: option- and data-dependent, but routine where column 4 contains code,motif,offset labels or where one requested input/index is missing or unreadable.
• Detectability or workaround: preflight every input and index externally, retain and group by the exact column-4 name, and independently count distinct contributing files rather than rows.

Affected versions and environment

• Released version: modkit 0.6.4
• Development revision: 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c
• Operating system and architecture: macOS 26.6, arm64
• Input format and index: bgzip-compressed bedMethyl with Tabix .tbi indexes
• Related tools: bgzip/tabix from HTSlib 1.23.1

Steps to reproduce

These fixtures are synthetic and contain no private data.

Create sizes.tsv:

chr1	100

Create input-a.bed:

chr1	10	11	m,CG,0	4	+	10	11	255,0,0	4	25.00	1	2	1	1	1	1	1
chr1	10	11	m,DRACH,2	5	+	10	11	255,0,0	5	40.00	2	2	1	2	0	1	0
chr1	10	11	m,CG,0	6	+	10	11	255,0,0	6	33.33	2	3	1	2	2	2	2
chr1	10	11	h,CG,0	3	+	10	11	255,0,0	3	33.33	1	2	0	3	3	3	3
chr1	10	11	h,CG,0	4	+	10	11	255,0,0	4	50.00	2	2	0	4	4	4	4

Create input-b.bed:

chr1	10	11	m,CG,0	5	+	10	11	255,0,0	5	80.00	4	1	0	4	4	4	4
chr1	10	11	m,DRACH,2	5	+	10	11	255,0,0	5	60.00	3	2	0	3	3	3	3

Compress, index, and merge:

cp input-a.bed input-a.input.bed
cp input-b.bed input-b.input.bed
bgzip -f input-a.input.bed
bgzip -f input-b.input.bed
tabix -f -p bed input-a.input.bed.gz
tabix -f -p bed input-b.input.bed.gz

modkit bedmethyl merge \
  input-a.input.bed.gz input-b.input.bed.gz \
  --genome-sizes sizes.tsv \
  --out-bed merged.bed \
  --threads 1 --io-threads 1 \
  --min-samples all

Then test a missing requested input while protecting an existing output with a sentinel:

printf 'SENTINEL\n' > protected-output.bed

modkit bedmethyl merge \
  input-a.input.bed.gz absent.input.bed.gz \
  --genome-sizes sizes.tsv \
  --out-bed protected-output.bed --force \
  --threads 1 --io-threads 1 \
  --min-samples all

printf 'exit=%s\n' "$?"
cat protected-output.bed

Control or independent oracle

The bedMethyl name is the complete fourth column, so m,CG,0 and m,DRACH,2 are distinct merge identities. A row is supported by --min-samples all only if each requested input contributes at least one qualifying record for that same identity; duplicate rows within one input may add counts but must not add samples.

Observed behavior

The first command exits 0 and emits:

chr1	10	11	h	7	+	10	11	255,0,0	7	42.86	3	4	0	7	7	7	7
chr1	10	11	m	25	+	10	11	255,0,0	25	48.00	12	10	3	12	10	11	10

The two m motif identities are collapsed into one row. The h row is retained even though both supporting records occur only in input-a.

The missing-input command also exits 0, silently skips absent.input.bed.gz, overwrites the sentinel, and emits rows derived only from input-a. Duplicate rows in that one surviving file satisfy the requested two-sample threshold.

The behavior was deterministic in repeated one-thread runs. The first output had SHA-256 60ee4d996c0d5cb7d8218f77e383f58213567860268168f1245daf1a2888cb42.

Expected behavior

The valid two-input merge should exit 0 and emit distinct motif rows, excluding h,CG,0 because it is absent from input-b:

chr1	10	11	m,CG,0	15	+	10	11	255,0,0	15	46.67	7	6	2	7	7	7	7
chr1	10	11	m,DRACH,2	10	+	10	11	255,0,0	10	50.00	5	4	1	5	3	4	3

If any requested input data file or supported Tabix index cannot be opened, merge should exit nonzero, identify that input, and fail before creating or truncating the output. --min-samples must continue to refer to the requested input cardinality, and duplicate rows from one input must count as one contributing sample for a given identity.

Root-cause evidence

• BedMethylLine parses only the first comma-delimited component of column 4, and its generic serializer consequently writes only that modification code at tabix.rs:25-60.
• The merge key omits the complete column-4 name and increments n_samples once per row rather than once per input at subcommands.rs:192-258.
• The output is opened before input validation, and filter_map silently removes every failed Tabix reader at subcommands.rs:289-326.
• Parent/release regression: the public fixture above exits 0 with collapsed names and a false two-sample h row; the missing-input variant exits 0 and overwrites the sentinel.

Proposed fix scope

• Validate and collect every requested input reader before opening the output; return contextual data/index errors without shrinking the reader set.
• Preserve the exact column-4 name in a merge-specific parsed record and include it in the merge identity.
• Track keys contributed by each reader so duplicate rows still sum their numeric fields but increment sample support only once per requested input.
• Preserve existing ordering and valid merge bytes where the complete name already equals the modification code.

Non-goals

• No full pre-read of each BGZF stream; decode failures that arise only after successful reader/index construction remain a separate late-pipeline issue.
• No writer/channel cancellation or compressed-output finalization redesign.
• No change to the existing CLI contract that the same path requested twice represents two requested samples.
• No performance or public schema change.

Acceptance criteria

• The two-input reproducer emits exactly the two expected motif-specific rows and excludes h,CG,0.
• Missing data, missing index, corrupt index, and mixed valid/invalid --min-samples all inputs return nonzero and identify the invalid path before output truncation.
• Valid .tbi and .csi inputs remain accepted.
• Each requested input contributes at most one sample count per complete merge identity, while all eight numeric count fields from duplicate rows remain summed.
• Bare modification names, motif-bearing names, and numeric ChEBI names are preserved exactly and output is invariant to input/row order.
• Existing valid merge output outside the affected identities remains byte-identical.
• Focused, affected-package, and full workspace tests pass.

Reproduction artifacts

Artifact	Size	SHA-256	Notes
sizes.tsv	9 B	b13960b8ab804bd52e9b12a9bed1e6348c898dc91b1e073ce1cf829741d591e7	One synthetic contig
input-a.input.bed.gz	167 B	3c1038492a0f27a10ec5d88673a79d6a9abf2225651911650fb426d5c1420ebb	Duplicate and motif-specific rows
input-a.input.bed.gz.tbi	108 B	deb8187d4eb4eb6b605b1becf750ea1daeadecc944941a77d8f83d3072f63e1a	Tabix index
input-b.input.bed.gz	123 B	ad59847dcfa5887c84fb4bbf89776aba449e0cbce7e565a203a193cf042a3728	Independent second input
input-b.input.bed.gz.tbi	107 B	cf184a9c3e0d93548d53a525466a0f8ae721109c164b0c85ffad784eeba63b89	Tabix index

Related work

• Related issue: #625 defines the intended distinct-input semantics for --min-samples; #334 reported a successful merge with missing output contig data, but does not establish this reproducer's input-open or row-identity mechanisms.
• Proposed PR: #664.
• Deferred issue: late BGZF/decode and downstream merge-pipeline errors remain separate from this bounded preflight and identity fix.