In a comments on a recent post,
Peter Beerli keenly
noticed that the optimization of the alignment of a phylogeny to a plotted
map produced by phytools::phylo.to.map was not very good & could
be improved considerably by eye:
At first, I assumed that this was just because phylo.to.map uses
a very simple & greedy
optimization routine. This is true, but it actually turns out that I was missing something
else as well. I forgot then when making optimizing the "phylo.to.map" object,
not just when plotting it, I needed to specify a plot direction for the tree. That's
because the optimal node rotation for a downward-facing tree (which is assumed by
default) can obviously be quite different than a rightward-facing tree.
Here's what I mean, using similar data to last time:
library(mapdata)
library(viridis)
Argentina.tree
##
## Phylogenetic tree with 12 tips and 11 internal nodes.
##
## Tip labels:
## T.plfuvkcyo, O.wcsuker, T.pchgeav, W.wjirnahuvq, S.dzfigk, R.ojnxqzgdi, ...
##
## Rooted; includes branch lengths.
Argentina
## lat long
## T.plfuvkcyo -24.42307 -65.97496
## O.wcsuker -25.35801 -61.96024
## T.pchgeav -28.65778 -58.49548
## W.wjirnahuvq -28.71277 -66.08495
## S.dzfigk -31.13261 -61.30028
## R.ojnxqzgdi -33.16747 -65.53499
## H.abxnr -35.64229 -67.89983
## G.twnhocpf -37.01720 -59.43041
## I.hfga -38.88707 -67.12988
## V.zqnkawixo -41.19691 -68.61478
## B.zompicgsvl -43.45175 -67.56985
## I.lcznqiy -47.52147 -69.76970
## T.plfuvkcyo -25.58899 -65.80964
## T.plfuvkcyo -25.71640 -65.91033
## O.wcsuker -24.52574 -62.78101
## O.wcsuker -25.17543 -60.74765
## W.wjirnahuvq -28.38156 -66.87795
## S.dzfigk -30.26453 -60.03134
## S.dzfigk -31.27273 -60.83937
## R.ojnxqzgdi -34.26012 -63.75400
## H.abxnr -35.54303 -67.84130
## H.abxnr -36.19707 -68.24773
## I.hfga -39.63231 -68.64393
## V.zqnkawixo -41.89663 -68.27408
## B.zompicgsvl -42.76742 -68.35690
## I.lcznqiy -46.06199 -69.94464
## first UNROTATED
obj<-phylo.to.map(Argentina.tree,Argentina,database="worldHires",
regions="Argentina",plot=FALSE,rotate=FALSE)
cols<-setNames(sample(viridis(n=Ntip(Argentina.tree))),
Argentina.tree$tip.label)
plot(obj,direction="rightwards",colors=cols,ftype="off",cex.points=c(0,1.5),
pts=FALSE)
## next ROTATED USING DEFAULT BUT PLOTTED 'rightwards'
obj<-phylo.to.map(Argentina.tree,Argentina,database="worldHires",
regions="Argentina",plot=FALSE)
## objective: 32
## objective: 30
## objective: 30
## objective: 28
## objective: 28
## objective: 26
## objective: 26
## objective: 26
## objective: 24
## objective: 22
## objective: 20
cols<-setNames(sample(viridis(n=Ntip(Argentina.tree))),
Argentina.tree$tip.label)
plot(obj,direction="rightwards",colors=cols,ftype="off",cex.points=c(0,1.5),
pts=FALSE)
## finally ROTATED 'rightwards' AND PLOTTED 'rightwards'
obj<-phylo.to.map(Argentina.tree,Argentina,database="worldHires",
regions="Argentina",plot=FALSE,direction="rightwards")
## objective: 40
## objective: 36
## objective: 34
## objective: 32
## objective: 26
## objective: 20
## objective: 14
## objective: 6
## objective: 4
## objective: 2
## objective: 2
cols<-setNames(sample(viridis(n=Ntip(Argentina.tree))),
Argentina.tree$tip.label)
plot(obj,direction="rightwards",colors=cols,ftype="off",cex.points=c(0,1.5),
pts=FALSE)
What a difference!
Hi Liam, I'm having problems with the colors portion of the plot function. My code is identical to yours, but I am getting the following traceback on an error that colors[cw$tip.label, ] has the incorrect number of dimensions
ReplyDelete5.
plot.xy(xy.coords(x, y), type = type, ...)
4.
lines.default(c(X[which(cw$edge[, 2] == i), 2] + if (from.tip) 0 else sh[i], coords[i, 1]), c(y[i], coords[i, 2]), col = colors[cw$tip.label, ][i, 1], lty = lty, lwd = lwd[2])
3.
lines(c(X[which(cw$edge[, 2] == i), 2] + if (from.tip) 0 else sh[i], coords[i, 1]), c(y[i], coords[i, 2]), col = colors[cw$tip.label, ][i, 1], lty = lty, lwd = lwd[2])
2.
plot.phylo.to.map(obj, direction = "rightwards", colors = cols, ftype = "off", cex.points = c(0, 0.5), pts = FALSE)
1.
plot(obj, direction = "rightwards", colors = cols, ftype = "off", cex.points = c(0, 0.5), pts = FALSE)
Any idea what might be going on? This is using phytools 0.6.72.
Cheers,
Mark
I should mention that without the colors argument the function plots beautifully.
DeleteHi Mr. Liam Revell, the same error occurs when I ran with my data. In addition, I was able to include only one record per terminal, even creating a matrix with some duplicated row names and the first and second colums as 'lat' and 'lon' coordinates. The "Objetive:" warning message load the data, but when I plot the phylo.to.map object, only one record is showed.
DeleteI got the error "Error in colors[cw$tip.label, ] : subscript out of bounds" when plotting the map with rotate=TRUE (but it worked fine with rotate=FALSE). I was able to solve it by replacing spaces with underscores in the tip labels and rownames of the coordinate matrix. After some sleuthing it looks like the minRotate function converts spaces to underscores in the tip names, which then causes a mismatch when subsetting colors (named with spaces) by the tip labels (now named with underscores).
DeleteNevermind, it was not minRotate that converted them to underscores. Must be something else.
DeleteHi, do you know where I could get the maps and trees datasets?
ReplyDeleteRegards
These are simulated data - so exactly the same data can probably not be obtained. The code I used to simulate the data is on my blog here.
DeleteHi, thank you for an amazing package and tutorial - this is exactly what I was looking for. Do you know if there is a way to only plot a subregion of a worldHires map? Fore example I am plotting the UK but only want the subregion "Great Britain"? Thank you in advance.
ReplyDelete