bedmethyl map-to-genome emits invalid intervals and transcript-relative strands

[SuhasSrinivasan]

Summary

modkit bedmethyl map-to-genome uses the input bedMethyl strand to choose the genomic interval endpoint. A negative-strand row therefore becomes [genome_start, genome_start - 1), producing a reversed or zero-length BED interval. The command also retains transcript-relative row strands when the selected transcript is on the genomic negative strand instead of composing them into genomic orientation.

Both behaviors reproduce on accepted indexed bedMethyl input with modkit 0.6.4 and current upstream revision 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c.

Severity

Severity: High — scientific correctness and output validity

Rationale: The command exits successfully while emitting invalid coordinates and incorrect genomic strands. Downstream interval tools may reject, drop, or misinterpret the affected rows, and a retained transcript-relative strand can invert biological strand interpretation.

User and scientific impact

• Affected result or workflow: transcriptome-to-genome conversion of bedMethyl rows, especially negative-strand calls and negative-strand transcripts.
• Direction of error: reversed/zero-length intervals and wrong genomic strand labels.
• Likely exposure: routine for transcriptome-aligned direct-RNA data containing negative-strand rows or transcripts.
• Detectability or workaround: independently remap every transcript position and compose transcript/call strand externally; validate that every output interval satisfies end = start + 1.

Affected versions and environment

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

Steps to reproduce

These fixtures are synthetic and contain no private data.

Create models.gtf:

chrP	test	exon	1	2	.	+	.	gene_id "g_pos"; transcript_id "tx_pos";
chrP	test	exon	11	12	.	+	.	gene_id "g_pos"; transcript_id "tx_pos";
chrN	test	exon	1	2	.	-	.	gene_id "g_neg"; transcript_id "tx_neg";
chrN	test	exon	11	12	.	-	.	gene_id "g_neg"; transcript_id "tx_neg";

Create tx_pos.bed and tx_neg.bed:

# tx_pos.bed
tx_pos	1	2	m	1	-	1	2	255,0,0	1	100.00	1	0	0	0	0	0	0

# tx_neg.bed
tx_neg	0	1	m	1	+	0	1	255,0,0	1	100.00	1	0	0	0	0	0	0
tx_neg	1	2	m	1	-	1	2	255,0,0	1	100.00	1	0	0	0	0	0	0

Compress, index, and map both transcripts:

cp tx_pos.bed tx_pos.input.bed
cp tx_neg.bed tx_neg.input.bed
bgzip -f tx_pos.input.bed
bgzip -f tx_neg.input.bed
tabix -f -p bed tx_pos.input.bed.gz
tabix -f -p bed tx_neg.input.bed.gz

modkit bedmethyl map-to-genome \
  tx_pos.input.bed.gz mapped-pos.bed \
  --transcript-id tx_pos --ignore-version --gtf models.gtf

modkit bedmethyl map-to-genome \
  tx_neg.input.bed.gz mapped-neg.bed \
  --transcript-id tx_neg --ignore-version --gtf models.gtf

cat mapped-pos.bed mapped-neg.bed

Observed behavior

modkit 0.6.4 exits 0 and reports all three rows processed, but emits:

chrP	1	0	m	1	-	1	0	255,0,0	1	100.00	1	0	0	0	0	0	0
chrN	11	12	m	1	+	11	12	255,0,0	1	100.00	1	0	0	0	0	0	0
chrN	10	9	m	1	-	10	9	255,0,0	1	100.00	1	0	0	0	0	0	0

The first and third rows have end < start; both negative-transcript rows retain transcript-relative rather than genomic strand.

Expected behavior

Every mapped base is the valid half-open single-base interval [genome_start, genome_start + 1). Row strand is composed with transcript strand: positive transcripts preserve +/-, negative transcripts flip + and -, and unstranded . remains ..

The exact expected output is:

chrP	1	2	m	1	-	1	2	255,0,0	1	100.00	1	0	0	0	0	0	0
chrN	11	12	m	1	-	11	12	255,0,0	1	100.00	1	0	0	0	0	0	0
chrN	10	11	m	1	+	10	11	255,0,0	1	100.00	1	0	0	0	0	0	0

Root-cause evidence

• genome_stop is selected from the input row strand, subtracting one for negative rows, at subcommands.rs:764-779.
• The parsed transcript model already retains genomic strand at isoform/mod.rs:146-155, but the map-to-genome path does not use it when transforming the bedMethyl strand.

Proposed fix scope

• Map the transcript coordinate once, then always emit [genome_start, genome_start + 1) with checked addition.
• Compose the bedMethyl row strand with the transcript model strand in one small helper.
• Expose a crate-private transcript-strand accessor rather than widening the model fields.
• Add a compact two-exon matrix covering both transcript orientations, both row strands, unstranded rows, coordinate zero, and exon boundaries.

Non-goals

• No change to transcript-ID prefix matching, GTF model construction, or version-selection policy.
• No genomic resorting of negative-transcript output; rows retain existing input/transcript order.
• No change to the existing policy for malformed bedMethyl rows or non-unit input intervals.
• No merge, DMR, sampling, or performance changes.

Acceptance criteria

• The minimal reproducer emits exactly the expected three rows.
• Every mapped interval satisfies end = start + 1, including genomic coordinate zero and exon boundaries.
• Positive-transcript strands are preserved; negative-transcript +/- are flipped; . remains ..
• Existing successful behavior outside the stated scope remains unchanged.
• Focused, affected-package, and full workspace tests pass.

Reproduction artifacts

Artifact	Size	SHA-256	Notes
models.gtf	268 B	def6fd8a5b816e5cc85c28419d36d9b586aa7d22e35df9ddb806216d0d423a86	Positive and negative two-exon transcripts
tx_pos.bed	52 B	e5b15a00c2df80f8e203db32ca7f84be57f49f72d2bbbb82baa9009004010c4c	Negative row on a positive transcript
tx_neg.bed	104 B	924c1d30fd8b1caa127c57f539151a4564551ef09f5b833d66ecfdb187c3ca4b	Positive and negative rows on a negative transcript

Related work

• Related issue: none.
• Proposed PR: Map bedMethyl intervals and strands into genomic space #662
• Deferred policy or upstream dependency: transcript-contig prefix ambiguity and genomic output sorting remain separate questions.

[SuhasSrinivasan]

The complete code change is in ready-for-review PR #662. Closing this standalone catalog issue in line with the maintainer's request that these submissions be represented by PRs rather than separate AI-generated issues.