Skip to content

gen-mut-model: refuse REF mismatches as NEAT2 did, with a bcftools fix in the error - #826

Open
joshfactorial wants to merge 5 commits into
developfrom
fix/gen_mut_model_ref_checks
Open

joshfactorial wants to merge 5 commits into
developfrom
fix/gen_mut_model_ref_checks

Conversation

@joshfactorial

@joshfactorial joshfactorial commented Oct 7, 2026 •

Copy link
Copy Markdown
Collaborator

Behavior decisions from the #819 audit, as chosen by the maintainer. Builds on #823, which edits the same test; merge #823 first, and this PR's diff then shrinks to its own commit.

Versioning: a VCF with any REF mismatch used to build a model and is now refused, and the fitted rate changes for any VCF with edge SNPs. Simulated content is not public API under versioning.md, and this lands in v4.0.0, so no separate bump question arises.

What NEAT2 did (~/code/neat2/utilities/genMutModel.py)

  • SNP REF ≠ reference: hard error, exit(1).
  • SNP whose context has an N: skipped, and counted nowhere; SNP_COUNT += 1 runs only after every check (line 329).
  • Indel REF: never checked.

What eidolon did

  • Mismatched SNPs were warned about and skipped, but counted first: snp_count += 1 came before the edge and REF checks, so a skipped SNP still raised snp_freq and mutation_rate. Measured: one good SNP plus one mismatched gave rate 2/13114.
  • Indel REFs were never checked.

Now (NEAT2's rule, with a better error)

  • Any SNP or indel whose REF disagrees with the reference refuses the fit. NEAT2 exited at the first mismatch; eidolon scans the whole VCF first, so the error gives the SNP and indel counts and the first five positions (H1N1_HA:1000 REF G, reference T).
  • The error says how to proceed deliberately: bcftools norm -f reference.fa --check-ref x in.vcf.gz -Oz -o checked.vcf.gz drops the mismatching records. Checked against bcftools 1.19: it kept the matching record and dropped the mismatched one (skipped: 1).
  • A SNP at a contig edge has no context. It is left out and no longer counted, matching NEAT2, with a warning. An all-edge VCF still errors (NoUsableSnps).
  • A BED that names no reference contig now fails with both contig lists and a chr1 vs 1 hint, instead of Trinuc counts are empty. Unknown error. This only fires when none of the BED's contigs is in the reference. A BED whose regions are all too short keeps the old path, where "check the names" would mislead.
  • docs-site/src/models/mutation-model.md describes all of it.

An earlier revision of this PR dropped mismatches with a warning and refused only above 1%. Review replaced that with this rule: it removes an unmeasured constant, and the user drops records explicitly rather than a threshold deciding silently.

Evidence

  • Known-answer tests, built from the H1N1 reference itself so every REF is exact:
    • 100 good SNPs give rate exactly 100/13114 (must-not-fire);
    • 100 good + 1 SNP with a wrong REF is refused, naming 1 SNP(s) and 0 indel(s), the position with both bases, and bcftools norm, and no model is written;
    • 100 good + 1 indel with a wrong REF is refused, naming H1N1_HA:80 REF ACG, reference ACA;
    • 100 good + 1 SNP at a contig edge builds with rate still 100/13114 (must-not-fire for the mismatch error);
    • two edge-only SNPs give NoUsableSnps { counted: 2, edge: 2 };
    • an unknown BED contig names chrZ_nonexistent and H1N1_HA.
  • The indel-only test's fixture had wrong REFs (H1N1_HA POS 50 is C, 80–83 ACAA). The new check refused it, and the fixture is corrected.
  • Mutations, each confirmed applied and caught by exactly the targeted test:
    • removing the mismatch error (warn-and-skip again) fails the SNP and indel refusal tests;
    • never checking indel REFs fails the indel test;
    • moving snp_count += 1 back before the checks fails the edge test;
    • removing the BED guard fails the BED test.
  • Full workspace suite (1,100 tests), fmt --check and clippy -D warnings pass locally. model_parity's mutation-model baseline is unchanged.

Not verified

  • Not run on a real VCF/reference pair, so how often real VCFs trip the refusal has not been measured.

Refs #819.

🤖 Generated with Claude Code

joshfactorial and others added 4 commits October 6, 2026 21:00
- test_transition_matrix_from_tsv asserts all four rows of the TSV; it
  passed with the TSV ignored.
- balanced_chimeric_offset's short-fragment test asserts the floor on
  both sides; its premise that the floor cannot hold at frag=read was
  wrong, and 'off >= 1' passed with the upper bound removed.
- Drop two max_reads smoke tests that could not see max_reads; three
  tests already pin it exactly.
- The degradation-floor config test checks what was stored.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
- test_runner_with_indels: assert variant-type weights (1/3 each),
  mutation_rate = 3/13114, and the fitted insertion (2) and deletion (3)
  lengths read from the written model. Fixture REFs now match H1N1.
- test_runner_skips_reference_mismatch_variant: assert the skipped SNP's
  context (GCT) carries no weight and the usable one (TCT) carries it all.
- config: assert exact r1/r2, VCF and BAM paths; rename
  test_overwrite_warn to test_overwrite_output_is_accepted, since no log
  capture exists to check the warning.
- generate_fragments: assert placed count and average depth for the
  three print-only depth tests.

Each was checked by mutating the production code it covers.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Each test now asserts a hand-derived expected value instead of existence
or is_ok(), and each was checked by mutating the covered production code
and watching it fail:

- fastq_tools test_apply_variants: fixture REF now matches the sequence,
  both variants homozygous, Q60; asserts read ATGTATGA and MMMMDMMMM.
- file_io test_create_output_file_new: reads the written bytes back.
- folder_tools test_check_parent: asserts the returned path; new
  should_panic test for a missing parent with create=false.
- sv_model extreme-lambda Poisson: +/-5 sigma bound instead of n > 0.
- filter_reads config: replaces the empty stub with known-answer tests of
  create_map_item and RunConfiguration::from.
- filter_lib test_prep_file_for_filtering: reads back through the
  returned reader and writer, both plain and gzipped input.
- frag length min_reads=0: a stray at 100,000 is kept at min_reads=0 and
  removed at min_reads=2; asserts both fitted Normal models.
- gc bias all-N contig: chr1 is now 200 bp so it reaches the N handling;
  asserts from the bin report that only chr2's two windows count.
- gen_mut_model empty transition_matrix_file: asserts the other parsed
  fields.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…819)

A SNP whose REF disagreed with the reference was skipped but still
counted toward snp_freq and mutation_rate; NEAT2 counted a SNP only once
it was usable. Count it nowhere, check indel REFs the same way, and
refuse the fit when more than 1% of checked records mismatch. A BED that
names no reference contig now says so instead of 'Unknown error'.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This was referenced Oct 7, 2026
Replace the 1% threshold with NEAT2's rule: a SNP or indel whose REF
disagrees with the reference refuses the fit. The whole VCF is scanned
first, so the error gives the SNP and indel counts and the first five
positions, plus the bcftools norm --check-ref x command that drops such
records deliberately. Removes the unmeasured constant and the counter
hidden in a match guard. Edge SNPs are still left out uncounted.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@joshfactorial joshfactorial changed the title gen-mut-model: drop REF-mismatched SNPs and indels, refuse above 1% gen-mut-model: refuse REF mismatches as NEAT2 did, with a bcftools fix in the error Oct 7, 2026

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant