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!

7 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
    3. I 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).

      Delete
    4. Nevermind, it was not minRotate that converted them to underscores. Must be something else.

      Delete
  2. Hi, do you know where I could get the maps and trees datasets?
    Regards

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

      Delete

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