Wednesday, March 12, 2014

Plotting a right-facing round phylogram

Don't ask me why I'm working on this. Here's how to do it:

roundPhylogram<-function(tree){
  n<-length(tree$tip.label)
  # reorder cladewise to assign tip positions
  cw<-reorder(tree,"cladewise")
  y<-vector(length=n+cw$Nnode)
  y[cw$edge[cw$edge[,2]<=n,2]]<-1:n
  # reorder pruningwise for post-order traversal
  pw<-reorder(tree,"pruningwise")
  nn<-unique(pw$edge[,1])
  # compute vertical position of each edge
  for(i in 1:length(nn)){
    yy<-y[pw$edge[which(pw$edge[,1]==nn[i]),2]]
    y[nn[i]]<-mean(range(yy))
  }
  # compute start & end points of each edge
  X<-nodeHeights(cw)
  ## end preliminaries
  # open & size a new plot
  plot.new(); par(mar=rep(1.1,4))
  plot.window(xlim=c(0,1.05*max(X)),ylim=c(1,max(y)))
  # plot edges
  for(i in 1:nrow(X)){
    b<-y[cw$edge[i,1]]
    c<-X[i,1]
    d<-if(y[cw$edge[i,2]]>y[cw$edge[i,1]]) 1 else -1
    xx<-X[i,2]
    yy<-y[cw$edge[i,2]]
    a<-(xx-c)/(yy-b)^2
    curve(d*sqrt((x-c)/a)+b,from=X[i,1],to=X[i,2],add=TRUE,
      lwd=2)
  }
  # plot tip labels
  for(i in 1:n)
    text(X[which(cw$edge[,2]==i),2],y[i],tree$tip.label[i],
      pos=4,offset=0.3,font=2)
}

Now let's see how it works:

> library(phytools)
> source("roundPhylogram.R")
> tree<-pbtree(n=26,tip.label=LETTERS[26:1])
> roundPhylogram(tree)

1 comment: