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)

}

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)

> source("roundPhylogram.R")

> tree<-pbtree(n=26,tip.label=LETTERS[26:1])

> roundPhylogram(tree)

Gosh, that makes me want to ask.

ReplyDelete