extract can omit documented motif metadata and overreport emitted rows

[SuhasSrinivasan]

Summary

modkit extract full and extract calls omit or misreport motif metadata in three accepted cases. A valid motif at 0-based reference coordinate zero is rendered as ., --cpg omits the documented motifs column and cannot be combined with its documented --annotate-motifs or --mask options, and extract calls --pass-only counts filtered-out candidates as rows written.

The first two behaviors make --cpg differ from its documented --motif CG 0 shorthand and can affect downstream motif-stratified analysis or schema parsing. The row-accounting defect changes only progress/log metadata; the emitted call rows themselves are already correct.

Severity

Severity: Medium — scientific metadata, output schema, and reproducibility

Rationale: Reference coordinate zero is a valid aligned position and can be common across many short transcriptome contigs, so labeling its motif as absent is scientifically wrong. The --cpg schema also contradicts the documented 22-column full and 24-column calls contracts and rejects documented option combinations. The accompanying row-counter defect is Low severity by itself because it does not alter emitted rows.

User and scientific impact

• Affected result or workflow: extract full and extract calls with --motif, --cpg, --annotate-motifs, --mask, or --pass-only.
• Direction of error: missing/misclassified motif labels, a missing final schema column, parser rejection of valid documented combinations, and inflated “rows” run metadata.
• Likely exposure: option-specific but routine for CpG extraction; coordinate-zero exposure is especially relevant to transcriptome references with many short contigs.
• Detectability or workaround: column counts and parser errors are visible. --motif CG 0 restores the motif column for most sites but still mislabels coordinate zero. The progress count can be checked against data-row counts. Recomputing motif overlap downstream is the only complete released-version workaround.

Affected versions and environment

• Released version: modkit 0.6.4
• Development revision: 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c
• Operating system and architecture: macOS 26.6, arm64
• Input: coordinate-sorted BAM with BAI and indexed FASTA
• Related tool: samtools 1.23.1 for the synthetic fixture

Steps to reproduce

Minimal input

Create ref.fa:

>chr1
CGCA

Create input.sam:

@HD	VN:1.6	SO:coordinate
@SQ	SN:chr1	LN:4
r1	0	chr1	1	60	4M	*	0	0	CGCA	IIII	MM:Z:C+m?,0,0;	ML:B:C,255,128	MN:i:4

The two explicit 5mC probabilities are at 0-based reference positions 0 and 2. Only position 0 overlaps CG,0.

Build the fixture:

samtools faidx ref.fa
samtools view -b -o input.bam input.sam
samtools index input.bam

Commands

First, annotate all explicit-motif output and inspect the final column:

modkit extract full input.bam explicit-full.tsv \
  --reference ref.fa --motif CG 0 --annotate-motifs \
  --threads 1 --io-threads 1 --suppress-progress --force

modkit extract calls input.bam explicit-calls.tsv \
  --reference ref.fa --motif CG 0 --annotate-motifs \
  --filter-threshold 0.9 --pass-only \
  --threads 1 --io-threads 1 --suppress-progress --force

awk -F '\t' 'FNR <= 2 { print $3, $NF }' \
  explicit-full.tsv explicit-calls.tsv

Next, use the documented CpG shorthand and count fields:

modkit extract full input.bam cpg-full.tsv \
  --reference ref.fa --cpg \
  --threads 1 --io-threads 1 --suppress-progress --force

modkit extract calls input.bam cpg-calls.tsv \
  --reference ref.fa --cpg --no-filtering \
  --threads 1 --io-threads 1 --suppress-progress --force

awk -F '\t' 'FNR <= 2 { print FILENAME, NF, $NF }' \
  cpg-full.tsv cpg-calls.tsv

Then exercise the documented dependent options:

modkit extract full input.bam cpg-annotated.tsv \
  --reference ref.fa --cpg --annotate-motifs --force

modkit extract full input.bam cpg-masked.tsv \
  --reference ref.fa --cpg --mask --force

Finally, capture the pass-only run summary:

modkit extract calls input.bam pass-only.tsv \
  --reference ref.fa --filter-threshold 0.9 --pass-only \
  --threads 1 --io-threads 1 --suppress-progress \
  --log-filepath pass-only.log --force

grep 'processed .* rows' pass-only.log
awk 'END { print NR - 1 }' pass-only.tsv

Control or independent oracle

• The extract schema defines ref_position as 0-based, so position 0 is valid and overlaps the first CG,0 occurrence in CGCA.
• The CLI describes --cpg as shorthand for --motif CG 0, and the published schemas state that motifs is present when either --motif or --cpg is used.
• With threshold 0.9 and --pass-only, ML 255 passes and ML 128 fails. Exactly one data row is emitted, so the run summary must report one row.

Removing only the repaired final motif field from explicit-motif output gives byte-identical installed/fixed rows. An independent rerun of the pass-only commands above retained and compared both 371-byte TSVs; they are byte-identical before and after the accounting repair, SHA-256 47b6836e78d95324f4a74a204f69bf2bb5a3a3ba6ca549291f94bae636cd913e.

Observed behavior

On modkit 0.6.4 and current upstream:

Explicit motif, ref_position 0: motifs = .       (expected CG,0)
--cpg full:                    21 columns        (expected 22)
--cpg calls:                   23 columns        (expected 24)
--cpg --annotate-motifs:       exit 2, demands --motif
--cpg --mask:                  exit 2, demands --motif
Pass-only output data rows:    1
Pass-only log:                 processed 1 reads, 2 rows, ...

The same coordinate-zero omission occurs in both the 22-column full renderer and the 24-column calls renderer when an explicit motif enables the field.

Expected behavior

• Reference coordinate zero is queried like every other non-negative aligned reference coordinate and is labeled CG,0 in this fixture.
• --cpg and --motif CG 0 produce the same motif-selection and annotation semantics.
• --cpg output has the documented final motifs column: 22 columns for full and 24 for calls.
• --annotate-motifs and --mask accept either --motif or --cpg, while still requiring at least one motif source.
• The “rows” progress/log count increases only after a row is actually written.
• Call probabilities, call codes, fail values, filtering membership, row ordering, and base schemas without a motif source remain unchanged.

Root-cause evidence

• Both renderers exclude reference coordinate zero with i > 0 at extract/writer.rs:68-81 and read_ids_to_base_mod_probs.rs:499-510.
• Both command paths derive the output schema from motif.is_some() and omit cpg at extract/subcommand.rs:279-309 and extract/subcommand.rs:719-754.
• --annotate-motifs and --mask require only the raw motif argument even though their help and the cpg help define the shorthand at extract/args.rs:99-132.
• The calls writer increments rows_written after every candidate even when to_row returns None under --pass-only at extract/writer.rs:280-308.
• The published schema requires the final field for either motif source at intro_extract.md:30-40 and intro_extract.md:49-74.

Proposed fix scope

• Query motif annotations at every non-negative reference position in both row renderers.
• Derive one motif.is_some() || cpg schema predicate for both extract subcommands.
• Make --annotate-motifs and --mask depend on an argument group satisfied by either --motif or --cpg.
• Increment rows_written only when to_row returns a row and that row is successfully written.
• Add exact coordinate-zero, CpG schema, parser, masking, and pass-only accounting regressions.

Non-goals

• No modification-calling, threshold-estimation, probability, fail, pass-only membership, or output-order change.
• No motif-discovery, motif-coordinate, reference mapping, or strand-semantics redesign.
• No change to the 21-column full or 23-column calls base schemas when no motif source is requested.
• No BGZF writer/finalization or performance change.
• No change to the existing empty-target error when masking removes every motif from a reference.

Acceptance criteria

• The synthetic coordinate-zero row ends in CG,0 for both full and calls.
• --cpg produces 22-column full and 24-column calls output with a final motifs header and CG,0 value.
• --cpg and --motif CG 0 are byte-identical across full/calls, restriction/annotation, and masked/unmasked controls.
• --annotate-motifs and --mask accept either motif source and reject use without either source.
• The pass-only fixture emits one data row, reports one row, and remains byte-identical apart from the corrected log metadata.
• Exact tests cover coordinate zero, nonzero matched/unmatched sites, soft masking, both schemas, and no-motif base schemas.
• Focused tests and the applicable full workspace suite pass.

Reproduction artifacts

Artifact	Size	SHA-256	Notes
ref.fa	11 B	273213eaf04f1e15a4ffdf3079e7eb331fa351a62caad2e58fbec6ef02f9faca	Four-base synthetic reference
ref.fa.fai	13 B	3c5c82b3f0d5e9f0550836f5c3eb5fb9928db3b10c8b426c83461b358419cb10	samtools FASTA index
input.sam	113 B	a9e86621433888a76108e662ee00c0c1554be27aa05abcb81bdcf1ff1b4b47d7	One read with two C+m calls
input.bam	287 B	d5a476b721fb7ce9449fa6339a3195f20501a3285c7dcdce33e6b888a433f533	Coordinate-sorted fixture data
input.bam.bai	96 B	82ca96e37b2227a2e0da2a1f67884620e4466b9ee18a6bf8163ce04cfca23434	BAM index
pass-only.tsv	371 B	47b6836e78d95324f4a74a204f69bf2bb5a3a3ba6ca549291f94bae636cd913e	Installed and repaired outputs are byte-identical; only logged row accounting differs

Related work

• Fresh dependency resolution selects incompatible hts-sys bindings #665 / Pin hts-sys to bindings compatible with rust-htslib 0.46 #666 address the separate fresh-build dependency incompatibility and do not change extract metadata or accounting.
• Proposed PR: to be linked only after the implementation has been rebased onto the accepted dependency resolution and the complete test gate passes.
• Downstream compatibility note: the repair intentionally appends the already documented motifs column to --cpg output; parsers built around the incorrect released 21/23-column behavior will need to accept the documented 22/24-column schema.