5
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ If the phylogenetic trees can be represented as strings then it is basically just string matching. Anyway, I think your question would be more suited for Bioinformatics as it is mostly about the programming aspect. $\endgroup$
    – WYSIWYG
    Feb 5 '19 at 10:36
  • $\begingroup$ Could you explain what is a "unique phylogenetic tree"? Perhaps if you upload an image it would be clearer what do you mean. $\endgroup$
    – llrs
    Feb 5 '19 at 15:51
  • $\begingroup$ @llrs When I mention a 'unique phylogenetic tree', I am referring to the topology of the tree. I would like to extract all the unique topologies from my BEAST output, and calculate how many times each of these topologies appears in my BEAST output. This would allow me to establish a percentage for each topology. The idea is to know which topologies are the most represented in addition to the consensus tree. $\endgroup$
    – jvddorpe
    Feb 6 '19 at 15:45
4
$\begingroup$

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)
$\endgroup$
1
  • 1
    $\begingroup$ This is fantastics, thanks for posting back $\endgroup$
    – M__
    Feb 11 '19 at 13:46
2
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ Many thanks for your answer! However, I already have a consensus tree and I am now trying to get the other most represented topologies from my BEAST output. $\endgroup$
    – jvddorpe
    Feb 6 '19 at 15:49
  • 1
    $\begingroup$ Ok. As far as I'm aware that program has not been written, at least a couple of years ago it wasn't, not sure about this now. Using R "ape" is your best bet. On this question I definitely can't help. $\endgroup$
    – M__
    Feb 6 '19 at 16:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.