localize --min-coverage does not filter bedMethyl records

[SuhasSrinivasan]

Summary

modkit localize --min-coverage is parsed and logged but is not applied to bedMethyl records. Changing the threshold therefore does not change the localized aggregate.

Severity

Severity: High — scientific correctness

Rationale: The command accepts the option and exits successfully but silently includes observations the user explicitly requested to exclude, changing localized coverage, modification counts, and percent-modified profiles.

User and scientific impact

• Affected workflow: modkit localize profiles around supplied regions.
• Direction of error: Low-coverage records inflate n_valid and can increase or decrease the aggregate modification percentage.
• Exposure: Routine whenever --min-coverage is used and the input contains records below the threshold; sparse direct-RNA profiles can be especially sensitive.
• Detectability: The option is logged without a warning that it has no effect.

Affected versions and environment

• modkit 0.6.4.
• Assessed source revision 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c.
• macOS on Apple silicon.

Steps to reproduce

Minimal input

From a modkit checkout, create two one-base focus regions and a genome-sizes file:

mkdir -p /tmp/modkit-localize-coverage
printf 'chr20\t9681998\t9681999\nchr20\t9838537\t9838538\n' \
  > /tmp/modkit-localize-coverage/regions.bed
printf 'chr20\t64444167\n' \
  > /tmp/modkit-localize-coverage/genome-sizes.tsv

The corresponding records in the existing indexed bedMethyl resource have (valid coverage, modified count) values (1,1) and (23,2):

tabix tests/resources/lung_00733-m_adjacent-normal_5mc-5hmc_chr20_cpg_pileup.bed.gz chr20:9681999-9681999
tabix tests/resources/lung_00733-m_adjacent-normal_5mc-5hmc_chr20_cpg_pileup.bed.gz chr20:9838538-9838538

Commands

modkit localize \
  tests/resources/lung_00733-m_adjacent-normal_5mc-5hmc_chr20_cpg_pileup.bed.gz \
  --regions /tmp/modkit-localize-coverage/regions.bed \
  --genome-sizes /tmp/modkit-localize-coverage/genome-sizes.tsv \
  --window 1 \
  --min-coverage 1 \
  --threads 1 \
  --io-threads 1 \
  > /tmp/modkit-localize-coverage/min1.tsv
printf 'min1_exit=%s\n' "$?"

modkit localize \
  tests/resources/lung_00733-m_adjacent-normal_5mc-5hmc_chr20_cpg_pileup.bed.gz \
  --regions /tmp/modkit-localize-coverage/regions.bed \
  --genome-sizes /tmp/modkit-localize-coverage/genome-sizes.tsv \
  --window 1 \
  --min-coverage 3 \
  --threads 1 \
  --io-threads 1 \
  > /tmp/modkit-localize-coverage/min3.tsv
printf 'min3_exit=%s\n' "$?"

Control or independent oracle

The option defines an inclusive per-record predicate applied before aggregation:

keep record iff valid_coverage >= min_coverage

At threshold 1, both records contribute: n_valid=24, n_mod=3, percent_modified=12.5. At threshold 3, only the (23,2) record contributes.

Observed behavior

Both commands exit 0 and create output, but both thresholds emit the same aggregate:

mod_code	offset	n_valid	n_mod	percent_modified
C	-1	24	3	12.5

Expected behavior

Both commands should exit 0. Threshold 1 should retain the 24/3 aggregate, while threshold 3 should exclude the coverage-1 record and emit:

mod_code	offset	n_valid	n_mod	percent_modified
C	-1	23	2	8.695652

Filtering after aggregation would be incorrect because multiple individually under-covered loci could combine and pass the threshold.

Root-cause evidence

• EntryLocalize::run assigns and logs self.min_coverage but does not pass it into GenomeRegion::into_localized_mod_counts.
• That function consequently folds every fetched bedMethyl record into the relative-offset aggregate.
• Passing the threshold into the function and filtering each record before the fold produces the exact oracle above.

Proposed fix scope

Pass the threshold into localized counting and retain a bedMethyl record only when valid_coverage >= min_coverage, before aggregation. Add an end-to-end regression covering below, equal-to, and above-threshold records.

Non-goals

• Does not change localize anchor, window, or offset conventions.
• Does not define or implement --batch-size behavior.
• Does not change chart rendering or broader worker-error handling.

Acceptance criteria

• The focused regression is red on the assessed parent and green on the proposed fix.
• Below/equal/above-threshold cases establish inclusive per-record filtering.
• Threshold 1 emits 24/3/12.5%; threshold 3 emits 23/2/8.695652%.
• Focused, affected integration, and serial workspace tests pass.

Related work

• Proposed fix: Apply localize minimum coverage before aggregation #652