Saturday, July 27, 2013

New phytools version with some phylogeny & geographic mapping options

I'm new to geographic maps (see 1, 2), but they are fun to play with so, after a little work, I've added the function phylo.to.map. It now creates an object of class "phylo.to.map" which can be used with the S3 plotting method plot.phylo.to.map. The function also allows optional user control of a bunch of different plotting options; however these are not yet thoroughly documented.

The new version of phytools with these updates (phytools 0.3-19) also depends on the R package maps.

Here's a quick demo:

> require(phytools)
Loading required package: phytools
> packageVersion("phytools")
[1] ‘0.3.19’
> # simulate tree & data
> tree<-pbtree(n=26,scale=100)
> tree$tip.label<-LETTERS[26:1]
> lat<-fastBM(tree,sig2=10,bounds=c(-90,90))
> long<-fastBM(tree,sig2=80,bounds=c(-180,180))
> # now plot projection
> phylo.to.map(tree,cbind(lat,long),asp=1.3)
objective: 168
...
objective: 78

We can also zoom in on particular regions of the map, although that is not well developed. So, for instance:

> tree<-pbtree(n=26,scale=100)
> tree$tip.label<-LETTERS[26:1]
> lat<-fastBM(tree,sig2=10,bounds=c(-36,-16),a=-20)
> long<-fastBM(tree,sig2=80,bounds=c(115,150),a=130)
> # now plot projection
> phylo.to.map(tree,cbind(lat,long),ylim=c(-40,-10), xlim=c(110,155))
objective: 166
...
objective: 114

It's also now possible to project the tree directly onto a geographic map. Users of this method, though, should keep in mind that "ancestral nodes" are just a projection of the tree - and are not intended as a reconstruction of ancestral areas:

> tree<-pbtree(n=26,scale=100)
> tree$tip.label<-LETTERS[26:1]
> tree<-paintSubTree(tree,length(tree$tip)+1,"1")
> lat<-fastBM(tree,sig2=10,bounds=c(-90,90))
> long<-fastBM(tree,sig2=80,bounds=c(-180,180))
> plot.new()
> plot.window(xlim=c(-180,180),ylim=c(-90,90),asp=1.3)
> map("world",add=TRUE)
> phylomorphospace(tree,cbind(long,lat), colors=setNames("red",1),node.by.map=TRUE,add=TRUE, label="horizontal")

4 comments:

  1. Would it be possible to Plot the base map at an angle, as if you were viewing the map from space, and plot a cladogram where the tips are anchored to the map?

    ReplyDelete
  2. Liam, is it possible to add latitude and longitude tick marks with this function phylo.to.map?

    Thanks, Lina

    ReplyDelete
    Replies
    1. Hi Lina. I just wrote a post about this here. Let me know if this is what you had in mind. - Liam

      Delete
  3. Thank you so much for this !
    It seems that this function doesn't work with any kind of tree.
    I tried this one, with success:

    (Copte:8.0,((Compostelle:2.0,(Byzance:1.0,Latin:1.0):1.0):3.0,(Roumanie:2.0,Russie:2.0,Bucovine:2.0):3.0):3.0,(Geez:1.0,Arabe:1.0):7.0);

    But this one failed, and I had a message "Error in xy.coords(x, y) : 'x' and 'y' lengths differ":

    (Copte:8.0,((Compostelle:2.0,(Byzance:1.0,Latin:1.0):1.0):3.0,(Roumanie:2.0,Russie:2.0,Bucovine:2.0):3.0):3.0,(Geez:1.0,Arabe:1.0):7.0);

    Is there a way to change this ?

    And is there a possibility to rotate the tree above the map?

    Ex.:

    I used this set of coordinates (named "coords.csv"):

    X lat long
    Copte 31.2 29.90
    Geez 9.0 38.70
    Arabe 30.0 31.00
    Compostelle 42.8 -8.50
    Byzance 41.0 28.90
    Latin 51.5 0.07
    Roumanie 47.0 25.50
    Russie 55.7 37.60
    Bucovine 48.0 26.00

    Then I wrote:

    coords <-read.csv("coords.csv", header = TRUE, sep= ";")
    row.names(coords)=coords$X
    coords=coords[,c(2,3)]
    phy <- "((Copte:2.0,(Geez:1.0,Arabe:1.0):1.0):2.0,((Bucovine:2.0,(Roumanie:1.0,Russie:1.0):1.0):1.0,(Compostelle:2.0,(Byzance:1.0,Latin:1.0):1.0):1.0):1.0);"
    tree <-read.newick(text=phy)
    plot(tree,use.edge.length=FALSE, edge.width=4,edge.color="brown", label.offset=0.15)
    obj<-phylo.to.map(tree,coords, fsize=0.9, psize=2, asp=1.2,split=c(0.3,0.5), mar=c(2.1,8.1,2.1,2.1), ylim=c(-10,75), xlim=c(-30,90))

    The result is OK, but I'd prefer to rotate the tree in order to have the left branch (Copte, Geez, Arabe) on the right site.

    Is this possible ?

    Thank you for your time !
    My best,

    JLLQ

    ReplyDelete