Repository navigation
Simulate ancestral alignments under the TKF and substitutions models - #140
Conversation
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## develop #140 +/- ##
===========================================
+ Coverage 96.32% 96.76% +0.43%
===========================================
Files 32 52 +20
Lines 4355 6486 +2131
===========================================
+ Hits 4195 6276 +2081
- Misses 160 210 +50 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
- Replace local WILDCARD_CHAR constant (b'N') with the shared AMB_CHAR from crate::alphabets for consistency - Move TKF92MSASimulationResult struct definition after the link/fate enums so related types are grouped more logically - Add clarifying doc comments on Seqs type alias and the indel_model field
- Remove the local tkf92_fixed test helper and naive_merge utility, which duplicated indel cost logic already encapsulated in TKF92FixedIndelCostBuilder - Update the simulate test to build and invoke the cost via TKF92FixedIndelCostBuilder, keeping the assertion equivalent - Clean up unused imports (AncestralAlignment, h1, log_i1, log_n1, get_mapping_for_any_node)
- Add FragmentSampler trait with sample_fragment_length<R: Rng> -> (usize, f64), returning the sampled length and its log-probability - Implement FragmentSampler for TKF91IndelModel (always returns length 1, log-prob 0.0) - Implement FragmentSampler for TKF92IndelModel (geometric draw using r parameter) - Rename TKF92MSASimulator -> TKFMSASimulator<T: TKFModel + FragmentSampler, Q, R>, replacing the hardcoded TKF92IndelModel field with the generic T - Delegate sample_fragment_length on the simulator to the trait, eliminating all model-specific logic from the simulator itself - Rename simulate test to tkf92_simulate; add tkf91_simulate test that asserts every fragment has length 1 and logl matches TKF91IndelCostBuilder
- Replace root_residue_count accumulator with fragmentation.last().unwrap_or(&0) + current_link.length - Drop the root_residue_count: usize parameter from append_link_to_msa and links_to_msa - Boundary computation is now self-contained: each entry is the previous boundary plus the current fragment's length, which is correct because fragmentation is a monotonically increasing cumulative sequence
- Extract dispatch_immortal_children to handle immortal link traversal - Extract process_link to handle non-immortal link MSA writes - Extract apply_fate to apply a single branch fate (homolog/deletion/non-homolog) - Add explicit 'a lifetime annotations to process_link, apply_fate and dispatch_immortal_children so that references pushed onto tree_stack and insertions are tied to the lifetime of the source TKFLink data
- Rename TKF92MSASimulationResult to TKFMSASimulationResult<AA> - Make simulate_msa() return typed struct instead of tuple - Add msa_to_alignment_with_non_emitting_cols() method - Fix typos: didnt→didn't, comulative→cumulative, bc→because - Rename test tkf_homlog_probs to tkf_homolog_probs
- Move trait bounds into a clause for readability in
- it seems to be simpler to write simulation logl match in tests only
- Add SubstitutionSimulator and SubstitutionSimulatorBuilder in src/substitution_models/simulate.rs - Precompute transition matrices P(edge) once and reuse for all sites - Fixed alignment length provided via builder (set at construction) - Root sampled from model equilibrium freqs; traverse preorder to sample children - Export new module in src/substitution_models/mod.rs - Add unit tests for correctness and reproducibility
Replace EvoModel with SubstModel<Q> in substitution simulator - Make builder and simulator generic over Q: QMatrix and store SubstModel<Q> - Reimplement p(time) locally when precomputing P = exp(Q * t) to avoid EvoModel dependency - Access frequencies and alphabet via the underlying QMatrix - Update tests to construct SubstModel<JC69> (no API changes to tests) Files changed: - src/substitution_models/simulate.rs Rationale: keeps simulator tied to concrete substitution model representation and avoids using the dynamic EvoModel trait where unnecessary. This simplifies access to Q and its frequencies while preserving behavior and reproducibility.
Refactor generics in substitution simulator to use where clauses - Move generic bounds into clauses for builder, simulator and impl blocks - Improves readability by reducing line length - No functional changes File changed: - src/substitution_models/simulate.rs
- Replace JC69 with GTR in substitution simulator private tests - Provide concrete GTR frequencies and rate parameters for deterministic behavior - Keep existing tree! macro usage unchanged This makes tests use a realistic, parameterized DNA model rather than JC69, improving test coverage for general substitution models.
- Rename substitution simulator module file to - Update module export in to use - Preserve tests and functionality while switching filename to better reflect purpose
d4b2a1d to
1e31b3a
Compare
Add default method to trait and support methods in MASA. - Detect columns where all leaf maps are and remove them - Update leaf and ancestral mappings and sequences accordingly - Add helpers on MASA to efficiently update internal state - Note: update deferred (see issue acg-team#150)
- Add complex test `into_alignment_masa_to_msa` in `src/alignment/tests.rs
| } | ||
|
|
||
| /// Rebuilds a clean tree from a (possibly mutated) tree, dropping any nodes that are no | ||
| /// longer connected to the root. Serialization to Newick and re-parsing re-computes all |
There was a problem hiding this comment.
Uhhhh, I understand why this is the easiest way to handle tree creation in this scenario, but this seems wrong somehow... using serialisation as a tool for anything other than output.
Would there be an easy way to create the tree without resorting to newick that we could add?
There was a problem hiding this comment.
Lets discuss in person if we even want the prune_empty_leaves method (its the only place where this rebuild_tree is used.
There was a problem hiding this comment.
If I remember correctly we decided to keep it for now.
|
|
||
| /// Returns all sequences (leaf and ancestral) in aligned form, mirroring | ||
| /// [`Display`](`std::fmt::Display`) for MASA. Node ids are resolved through the tree, | ||
| /// since `idx_to_id` is private to the alignment module. |
There was a problem hiding this comment.
This also somehow feels wrong that you have to resolve node ids manually. Could we implement a helper that isn't relying on private indices, or work around this a different way?
|
(An alternative would be to have this code as a separate crate that uses the lib via its API) |
- removed into_alignment, its now discard ancestral seqs in phyloinfo - use our rng in fragment sampler - keep col mask is now bitverctor - sim subst msa split into multiple methods
|
|
||
| let aligned = self.aligned_seqs(); | ||
| let pruned_tree = prune_empty_leaves_from_tree(&self.tree, &leaves_to_remove) | ||
| .unwrap_or_else(|e| { |
There was a problem hiding this comment.
wouldn't it make sense to call it with a ? instead, since prune_empty_leaves_from_tree also returns a result?
There was a problem hiding this comment.
changed the method to panic on bug.
| } | ||
|
|
||
| impl TKFLink { | ||
| fn new_immortal(node: NodeIdx) -> Self { |
There was a problem hiding this comment.
could TKFLink be an enum type instead? So that instead of checking if link.is_immortal you could do a match?
There was a problem hiding this comment.
I think the refactor is not worth the hustle.
There was a problem hiding this comment.
Only one comment (phyne/src/alignment/mod.rs Line 140) is really a potential issue here, the rest are nitpicks in reality.
Good job!
This PR implements an
AlignmentSimulationtrait which simulates an ancestral alignment and implements it for the TKF91 and TKF92 models and substitution models. It also introduces a utility for removing non surviving homology paths and ancestral alignment transformation into standard leaf-only alignments.Key Changes