R code for working with Buddha cladistics data, part of the analysis documented here: http://www.academia.edu/1773421/A_Cladistic_Evaluation_of_Ancient_Thai_Bronze_Buddha_Images_Six_Tests_for_a_Phylogenetic_Signal_in_the_Griswold_Collection See here for the data that accompanies this code: https://gist.github.com/benmarwick/6081694

buddha_cladistics.R

rm(list = ls()) # clear workspace
# see https://gist.github.com/benmarwick/6081694 for data file to use
# with this code
library(phangorn) # load library phangorn 'Phylogenetic analysis in R'
# change the following line to match your computer
setwd("E:\\My Documents\\My Papers\\EurASEAA2010")
b <- read.nexus.data("infileb12.txt.nex") #make sure this nex file is in wd
b.pd <- phyDat(b, type="USER", levels=c(0:2)) # convert nex to phyDat format
b.pd #check what went in....
b.pd.pscore <- pratchet(b.pd, method="fitch")
#  maximum parsimony score (pscore) = 36
#
# returns all most parsimonious trees in an object of class multiPhylo
b.pd.bab <- bab(b.pd)
#
# plots all most parsimonious trees
library(ape)
plot(b.pd.bab, type="phylogram")
#
# find all trees, exhaustive and branch & bound MP searches
# and returns parsimony scores
library(phytools)
b.pd.e <- exhaustiveMP(b.pd, method="exhaustive")
# searches 945 trees, keeps 3
# find length of the kept trees
b.pd.e.pscores <- NULL
for(i in 1:length(b.pd.e))
{
  b.pd.e.pscores[i] <- attr(b.pd.e[[i]],"pscore")
}
b.pd.e.pscores
# all are 36 long
#
# find all trees, another method
b.pd.bb <- exhaustiveMP(b.pd, method="branch.and.bound")
# returns 3 phylogenetic trees, how long are they?
b.pd.bb.pscores <- NULL
for(i in 1:length(b.pd.bb))
{
  b.pd.bb.pscores[i] <- attr(b.pd.bb[[i]],"pscore")
}
b.pd.bb.pscores
# the same, 36
#
#  have a look a the trees resulting from the exhaustive search... they look good, best one is 3 because it appears to be rooted consistent with stratigraphic data
plot(b.pd.e, type="phylogram", layout=4)