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.

1 comment:

  1. Dear Liam, how are you?
    I have obtained a boring error message even having named the rows of my matrix of coordinates. The message is as follows:
    Error in res[edge[i, 2]] <- res[edge[i, 1]] + el[i] :
    replacement has length zero
    My coordinates are as you recommend: latitude first and longitude after. Could you help me? My database is South America in phylo.to.map function.
    Here part of my dataset and the phylogeny:
    lat long
    E_L_machadus -29.4500 -50.5833
    E_L_machadus -22.6333 -44.5833
    E_L_machadus -22.7333 -45.5833
    E_L_longicornis -13.9833 -56.4500
    E_L_longicornis -4.2917 -44.7917
    E_L_longicornis -23.0333 -45.5500
    E_L_imitator -28.5500 -56.0500
    E_L_imitator -34.5667 -58.4667
    E_L_imitator -34.5875 -58.6725
    E_L_imitator -32.0000 -60.2500
    E_L_imitator -33.0103 -58.5142
    E_L_imitator -34.5833 -58.4167
    E_L_imitator -32.7655 -52.6040
    E_L_imitator -29.3500 -49.7333
    E_L_imitator -29.9756 -50.1281
    E_L_imitator -27.6333 -52.2833
    E_L_illotus -24.3333 -67.0167
    E_L_illotus -23.1333 -64.3333
    E_L_illotus -23.2101 -64.0972
    E_L_illotus -26.8833 -65.4166
    E_L_illotus -24.3000 -65.9500
    E_L_illotus -23.3000 -51.3667
    E_L_illotus -22.9000 -47.0833
    E_L_illotus -22.8667 -48.4333
    E_L_illotus -22.8177 -47.0683
    E_L_illotus -23.1958 -51.1997
    E_L_illotus -22.4833 -48.5333
    E_L_illotus -22.7403 -47.3344
    E_L_illotus -22.1833 -51.2333
    E_L_illotus -22.7167 -47.6333
    E_L_illotus -20.8667 -51.4833
    E_L_illotus -23.9658 -45.1367
    E_L_illotus -21.1667 -47.8000
    E_L_illotus -24.6500 -56.4333
    E_L_illotus -17.3667 -63.3000
    E_L_illotus -22.1667 -64.7000
    E_L_sharpi -31.7000 -65.1833
    E_L_sharpi -34.9000 -54.9500
    E_L_sharpi -23.6500 -46.7000
    E_L_sharpi -32.7655 -52.6040
    E_L_sharpi -31.7667 -52.3333
    E_L_sharpi -30.3833 -56.4500

    #######
    tree:
    tree = (ROOT,(E_L_longicornis,(((E_L_picticornis,E_L_triangulator),((E_L_cornutus,(E_L_machadus,E_L_riograndensis)),(E_L_illotus,(E_L_imitator,(E_L_circumfusus,E_L_sharpi))))),(E_L_aceratos,E_L_cribrarius))));

    My best regards, José.

    ReplyDelete

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