Calculate the percentage of each unique phylogenetic tree in a BEAST output

Question

I have a nexus formatted BEAST output containing 20,000 phylogenetic trees of seven taxa. Is there any way to get the percentage of each unique phylogenetic tree contained in this output?

I already made an unsuccessful attempt with R.

Answer 1

I finally managed to do it in R. Here is my code:

install.packages('devtools')
library(devtools)
install_github('santiagosnchez/rBt')
library(rBt)

beast_output <- read.annot.beast('beast_output.trees')
beast_output_rooted <- root.multiPhylo(beast_output, c('taxon_A', 'taxon_B'))
unique_topologies <- unique.multiPhylo(beast_output_rooted)

count <- function(item, list) {
  total = 0
  for (i in 1:length(list)) {
    if (all.equal.phylo(item, list[[i]], use.edge.length = FALSE)) {
      total = total + 1
    }
  }
  return(total)
}

result <- data.frame(unique_topology = rep(0, length(unique_topologies)),
                     count = rep(0, length(unique_topologies)))
for (i in 1:length(unique_topologies)) {
  result[i, ] <- c(i, count(unique_topologies[[i]], beast_output_rooted))
}

result$percentage <- ((result$count/length(beast_output_rooted))*100)

Answer 2

You require a consensus phylogeny.

This is available in RAxML, PAUP, the old Phylip package and possibly the stats companion to Beast, Tracer. I'd need to dig up the code for RAxML and PAUP. The easiest is Consense here because you just download it and input the sequence. NOTE their is possibly a format difference (RAxML and Consense use Phylip format), I suspect you've got nexus, or Newick format.

Beast is likely to have its own consensus formating OR the Beast companion "Tracer" may have this (I can't remember). Tracer is here and would be cool if it did consensus phylogeny because there will be no conversion problems (it doesn't use Phylip format). The alternative would be to use PAUP (no conversion problems), which is currently in a state of flux from various stages of freeware and updates. PAUP is contained within the Genious package, which is available as a 1 month free-ware.

Have a look and let me know: because if you want to use RAxML or Consense you'll likely need to remind me of the Beast output, if it is what I think it is a format conversion app is needed. The calculation is easy, its easy to script a conversion program, but I think tree formating will be an issue.

You could use Ape, within R, but as you have gathered there are alot of other solutions. There are more recent, impressive, solutions as stand-alones which I've not used yet.