#Lewton et al., "The effects of phylogeny, body size, and locomotor behavior on the three-dimensional shape of the pelvis in extant carnivorans", submitted to PeerJ, 2020

library(geomorph)
library(geiger)

filelist <- list.files(pattern = "*.txt") ## list all morphologika files
carnivore <- read.morphologika(filelist)
dimnames(carnivore)[[3]] <- gsub(".txt","",filelist)

classifier<-read.csv("carnivore_group_matrix.csv", header=T, row.names=1)
loco<-classifier$loco
names(loco)<-dimnames(carnivore)[[3]]

taxon<-classifier$taxon
names(taxon)<-dimnames(carnivore)[[3]]

#estimate missing data
any(is.na(carnivore))
carnivore[which(carnivore == 9999)]<-NA
carnivore <- estimate.missing(carnivore,method="Reg")

#read in the curve sliders
curvepts<-as.matrix(read.csv("carnivore_curvepts.csv", header=T))

#read in the wireframe links
links<-as.matrix(read.csv("links.csv"))

carnivore.gpa <- gpagen(carnivore, curves=curvepts, PrinAxes=TRUE, max.iter=NULL, Proj=TRUE, print.progress=TRUE)

#Phylogenetic analyses
tree<-read.nexus("Nyakatura_Bininda-Emonds_2012_carn_tree_33sp_bestest.nex")
data<-read.delim("species mean pc scores.csv",sep=",",header=TRUE, row.names=1)
td <- treedata(tree,data, sort=TRUE, warnings=TRUE)
loco<-factor(td$data[,'loco'])

PC1<-td$data[,3]
PC1<-as.numeric(PC1)
names(PC1)<-rownames(td$data)

PC2<-td$data[,4]
PC2<-as.numeric(PC2)
names(PC2)<-rownames(td$data)

PC3<-td$data[,5]
PC3<-as.numeric(PC3)
names(PC3)<-rownames(td$data)

PC4<-td$data[,6]
PC4<-as.numeric(PC4)
names(PC4)<-rownames(td$data)

CS<-td$data[,2]
CS<-as.numeric(CS)
names(CS)<-rownames(td$data)

fit1 <- procD.pgls(PC1 ~ loco2 + log(CS), phy=td$phy, iter=9999)
fit2 <- procD.pgls(PC2 ~ loco2 + log(CS), phy=td$phy, iter=9999)
fit3 <- procD.pgls(PC3 ~ loco2 + log(CS), phy=td$phy, iter=9999)
fit4 <- procD.pgls(PC4 ~ loco2 + log(CS), phy=td$phy, iter=9999)

#make species means dataframe for phylomorphospace plot
carnivore.gpa.df <- geomorph.data.frame(carnivore.gpa, taxon=classifier$taxon)
new.coords <- coords.subset(A = carnivore.gpa.df$coords, group = carnivore.gpa.df$taxon)
names(new.coords)
spmean <- lapply(new.coords, mshape)
spmean <- simplify2array(spmean, higher=TRUE)

#test for phylogenetic signal
physignal(A=spmean, phy=td$phy)

#Plot phylomorphospace
carnPCA.phylo <- gm.prcomp(spmean, phy=td$phy)

plot(carnPCA.phylo, phylo = TRUE, phylo.par = list(edge.color = "black", edge.width = 1, edge.lty = 1,
                                                   node.bg = "black", node.pch = 21, node.cex = 1))
gps <- loco
  
plot(carnPCA.phylo, pch=21, bg=gps, cex=1.5, phylo = TRUE,
     phylo.par = list(edge.color="grey", node.cex=0)); title(main="phylomorphospace")

text(carnPCA.phylo$x, labels = labels(carnPCA.phylo$x)[[1]],
     pos = 2, font = 3, cex=0.6) 

#Make warped meshes for figures
ref<-mshape(spmean)
mesh <- read.ply("LACM_54170_Vulpes_Pelvis.ply")
mesh.coord <- read.morphologika("LACM_54170 Vulpes.txt")
meanmesh <- warpRefMesh(mesh=mesh, mesh.coord=mesh.coord[,,1], ref=ref, centered=FALSE)

PC1minmesh <- plotRefToTarget(ref,carnPCA.phylo$shapes$shapes.PC1$min, mesh=meanmesh, method="surface")
PC1maxmesh <- plotRefToTarget(ref,carnPCA.phylo$shapes$shapes.PC1$max, mesh=meanmesh, method="surface")
PC2minmesh <- plotRefToTarget(ref,carnPCA.phylo$shapes$shapes.PC2$min, mesh=meanmesh, method="surface")
PC2maxmesh <- plotRefToTarget(ref,carnPCA.phylo$shapes$shapes.PC2$max, mesh=meanmesh, method="surface")
PC3minmesh <- plotRefToTarget(ref,carnPCA.phylo$shapes$shapes.PC3$min, mesh=meanmesh, method="surface")
PC3maxmesh <- plotRefToTarget(ref,carnPCA.phylo$shapes$shapes.PC3$max, mesh=meanmesh, method="surface")
PC4minmesh <- plotRefToTarget(ref,carnPCA.phylo$shapes$shapes.PC4$min, mesh=meanmesh, method="surface")
PC4maxmesh <- plotRefToTarget(ref,carnPCA.phylo$shapes$shapes.PC4$max, mesh=meanmesh, method="surface")