feat: detect putative recombinant regions - #1774
Open
ivan-aksamentov wants to merge 44 commits into
Open
Conversation
Two-state HMM (wildtype vs recombinant) decoded by Viterbi in log-space. Three-valued observations (Ref, Mut, Missing) where Missing emits with probability 1 in both states, contributing no evidence while allowing transitions across uncovered stretches. Adds pathogen.json config schema with gamma, muW, muR parameters.
- gamma = 1/L: one expected state switch per genome - muW = mean terminal branch length / L: typical divergence of a new sequence - muR = median pairwise inter-clade leaf distance / L: expected mutation density in a region from a different clade
Run detection per-sequence in the parallel phase after private mutation calling. Surface skip reasons through the warnings output when detection cannot run (no tree, too few clades, no branch mutations).
Purple markers on the sequence view with tooltip showing "Putative recombinant", nucleotide range, and length.
Enforce invariants on construction: parameters in (0,1), gamma < 0.5 (state switching rarer than staying), muR > muW (recombinant state carries elevated divergence). Trim leading/trailing Missing positions from decoded intervals so reported ranges start and end at covered positions.
Use pairwise leaf-to-leaf distances via MRCA path lengths instead of founder-to-founder distances. Founder distances underestimate inter-clade divergence when founder branches are short relative to the tips.
… as missing Exclude deletions from per-branch mutation counts: a deletion is a single event, not a run of per-site mutations. Treat placement-masked positions as Missing to prevent false recombinant calls at homoplasic sites.
Nested RecombinationRegion with range, length, and optional confidence. Dot-notation TSV columns in a dedicated Recombination output category. Rec. column in the results table with sortable total length and tooltip.
Forward-backward algorithm computes per-site posterior marginals in log-space. Each Viterbi-decoded interval receives a confidence score equal to the mean posterior marginal within it. Higher confidence produces more opaque markers in the viewer.
Show markers in absolute and all relative views. Add per-region confidence to the Rec. column tooltip. Cap the tooltip region list before building React elements.
Resolve each HMM parameter independently: explicit pathogen.json values bypass estimation for that parameter. Validate RecombinationHmmParams on deserialization via TryFrom. Parse branch mutation tokens structurally to classify deletions by gap base. Add minPrivateSubsToRun gate (default 1).
Introduce shared test macros (pretty_assert_ulps_eq!, pretty_assert_abs_diff_eq!, assert_error!). Derive PartialEq on result types and compare whole values.
Add algorithm page to user docs. Document pathogen.json recombination config, output columns, and parameter bounds. Fix broken pipe escaping in the QC-status TSV table.
Decision record covering model specification, parameter estimation, observation model, pipeline integration, configuration, output format, and viewer behavior. Includes follow-up proposals and minimizer scoring notes.
Add explicit lifecycle rules: a resolved issue or an implemented/resolved proposal no longer belongs to any KB category and must be deleted (or, for proposals, moved to `decisions/`), never marked "resolved" in place. Git history preserves the record. Remove the mu_r estimator proposal, which had been left in `proposals/` with a "resolved" status after its design shipped.
Point the gamma and muW references at the current `estimate_gamma` and `estimate_mu_w` functions in `recombination_estimate.rs`; the old line ranges had drifted onto unrelated code as helpers were added above them.
2 tasks
- Separate from the Rust/web builder (docker/dev); provides numpy, scipy, pandas, matplotlib, JupyterLab via micromamba for running capture scripts and notebooks against the project
- Pin `viterbi_decode` and `forward_backward_marginals` against the `recomb_inference` Python prototype at commit 944ec48 - Two tolerance regimes for forward-backward: 1e-9 (ordinary scale) and 1e-7 (genome scale, absorbs the prototype's `eps = 1e-10` log-smoothing)
…ummaries - Viterbi optimality vs brute-force over all hidden paths (random obs and params) - log_sum_exp_2 matches naive, commutative, lower-bounded - build_observations length preservation and missing-precedence - RecombinationResult summary consistency; forward-backward all-Missing marginals - interval-confidence bounds; params JSON round-trip (ULP-bounded) - estimator median invariants and verbatim parameter overrides
- Guard the 1-based-closed range rendering and 3-decimal confidence formatting of the recombination CSV columns against the documented output contract. - Pin the recombination column category to its default-on toggle and its six-column map, so a change to the selection wiring surfaces in tests.
…teriors Extract the forward-backward alpha/beta recurrence into forward_backward_posteriors, returning the full per-site [P(wildtype), P(recombinant)] distribution, and project forward_backward_marginals from it. Keeps the recurrence in one place and exposes both states so the posteriors can be checked against an independent oracle.
Replace the tautological forward-backward sum-to-one test, which re-derived the alpha/beta recurrence in the test body and asserted an identity guaranteed by the normalizer, with a brute-force marginal oracle: enumerate all hidden-state paths for short vectors and compare production posteriors against the independent marginalization. Parameterize the log-sum-exp naive-oracle cross-check as rstest cases so each operand pair reports independently. Assert the serde rejection message with assert_error!. Compare serialized RecombinationResult JSON against exact indoc literals instead of substring checks. Drop dead #[allow(clippy::float_cmp)] on ulps-based estimator tests. Annotate the qualitative forward-backward thresholds as loose sanity bands.
- Move the two large inline `mod tests` blocks out of `recombination.rs` and `recombination_estimate.rs` into topic-scoped files under `analyze/__tests__/`, sharing observation builders, tree fixtures, and the brute-force oracles through one helper module instead of duplicating them per file. - Widen a few internals (`WILDTYPE`, `RECOMBINANT`, `intervals_sorted_disjoint_nonempty`, `log_sum_exp_2`, `log_emission`, `median`, `as_probability`) to `pub(crate)` so the relocated white-box tests reach them; the validation-gated `RecombinationHmmParams` fields and the `RecombinationHmmParamsRaw` wire type stay private, with `new` asserted through its getters.
- Separate the single detection file into submodules by concern (parameters, observations, Viterbi decoding, forward-backward scoring, config, result), so each concern is independently readable and testable instead of interleaved in one file. - Group the tree-based parameter estimator and the co-located tests under the same `recombination/` directory as its family rather than as loose analyze-level siblings.
- Move the detection pipeline out of the analysis caller into `run_recombination`, so the module owns its step order (gate, observations, decode, scoring) instead of exposing loose building blocks the caller must sequence correctly. - Demote those building blocks to crate-internal now that only the façade is called from outside the module.
- Noun-phrase function names renamed to verb-led forms across the module (e.g. `recombination_missing_ranges` → `collect_missing_ranges`, `median` → `compute_median`, `log_emission` → `compute_log_emission`) - All call sites, imports, and doc-comment cross-references updated in source and test files
The pattern used a trailing slash, so it matched only a real node_modules directory. Worktrees symlink node_modules into the checkout, and a symlink is not a directory, so `git add -A` staged it. Dropping the slash ignores both forms.
…documentation - Add user-facing algorithm description and pathogen-config documentation - File KB issue for insertion-calibration mismatch between estimator and decoder - Remove resolved test-coverage tracking proposal - Condense verbose doc comments, schema descriptions, and KB entries
- Reuse inclusive ancestor iteration and cached leaf paths to make pairwise distance calculation explicit and avoid repeated tree walks. - Use established statistics and iterator APIs while preserving fallback behavior with focused boundary coverage.
Replace the hand-rolled `compute_median` with `statrs`'s `OrderStatistics::median`, the same crate already used for the muW mean in this module. The estimate now returns early on an empty distance set, keeping the undefined-rate contract without a separate helper. Drop the helper's dedicated unit and property tests: they exercised third-party median behavior. Odd, even, and empty cases remain covered end-to-end through the muR resolution tests.
The branch audit surfaced hot-path costs on the recombination decoder: a per-site `ln()` recompute, per-sequence scratch allocation with oversize backpointers, and quadratic muR inter-clade pairing. The obvious fixes are un-modeled (a positional `[[f64;2];3]` emission table, `thread_local!` scratch pooling), so record properly typed, exactness-preserving designs as issues instead of applying crude fixes. The muR pairing issue gains an exact-fast option (LCA preprocessing removes the per-pair root walk) and a bounded-memory option (integer-histogram streaming median) alongside the existing subsample tradeoff, since inter-clade distances are integer substitution counts.
Pair each decoded region with its own optional confidence at the `from_ranges` input, so a per-region confidence-count mismatch cannot be constructed. This removes the release-mode index panic (short slice) and silent-accept (long slice) that the previous `debug_assert`-only length guard left open. Also replace the unreachable `count == 0` confidence branch (it returned a plausible `0.0` rather than surfacing the violation) with a `debug_assert` on the decoder's non-empty-interval postcondition, and document the longest-region invariant with `expect` instead of a bare `unwrap`.
The `gamma` schemars annotation advertised `maximum: 1.0`, but the HMM validator requires `gamma < 0.5`. A `pathogen.json` with `gamma` between 0.5 and 1.0 passed generated-schema validation and then failed at dataset load. Set the schema bound to `0.5` and regenerate the input-pathogen and auspice-extensions schemas. `muW`/`muR` open-interval endpoints stay documented in the field descriptions: schemars 0.8 cannot express exclusive bounds in the derive.
Adding the `includeRecombination` toggle without a serde default would reject an older serialized `CsvColumnConfig` that predates the field. Deserialize missing fields from the existing `Default` impl (the single source of truth), so shipped configs keep loading and each toggle takes its established value.
- Add a direct `compute_log_emission` test with hand-derived `ln` values, independent of the brute-force oracles (which reuse the same emission function) so a Ref/Mut or wildtype/recombinant column swap is caught. - Collapse per-element assertion loops: the brute-force marginal test compares each state column as a whole slice (`approx` implements `AbsDiffEq` for `[f64]`), and the qualitative-band and confidence loops become single min/max assertions that report the offending extreme. - Add boundary cases: zero reference length for `build_observations`, a per-field `muW` out-of-range override, and full `RecombinationConfig` camelCase serde round-trip. - Use `json_parse` for golden-master fixtures; drop assertions entailed by a preceding whole-value equality; pair range and confidence at the `from_ranges` call sites. The serde round-trip keeps a 2-ULP tolerance: serde_json's default float parser is best-effort, not bit-exact (its `float_roundtrip` feature would guarantee exactness), so a bit-exact assertion flakes by ~1 ULP. The genome-scale golden-master tolerance stays at 1e-7: the measured max marginal divergence from the prototype's `log(0)` smoothing is 3.19e-8, so 1e-8 does not pass.
Repoint KB issue and proposal links, and the golden-master capture-script header, to the split `recombination/*.rs` files: the monolithic `recombination.rs` and `recombination_estimate.rs` no longer exist. Stale line anchors are dropped, since the split changed line numbers; the links now target the correct concern file. Also convert the algorithm page's definition lists from `term -- description` to `term: description`.
`ColumnRecombination` and `SequenceMarkerRecombination` each reimplemented percentage rendering. Extract `formatPercentage` and `formatPercentageOfTotal` so the precision and the zero-denominator policy live in one place.
- `ln(0.1)` in the emission test equals `-LN_10`; the hand-derived `ln` expected values are a deliberate oracle, not an approximation of a named constant.
The muW/muR tree-based estimates counted insertion tokens (gap reference, e.g. "-10A") as branch substitutions, but the decoder never observes insertions: its `Mut` emissions come from `private_substitutions` only and `build_observations` has no insertion channel. Counting them calibrated the emission rates against events the model cannot score, inflating both rates on externally produced trees (augur/TreeTime) that carry insertion tokens. Count a branch mutation only when both reference and query bases are present (non-gap), so calibration matches the substitution-only event set the decoder sees. Nextclade-built trees carry substitutions only and are unaffected.
Reconcile the decision record with the substitution-only calibration: both muW and muR count a branch mutation only when reference and query bases are present, so neither insertions nor deletions contribute. Update the nextclade_data mirror description accordingly and remove the now-resolved known issue.
- Function-level `#[allow]` does not suppress warnings from `#[case]` attribute arguments; module-level does
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.
nextclade_data: #464Try it
Provenance
Design and function contract specified in the Nextclade multi-reference meeting (2026-06-26). Implementation follows a Python prototype and its algorithm notes.
Motivation
A recombinant sequence -- one whose genome is a mosaic of segments from different lineages -- produces a localized cluster of private mutations when placed on the reference tree. Nextclade currently has no way to distinguish this pattern from a poorly sequenced region, a novel variant, or contamination. The private-mutation QC rule penalizes all three equally.
This PR adds a per-sequence detector that identifies where the breakpoints are, without naming the donor lineage. Downstream consumers can mask the flagged regions and re-place if needed.
Design
The full specification is in the decision record, which covers the model, estimation, pipeline integration, configuration, output, viewer, and deferred work. The summary below is for reviewers who want to evaluate the approach without reading the full document.
Model
A two-state Hidden Markov Model walks the alignment with states wildtype ($w$ ) and recombinant ($r$ ), decoded by Viterbi in log-space [src]. At each reference position the model sees one of three observations relative to the sequence's inferred parent:
RefMutMissingMissingpositions contribute no evidence but allow transitions to persist across uncovered stretches, so a recombinant region bridging an N run is not split. Deletions areMissingbecause a large deletion is one mutational event, not a cluster of per-site mutations that would inflate local density and trigger false calls. Placement-masked sites areMissingto prevent false calls at homoplasic positions.The joint likelihood for observation sequence$s$ given hidden states $h$ is:
where$s = (s_1, \ldots, s_L)$ is the observation sequence, $h = (h_1, \ldots, h_L)$ is the hidden state sequence with $h_l \in {w, r}$ , $P(s_l \mid h_l)$ is the emission probability, and $T$ is a symmetric transition matrix with off-diagonal $\gamma$ .
Parameters
Three parameters govern the model [src]:
Parameters are resolved once per dataset, not per sequence. Each parameter resolves independently: an explicit$\mu_r$ uses leaf-to-leaf MRCA path distances, not founder-to-founder distances, because short founder branches underestimate inter-clade divergence.
pathogen.jsonvalue bypasses estimation for that parameter. Estimation counts substitutions and insertions but excludes deletions (a deletion token is identified structurally by a gap query base, not by string suffix).Each state transition costs$c = \ln((1-\gamma)/\gamma)$ nats; the decoder switches only when accumulated emission evidence exceeds $c$ . For a 30 kb genome, $c \approx 10.3$ nats.
Three invariants are enforced on every construction and deserialization path: all parameters in$(0,1)$ ; $\gamma < 0.5$ ; $\mu_r > \mu_w$ .
Forward-backward confidence
After Viterbi identifies intervals, the forward-backward algorithm [src] computes per-site posterior marginals$P(h_l = r \mid s)$ in log-space. Each interval receives a confidence score equal to the mean posterior marginal within it -- a per-call reliability measure that Viterbi's single hard path cannot express. Values near 1 indicate high posterior certainty; values near 0.5 indicate ambiguous evidence.
Forward-backward runs only when Viterbi finds at least one interval, avoiding the$O(L)$ memory cost for the common case.
Pipeline integration
Detection runs per-sequence in the parallel phase, after private mutation calling [src]:
minPrivateSubsToRun(default 1) -- a sequence with no private substitutions has an all-Ref/Missingobservation vector, so the recombinant state can never outscore wildtypeReffor aligned positions,Mutfor private substitutions, then non-comparable ranges override toMissingMissingConfiguration
pathogen.jsonsnippet (all fields optional):{ "recombination": { "enabled": true, "minPrivateSubsToRun": 1, "gamma": 3.3e-5, "muW": 0.001, "muR": 0.01 } }enabled: absent = default-on,true= explicit opt-in,false= disabledminPrivateSubsToRun: skip sequences below this count (default 1)gamma,muW,muR: when set, bypass tree-based estimation for that parameter. Must satisfyThe companion
scripts/recombination_paramsinnextclade_datacomputes these values from a dataset's tree and reference, for freezing reviewed parameters intopathogen.json.When a dataset cannot support detection (no tree, fewer than two clades, no branch mutations):
enabledabsent): silently skippedenabled: true: dataset-level error naming the causeInvalid explicit parameters are always errors. The skip reasons are enumerated by
RecombinationSkipReason[src].Output
Recombination is an observation about sequence origin, not a quality metric, so results go into dedicated output columns rather than the QC system.
JSON:
recombinationfield on each sequence result, omitted when no intervals found:{ "regions": [{ "range": { "begin": 100, "end": 200 }, "length": 100, "confidence": 0.95 }], "totalRegions": 1, "totalLength": 100, "longestRegion": { "range": { "begin": 100, "end": 200 }, "length": 100, "confidence": 0.95 } }TSV: dot-notation columns in a dedicated
Recombinationcategory (recombination.regions,recombination.regionConfidences,recombination.totalRegions,recombination.totalLength,recombination.longestRegion.range,recombination.longestRegion.length).Viewer
#8f2fd4fill,#e0b3ffborder) in all nucleotide views (absolute and relative)Rec.column in the results table, sortable by total recombinant length, with a tooltip showing region count, total length, longest region, and per-region confidenceDifferences from the prototype
Missingwith emission probability 1 instead of binaryMissingto prevent false calls at homoplasic positionsMissingstripped so reported ranges start and end at covered positionsTryFromin addition to range checksWork items
RecombinationRegionwith dot-notation TSV columns andRec.table columnminPrivateSubsToRungate, deserialization validation, structural branch-mutation parsingpretty_assert_ulps_eq!,assert_error!) and whole-value comparisonKnown issues
Possible improvements