Repository navigation
Name the cause when gen-seq-error-model has no weight to fit (#767) - #817
Merged
Merged
Conversation
Both paths the #765 audit traced to a bare InvalidWeights are reachable, and both now fail, or fit, for a stated reason. 1. MD tags that disagree with SEQ. TransitionObserver counted an MD mismatch position without checking that the read base differs from the reference base, so a stale tag recorded a self-transition; the diagonal is zeroed when the matrix is built, and a row holding only those had no weight. Such positions are now left out and tallied. The runner warns with the count, and when they are the only evidence it says the MD tags disagree with the reads, how many positions, and to rewrite them with samtools calmd. A BAM with some stale positions now fits from its genuine mismatches instead of failing. 2. A population with no bases. The coverage check counted transition positions only, which a 1 bp model has none of, so it could not fire: 1 bp R1 reads with an all-empty R2 reached the quality fit with no seed and failed as InvalidWeights([0.0]). The check now counts the seed: a population with none covers 0 bp and is refused by name. A pair of 1 bp files still fits. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Plain wording for comments that described populations whose reads stop at different lengths as 'ragged'. Test and fixture names are unchanged. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
A self-transition is the only detectable symptom of a stale MD tag; the same tag can mis-state the reference base at a real mismatch, and that counts as ordinary evidence. Fitting from the remaining mismatches would rest on tags already shown not to describe the reads, so any disagreement is now a hard error naming both counts and calmd. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
The #765 audit traced two inputs to
gen-seq-error-modelthat could still fail as a bareInvalidWeightserror. Both are reachable, and each now fails with a message naming the cause.Versioning: PATCH-level under
versioning.md. No config key, output format or model format changes. Inputs that failed still fail, now with a message naming the cause.1. MD tags that disagree with the reads
Reproduced: a BAM whose MD tag names a mismatch where SEQ has the reference base (MD
0A3on readACGT) failed withInvalidWeights([0.0, 0.0, 0.0, 0.0]).Cause:
TransitionObservercounted every MD mismatch position as a (reference → read) substitution without checking that the bases differ, so a stale tag recorded a self-transition. The matrix builder zeroes the diagonal, so a row holding only those had no weight left. One stale row was enough to fail an otherwise good fit.Fix: such positions are tallied (
md_seq_disagreements), and any disagreement is a hard error. The error gives the count of stale positions and of consistent mismatches, says the tags cannot be trusted, and points tosamtools calmd -b. It is separate from the existing "no MD tags" message.Why refuse rather than fit the rest: a self-transition is the only symptom of a stale tag that can be detected. The same tag can also name the wrong reference base at a real mismatch, or miss one entirely, and both look like ordinary evidence. Dropping the visible cases would leave a count built on tags already shown not to match the reads. An earlier revision of this PR fitted from the remaining mismatches; review caught it, and it now refuses the BAM.
read_bam_transition_reportreturns counts, masked and disagreement tallies together;read_bam_transitionsandread_bam_transitions_maskedare thin wrappers, so their tests are unchanged.2. A population with no bases
Reproduced: 1 bp R1 reads with an R2 file of empty records failed with
InvalidWeights([0.0]). With 2 bp R1 reads the same input was already refused clearly.Cause: the population coverage check counted positions after the first. A 1 bp model has none, so the check could not fire, and an R2 with no seed reached the quality fit.
Fix: the check counts the seed. A population with no seed covers 0 bp and is refused as
the R2 population covers 0 bp, but the model covers 1 bp …. This is the same rule the check already applied, now extended to position 1.How much this matters: the crash needs a contrived input (every R1 read 1 bp, every R2 record empty), which real data does not produce. The change is kept because the check was miscounting in ordinary cases too. It assumed every population had a first base, so a population with no data was reported as covering 1 bp. With 2 bp R1 reads and an empty R2, the input was correctly refused, but the message said "covers 1 bp"; it now says 0 bp.
Single-file cases did not reach
InvalidWeights, checked with the binary:Evidence
fmt --checkandclippy -D warningspass locally.Not verified
Closes #767.
🤖 Generated with Claude Code