Skip to content

summary output can be nondeterministic and omit filtered-only categories #672

Description

@SuhasSrinivasan

summary output can be nondeterministic and omit filtered-only categories

Summary

modkit summary writes category rows by iterating unordered maps and derives displayed states primarily from passing calls. Equivalent runs can therefore emit TSV/table rows in different orders, while modification states represented only among filtered calls can be absent or printed with zero counts even though they contribute to the reported totals.

The same command logs --sampling-frac 0.1 as 0.1%, although the actual requested fraction is 10%.

Severity

Severity: High — scientific reporting and reproducibility

Rationale: The filtered-only defect can make category rows fail to reconcile with the command's own pass/fail/all totals and can omit biologically relevant canonical or modified states from downstream parsing. The nondeterministic row order and incorrect sampling log are Medium reporting/reproducibility defects on their own. Underlying sampling decisions and aggregate call totals are not changed by this fix.

User and scientific impact

  • Affected result or workflow: summary TSV/table output, category-based QC or downstream parsing, byte-level reproducibility, and sampling provenance logs.
  • Direction of error: omitted or falsely zero filtered categories; non-reconciling category totals; row-order nondeterminism; a 100-fold logging-unit error.
  • Likely exposure: filtered-only states are threshold/data dependent; nondeterministic ordering and the sampling log affect routine invocations of the corresponding output/options.
  • Detectability or workaround: users can compare category sums with total fields and parse rows by key rather than position. The missing canonical count in the CLI example is derivable by subtraction, but a TSV-omitted filtered-only modified state cannot in general be reconstructed reliably from category totals alone.

Affected versions and environment

  • Released version: reproduced with modkit 0.6.4.
  • Development revision: 5cecc3fb3a9336068d9e3c68d5c08d678153dd2c.
  • Operating system and architecture: macOS 26.6, Apple Silicon.
  • Input format: modBAM.
  • Rust toolchain used for source-level regressions: rustc/cargo 1.90.0.

Steps to reproduce

From a modkit checkout containing the repository test resources:

Filtered-category reconciliation

modkit summary \
  --ignore-index \
  --tsv \
  --filter-threshold 1 \
  tests/resources/bc_anchored_10_reads.sorted.bam \
  > summary.tsv

grep '^C_.*fail' summary.tsv

Observed behavior

C_fail_calls_modified_m	53
C_fail_calls_modified_h	44
C_total_fail_mod_calls	107

The emitted category rows account for only 97 failed calls, while the total reports 107. The ten filtered canonical-C calls have no corresponding C_fail_calls_unmodified row. A separate exact source-level regression with an empty pass map and filtered counts {canonical: 1, m: 3, ChEbi(123): 2} shows that the affected TSV can omit the entire represented base block, while the table prints the filtered-only states as zero; the corrected output reports all three categories and a failed total of six.

Nondeterministic row order

for run_id in 1 2 3 4 5 6 7 8 9 10 11 12; do
  modkit summary \
    --ignore-index \
    --tsv \
    --filter-threshold 1 \
    tests/resources/bc_anchored_10_reads.sorted.bam \
    2>/dev/null | shasum -a 256
done

In one twelve-process reproduction, equivalent results produced two stdout hashes:

8bcdff0b50590a7efb3cfed0629b81a8c45d9d90571f487f1d983e8e8b77096e
f28c75cfbc14a687f22929b6e500302c271e03a821b4a101e191e9b029de47ec

The key/value content is equivalent; modification-category row order differs.

Sampling log units

modkit summary \
  --ignore-index \
  --sampling-frac 0.1 \
  --seed 42 \
  --filter-threshold 1 \
  tests/resources/bc_anchored_10_reads.sorted.bam \
  >/dev/null

Observed log:

sampling 0.1% of reads

Control or independent oracle

  • 53 + 44 = 97, not the reported failed total of 107; an additional canonical filtered category of 10 is required for reconciliation.
  • A fraction of 0.1 is 10%, not 0.1%.
  • Equivalent summary objects constructed with opposite map/set insertion orders must produce identical output bytes.

Expected behavior

  • For every represented canonical base, the output contains its canonical state plus the union of observed, passing, and filtered modification codes.
  • Canonical state is printed first, followed by letter modification codes lexicographically and numeric ChEBI identifiers numerically.
  • Filtered-only rows have pass count/fraction zero and their exact fail/all count and fraction.
  • Category sums reconcile exactly with the corresponding pass/fail/all totals; zero denominators produce finite zero fractions rather than NaN.
  • Equivalent runs and insertion orders produce byte-identical TSV and table output.
  • --sampling-frac 0.1 is logged as 10%, without changing the sampled population.

Root-cause evidence

  • modkit-core/src/writers.rs: Summary output iterates nested HashMap/HashSet values directly and derives displayed rows from passing counts plus observed code metadata rather than the union of pass and filtered states.
  • modkit-core/src/modbam_util/subcommands.rs: the Summary sampling log interpolates the raw fraction while appending a percent sign.
  • Parent-revision regressions reproduce non-reconciling filtered counts, opposite-insertion output differences, and the incorrect log unit.

Proposed fix scope

  • Define one writer-local output order: A/C/G/T, canonical first, then letter codes, then numeric ChEBI identifiers.
  • Emit the union of passing states, filtered states, and observed modification codes with zero-safe fractions.
  • Scale the sampling fraction by 100 only in the diagnostic message.
  • Keep deterministic-output and filtered-only reconciliation as separate commits within one focused Summary PR.

Non-goals

  • No change to sampling decisions, thresholds, eligibility, mapped/unmapped handling, modification classification, or aggregate counts.
  • No global chemical/ontology ranking for modification codes; the order is a presentation contract.
  • No inclusion of separate indexed sampling, unmapped-classification, coordinate-selection, or read-metadata issues.
  • No output-system or performance refactor.

Acceptance criteria

  • The filtered-only reproducer emits the missing canonical category and all category sums reconcile with the reported totals.
  • Opposite insertion orders, repeated processes, and supported thread counts produce byte-identical TSV and table output.
  • Canonical, letter-code, and numeric-ChEBI rows follow the documented deterministic order.
  • Fractions remain finite when the pass population is empty.
  • A requested sampling fraction of 0.1 logs as 10%, while selected reads and scientific counts remain unchanged.
  • Focused writer/log/integration tests and the applicable full workspace test suite pass.

Reproduction artifacts

Artifact Size SHA-256 Notes
tests/resources/bc_anchored_10_reads.sorted.bam 5,527 bytes 4441acbe1ad59caf6ab5d56bf18097d8999c08351395219c919cb863652d85f0 Existing modkit test fixture used for CLI reproductions.

Related work

  • Related issues: none found.
  • Proposed PR: to be linked once the completed branch has been rebased and regated.
  • Deferred policy question or upstream dependency: the separate ModCodeRepr ordering-law issue preserves the same presentation order but changes a global Rust trait contract; dependency-resolution PR Pin hts-sys to bindings compatible with rust-htslib 0.46 #666 affects publication timing only.

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