Skip to content

summary --cpg includes unmapped modification calls #649

Description

@SuhasSrinivasan

Summary

For an indexed BAM, modkit summary --cpg includes modification calls from unmapped records, while the documented equivalent --motif CG 0 excludes them. This changes the summary denominator and modification counts even though an unmapped record cannot be assigned to a reference CpG.

Severity

Severity: High — scientific correctness

Rationale: A routine documented shorthand can silently include a biologically ineligible population and alter read counts and modification statistics while the command exits successfully.

User and scientific impact

  • Affected workflow: Indexed modkit summary --cpg.
  • Direction of error: Read and call totals can be inflated by unmapped records.
  • Exposure: Proportional to the unmapped-read fraction and the MM/ML calls carried by those records.
  • Detectability: The output contains no warning; comparison with explicit --motif CG 0 or a mapped-only control is required.

Affected versions and environment

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

Steps to reproduce

Minimal input

From a modkit checkout, combine the existing mapped and unmapped test resources and create an index:

mkdir -p /tmp/modkit-summary-cpg
samtools merge -f \
  -o /tmp/modkit-summary-cpg/mixed.bam \
  tests/resources/bc_anchored_10_reads.sorted.bam \
  tests/resources/bc_anchored_10_reads.unmapped.bam
samtools index /tmp/modkit-summary-cpg/mixed.bam

Commands

modkit summary /tmp/modkit-summary-cpg/mixed.bam \
  --reference tests/resources/CGI_ladder_3.6kb_ref.fa \
  --cpg \
  --filter-threshold 0 \
  --tsv \
  --suppress-progress \
  --threads 1 \
  --io-threads 1 \
  > /tmp/modkit-summary-cpg/cpg.tsv
printf 'cpg_exit=%s\n' "$?"

modkit summary /tmp/modkit-summary-cpg/mixed.bam \
  --reference tests/resources/CGI_ladder_3.6kb_ref.fa \
  --motif CG 0 \
  --filter-threshold 0 \
  --tsv \
  --suppress-progress \
  --threads 1 \
  --io-threads 1 \
  > /tmp/modkit-summary-cpg/motif.tsv
printf 'motif_exit=%s\n' "$?"

Control or independent oracle

--cpg is documented as shorthand for --motif CG 0. The explicit motif invocation and the mapped-only input are independent controls: all three normalized summaries must agree because reference-motif membership is undefined for unmapped records.

Observed behavior

Both commands exit 0 and create output, but --cpg includes all 20 records:

count_reads_C	20
C_total_mod_calls	186
total_reads_used	20

The explicit --motif CG 0 control excludes the 10 unmapped records:

count_reads_C	10
C_total_mod_calls	77
total_reads_used	10

Expected behavior

Both commands should exit 0 and produce identical normalized summary values. The expected values on this fixture are total_reads_used=10, count_reads_C=10, and C_total_mod_calls=77; no unmapped call should change any canonical or modified category.

Root-cause evidence

  • The indexed Summary path constructs the same CG motif for --cpg.
  • The later guard that prevents unmapped collection checks mapped_only, matched_only, and an explicit motif, but omits the cpg flag.
  • Consequently, the path fetches FetchDefinition::Unmapped and merges unmapped calls into the reference-CpG histogram.
  • Including the CpG shorthand in that guard makes the two commands agree.

Proposed fix scope

Treat --cpg as a resolved motif selection when deciding whether indexed Summary may collect unmapped records, and add an end-to-end regression with a valid unmapped MM/ML/MN record.

Non-goals

  • Does not change mapped motif matching.
  • Does not change sampling or automatic-threshold calculation.
  • Does not redefine modification classes or Summary output schema.

Acceptance criteria

  • The regression is red on the assessed parent and green on the proposed fix.
  • --cpg, --motif CG 0, and the mapped-only control have identical normalized scientific counts.
  • The unmapped record changes neither read totals nor canonical/modified categories.
  • Focused, affected integration, and serial workspace tests pass.

Related work

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions