Skip to content
Snippets Groups Projects
Commit fc6ef36c authored by iago-lito's avatar iago-lito
Browse files

Trim branches lengths to not count root lengths.

parent 6ecf34b0
Branches
Tags
No related merge requests found
......@@ -36,20 +36,65 @@ impl<'i> GeneTree<'i> {
// Iterate over terminal nodes names
// and their corresponding branches length from root.
pub(crate) fn branches_lengths(
&self,
) -> impl Iterator<Item = (usize, (&'_ str, BranchLength))> {
// The length is "trimmed" because root branch length is *not* taken into account.
fn trimmed_branches_lengths(&self) -> impl Iterator<Item = (usize, (&'_ str, BranchLength))> {
// Accumulate "distance" with branches lengths from root to tips.
self.tree.root_to_leaves(
0.,
|distance, (_, inode)| distance + inode.branch,
|distance, (_, tnode)| (tnode.name, distance + tnode.branch),
|distance, (i, inode)| {
let is_root = inode.parent == i;
if is_root {
distance
} else {
distance + inode.branch
}
},
|distance, (i, tnode)| {
let is_root = tnode.parent == i;
let total = if is_root { distance } else { distance + tnode.branch };
(tnode.name, total)
},
)
}
pub fn mean_branch_length(&self) -> BranchLength {
self.branches_lengths()
self.trimmed_branches_lengths()
.map(|(_, (_, length))| length)
.summed_mean::<usize>()
}
}
#[cfg(test)]
mod tests {
use crate::lexer::Lexer;
#[test]
fn trimmed_lengths() {
// Ignore the root branch length.
let input = "(A:1,(B:2,(((C:3,D:4):5,E:6):7,F:8):9):10):11; 12";
//
// (branches lengths) (nodes indexes)
// 1
// ┌───────────── A ┌────────────1 A
// 11│ 2 │
// ──┤ ┌──────── B ─0 ┌───────3 B
// │ │ 3 │ │
// │ 10 │ 5┌── C │ │ ┌─7 C
// └────┤ 7┌─┤ 4 └────2 ┌─6
// │ ┌─┤ └── D │ ┌─5 └─8 D
// │ │ │ 6 │ │ │
// │9│ └──── E │ │ └───9 E
// └─┤ 8 └─4
// └────── F └────10 F
//
let tree = Lexer::new(input).read_gene_tree().unwrap();
let mut it = tree.trimmed_branches_lengths();
assert_eq!(it.next(), Some((1, ("A", f64::from(1)))));
assert_eq!(it.next(), Some((3, ("B", f64::from(10 + 2)))));
assert_eq!(it.next(), Some((7, ("C", f64::from(10 + 9 + 7 + 5 + 3)))));
assert_eq!(it.next(), Some((8, ("D", f64::from(10 + 9 + 7 + 5 + 4)))));
assert_eq!(it.next(), Some((9, ("E", f64::from(10 + 9 + 7 + 6)))));
assert_eq!(it.next(), Some((10, ("F", f64::from(10 + 9 + 8)))));
assert_eq!(it.next(), None);
}
}
......@@ -38,7 +38,7 @@ pub(crate) enum Node<'i, B> {
#[cfg_attr(test, derive(Debug, PartialEq))]
pub(crate) struct InternalNode<'i, B> {
parent: usize,
pub(crate) parent: usize,
children: [usize; 2],
pub(crate) branch: B,
name: Option<&'i str>,
......@@ -46,7 +46,7 @@ pub(crate) struct InternalNode<'i, B> {
#[cfg_attr(test, derive(Debug, PartialEq))]
pub struct TerminalNode<'i, B> {
parent: usize,
pub(crate) parent: usize,
pub(crate) branch: B,
pub name: &'i str,
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment