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!

Posted by at

10 comments:

  1. Mark ScherzMarch 19, 2019 at 12:59 PM

    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. Mark ScherzMarch 19, 2019 at 12:59 PM

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

      Delete
    2. Lucas JordãoJune 5, 2019 at 9:09 PM

      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. John PowersFebruary 12, 2020 at 3:47 PM

      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. John PowersFebruary 16, 2020 at 10:12 PM

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

      Delete
  2. apardoMay 29, 2020 at 5:34 PM

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

    ReplyDelete
    Replies
    1. Liam RevellMay 29, 2020 at 7:56 PM

      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
  3. QPaps04May 19, 2021 at 2:37 PM

    Hi, 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
  4. DiegoJuly 11, 2022 at 1:06 AM

    Hi Liam Revell

    Object of class "phylo.to.map" containing:

    (1) A phylogenetic tree with 11 tips and 10 internal nodes.

    (2) A geographic map with range:
    -4.22N, 13.38N
    -81.72W, -66.87W.

    (3) A table containing 23 geographic coordinates (may include
    more than one set per species).

    The nodes of the tree may not have yet been rotated to maximize alignment between the phylogeny & the map.

    The map of Colombia appears but the points do not appear, nor the phylogenetic tree, I have tried to search by all means but I have not been able to visualize it completely

    ReplyDelete
  5. DiegoJuly 11, 2022 at 1:07 AM

    This comment has been removed by the author.

    ReplyDelete

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

Subscribe to: Post Comments (Atom)