Correct sequence extraction, contig-boundary handling and pairing orientation; migrate docs to Zensical - #82
Merged
Merged
Conversation
ruff v0.14.11 -> v0.16.1, numpydoc v1.10.0 -> v1.11.0rc0, actionlint v1.7.8 -> v1.7.12, yamlfmt v0.20.0 -> v0.21.0. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_014Vs6ZCghyMszqBMZ2jKJa3
Sequence extraction from an indexed FASTA (--genome) and from a BLAST database (--blastdb) are presented as interchangeable, but produced different sequences for the same coordinates. Each of ~10 call sites reimplemented coordinate conversion, clamping and strand handling in its own `if blastdb:` branch, and the two branches had drifted. No test exercised blastdbcmd at all. Add src/tirmite/utils/extract.py as the single extraction primitive: 1-based inclusive plus-strand coordinates in, uppercase plus-strand sequence out, clamped and length-verified. All call sites route through it, removing the duplication that allowed the drift. Backend divergences fixed: * With --padlen on a minus-strand hit, the lowercase pad markers were computed in forward-genomic space but applied to a string blastdbcmd had already reverse-complemented, so they landed on the wrong end and mislabelled which bases were the hit. Case marking now happens on the plus strand, before reverse complement. * Elements and left TIRs came back reverse-complemented from the BLAST path but plus-strand from the FASTA path. Plus strand is now canonical, matching the GFF coordinates and the flank/TSD convention. * blastdbcmd silently returns the ENTIRE sequence, exit 0, when a range starts past the contig end. Clamping plus a returned-length check now rejects that rather than substituting a whole contig for a region. * blastdbcmd discards soft-masking while pyfaidx preserves it; both are now uppercased, so masked genomes cannot make the backends disagree. * Record IDs, descriptions and failure modes are now identical: '_rc' suffixes, genome descriptions (read from the BLAST title when needed), and a missing contig returning None from both. Contig-boundary handling, previously silent truncation: * Add fetch_region_padded: partial overlap is padded with N to the requested width and reports the pad counts; a region with no overlap still returns None rather than an all-N record. Applied to flanks, in-model TSDs and the --padlen window; padded records carry padded:L,R in their description. --no-pad-flanks restores truncation. * Deliberately NOT applied to hmm_build seed extraction, whose output becomes HMM training data where Ns would become model columns. * A TSD truncated at a contig end came back short, the length mismatch made the comparison be skipped, and the result was written as tsd_hamming=0 - indistinguishable from a verified duplication. Add compare_tsds, which excludes padded positions and reports NA when there is no informative position. * Clamp model-coverage deficits to >= 0. A model_len below hmm_end (wrong --lengths-file, mismatched HMM) made the offset negative and shifted the flank window into the hit, extracting element sequence as flanking; the flank_max_offset guard could not catch it. * Reject inverted hit windows in fetch_padded_hit instead of emitting a nonsense arrangement of flanks around nothing. * Fix 1-based nhmmer coordinates being fed into a 0-based pyfaidx slice in hmm_build, which dropped the first base of every extracted hit. Add --extend-hits-to-model to emit partial hits at full model length, using the deficits already computed for flank projection. The offset correction for hits that only partially cover the model was audited and is correct on the appropriate side in all four strand/terminus combinations; tests now pin that, including hits incomplete at both ends. Verified with 405 tests, including a differential suite that builds a FASTA and a BLAST database from the same file and asserts byte-identical results, and end-to-end runs of the pipeline both ways whose outputs now diff clean. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_014Vs6ZCghyMszqBMZ2jKJa3
Material for MkDocs has entered maintenance mode. Zensical is the successor from the same team, renders the same HTML, and reads mkdocs.yml natively - so the Markdown sources needed no changes. Replace mkdocs.yml with zensical.toml. Everything carries over: the teal/cyan palette in both schemes, all eight theme features, logo, favicon, nav, repo and edit links, and the social icon. theme.variant is set to "classic" so the site keeps the familiar Material appearance rather than silently switching readers to Zensical's newer default look. Two simplifications fall out of the move: search is built in, so the `plugins: [search]` entry is gone, and mermaid's custom fence takes the format as a plain string instead of the `!!python/name:` YAML tag. The deploy workflow now builds with `zensical build --clean --strict` and publishes via the Pages artifact rather than pushing a gh-pages branch, so it no longer needs write access to the repository. This requires the repository's Pages source to be set to "GitHub Actions". Also correct usage examples in the tutorials, which had drifted from the CLI. Every command in docs/ and README.md was checked against the actual argparse definitions; these two files carried the bulk of the errors: * `tirmite search` was documented with an older interface throughout: --hmm-file is --hmm, --blast-query is --fasta, --nhmmer-file is --nhmmer-results, --maxeval is --max-evalue/--hmm-max-evalue/--blast-max-evalue, and --query-len (an int) is --lengths-file (a file). --mincov does not exist on search at all. * `tirmite seed` was shown with --max-gap, removed in 1.4.0, and --blast-file, which is --blast-hits. The pre-computed-hits example piped raw blastn output into it, but that option takes TIRmite's own format from --save-blast-hits. * The HMM-update example omitted --update, without which --hmm-file does nothing. * Asymmetric output files are named `_left.hmm`/`_right.hmm`, not `_LEFT`/`_RIGHT`. * The search-to-pair handoff examples omitted --lengths-file, which --blast-file and asymmetric BLAST pairing both require. Verified by running the corrected commands: they succeed and produce exactly the filenames now documented. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_014Vs6ZCghyMszqBMZ2jKJa3
Audit of the pairing logic across all four --orientation settings, for both symmetric and asymmetric elements, checking that the external side of an element is correctly identified before flanking regions are extracted. parseHitsGeneral enumerates two strand combinations per orientation - the canonical one and its reverse complement - so a single run legitimately pairs elements inserted in either genomic direction. For a reverse-inserted asymmetric element the config's left model then sits at the HIGHER coordinate. Two subsystems disagreed about what "left terminus" meant: flipTIRs is purely positional and drives every paired consumer, while _determine_terminus_type and filter_hits_by_anchor derive it from model identity, ignoring strand. They coincide only for forward insertions, and nothing reconciled them. Paired flank and TSD geometry was already correct - is_left is positional and combines with the hit's own strand in compute_flank_coordinates - and is left unchanged, now with tests pinning it. Add resolve_terminus(), which separates terminus ROLE (which end of the element, from the model) from is_lower (which genomic side its outer edge faces, from the strand), and route the fixes through it: * filter_hits_by_anchor fed a role where a genomic side was expected, so for a reverse insertion --max-offset was measured against the hit's INNER edge and valid hits were silently discarded before pairing. Highest impact, since the loss is invisible downstream. * Unpaired asymmetric flanks were taken from the lower-coordinate side regardless of strand, i.e. from inside the element for a reverse insertion. Wrong sequence, not just a wrong label. * A pairing-map row naming the same feature twice had 'left' immediately overwritten by 'right', so symmetric models never reached the both-ends anchor test. * Asymmetric output now routes by terminus role, so each model's termini stay in that model's file; previously a reverse insertion wrote the right model's terminus into the left model's file. Sequence and coordinates are unchanged. Reconstructed target sites gain an element_orientation field. * --orientation is parsed and validated once, in PairingConfig, and reused by filter_hits_by_anchor. 'f,r' previously set both strands to '-' in pairing while the anchor filter read it correctly; a malformed value left the strand table empty and produced zero pairs with no error. Two further defects found while investigating, both confirmed by test first: * Symmetric F,F never paired at all. The if/elif meant every hit took only the left branch when both strands match, so no hit collected an upstream candidate and reciprocity could never be established. F,F is the documented LTR case and returned zero pairs; the existing F,F tests pass because they drive the asymmetric path via two model names. * Symmetric R,R paired hits with themselves. _find_candidates never excluded the reference hit, and on the minus strand the 5'/3' swap makes the self-distance positive, so a lone hit produced a "pair" of one. Align the two meanings of --maxdist. The legacy path measured the literal gap between hits; the current path measured to the partner's strand-relative 5' end, which for a minus-strand partner is its far edge, adding the whole length of that terminus. Both now use one shared inter_hit_distance(): the gap between the facing inner edges, independent of strand and model length. Candidate ranking uses the same measure as the filter. This also rejects overlapping termini, which the old measure accepted. Runs that leave --maxdist unset are unaffected; real Sanctuary output is byte-identical. An existing --maxdist admits a slightly narrower set of pairs, short by roughly one terminus length. Every fix is mutation-checked. Verified end-to-end on a genome carrying the same asymmetric element in both orientations: both pair, all four flanks come from outside the element, each model lands in its own file, and --genome and --blastdb outputs stay byte-identical. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_014Vs6ZCghyMszqBMZ2jKJa3
Coverage was enabled by addopts in pyproject.toml, so every matrix job instrumented the code and wrote coverage.xml even though only the 3.14 job uploads to Codecov. Pass --no-cov on the other versions. The extraction equivalence tests shell out to makeblastdb and blastdbcmd and are skipif-gated on those binaries, which CI never installed. That meant 69 tests - the whole FASTA-vs-BLAST-database equivalence suite - skipped silently while the run reported green. Install BLAST+ for the 3.14 job only, so the other five stay fast, and verify the binaries in a separate step so a broken install fails loudly instead of quietly reverting to skips. Locally, without BLAST+ on PATH: 386 passed, 69 skipped. With it: 455 passed. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_014Vs6ZCghyMszqBMZ2jKJa3
Covers the sequence-extraction unification, contig-boundary padding, pairing orientation fixes, the --maxdist alignment, and the move to Zensical. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_014Vs6ZCghyMszqBMZ2jKJa3
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #82 +/- ##
============================================
- Coverage 100.00% 41.75% -58.25%
============================================
Files 1 11 +10
Lines 5 5276 +5271
Branches 0 1095 +1095
============================================
+ Hits 5 2203 +2198
- Misses 0 2866 +2866
- Partials 0 207 +207 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
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
1 participant
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.
Summary
Three related pieces of work, plus infrastructure updates. The headline items are correctness fixes in sequence extraction and terminus pairing, several of which were producing wrong sequences or silently discarding valid hits.
See the 1.5.0 CHANGELOG section for the full list.
Sequence extraction was backend-dependent
--genomeand--blastdbare presented as interchangeable, but returned different sequences for the same coordinates. Each of ~10 call sites reimplemented coordinate conversion, clamping and strand handling in its ownif blastdb:branch, and the two had drifted.All sites now route through a single primitive (
tirmite/utils/extract.py) with one documented contract. Notable fixes:--padlenon a minus-strand hit, lowercase pad markers were applied to a sequenceblastdbcmdhad already reverse-complemented, marking the wrong end and mislabelling which bases were the hit.blastdbcmdsilently returns the entire sequence, exit 0, when a range starts past a contig end. Now rejected by a returned-length check.blastdbcmdbut preserved bypyfaidx; both are now uppercased.pyfaidxslice intirmite seed, dropping the first base of every hit.Before this change 8 of 9 output files differed between backends; they are now byte-identical, proven by a differential test suite that builds a FASTA and a BLAST database from the same file.
Contig boundaries
Regions overrunning a contig were silently truncated. They are now N-padded to the requested width and marked
padded:<l>,<r>;--no-pad-flanksrestores the old behaviour.The consequential fix here: a TSD truncated at a contig end was reported as
tsd_hamming=0— the length mismatch made the comparison be skipped, leaving a truncated, unverifiable TSD indistinguishable from a verified perfect duplication. This silently corrupted the signaltirmite validateexists to check. Unverifiable TSDs now reportNA.Pairing orientation
parseHitsGeneralaccepts both strand combinations per orientation, so a single run pairs elements inserted either way round. For a reverse-inserted asymmetric element the config's left model then sits at the higher coordinate, and two subsystems disagreed about what "left terminus" meant.F,Fnever produced any pairs — the documented LTR use case returned nothing.R,Rpaired hits with themselves; a lone hit produced a "pair" of one.--max-offsetmeasured the wrong model edge for reverse insertions, silently discarding valid hits before pairing.Paired flank geometry was already correct and is unchanged, now with tests pinning it.
Behaviour changes to review
--maxdistnow measures the gap between the facing inner edges of the two hits. The pairing path previously measured to the partner's far edge, adding a whole terminus length, and disagreed with the legacy path. Runs leaving it unset are unaffected; an explicit value now admits a slightly narrower set of pairs._L/_Rsuffixes follow terminus role rather than genomic order.reconstruct_target_sitereturnsNoneinstead of0for an uncomparable TSD.Docs
Migrated from Material for MkDocs (now in maintenance mode) to Zensical.
mkdocs.yml→zensical.toml; Markdown sources unchanged, site appearance unchanged.Also corrected the tutorial usage examples, which had drifted from the CLI —
tirmite searchwas documented with an older interface throughout, and several--insertion-siteexamples omitted the required--flanks/--flanks-pairedand would have failed at runtime.Testing
455 tests pass (320 before this branch). Every fix is mutation-checked: reintroducing the defect fails the suite.
Verified end-to-end by running the pipeline both ways and diffing outputs, on real Sanctuary/BMP2768_Av data and on synthetic genomes covering contig edges, soft-masking, and the same asymmetric element in both orientations.
Action required before merge
gh-pagesbranch; until this is changed the site will keep serving the stale branch.🤖 Generated with Claude Code
https://claude.ai/code/session_014Vs6ZCghyMszqBMZ2jKJa3