Create a gist now

Instantly share code, notes, and snippets.

@cmtucker /caRd2015.R Secret
Last active Dec 17, 2015

What would you like to do?
The EEB & Flow's annual holiday caRd, 2015
##2015 Holiday caRd from the EEB & Flow, by C. Tucker
##If you don't have these libraries, please install before running...
#I didn't test this on different OS/versions of R this year, so if you run into problems refer to the blog post instead :-)
##Req'd libraries
require(ape)
require(rworldmap)
require(FD)
# function modified from on phytools (& dependencies) and rworldmap (& dependencies) (written by Liam J. Revell 2013)
phylo.to.map2<-function(tree,coords){
# open & size a new plot
par(mai=c(0.1,0.1,0.1,0.1))
map<-getMap(resolution="low")
plot(map,ylim=c(-30,100))
# rescale tree so it fits in the upper half of the plot
# with enough space for labels
sh<-max(strwidth(tree$tip.label))/((par()$usr[2]-par()$usr[1])*(par()$usr[4]-par()$usr[3]))
tree$edge.length<-tree$edge.length/max(nodeHeights(tree))*(90-sh)
n<-length(tree$tip.label)
# reorder cladewise to assign tip positions
cw<-reorder(tree,"cladewise")
x<-vector(length=n+cw$Nnode)
x[cw$edge[cw$edge[,2]<=n,2]]<-0:(n-1)/(n-1)*360-180
# reorder pruningwise for post-order traversal
pw<-reorder(tree,"pruningwise")
nn<-unique(pw$edge[,1])
# compute horizontal position of each edge
for(i in 1:length(nn)){
xx<-x[pw$edge[which(pw$edge[,1]==nn[i]),2]]
x[nn[i]]<-mean(range(xx))
}
# compute start & end points of each edge
Y<-180-nodeHeights(cw)
# plot vertical edges
for(i in 1:nrow(Y))
lines((rep(x[cw$edge[i,2]],2)+10),(Y[i,]+10),lwd=2,lend=2, col="forestgreen")
# plot horizontal relationships
for(i in 1:tree$Nnode+n)
lines((range(x[cw$edge[which(cw$edge[,1]==i),2]])+10),(Y[which(cw$edge[,1]==i),1]+10),lwd=2,lend=2, col="forestgreen")
# plot tip labels
lab_short <- c("Coca \ncola Santa","Mall \nSanta","Salvation \nArmy.Santa","Kris.Kringle", "Santa \nClaus","Pere.Noel","Sinterklaas","Saint \nNicholas", "Father \nChristmas", "Ded \nMoroz")
for(i in 1:n)
text(x[i],Y[which(cw$edge[,2]==i),2],lab_short[i],srt=0, col="red", cex=0.5)
coords<-coords[tree$tip.label,2:1]
points(coords,pch=16,cex=1,col="red")
for(i in 1:n) lines(c(x[i],coords[i,1]),c(Y[which(cw$edge[,2]==i),2],coords[i,2]),col="red",lty="dashed")
}
#also from phytools
nodeHeights <- function (tree)
{
if (attr(tree, "order") != "cladewise" || is.null(attr(tree,
"order")))
t <- reorder(tree)
else t <- tree
root <- length(t$tip.label) + 1
X <- matrix(NA, nrow(t$edge), 2)
for (i in 1:nrow(t$edge)) {
if (t$edge[i, 1] == root) {
X[i, 1] <- 0
X[i, 2] <- t$edge.length[i]
}
else {
X[i, 1] <- X[match(t$edge[i, 1], t$edge[, 2]), 2]
X[i, 2] <- X[i, 1] + t$edge.length[i]
}
}
if (attr(tree, "order") != "cladewise" || is.null(attr(tree,
"order")))
o <- apply(matrix(tree$edge[, 2]), 1, function(x, y) which(x ==
y), y = t$edge[, 2])
else o <- 1:nrow(t$edge)
return(X[o, ])
}
########################
##Make tree
edge <- matrix(c(11,12,12,13,13,14,14,9,14,15,15,6,15,7,13,8, 12, 10,11,16,16,17,17,18,18,1,18,2,17,3,16,19,19,4,19,5), nrow=18, ncol=2, byrow=TRUE)
edge.length <- c(0.1, 0.4, 0.1, 0.5, 0.3, 0.2, 0.2, 0.6, 1.0, 0.2, 0.3, 0.3, 0.3, 0.3, 0.6, 0.1, 0.8, 0.8)
tip.label <- c("Coca.cola.Santa", "Mall.Santa", "Salvation.Army.Santa", "Kriss.Kringle", "Santa.Claus", "Pere.Noel", "Father.Christmas", "Sinterklaas", "Saint.Nicholas", "Ded.Moroz")
Nnode <- 9
xmas.tree <- list(edge, edge.length, Nnode, tip.label)
class(xmas.tree) <- "phylo"
names(xmas.tree) <- c("edge", "edge.length", "Nnode", "tip.label")
####Start
#taxa exploration
#functional ecology
heft <- c("fat", "fat", "fat", "fat", "fat", "thin", "fat","thin", "thin", "thin")
transport <- c("reindeer", "reindeer", "reindeer", "reindeer", "reindeer", "donkey", "foot", "horse", "foot", "foot")
first.appearance <- c(1900, 1900, 1890, 1800, 1700, 1400, 1400, 400, 300, 1937)
santatraits <- cbind.data.frame(heft, transport, first.appearance)
rownames(santatraits) <- xmas.tree$tip.label
###Card plotting
#Make device with set size, works on multiple OS
dev.new <- function(width = 5, height = 5){
platform <- sessionInfo()$platform
if (grepl("linux",platform)) {
x11(width=width, height=height)
}else{
if (grepl("pc",platform)) {
windows(width=width, height=height)
}else{
if (grepl("w32",platform)) {
windows(width=width, height=height)
}else{
if (grepl("apple", platform)){
quartz(width=width, height=height)} }}}}
dev.new(10, 6)
par(mfrow=c(1,2), bg="darkseagreen1", mai=c(0.1,0.1,0.1,0.1), oma=c(0.1,0.1,0.1,0.1))
plot(xmas.tree, type = "c", FALSE, edge.color="darkgreen", edge.lty=1, edge.width=18, label.offset = 1, direction="downward", font=3, tip.color="darkred")
#reconstruct ancestral state
cc <- ace(transport, xmas.tree, type="discrete")
co2 <- c("yellow", "gold", "darkorange", "red")
nodelabels(pie = cc$lik.anc, piecol = co2, cex = c(1.5, rep(1, 8)))
#Plot traits (fatness) against Santa
co1 <- c("blue", "purple")
tiplabels(pch = 12, col = co1[as.factor(heft)], cex = 3.5, adj=c(0.5, 0), lwd=2)
#Let's see the transportation mode trait too:
tiplabels(pch = 8, col = co2[as.factor(transport)], cex = 2, adj=c(0.5, 0), lwd=2)
#where do santa taxa live?
locales <- cbind(c(40.7, 40.7, 40.7, 70.93, 70.93, 47, 51.4, 52.2, 39.8, 55.7), c( -73.93, -73.93, -73.93, -95.4, -95.4, 2.7, -0.85, 5.7, 33.3, 37.6))
rownames(locales) <- xmas.tree$tip.label
phylo.to.map2(xmas.tree, locales)
text(x=-20, y=-120, "Happy Holidays from \nthe EEB & Flow", col="darkred", cex=2, font=4)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment