This is an archived static version of the original phylobabble.org discussion site.

Balanced minimum evolution computation in R

ematsen

Hello friends–

I’m wondering if there is some way to calculate the BME score in an existing R package. By this score, I mean l in the sense of

ncbi.nlm.nih.gov

Neighbor-joining revealed.

O Gascuel and M Steel, Molecular biology and evolution, Nov 2006

It is nearly 20 years since the landmark paper (Saitou and Nei 1987) in Molecular Biology and Evolution introducing Neighbor-Joining (NJ). The method has become the most widely used method for building phylogenetic trees from distances, and the original paper has been cited about 13,000 times (Science Citation Index). Yet the question "what does the NJ method seek to do?" has until recently proved somewhat elusive, leading to some imprecise claims and misunderstanding. However, a rigorous answer to this question has recently been provided by further mathematical investigation, and the purpose of this note is to highlight these results and their significance for interpreting NJ. The origins of this story lie in a paper by Pauplin (2000) though its continuation has unfolded in more mathematically inclined literature. Our aim here is to make these findings more widely accessible.

There are heaps of functions to infer trees using distance matrices, but here we have a distance matrix and a tree and want to get a number.

@klaus_schliep?

Thanks!

Erick

klaus_schliep

Hi @ematsen, if your tree is binary this is trivial, as all the necessary functions exist:

library(ape)
library(phangorn)
library(magrittr)

# a little function for you
get_l <- function(tree, dm){
  # make sure the w and d will be in the same order
  label <- tree$tip.label  
  dm <- as.matrix(dm)[label, label]
  dm <- dm[lower.tri(dm)]
  tmp <- designTree(tree)
  w <- 2^(rowSums(tmp)-2)
  sum(dm / w)
}

# a little example
data(yeast)
dm <- dist.ml(yeast)

fastme_tree <- fastme.bal(dm)
get_l(fastme_tree, dm)

# compare with brute force
trees <- allTrees(8, tip.label = names(yeast)) # compute all possible trees
trees %<>% .uncompressTipLabel() %>% unclass   # this speeds up things a bit

bf_l <- sapply(trees, get_l, dm=dm)
bf_tree <- trees[[which.min(bf_l)]]
RF.dist(bf_tree, fastme_tree)

I have no code for the general case with multifurcations (yet). Hope this gets you started.

Cheers, Klaus

ematsen

Very slick! Thanks @klaus_schliep! This will do just fine for us.