Tuesday, January 31, 2017

New tree drawing method to extract phylogeny from a plotted figure

After yesterday posting a video illustrating how to build a tree by hand using bind.tip in the phytools package, it occurred to me that the next logical step would be a function that would allow the user to build a tree completely free-hand, including with edge lengths, on top of an image of a tree.

The following is some code that will do this, below which I've posted another video demo. For the example I have used the tree from Amemiya et al. (2013).

library(phytools)
library(jpeg)
get.treepos<-phytools:::get.treepos

tree.drawer<-function(img){
    par(fg=make.transparent("grey",0.8))
    img<-readJPEG(img)
    plot.new()
    par(mar=rep(0.1,4))
    plot.window(xlim=c(0,10),ylim=c(0,10))
    rasterImage(img,0,0,10,10)
    cat("  Click the position of the GLOBAL ROOT.\n")
    flush.console()
    root<-unlist(locator(1))
    cat("  Enter the name of a tip RIGHT of the root. > ")
    flush.console()
    right<-readLines(n=1)
    cat(paste("  Click on the position of ",right,".\n",sep=""))
    flush.console()
    right.xy<-unlist(locator(1))
    cat("  Enter the name of a tip LEFT of the root. > ")   
    flush.console()
    left<-readLines(n=1)
    cat(paste("  Click on the position of ",left,".\n",sep=""))
    flush.console()
    left.xy<-unlist(locator(n=1))
    left.xy
    tree<-list(edge=matrix(c(3,3,1,2),2,2),
        edge.length=c(right.xy[1]-root[1],left.xy[1]-root[1]),
        Nnode=1,tip.label=c(right,left))
    class(tree)<-"phylo"
    tips<-setNames(c(right.xy[2],left.xy[2]),tree$tip.label)
    names(tips)<-gsub(" ","_",names(tips))
    plotTree(tree,add=TRUE,tips=tips,xlim=c(0,10)-root[1],ylim=c(0,10),
        color=make.transparent("blue",0.4),lwd=4)
    tip<-0
    cat("  Enter the name of tip to add (or press ENTER). > ")
    flush.console()
    tip<-readLines(n=1)
    while(tip!=""){
        cat(paste("  Click on the position of ",tip,".\n",sep=""))
        flush.console()
        xy<-unlist(locator(1))
        cat("  Click on the position of its MRCA in the built tree.\n")
        flush.console()
        obj<-get.treepos(message=FALSE)
        tree<-bind.tip(tree,tip,edge.length=xy[1]-(nodeheight(tree,obj$where)-
        obj$pos),where=obj$where,position=obj$pos)
        tips<-c(tips,setNames(xy[2],tip))
        names(tips)<-gsub(" ","_",names(tips))
        plot.new()
        par(mar=rep(0.1,4))
        plot.window(xlim=c(0,10),ylim=c(0,10))
        rasterImage(img,0,0,10,10)
        plotTree(tree,add=TRUE,tips=tips,xlim=c(0,10)-root[1],ylim=c(0,10),
            color=make.transparent("blue",0.4),lwd=4)
        old<-tip
        cat("  Enter the name of tip to add (or press ENTER). > ")
        flush.console()
        tip<-readLines(n=1)
        while(tip=="GOBACK"){
            cat(paste("  Dropping ",old,".\n",sep=""))
            tree<-drop.tip(tree,old)
            plot.new()
            par(mar=rep(0.1,4))
            plot.window(xlim=c(0,10),ylim=c(0,10))
            rasterImage(img,0,0,10,10)
            plotTree(tree,add=TRUE,tips=tips,xlim=c(0,10)-root[1],ylim=c(0,10),
            color=make.transparent("blue",0.4),lwd=4)
            cat("  Enter the name of tip to add (or press ENTER). > ")
            flush.console()
            tip<-readLines(n=1)
        }   
    }
    par(fg="black")
    tree
}

tree<-tree.drawer(img="Amemiya_etal_2013-tree.jpg")

Here is our final tree:

plotTree(tree)

plot of chunk unnamed-chunk-2

The result, of course, is an object of class "phylo" that we can write to file or use in any other analysis that we might be interested in.

So far this only works for right-facing phylograms, but in principle one could extend to any plotting style - with a little work.

Astute readers may notice that the node rotations of the final tree are different than in the original tree - although the topology & branch lengths are correct. Well, should we be so inclined, we can 'fix' this using tipRotate as follows:

tip.order<-c("Elephant_shark","Little_skate","Spotted_catshark","Zebrafish",
    "Pufferfish","Tilapia","Coelecanth","Lungfish","Chinese_brown_frog",
    "Western_clawed_frog","Lizard","Zebra_finch","Turkey","Chicken","Platypus",
    "Opossum","Tammar_wallaby","Armadillo","Elephant","Mouse","Human","Dog")
tree<-tipRotate(tree,setNames(1:Ntip(tree),tip.order))
plotTree(tree)

plot of chunk unnamed-chunk-3

Neat.

1 comment:

  1. I really like your blog, it looks very nice, I'm happy to visit again to see your blog because it's very good indeed, thanks’ for all. You have such a great sense of style. I love reading your blog posts because you always give me inspiration to try something new and different. Last time I’m visit write my essay helps college students to enhance traditional essay writing abilities. Other than the college syllabus, writing ability is likewise needed to broaden self-development.

    ReplyDelete