Saturday, June 29, 2013

phytools page down again - new version of phytools submitted to CRAN

The phytools page is having problems again. This is a UMass-Boston server-wide issue, so there is not much I can do about it unfortunately (aside from complain to the system administrator). To make all the recent updates to the phytools package available to all users, however, I have just submitted a new version of phytools to CRAN. Hopefully this version is accepted and made available via the CRAN mirrors as well as in binary form soon.

Some updates in this version over the previous CRAN release (phytools 0.2-80) include the following:

1. A new function to add a legend to a plotted stochastic map tree (add.simmap.legend).

2. A new version of plotSimmap that sets the environmental variable "lastplot.phylo" to be compatible with ape functions nodelabels, tiplabels, and edgelabels.

3. A new version of findMRCA that can also (optionally) return the height above the root of the MRCA of a pair of taxa.

4. A robust Newick tree reader that can handle singleton nodes, if present.

5. A bug fix in brownieREML.

6. A new function (collapse.to.star) that collapses specified subtrees into a star tree, while retaining the same total height above the root for each tip.

7. New versions of the phytools 3D functions phylomorphospace3d and fancyTree(...,type="traitgram3d") that work a little better with the functions of rgl (described here).

8. User control of tip & internal node sizes in phylomorphospace..

9. A hacky work-around fix for the problem of attaching a new tip separated from an existing tip by zero lenth using bind.tip.

10. Finally, today, a fix to the broken function evol.vcv, which is the function implementing the method of Revell & Collar (2009).

Hopefully this new phytools update is accepted and posted to CRAN soon.

Friday, June 28, 2013

Upper triangle of a matrix to a vector by row

Say we want to get a square matrix into a vector. We can do:

to.upper<-function(X) X[upper.tri(X,diag=TRUE)]
but this will give us our upper triangular matrix as a vector by column, i.e.:
> X
       [,1]   [,2]   [,3]   [,4]   [,5]
[1,] -0.329 -0.825 -0.906  0.834  0.702
[2,]  1.195 -0.237 -1.132 -0.827  0.150
[3,] -0.097 -1.270 -1.397  0.450  0.791
[4,]  0.883 -0.574  1.538 -2.632  0.135
[5,]  1.993 -0.520 -0.071  0.094 -0.064
> to.upper(X)
[1] -0.329 -0.825 -0.237 -0.906 -1.132 -1.397  0.834 -0.827  0.450 -2.632  0.702  0.150  0.791  0.135 -0.064

If we want to get our upper triangular matrix as a vector by row (as I did), we can use:

to.upper<-function(X) t(X)[lower.tri(X,diag=TRUE)]
which works just as we'd hoped, i.e.:
> X
       [,1]   [,2]   [,3]   [,4]   [,5]
[1,] -0.329 -0.825 -0.906  0.834  0.702
[2,]  1.195 -0.237 -1.132 -0.827  0.150
[3,] -0.097 -1.270 -1.397  0.450  0.791
[4,]  0.883 -0.574  1.538 -2.632  0.135
[5,]  1.993 -0.520 -0.071  0.094 -0.064
> to.upper(X)
[1] -0.329 -0.825 -0.906  0.834  0.702 -0.237 -1.132 -0.827  0.150 -1.397  0.450  0.791 -2.632  0.135 -0.064

(Note that if we have a symmetric matrix than the upper triangular matrix by row & the lower triangular matrix by column are the same. Here I used a non-symmetric square matrix so that we could tell these apart.)

phytools page temporarily down

The phytools page is temporarily down. This seems to be a server-wide issue and is thus not my fault - hopefully this will be fixed in the morning.

Thursday, June 27, 2013

Fix for bind.tip

At this year's Evolution meeting it was reported to me that bind.tip has some weird behavior if you try and attach a tip to the end of a terminal edge. One might want to do this if, say, identical haplotypes were removed for inference - and then one wanted to add them back in for plotting.

Here's a demo of the problem:

> tree<-pbtree(n=20)
> plotTree(tree,setEnv=TRUE,offset=0.7)
setEnv=TRUE is experimental. please be patient with bugs
> tiplabels()
> aa<-bind.tip(tree,"t21",where=1,position=0)
> plotTree(aa,setEnv=TRUE,offset=0.7)
setEnv=TRUE is experimental. please be patient with bugs
> tiplabels()

Holy cow! Weird, right? Instead of binding our new tip zero-length below the tip node, bind.tip is replacing the tip label of our target terminus! Note that this will not happen even if we add the the new tip just below the target tip.

In this case, bind.tip actually just wraps around the ape function bind.tree which attaches two trees together.

I have built a work-around into the latest version of bind.tip (source code now here). This is also in a new build of phytools (phytools 0.2-89). Let's verify that this works:

> require(phytools)
Loading required package: phytools
> packageVersion("phytools")
[1] ‘0.2.89’
> bb<-bind.tip(tree,"t21",where=1,position=0)
> plotTree(bb,setEnv=TRUE,offset=0.7)
setEnv=TRUE is experimental. please be patient with bugs
> tiplabels()
> bb$edge.length
[1] 0.49852021 0.90133001 0.63759396 0.36295139 0.00000000 0.00000000 ....

This fix is a bit of a hack - but at least it works.

Wednesday, June 26, 2013

User control of node sizes in phylomorphospace

A phytools user requests:

Is there a way to turn off the internal nodes in the phylomorphospace function in phytools?

I decided to interpret this more generally as a request for user control over the size of the plotted points at terminal and internal nodes in a phylomorphospace (that is, a bivariate projection of the tree into morphospace) plot.

This was easy enough to add and is now an optional argument in the phylomorphospace function. Point size is controlled using the cex argument in points - which adjusts the plotted point size relative to the default (1.0). Code for the new version is here and I also posted a new phytools build (phytools 0.2-88) containing this update.

Here's a demo, using the original request above:

> require(phytools)
Loading required package: phytools
> packageVersion("phytools")
[1] ‘0.2.88’
> tree<-pbtree(n=30)
> X<-fastBM(tree,nsim=2)
> layout(matrix(c(1,2),1,2))
> plotTree(tree,mar=c(3.2,0.1,1.4,0.1))
> par(mar=c(4.1,4.1,2.1,0.5))
> phylomorphospace(tree,X,node.size=c(0,1.3),lwd=2, xlab="x",ylab="y",xlim=c(-4.2,2.2))
(Full resolution.)

Monday, June 24, 2013

New versions of phytools 3D methods

In my Evolution 2013 meeting talk (which I will put online soon), I showed a couple of old 3D visualization methods: a three-dimensional projection of the tree into morphospace ("phylomorphospace"); and a 3D traitgram visualization method in which x & y are phenotypic trait axes, while z represents the time since the root.

To make the plots I made a few minor updates to the two functions in which these methods are implemented: phylomorphospace3d and fancyTree(...,type="traitgram3d"). These updates are in a new phytools version (phytools 0.2-87). Here's a quick demo of how to create these as .gifs:

> require(phytools)
Loading required package: phytools
> packageVersion("phytools")
[1] ‘0.2.87’
> # simulate tree & data
> tree<-pbtree(n=40)
> X<-fastBM(tree,nsim=3)
> # create 3D phylomorphospace
> xx<-phylomorphospace3d(tree,X)
> movie3d(xx,duration=20,dir=".",movie= "phylomorphospace3d")
Will create: ./phylomorphospace3d.gif
Executing: convert -delay 1x10 phylomorphospace3d*.png phylomorphospace3d.gif
Deleting frames.
> # create 3D traitgram
> xx<-fancyTree(tree,type="traitgram3d",X=X[,1:2])
> movie3d(xx,duration=20,dir=".",movie="traitgram3d")
Will create: ./traitgram3d.gif
Executing: convert -delay 1x10 traitgram3d*.png traitgram3d.gif
Deleting frames.

That's it for now.

Sunday, June 23, 2013

Function to collapse some subtrees into a star - while retaining the same total height of the tips

The following was recently requested via the R-sig-phylo mail list serve:

I'd like to collapse the descendants of a node, identified using something like node <- mrca(tree)["A", "B"]. I did not see a function in ape, geiger, phyloch, or picante to do something like collapse.descendants(node). Is there a package with a function like this?

This can be pretty easily done using the functions of phytools. Here's a little function to do this:

function(tree,node){
 tt<-splitTree(tree,split=list(node=node,bp=
  tree$edge.length[which(tree$edge[,2]==node)]))
 ss<-starTree(species=tt[[2]]$tip.label,branch.lengths=
  diag(vcv(tt[[2]])))
 ss$root.edge<-0
 tree<-paste.tree(tt[[1]],ss)
 return(tree)
}

And here is a quick demo:

> tree<-pbtree(n=50,scale=1)
> plotTree(tree,node.numbers=T,fsize=0.8)
> tree<-collapse.to.star(tree,fastMRCA(tree,"t7","t12"))
> tree<-collapse.to.star(tree,fastMRCA(tree,"t50","t18"))
> plotTree(tree,node.numbers=T,fsize=0.8)

Cool. That did exactly what we wanted tit to do. Note that every time we collapse a subtree, the node numbers of the tree will change - so we cannot use the node numbers from the original tree to collapse multiple subtrees (we need to recompute the target node each time).