Skip to content

Feature/arpip impl - #211

Open
junniest wants to merge 48 commits into
developfrom
feature/ARPIP_impl
Open

junniest wants to merge 48 commits into
developfrom
feature/ARPIP_impl

Conversation

@junniest

@junniest junniest commented Nov 4, 2026

Copy link
Copy Markdown
Contributor

Implementation of indel event reconstruction under PIP and refactoring.

ARPIP implementation:

  • Added the asr module to pip_model with a PIPIndels struct that implements AncestralSequenceReconstruction.

API Improvement:

  • Added the from_msa_with_ancestral and from_msa_with_ancestral_unchecked methods to the AncestralAlignment trait, to construct a PhyloInfo with ancestral sequences from an existing alignment and a set of ancestral sequences.
    from_msa_with_ancestral includes default implementation with validation checks for matching alphabets and sequence lengths between MSA and ancestral sequences and matching node IDs.

Utility Functions and Refactoring:

  • Introduced the mapping_to_bits utility function, which converts a sequence mapping to a FixedBitSet (will be potentially refactored away when the Mapping type changes to FixedBitSet from Vec<Option<usize>>).
  • Refactored the alphabet method in the Alignment trait and its implementations to always return a reference to a static Alphabet.
  • Moved validation logic for ancestor IDs to a dedicated function (validate_ancestor_ids) and adjusted imports accordingly.
  • Added the warn_on_extinct_columns function to provide users with warnings when columns in the alignment go extinct in all leaf sequences, suggesting the use of remove_extinct_columns.

Initial version of the ARPIP algo with homolgy path recpntruction
Integration test for PIP ASR using data from the ARPIP paper.
Checks that homology paths are reconstructed correctly.
Added function to transform Mapping to FixedBitSet to alignment mod
Fixed test scenarios for ARPIP as they were checking the wrong tree node ids
Refactored PIP ASR to separate struct to match e.g. `ParsimonyPresenceAbsence`
Refactor to have `AncestralSequenceReconstruction` trait impl first for PIP asr and
helpers second
Rename deletion_probabilities_at_nodes -> deletion_probabilities
Correct implementation of ARPIP indel inference with proper insertion probability computation
Added tests for when the most likely insertion location neitger at first node (MRCA) or root
Changed `Alignment::alphabet()` to return static ref to the alphabet to allow creating an MASA
from an MSA.
Split `validate_ids_with_ancestors` in two methods, `validate_taxa_ids` and `validate_ancestor_ids`
to allow validating them independently.
Added method to create an MASA from an Alignment and provided ancestral sequences in a `Sequence`
object (`from_msa_with_ancestral`).
Added a function `warn_on_extinct_columns` to reduce duplication for the check.
Use `warn_on_extinct_columns` in AncestralAlignment creation from all sequences
Added tests for `MASA::from_msa_with_ancestral` including correct creation and different
types of possible errors
Cleaned up the tests for pip asr intermediate (tmp) values
Use `from_msa_with_ancestral` in PIP ASR to reuse MSA sequences rather than recreate a
full `Sequences` object and use `from_aligned_with_ancestral_unchecked`.
Split `from_msa_with_ancestral` into two methods, `from_msa_with_ancestral` and
`from_msa_with_ancestral_unchecked`.
`from_msa_with_ancestral` has a default impl with the necessary checks and calls
`from_msa_with_ancestral_unchecked`.
Add doctests for `from_msa_with_ancestral` and `from_msa_with_ancestral_unchecked`
Rename `set_internal` to `compute_internal_probabilities` and `set_leaf` to
`compute_leaf_probabilities` to distinguish between the upward pass and the downward
pass through the tree.
Refactored internal character presence determination into separate method and
include presence in the tmp storage
Refactor to extract deciding deletion location (node or both children) to a separate method
@junniest
junniest requested review from MattesMrzik and a balanced review from Copilot November 4, 2026 13:56
@github-actions

github-actions Bot commented Nov 4, 2026 •

Copy link
Copy Markdown
Contributor

Test Results

1 184 tests  +55   1 183 ✅ +55   15m 35s ⏱️ - 12m 16s
    6 suites + 1       1 💤 ± 0 
    1 files   ± 0       0 ❌ ± 0 

Results for commit 5eaf7e0. ± Comparison against base commit 71d7ce2.

This pull request removes 8 and adds 63 tests. Note that renamed tests count towards both.
src/alignment/mod.rs - alignment ‑ surviving_columns_mask (line 167)
src/alignment/mod.rs - alignment::MASA ‑ from_aligned (line 480)
src/alignment/mod.rs - alignment::MASA ‑ from_aligned_with_ancestral_unchecked (line 588)
src/alignment/mod.rs - alignment::MASA ‑ remove_extinct_columns (line 662)
src/alignment/mod.rs - alignment::MSA ‑ alphabet (line 257)
src/alignment/mod.rs - alignment::MSA ‑ from_aligned_unchecked (line 349)
src/alignment/mod.rs - alignment::MSA ‑ len (line 285)
src/alignment/mod.rs - alignment::MSA ‑ seq_count (line 313)
alignment::tests ‑ masa_from_msa_and_ancestral
alignment::tests ‑ masa_from_msa_and_ancestral_unaligned_error
alignment::tests ‑ masa_from_msa_and_ancestral_wrong_length_error
alignment::tests ‑ masa_from_msa_with_ancestral_extra_ids
alignment::tests ‑ masa_from_msa_with_ancestral_missing_seqs
alignment::tests ‑ masa_from_msa_with_ancestral_wrong_alphabet
alignment::tests ‑ masa_from_msa_with_ancestral_wrong_seq_ids
alignment::tests::test_mapping_to_bits ‑ case_1
alignment::tests::test_mapping_to_bits ‑ case_2
alignment::tests::test_mapping_to_bits ‑ case_3
…

♻️ This comment has been updated with latest results.

Copilot AI left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🟡 Changes recommended

Unresolved critical defects can select incorrect ancestors, require quadratic memory, and admit invalid or ambiguous sequence IDs that distort alignment statistics and ancestral updates.

5 open issues
What changed in this PR

Adds PIP-based ancestral indel reconstruction to the existing ancestral sequence reconstruction interface, with alignment-construction APIs and ARPIP comparison tests.

Changes:

  • Introduces PIPIndels to infer ancestral character presence from insertion and deletion probabilities.
  • Adds checked and unchecked ancestral-alignment constructors, refactors ID validation, and standardizes static alphabet references.
  • Adds reconstruction tests and ARPIP reference data.
File Description
phyne/​tests/​pip_asr.rs Compares reconstructed ancestral mappings with ARPIP results.
phyne/​src/​pip_model/​tests.rs Adds a reconstruction integration test.
phyne/​src/​pip_model/​mod.rs Exports the ASR module and PIPIndels.
phyne/​src/​pip_model/​asr.rs Implements indel reconstruction and focused unit tests.
phyne/​src/​phylo_info/​phyloinfo_builder.rs Refactors leaf and ancestor ID validation with tests.
phyne/​src/​alignment/​tests.rs Tests mapping conversion and ancestral-alignment constructors.
phyne/​src/​alignment/​mod.rs Adds constructors and mapping utilities; refactors alphabet access and extinct-column warnings.
phyne/​data/​pip/​arpip/​sim-0_new.newick Supplies the ARPIP comparison tree.
phyne/​data/​pip/​arpip/​sim-0_msa_new.fasta Supplies the observed leaf alignment.
phyne/​data/​pip/​arpip/​sim-0_inferred.fasta Supplies reference ancestral and leaf sequences.

Review effort: Balanced

"ancestral sequences must have the same length as the provided alignment"
);
}
validate_ancestor_ids(&phylo.tree, &ancestral_seqs)?;

@MattesMrzik MattesMrzik Nov 5, 2026 •

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

copilot added

     phylo.tree.node_ids_are_unique()?;

This seems like a valid suggestion to me

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Like with the validate_taxa_ids it think this sequences object could contain seqs with ids that are in the tree tips set. they would not be detected by validate_ancestor_ids, in from_msa_with_ancestral they are silently ignored. So not quite like the MSA::from_aligned which added the additional seqs to the .seqs anyway rather than silently ignoring it.

Comment on lines +457 to +461
let tip_sequence_ids = HashSet::from_iter(sequences.into_iter().filter_map(|rec| {
let id = rec.id().to_string();
if !internal_ids.contains(&id) {
Some(id)
} else {

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is true:

  let tree = tree!("((A2:1,B3:1)I1:1,(C5:1,D6:1)I4:1)I0:0;");
    let sequences = Sequences::with_alphabet_unchecked(
        vec![
            record_wo_desc!("A2", b"XXX--X"),
            record_wo_desc!("B3", b"X-X--X"),
            record_wo_desc!("C5", b"X---X-"),
            record_wo_desc!("D6", b"--XXX-"),
            record_wo_desc!("I4", b"--XXX-"),
        ],
        Alphabet::dna(),
    );
    let msa = MSA::from_aligned(sequences, tree.clone()).unwrap();
    for seq in msa.msa.seqs() {
        println!("{}: {}", seq.id(), String::from_utf8_lossy(seq.seq()));
    }
    for map in msa.msa.leaf_maps().values() {
        println!("leaf map: {map:?}");
    }

prints

A2: XXXX
B3: XXX
C5: XX
D6: XXX
I4: XXX
leaf map: [Some(0), None, None, None, Some(1), None]
leaf map: [Some(0), Some(1), Some(2), None, None, Some(3)]
leaf map: [Some(0), None, Some(1), None, None, Some(2)]
leaf map: [None, None, Some(0), Some(1), Some(2), None]

@MattesMrzik MattesMrzik Nov 5, 2026 •

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Either add a comment that set sizes should still be check or go back to 3 separate methods:

  • one that checks seq ids = tree tip ids
  • one that checks seq ids = tree internal ids
  • one that check seqs ids = all tree ids

This way there might be some duplicate code, but its more clear what each method does and when to apply it and it would catch incorrect seqs set sizes automatically.

Comment thread phyne/src/pip_model/asr.rs Resolved
Comment on lines +209 to +212
tmp.f[idx] = f_child1.clone().component_mul(&tmp.f[children[1]]);
tmp.f[idx] *= 1.0 - tmp.deletion_probs[idx];
tmp.p[idx] =
f_child1.component_mul(&tmp.f[children[1]]) * tmp.insert_and_survive_probs[idx];
Comment thread phyne/src/alignment/mod.rs Outdated
Comment on lines +177 to +178
/// Tree is consumed and stored together with the ancestral alignment, because the alignment is only
/// meaningful in the context of the tree that it was built from.

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this comment could be reworked to fit better: Its not the tree that is consumed, but the PhyloInfo.

Suggested change
/// Tree is consumed and stored together with the ancestral alignment, because the alignment is only
/// meaningful in the context of the tree that it was built from.
/// PhyloInfo (i.e., Tree and MSA) is consumed, the ancestral alignment and corresponding tree are returned together,
///because the alignment is only meaningful in the context of the tree that it was built from.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, definitely. I think I copied the initial comment from from_aligned_with_ancestral_unchecked, which only gets a tree. I will rewrite, but now I am also thinking if I should flip the order of arguments, since from_aligned_with_ancestral_unchecked takes all_seqs and tree, whereas from_msa_with_ancestral_unchecked takes a PhyloInfo and ancestral_seqs.

Or, potentially, should we flip the arguments for from_aligned_with_ancestral_unchecked? It makes more sense to have tree+msa first and sequences second for the from_msa_... method, but then it clashes with the sequences+tree for the other one in my mind. It did not matter while there was only one method, but now it bothers me...

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i guess either
tree, seqs & phyloInfo, seqs
or
seqs, tree & seqs, PhyloInfo

Comment on lines +457 to +461
let tip_sequence_ids = HashSet::from_iter(sequences.into_iter().filter_map(|rec| {
let id = rec.id().to_string();
if !internal_ids.contains(&id) {
Some(id)
} else {

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is true:

  let tree = tree!("((A2:1,B3:1)I1:1,(C5:1,D6:1)I4:1)I0:0;");
    let sequences = Sequences::with_alphabet_unchecked(
        vec![
            record_wo_desc!("A2", b"XXX--X"),
            record_wo_desc!("B3", b"X-X--X"),
            record_wo_desc!("C5", b"X---X-"),
            record_wo_desc!("D6", b"--XXX-"),
            record_wo_desc!("I4", b"--XXX-"),
        ],
        Alphabet::dna(),
    );
    let msa = MSA::from_aligned(sequences, tree.clone()).unwrap();
    for seq in msa.msa.seqs() {
        println!("{}: {}", seq.id(), String::from_utf8_lossy(seq.seq()));
    }
    for map in msa.msa.leaf_maps().values() {
        println!("leaf map: {map:?}");
    }

prints

A2: XXXX
B3: XXX
C5: XX
D6: XXX
I4: XXX
leaf map: [Some(0), None, None, None, Some(1), None]
leaf map: [Some(0), Some(1), Some(2), None, None, Some(3)]
leaf map: [Some(0), None, Some(1), None, None, Some(2)]
leaf map: [None, None, Some(0), Some(1), Some(2), None]

"ancestral sequences must have the same length as the provided alignment"
);
}
validate_ancestor_ids(&phylo.tree, &ancestral_seqs)?;

@MattesMrzik MattesMrzik Nov 5, 2026 •

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

copilot added

     phylo.tree.node_ids_are_unique()?;

This seems like a valid suggestion to me

"ancestral sequences must have the same length as the provided alignment"
);
}
validate_ancestor_ids(&phylo.tree, &ancestral_seqs)?;

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Like with the validate_taxa_ids it think this sequences object could contain seqs with ids that are in the tree tips set. they would not be detected by validate_ancestor_ids, in from_msa_with_ancestral they are silently ignored. So not quite like the MSA::from_aligned which added the additional seqs to the .seqs anyway rather than silently ignoring it.

Comment on lines +217 to +218
/// Tree is consumed and stored together with the ancestral alignment, because the alignment is only
/// meaningful in the context of the tree that it was built from.

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same remark as the comment I added to: from_msa_with_ancestral

Comment on lines +850 to +852
let ancestral_seqs =
Sequences::with_alphabet_unchecked(ancestral_records, phylo.msa.alphabet())
.into_gapless();

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess in this method:

    fn from_msa_with_ancestral_unchecked(
        phylo: PhyloInfo<impl Alignment>,
        ancestral_seqs: Sequences,
    ) -> PhyloInfo<Self> {

when we check that we have a correct number of sequences and that they match the internal node ids (this is currently not the case as seqs with leaf ids can be in the passed ancestral_seqs), then we dont need to create this vec of ancestral_records in the for loop but can straight up take the passed ancestral_seqs and save it to MASA. we just need to assert the seqs alphabet is the same as the msa alphabet.

Comment thread phyne/tests/pip_asr.rs
// sim-0_msa_new.fasta as the input msa
// sim-0_new.newick as the input tree (changed node ids from "x" -> "i_x")
// sim-0_inferred.fasta is the result for comparison (concatenated alignment from input and ARPIP reconstruction)
// Note: this dataset has no insertions that are more likely to be higher than the MRCA.

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perhaps say that you have tests for the case that insertion is higher than MRCA, so users reading that here know that you also test for that case that is not covered here.

Comment on lines +29 to +30
// Deletion probabilities for each node, indexed by node (single value per branch).
deletion_probs: Vec<f64>,

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I dont know how its for you in vscode but for me its nice to have these as /// comments since then I can get the hover text on that variable (also at every place where its used) and read the comment without having to jump to this line.

}

#[derive(Debug, Clone, PartialEq)]
pub(crate) struct PIPAsrInfo<Q: QMatrix> {

@MattesMrzik MattesMrzik Nov 5, 2026 •

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perhaps link the paper where the diagram that you made is in

Comment on lines +204 to +209
// the potential values for f and p are precomputed but can be reset by deletion and insertion processing
let f_child1 = tmp.f[children[0]].clone();
tmp.f[idx] = f_child1.clone().component_mul(&tmp.f[children[1]]);
tmp.f[idx] *= 1.0 - tmp.deletion_probs[idx];
tmp.p[idx] =
f_child1.component_mul(&tmp.f[children[1]]) * tmp.insert_and_survive_probs[idx];

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

By the comment you mean that those values:

            insert_and_survive_probs: insert_and_survive_probabilities(model, info),
            deletion_probs: deletion_probabilities(model, info),

are precomputed?
Or you are doing the precomputing here and some time later these (f and p) values might get set to a different value?

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think yes, perhaps change comment to:

/// precompute potential values for f and p. They might be overwritten by deletion and insertion processing

so that it is clear that the precomputing happens here and was not "precomputing performed somewhere else and the place here is where they are overwritten".

let ancestral_seqs =
Sequences::with_alphabet_unchecked(ancestral_records, phylo.msa.alphabet());

AA::from_msa_with_ancestral(phylo, ancestral_seqs).unwrap()

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

perhaps also call the unchecked version here?

Comment thread phyne/tests/pip_asr.rs
Comment on lines +12 to +16
// Data from the ARPIP paper, javapip simulation with mu = 0.1
// sim-0_msa_new.fasta as the input msa
// sim-0_new.newick as the input tree (changed node ids from "x" -> "i_x")
// sim-0_inferred.fasta is the result for comparison (concatenated alignment from input and ARPIP reconstruction)
// Note: this dataset has no insertions that are more likely to be higher than the MRCA.

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

perhaps also add this to a readme in the data/pip/arpip/readme

}
}

#[cfg(test)]

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
#[cfg(test)]
#[cfg(test)]
#[cfg_attr(coverage, coverage(off))]

i guess this is also missing in:

  • src/parsimony/matrices/mod.rs:502
  • src/parsimony_presence_absence/mod.rs:127
  • src/alignment/sequences.rs:412
  • src/tkf_model/tkf91.rs:221
  • src/tkf_model/mod.rs:196 (tkf_numerical_tests)

phylo: PhyloInfo<impl Alignment>,
ancestral_seqs: Sequences,
) -> Result<PhyloInfo<Self>> {
if !ancestral_seqs.alphabet().eq(phylo.msa.alphabet()) {

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
if !ancestral_seqs.alphabet().eq(phylo.msa.alphabet()) {
if ancestral_seqs.alphabet() != phylo.msa.alphabet() {

-(-pip.mu() * node.blen).exp_m1() / pip.mu()
}) * common_factor
})
.collect::<Vec<f64>>()

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
.collect::<Vec<f64>>()
.collect()

#[test]
fn insertion_at_i5() {
let tree =
tree!("(((L3:1.0,L4:1.0)I2:1.0,(L6:0.001,L7:1.0)I5:10000.0)I1:1.0,L8:1.0)I0:0.0;");

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In this test and similar ones we dont need the branches to be that short and that long, more moderate lens are sufficient, but could also keep it like this.

Simplified sequence/node match checks with helper function `check_node_sequence_ids_match`.
Include validation that internal node ids are actually set (not empty) for `validate_ids_with_ancestors` and
`validate_ancestor_ids`
Added tests for `ensure_ancestral_ids_set` and `check_node_sequence_ids_match` helpers
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.

3 participants