Sunday, March 10, 2019

Update to "phylo.to.map" plot method to prevent the object from being optimized & plotted for different directions

I recently added a new feature to the phytools function phylo.to.map such that more than one locality could be linked to each tip of the plotted tree. When I shared this on Facebook

it was quickly pointed out that the node rotations could be easily improved to result in fewer linking lines crossing. E.g.:

Now, the function already does an optimization of precisely this - but as I realized later, my example had been optimized as though the tree was to be plotted above the map facing downwards, even though I ended up plotting the tree to the left of the tree facing rightwards. This is because the function is written in a way such that the generation & optimization of the "phylo.to.map" is done as a separate step to the plotting.

To try to prevent phytools users from running into precisely this same issue I have now added an additional element, direction, to the "phylo.to.map" object. The plot method now checks if the object direction matches the direction of the plot, and, if it does not, re-optimizes the node rotations to match the plot direction. We can still plot a tree without optimizing the node rotation as direction can also take the value "unoptimized" in which case nothing will be done when the user goes to plot the object.

I just pushed this update to GitHub so we can already see how it works.

library(phytools)
library(mapdata)
library(viridis)
## our tree & data:
Colombia.tree
## 
## Phylogenetic tree with 12 tips and 11 internal nodes.
## 
## Tip labels:
##  G.cnripkth, H.vlyc, T.zhloa, B.gotufrp, S.iaoz, B.lhdkxofue, ...
## 
## Rooted; includes branch lengths.
Colombia
##                     lat      long
## B.gotufrp     6.8645111 -73.51948
## B.gotufrp     7.5847074 -73.04759
## B.lhdkxofue   3.8413048 -68.70622
## E.rxaspgo     1.5739001 -76.26423
## G.cnripkth   10.2457286 -73.83771
## G.cnripkth   10.3915751 -74.01978
## H.vlyc        8.8932416 -75.19020
## H.vlyc        8.7695991 -76.19651
## H.vlyc       10.0501305 -75.31938
## K.paljn       1.3750050 -71.72943
## L.wyue        4.8755596 -75.11064
## L.wyue        5.7093001 -74.01726
## O.cwzqnrmf   -1.2504110 -70.53605
## S.iaoz        5.5915821 -70.37694
## S.iaoz        5.7544784 -69.51145
## S.iaoz        4.6409290 -69.83929
## T.lzbgvsaix   3.1650613 -75.94600
## T.lzbgvsaix   4.1932495 -76.09112
## T.zhloa       7.3020804 -76.46313
## X.ybcluhasfx -0.5741675 -72.08744

We can already see that these data also contain repeating species. First, let's plot the data without optimizing the node rotation at all:

cols<-setNames(sample(viridis(n=Ntip(Colombia.tree))),
    Colombia.tree$tip.label)
obj<-phylo.to.map(Colombia.tree,Colombia,rotate=FALSE,database="worldHires",
    regions="Colombia",plot=FALSE)
obj
## Object of class "phylo.to.map" containing:
## 
## (1) A phylogenetic tree with 12 tips and 11 internal nodes.
## 
## (2) A geographic map with range:
##      -4.22N, 13.38N
##      -81.72W, -66.87W.
## 
## (3) A table containing 20 geographic coordinates (may include more than one set per species).
plot(obj,direction="rightwards",colors=cols,ftype="off",cex.points=c(0,1.5),
    pts=FALSE)

plot of chunk unnamed-chunk-2

Obviously, this is pretty much garbage. Now, I'm going to try to repeat the error I made the other day. First, I'll optimize as if the tree is to be plotted downwards & then go ahead & plot it downwards:

obj<-phylo.to.map(Colombia.tree,Colombia,rotate=TRUE,database="worldHires",
    regions="Colombia",plot=FALSE)
## objective: 44
## objective: 44
## objective: 42
## objective: 40
## objective: 38
## objective: 30
## objective: 30
## objective: 28
## objective: 26
## objective: 26
## objective: 24
obj
## Object of class "phylo.to.map" containing:
## 
## (1) A phylogenetic tree with 12 tips and 11 internal nodes.
## 
## (2) A geographic map with range:
##      -4.22N, 13.38N
##      -81.72W, -66.87W.
## 
## (3) A table containing 20 geographic coordinates (may include more than one set per species).
plot(obj,colors=cols,ftype="off",cex.points=c(0,1.5),pts=FALSE)

plot of chunk unnamed-chunk-3

Now I realize that I want to instead have the phylogeny to the left of the map & plot it facing rightwards, but I try to do that without re-optimizing the object:

plot(obj,direction="rightwards",colors=cols,ftype="off",cex.points=c(0,1.5),
    pts=FALSE)
## "phylo.to.map" direction is "downwards" but plot direction has been given as "rightwards".
## Re-optimizing object....
## objective: 50
## objective: 18
## objective: 18
## objective: 18
## objective: 6
## objective: 2
## objective: 2
## objective: 2
## objective: 2
## objective: 2
## objective: 2

plot of chunk unnamed-chunk-4

Neat.

We can see that with the update when I try to plot a "phylo.to.map" object that has been optimized in the incorrect direction the plot method essentially re-runs the optimization to try to obtain a better node rotation for the plotting direction. Note that the original object is unchanged so it would be more efficient to optimize it in the planned plotting direction in the first place!

That's all for now.

No comments:

Post a Comment

Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.