Monday, April 7, 2014

New function to paint individual branch or branches of the tree

A few days ago I received the following email:

"I have recently come across your excellent paper entitled "On the analysis of evolutionary change along single branches in a phylogeny". In that you employed a software "available on request", so I am writing to request permission to obtain and use your software. Alternatively, I wonder whether the functions from that software were incorporated into the Phytools R package. I am sorry if I am missing the obvious and it is indeed in the Phytools package - I am not familiar with all functions in that package yet."

(One interesting attribute of this message is that it harks back to a day in the not too distant past in which sharing data or software 'on request' seemed more than reasonable - how quaint & passé that seems today!)

The paper that's being referred to here is Revell (2008; Am. Nat.) in which I show that a (at that time) new likelihood method by O'Meara et al. (2006) could be used to test hypotheses about exceptional evolution along (one or more) single branch(es) of the tree. Previously McPeek (1995) had proposed a cruder approach based loosely on contrasts, which I showed had reasonable type I error; however I also showed that the likelihood method was probably preferable under most circumstances.

Indeed this method is available in phytools in the form of brownie.lite; although to 'paint' branches which we have a priori hypothesized to have a higher or lower evolutionary rate is a little bit tricky. I decided to make this a whole lot easier by adding the function paintBranches to the phytools package (new source build here).

Here is an example workflow using the function to test the hypothesis of a higher rate of evolution along a specific branch of the tree:

> library(phytools)
Loading required package: ape
Loading required package: maps
Loading required package: rgl
> packageVersion("phytools")
[1] ‘0.3.99’
> ## simulate data with a high rate on one branch
> tree<-pbtree(n=26,tip.label=LETTERS[26:1])
> plotTree(tree)
> nodelabels()
> ## stretch one of the branches by x 100
> tt<-tree
> tt$edge.length[which(tt$edge[,2]==29)]<-
 tt$edge.length[which(tt$edge[,2]==29)]*100
> ## simulate data on distorted tree
> x<-fastBM(tt)
> ## paint branch with second regime to fit 2-rate model
> tree<-paintBranches(tree,29,"2")
> plotSimmap(tree,lwd=4)
no colors provided. using the following legend:
      1       2
"black"   "red"
> nodelabels()
> fit<-brownie.lite(tree,x)
> fit
ML single-rate model:
        s^2     se      a       k  logL
value   1.41    0.39    1.35    2  -39.37 

ML multi-rate model:
        s^2(1)  se(1)   s^2(2)  se(2)   a       k  logL 
value   0.96    0.26    40.24   60.27   0.12    3  -35.82

P-value (based on X^2): 0.01

R thinks it has found the ML solution.

We can use the same approach if our hypothesized exceptional regime is found on several branches in the tree. In this case, the function argument edge is an integer vector instead of a single node number. Finally, if we have an a priori hypothesis of multiple classes of exceptional rate, we can test that too using multiple calls to paintBranches.

That's it.

2 comments:

  1. Thanks a lot for that! It is simpler than I thought. Ideal to test predictions about selection regime shifts, especially if you have a set of rates of nucleotide substitutions for the same phylogeny to compare this to, to discard the possibility that the identified shift in rate was (at least partly) a by-product of different generation times or different DNA copy fidelity/repair mechanisms, as these affect synonymous substitution rates too. Thanks again.

    ReplyDelete
  2. Dear Liam,

    Is there any way to test models and to reconstruct ancestral states based on "painted tree" for discret traits?

    Thanks in advance.
    Dimitri

    ReplyDelete