#Stratigraphic Tree for Goviraptor

library(strap)
library(Claddis)
library(phytools)


#Import data: tree and age table - for file.choose, pick biogeotree, for age table, pick ages
data = read.nexus.data(file.choose()) #use master matrix file
tree = read.tree(file.choose()) #use majrule .phy file, not the majrule outgroups!
outgrouptree = read.tree(file.choose()) #use outgroup .phy tree
strict = read.tree(file.choose()) #use strict .phy file
phylodata = read.csv(file.choose()) #use phylodataoutgroup file


#Modify tip label and ensure names match
tree$tip.label[42] = "Oksoko_avarsan"
outgrouptree$tip.label[48] = "Oksoko_avarsan"
rownames(phylodata) = outgrouptree$tip.label
phylodata = phylodata[,-1]

#Select variables from dataset
age = phylodata[,1:2]
rownames(age) = outgrouptree$tip.label
biogeo = phylodata[,5]
names(biogeo) = outgrouptree$tip.label
biogeo = biogeo[-c(1:9)]
digits = phylodata[,6]
names(digits) = outgrouptree$tip.label
digits2 = digits[!is.na(digits)]
forelimb = phylodata[,11]
names(forelimb) = outgrouptree$tip.label
forelimb2 = forelimb[!is.na(forelimb)]
hum = phylodata[,12]
names(hum) = outgrouptree$tip.label
hum2 = hum[!is.na(hum)]

#Time-calibrate the outgroup tree for forelimb evolution (two approaches: one with and one without tip lengths equal to age range)
straptree = DatePhylo(tree = outgrouptree, ages = age, rlen = 10, method = "equal", add.terminal = FALSE)
straptree2 = DatePhylo(tree = outgrouptree, ages = age, rlen = 10, method = "equal", add.terminal = TRUE)

#Time-calibrate the ingroup tree for biogeographic analyses
biogeotree = drop.tip(tree, c(1:3))
biogeostraptree = DatePhylo(tree = biogeotree, ages = age[10:nrow(age),1:2], rlen = 10, method = "equal", add.terminal = TRUE)

#Simulate biogeographic history
bioestsim = make.simmap(biogeostraptree, biogeo, model = "ER", nsim = 1000)
bioest = summary(bioestsim)
cols = c("purple","orange","black","white")
names(cols) = c("1","2","3","4")
plot(bioest, ftype = "i", colors = cols, cex = c(0.6,0.25), ylim = c(42,0))
add.simmap.legend(leg = c("Asia", "North America", "Southern China", "Western Gobi"), colors=cols,prompt=FALSE,x=4, y = 35, fsize=1, shape = "circle")

#Estimate ancestral states of digits (colour limits restricted to range exhibited by oviraptorosaurs)
digest = contMap(straptree2, digits2, ftype = "i", fsize = 0.7, lims = c(0, 0.962), lwd = 3, method = "anc.ML", direction = "upwards")
nodelabels(node = c(58:95), pie=bioest$ace, piecol = c("purple", "orange", "black", "white"), cex = 0.35)

#Estimate ancestral states of forelimbs (colour limits shifted for clarity -- difficult to distinguish greens on opposite sides of the spectrum)
forest = contMap(straptree, forelimb2, ftype = "i", fsize = 0.7, lims = c(1.0, 1.5), lwd = 3, method = "anc.ML", direction = "upwards")
nodelabels(node = c(58:95), pie=bioest$ace, piecol = c("purple", "orange", "black", "white"), cex = 0.35)

#Estimate ancestral states of Humerus/forelimb ratio (No pattern: forelimb reduction is from antebrachial elements, not overall forelimb)
humest = contMap(straptree, hum2, ftype = "i", fsize = 0.7, lims = c(1.0, 1.5), lwd = 3, method = "anc.ML", direction = "upwards")
nodelabels(node = c(58:95), pie=bioest$ace, piecol = c("purple", "orange", "black", "white"), cex = 0.35)