plot a radial ("fan") tree with a phenotypic trait as a bar plot

radial_phylo_bars.R

require(ape)
require(geiger)

# use Geospiza data for example
data(geospiza)

phylo <- drop.tip(geospiza$geospiza.tree, 'olivacea')

bardata <- setNames(geospiza$geospiza.data$beakD, rownames(geospiza$geospiza.data))

# create dummy categorical variable, highlighting ground finches

statedata <- setNames(c(1,0,0,0,1,1,0,0,0,0,0,0,0), rownames(geospiza$geospiza.data)) +1

# order data so bars align with tips
bar.ordered <- bardata[phylo$tip.label[phylo$edge[phylo$edge[, 2] <= Ntip(phylo), 2]]]
state.ordered <- statedata[phylo$tip.label[phylo$edge[phylo$edge[, 2] <= Ntip(phylo), 2]]]

# plot phylogeny

par(mar=c(0,0,0,0), oma=c(0,0,0,0))
plot(phylo, type='fan', show.tip.label=F, x.lim=c(-2.2,2.2), y.lim=c(-2.2,2.2), edge.width=3)

# set angles, convert to radians

angle <- seq(0,360, length.out=Ntip(phylo)+1)  # to rotate so it starts counting on right
angle <- angle[-length(angle)] # remove last one so it doesn't overlap with first one
angle <- angle*pi/180 # convert to radians

# set starting point for bars where tips end

start <- 0.65

# scale to apply so its not disproportional to plotting area

scaledat <- 5

# grid lines
# these have to be set by hand depending on range of data
# remember to scale the same way as the bars 


circle.lines <- seq(0,360, length.out=100)*pi/180

points(cos(circle.lines)*(start), 
  sin(circle.lines)*(start), type='l', col='grey')
	
points(cos(circle.lines)*(start+1/scaledat), 
	sin(circle.lines)*(start+1/scaledat), type='l', col='grey')

points(cos(circle.lines)*(start+2/scaledat), 
	sin(circle.lines)*(start+2/scaledat), type='l', col='grey')

points(cos(circle.lines)*(start+3/scaledat), 
	sin(circle.lines)*(start+3/scaledat), type='l', col='grey')

# "radius" determines size of bars; 
# also, add starting point

radius <- (bar.ordered/scaledat) + start

segments(cos(angle)*start, sin(angle)*start,cos(angle)*radius, sin(angle)*radius, lwd=10, 
	col=state.ordered, lend=1)

# if species labels are needed:
# condition rotates labels to the left so it is in reading direction
# NOTE adjustment difference

for(i in 1:length(angle)){

if(cos(angle[i]) > 0)
text(cos(angle[i])*(start+3/scaledat), sin(angle[i])*(start+3/scaledat), 
	labels=names(bar.ordered)[i], cex=0.5, srt=angle[i]*180/pi, adj=-0.1)

if(cos(angle[i]) < 0)
text(cos(angle[i])*(start+3/scaledat), sin(angle[i])*(start+3/scaledat), 
	labels=names(bar.ordered)[i], cex=0.5, srt=angle[i]*180/pi + 180, adj=1.1)

	
}