Friday, March 8, 2019

A follow-up comment on phylo.to.map for making & plotting rightward-facing trees

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)

plot of chunk unnamed-chunk-1

## 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)

plot of chunk unnamed-chunk-1

## 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)

plot of chunk unnamed-chunk-1

What a difference!

3 comments:

  1. 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

    5.
    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

    ReplyDelete
    Replies
    1. I should mention that without the colors argument the function plots beautifully.

      Delete
    2. Hi 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.

      Delete

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