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)
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)
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
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.
Dear Liam, how are you?
ReplyDeleteI 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é.