Friday, March 8, 2019

Projecting a phylogenetic tree onto a map with multiple geographic points per taxon

In response to a recent user request I recently added the capability to map more than one geographic location per taxon to the phytools plotting function phylo.to.map. To see the code click here. A phytools version containing this update can already be installed via GitHub using devtools.

In the following demo I'll use some fictional data to illustrate how it can be used - first on a world map projecting from the top, then on a single-country map projecting from the left:

library(phytools)
packageVersion("phytools")
## [1] '0.6.71'
library(mapdata)
## here are my data
World.tree
## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  B.ijyqhp, J.wjbgkp, I.zgm, U.otngy, E.bgzoej, L.cyl, ...
## 
## Rooted; includes branch lengths.
head(World,20)
##                      lat       long
## B.ijyqhp    -53.75945233   89.82698
## B.ijyqhp    -49.27526159   92.08920
## B.iownj      -1.78488475 -162.17161
## B.iownj      -5.28332361 -160.89541
## B.iownj      -1.15704991 -169.82635
## B.xiwt      -18.32161500  -39.76940
## B.xiwt      -21.81437265  -42.49659
## B.xiwt      -16.38054009  -42.98694
## E.bgzoej    -33.49914588  167.24997
## E.bgzoej    -34.83267341  168.81525
## E.bgzoej    -33.50117255  164.69732
## F.lbfqkh      6.04544745  -99.61004
## F.lbfqkh     -0.06604802  -95.04093
## F.lbfqkh      3.37645461 -101.15848
## G.jltpcivk  -16.65924503  -49.84077
## G.jltpcivk  -14.68628460  -48.19405
## G.thjdsean   10.95726052  -73.92540
## G.thjdsean   10.30423882  -72.88593
## G.thjdsean    9.79876619  -70.26079
## G.xqpudsjri   4.33975876   27.89086
obj<-phylo.to.map(World.tree,World,plot=FALSE)
## objective: 82
## objective: 82
## objective: 72
## objective: 72
## objective: 72
## objective: 72
## objective: 70
## objective: 70
## objective: 70
## objective: 70
## objective: 68
## objective: 46
## objective: 46
## objective: 46
## objective: 40
## objective: 40
## objective: 40
## objective: 38
## objective: 38
## objective: 38
## objective: 38
## objective: 38
## objective: 38
## objective: 36
## objective: 34
plot(obj,ftype="i")

plot of chunk unnamed-chunk-1

If we like, it is very easy to plot the points & linking lines to our species to be shown in different colors. Here's a simple example:

cols<-setNames(sample(rainbow(n=Ntip(World.tree))),
    World.tree$tip.label)
plot(obj,colors=cols,ftype="i",fsize=0.6,cex.points=c(0.7,1.2))

plot of chunk unnamed-chunk-2

This also should work for sideways plotted trees. For the following it's necessary to install both mapdata (for the higher resolution country map) and viridis (for the color palette that I used).

library(mapdata)
library(viridis)
Argentina.tree
## 
## Phylogenetic tree with 12 tips and 11 internal nodes.
## 
## Tip labels:
##  J.zwilcgka, Z.higpu, X.usrfnti, H.noswc, K.jgquihs, P.xdwzvfj, ...
## 
## Rooted; includes branch lengths.
Argentina
##                   lat      long
## J.zwilcgka  -49.22635 -70.70463
## Z.higpu     -48.56640 -68.50478
## X.usrfnti   -47.35648 -70.04468
## H.noswc     -44.05671 -70.53964
## K.jgquihs   -42.90179 -67.73484
## P.xdwzvfj   -40.31697 -65.47999
## G.xydpt     -36.19226 -67.29487
## N.yinzpjo   -33.44245 -68.22980
## H.hpsycu    -31.46258 -64.05009
## Q.fstnuprx  -29.70271 -60.42035
## B.pnhma     -27.17288 -62.01523
## L.lfbpxwutn -24.25808 -64.98503
## J.zwilcgka  -49.92165 -70.27778
## Z.higpu     -47.80319 -69.42856
## X.usrfnti   -48.19909 -71.25053
## X.usrfnti   -48.46782 -70.90932
## H.noswc     -44.78242 -70.02675
## H.noswc     -43.45621 -71.00061
## K.jgquihs   -44.22061 -65.83883
## K.jgquihs   -42.57006 -67.32508
## G.xydpt     -35.57724 -67.16105
## N.yinzpjo   -33.73010 -68.13654
## H.hpsycu    -30.38485 -61.98099
## Q.fstnuprx  -28.14004 -60.40019
## Q.fstnuprx  -29.56328 -60.08278
## B.pnhma     -26.15989 -62.03782
## B.pnhma     -25.90764 -62.64044
## L.lfbpxwutn -24.41996 -65.55994
obj<-phylo.to.map(Argentina.tree,Argentina,database="worldHires",
    regions="Argentina",plot=FALSE)
## objective: 22
## objective: 22
## objective: 22
## objective: 20
## objective: 18
## objective: 18
## objective: 16
## objective: 16
## objective: 16
## objective: 14
## objective: 12
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-3

This is brand new & involved me working on code that I haven't really touched in a few years so please let me know if you come across any bugs!

For those interested in such things, here's how the data for this example were generated:

World.tree<-pbtree(n=26,scale=100)
World.tree$tip.label<-replicate(Ntip(World.tree),
    paste(sample(LETTERS,1),".",
    paste(sample(letters,round(runif(n=1,min=3,max=10))),
    collapse=""),
    sep=""))
lat<-fastBM(World.tree,sig2=10,bounds=c(-90,90))
long<-fastBM(World.tree,sig2=80,bounds=c(-180,180))
World<-cbind(lat,long)
for(i in 1:Ntip(World.tree)){
    ni<-sample(0:2,1)
    for(j in 1:ni){ 
        World<-rbind(World,c(World[i,1]+rnorm(n=1,sd=4),
            World[i,2]+rnorm(n=1,sd=4)))
        rownames(World)[nrow(World)]<-rownames(World)[i]
    }
}
Argentina<-locator(n=12) ## with world map plotted
Argentina<-cbind(Argentina$y,Argentina$x)
colnames(Argentina)<-c("lat","long")
Argentina.tree<-pbtree(n=nrow(Argentina),scale=100)
Argentina.tree$tip.label<-replicate(Ntip(Argentina.tree),
    paste(sample(LETTERS,1),".",
    paste(sample(letters,round(runif(n=1,min=3,max=10))),
    collapse=""),
    sep=""))
rownames(Argentina)<-Argentina.tree$tip.label
for(i in 1:Ntip(Argentina.tree)){
    ni<-sample(0:2,1)
    if(ni>0){
        for(j in 1:ni){
            Argentina<-rbind(Argentina,c(Argentina[i,1]+rnorm(n=1),
                Argentina[i,2]+rnorm(n=1)))
            rownames(Argentina)[nrow(Argentina)]<-rownames(Argentina)[i]
        }
    }
}

5 comments:

  1. I am having trouble with your artificial data for Argentina
    "Argentina<-locator(n=12)"
    kicks back the following error: "Error in locator(n = 12) : plot.new has not been called yet"

    ReplyDelete
    Replies
    1. Update: I think this was an issue with my World.map fake data, but now that it works and I have built my map, running the argentina lines cause R to think way to hard about that first line "Argentina<-locator(n=12)"

      Delete
  2. Is there any way to color two points coming from the same tip differently?

    For instance, I would like to have one map point colored blue and another red, but both should connect to the same tip of the tree.

    ReplyDelete
    Replies
    1. Hi Alex. Yes. The easiest way to do this is by simply adding your colored points at the end. There coordinates on the plot are just their longitudes (x) and latitudes (y) on a decimal scale. - Liam

      Delete

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