# PACKAGES -----------------------
{
  install.packages("Morpho")
  install.packages("rgl")
  install.packages("geomorph")
  install.packages("Rvcg")
  install.packages("ape")
  install.packages("ggplot2")
  install.packages("dplyr")
  install.packages("tidyr")
  install.packages("BBmisc")
  install.packages("ggtree")
  install.packages("colorRamps")
  install.packages("vegan")
  install.packages("phytools")
}

{
  library(Morpho)
  library(rgl)
  library(geomorph)
  library(Rvcg)
  library(ape)
  library(ggplot2)
  library(dplyr)
  library(tidyr)
  library(BBmisc)
  library(ggtree)
  library(colorRamps)
  library(vegan)
  library(phytools)
}

# TEMPLATE CREATION -------------------------
template.lm <- read.pts("./Ceratotherium_simum_cottoni_M51854_Patella_LM.pts")
template.mesh <- ply2mesh(file="./Ceratotherium_simum_cottoni_M51854_Patella.ply") 

# Anatomical landmarks
fixed <- as.integer(c(1:6))

# Matrix with surface sliding semiLMs
patch <- as.matrix(template.lm[110:570,])

# Anatomical LM and curve sliding semiLMs
template.lm2 <- template.lm[1:109,]

# Vector containing the rowindices of all semiLMs to be slided
slidings <- as.integer(c(7:570))

# Vector containing the rowindices of surface semiLMs to be slided
surface <- as.integer(c(110:570))

# Curves
# Edge of the articular surface (clockwise rotation)
SC1d <- as.integer(c(7:92))
# Medial ridge
SC2d <- as.integer(c(93:109))

# Curves list
curves <- list(SC1d,SC2d)

# Atlas creation
atlas <- createAtlas(template.mesh, as.matrix(template.lm2), patch = patch, corrCurves = curves, patchCurves = NULL, keep.fix = fixed)
plotAtlas(atlas, pt.size = 0.1, render = "s", point="s", legend=F)

# SPECIMENS DEFINTION ------------------------------

{
  datadef <- list.files(full.names = TRUE, pattern=".pts")
  def <- lapply(datadef,read.pts) 
  array.lm <- list2array(def)
  lab<-as.factor(datadef)
  lab<-substr(lab, start=3, stop=80)
  lab<-gsub(".pts", "", lab)
  dimnames(array.lm)[[3]]<-lab
}

# LMs SLIDING -----------------------------

# Projection onto the template
patched <- placePatch(atlas, array.lm, inflate = 2, ray = FALSE, path="./")
checkLM(patched, path="./", atlas = atlas, render = "s")

# Relaxation of the LMs on the template
template.tot <- rbind(as.matrix(template.lm2),patch)
relaxed.template <- array(dim=dim(patched))
for (i in 1:dim(patched)[3])
{
  relaxed.template[,,i]<-relaxLM(patched[,,i], template.tot, 
                                 slidings, outlines = curves,
                                 surp = surface, sur.name= paste("./",dimnames(patched)[[3]][i],".ply",sep=""), 
                                 mesh = NULL, tol = 1e-09,
                                 deselect = FALSE, inc.check = TRUE, iterations = 5, fixRepro = TRUE,
                                 missing = NULL, bending = TRUE, stepsize = 1,
                                 use.lm = NULL)
}
dimnames(relaxed.template)[[3]]<-dimnames(patched)[[3]]
checkLM(relaxed.template, path="./", atlas=atlas, render="s",alpha=1)

# Sliding onto the specimens
slided<-slider3d(relaxed.template, slidings, outlines = curves, surp = surface,
                 sur.path = ".", sur.name = NULL, meshlist = NULL, ignore = NULL,
                 sur.type = "ply", tol = 1e-07, deselect = FALSE, inc.check = T,
                 recursive = TRUE, iterations = 2, initproc = TRUE,
                 pairedLM = 0,mc.cores = 2, bending = TRUE,
                 fixRepro = TRUE) 

checkLM(slided$dataslide, path="./", atlas=atlas, render="s")

# Subsample 1: adult rhinos only, creation of a dataslide called "slided1"
# Subsample 2: all perissodactyls, creation of a dataslide called "slided2"

#################################
########## SUBSAMPLE 1 ##########
#################################

pca1 <- procSym(slided1$dataslide, scale = TRUE)

# colors
couleurs_sub1 <- c("grey","grey","grey","grey","grey","grey","grey","grey","grey","grey","red","red","red","red","red","red","black","black","black","blue","blue","green","green","green","green","green","green","green")
sex_sub1 <- c("M","M","M","F","F","M","F","M",NA,NA,"M",NA,NA,NA,NA,"M",NA,"M",NA,"M",NA,"F","M","F",NA,"F","M","F")
circ_femur_sub1 <- c(240.504,227.902,228.179,233.056,218.029,236.607,227.007,223.436,231.777,209.875,168.993,174.063,179.173,163.460,170.566,159.646,NA,185.431,211.147,206.699,209.406,230.095,184.772,224.071,231.211,210.418,208.353,217.742)

# Plot with convex hulls
PCscores_sub1 <- as.data.frame(pca1$PCscores[,1:2])
PCscores_sub1 <- cbind(PCscores_sub1, couleurs_sub1)

varPC1_sub1 <- paste0('PC1 = ', round((pca1$Variance[1,2]), digits = 1), '%')
varPC2_sub1 <- paste0('PC2 = ', round((pca1$Variance[2,2]), digits = 1), '%')

hull_data_sub1 <-   
  PCscores_sub1 %>%
  drop_na() %>%
  group_by(couleurs_sub1) %>% 
  slice(chull(PC1, PC2))

ggplot(PCscores_sub1, aes(x = PC1, y = PC2)) +
  geom_point(aes(color=couleurs_sub1, size=pca1$size/20)) +
  geom_polygon(data = hull_data_sub1, color = NA, aes(fill = couleurs_sub1, color = couleurs_sub1), alpha = 0.3, show.legend = FALSE) +
  geom_text(label = labels_sub1, size = 3, nudge_x = 0.004) +
  scale_color_manual(values=c("black","blue","green","grey","red")) +
  scale_fill_manual(values=c("black","blue","green","grey","red")) +
  labs(x = varPC1_sub1, y = varPC2_sub1) +
  coord_fixed() +
  scale_size_identity()

# NEIGHBOUR JOINING TREE ----------------------------------
distances_sub1 = dist(pca1$PCscores)
arbre_nj_sub1 = nj(distances_sub1)
ggtree(arbre_nj_sub1, layout = "ape") +
  geom_tippoint(color=couleurs_sub1, size=pca1$size/20) +
  geom_tiplab(size=2, hjust = 0.5, angle = 0)

# CORRELATION TESTS --------------
#PC1
log_taille_sub1 <- as.vector(log(pca1$size))

PC1_sub1 <- as.vector(pca1$PCscores[,1])
summary(lm(PC1_sub1~log_taille_sub1))
correlation_PC1_sub1 <- cor.test(PC1_sub1,log_taille_sub1, method="pearson")
correlation_PC1_sub1

# PC2
PC2_sub1 <- as.vector(pca1$PCscores[,2])
summary(lm(PC2_sub1~log_taille_sub1))
correlation_PC2_sub1 <- cor.test(PC2_sub1,log_taille_sub1, method="pearson")
correlation_PC2_sub1

# PERMANOVA -------------------
adonis2(pca1$PCscores ~ as.factor(couleurs_sub1) * as.factor(sex_sub1), method = "euclidean", permutations = 999)

# THEORETiCAL SHAPES ---------------------
meanshape_rhinos_sub1 <- as.matrix(pca1$mshape)
findMeanSpec(pca1$rotated)
closest_lm_rhinos_sub1 <- as.matrix(dataslide_sub1[,,4]) # coordoonées slidées de Dicerorhinus_sumatrensis_M81892_Patella
# Deformed meanshape coordinates on the nearest mesh to create a meanshape mesh
closest_mesh_rhinos_sub1 <- ply2mesh(file="Ceratotherium_simum_cottoni_C37_Patella.ply")
mesh_meanshape_rhinos_sub1 <- tps3d(closest_mesh_rhinos_sub1, closest_lm_rhinos_sub1, meanshape_rhinos_sub1)
shade3d(mesh_meanshape_rhinos_sub1, col = "ivory")

# Deformation of the meanshape obtained on the minima and maxima of each axis
# PC1
max_PC1_rhinos_sub1 <- restoreShapes(max(pca1$PCscores), pca1$PCs[,1],pca1$mshape)
min_PC1_rhinos_sub1 <- restoreShapes(min(pca1$PCscores), pca1$PCs[,1],pca1$mshape)

Patella_rhinos_sub1_Min_PC1 <- tps3d(mesh_meanshape_rhinos_sub1, meanshape_rhinos_sub1, min_PC1_rhinos_sub1)
shade3d(Patella_rhinos_sub1_Min_PC1, col = "ivory")
Patella_rhinos_sub1_Max_PC1 <- tps3d(mesh_meanshape_rhinos_sub1, meanshape_rhinos_sub1, max_PC1_rhinos_sub1)
shade3d(Patella_rhinos_sub1_Max_PC1, col = "ivory")
spheres3d(max_PC1_rhinos_sub1, radius = 0.001)

# PC2
max_PC2_rhinos_sub1 <- restoreShapes(max(pca1$PCscores), pca1$PCs[,2],pca1$mshape)
min_PC2_rhinos_sub1 <- restoreShapes(min(pca1$PCscores), pca1$PCs[,2],pca1$mshape)

Patella_rhinos_sub1_Min_PC2 <- tps3d(mesh_meanshape_rhinos_sub1, meanshape_rhinos_sub1, min_PC2_rhinos_sub1)
shade3d(Patella_rhinos_sub1_Min_PC2, col = "ivory")
Patella_rhinos_sub1_Max_PC2 <- tps3d(mesh_meanshape_rhinos_sub1, meanshape_rhinos_sub1, max_PC2_rhinos_sub1)
shade3d(Patella_rhinos_sub1_Max_PC2, col = "ivory")

# Transparents meshs + vectors (adapted from Pintore et al. 2023)
##PC1
acp <- gm.prcomp(gpa_sub1$coords)
arot <- arrayspecs(t(acp[["rotation"]]),p=(dim(acp[["rotation"]])[1]/3),k=3) # pcag = gm.prcomp
contribute <- apply((arot^2),3,rowSums)
dit2 <- cut(as.numeric(contribute[,1]),breaks = 10) #in contribute[,x] x=PC axis
idcolor <- blue2red(10)
for (i in 1:length(contribute[,1])) #in contribute[,x] x=PC axis
{
  tmpseg<-rbind(min_PC1_rhinos_sub1[i,],max_PC1_rhinos_sub1[i,])
  segments3d(tmpseg, lwd=3 ,color=idcolor[as.numeric(dit2[i])])
}
plot3d(Patella_rhinos_sub1_Min_PC1, col="burlywood1", axes=FALSE, box=FALSE, alpha=0.5, add=T)
plot3d(Patella_rhinos_sub1_Max_PC1, col="grey95", axes=FALSE, box=FALSE, alpha=0.5, add=T)

##PC2
arot <- arrayspecs(t(acp[["rotation"]]),p=(dim(acp[["rotation"]])[1]/3),k=3) # pcag = gm.prcomp
contribute <- apply((arot^2),3,rowSums)
dit2 <- cut(as.numeric(contribute[,2]),breaks = 10) #in contribute[,x] x=PC axis
idcolor <- blue2red(10)
for (i in 1:length(contribute[,2])) #in contribute[,x] x=PC axis
{
  tmpseg<-rbind(min_PC2_rhinos_sub1[i,],max_PC2_rhinos_sub1[i,])
  segments3d(tmpseg, lwd=3 ,color=idcolor[as.numeric(dit2[i])])
}
plot3d(Patella_rhinos_sub1_Min_PC2, col="burlywood1", axes=FALSE, box=FALSE, alpha=0.5, add=T)
plot3d(Patella_rhinos_sub1_Max_PC2, col="grey95", axes=FALSE, box=FALSE, alpha=0.5, add=T)

# ALLOMETRY ------------------
# With Csize
gpa_sub1 <- gpagen(dataslide_sub1)
gdf_sub1 <- geomorph.data.frame(gpa_sub1, couleurs_sub1, circ_femur_sub1)
allometry_sub1  <- procD.lm(coords ~ log(Csize), data=gdf_sub1, RRPP = T)
summary(allometry_sub1)
plot_allom_csize <- plotAllometry(allometry_sub1, size = gdf_sub1$Csize, logsz = T, method = "RegScore", col = couleurs_sub1, pch = 19, cex = gdf_sub1$Csize/50)
abline(lm(plot_allom_csize$plot.args$y) ~ as.numeric(plot_allom_csize$plot.args$x))

# Additional test without Dicerorhinus and Rh. sondaicus
dataslide_sans_dicero <- dataslide_sub1[,,-c(11:16,20:21)]
couleurs_sans_dicero <- couleurs_sub1[-c(11:16,20:21)]
gpa_sans_dicero <- gpagen(dataslide_sans_dicero)
gdf_sans_dicero <- geomorph.data.frame(gpa_sans_dicero, couleurs_sans_dicero)
allometry_sans_dicero  <- procD.lm(coords ~ log(Csize), data=gdf_sans_dicero, RRPP = T)
summary(allometry_sans_dicero)
plot_allom_sans_dicero <- plotAllometry(allometry_sans_dicero, size = gdf_sans_dicero$Csize, logsz = T, method = "RegScore", col = couleurs_sans_dicero, pch = 19, cex = gdf_sans_dicero$Csize/50)
abline(lm(plot_allom_sans_dicero$plot.args$y ~ as.numeric(plot_allom_sans_dicero$plot.args$x)))

# With femoral circumference
allometry_fem_sub1  <- procD.lm(coords ~ log(circ_femur_sub1), data=gdf_sub1, RRPP = T)
summary(allometry_fem_sub1)
plot_allom_circ_fem <- plotAllometry(allometry_sub1, size = gdf_sub1$circ_femur_sub1, logsz = T, method = "RegScore", col = couleurs_sub1, pch = 19, cex = gdf_sub1$Csize/50)

# Correlation Csize / Femoral circumference
correlation_csize_circ <- cor.test(log(gdf_sub1$Csize), log(circ_femur_sub1), method="pearson")
correlation_csize_circ

# Allometric theoretical shapes with Csize
preds_csize <- shape.predictor(allometry_sub1$GM$fitted, x = plot_allom_csize$RegScore, 
                         predmin = min(plot_allom_csize$RegScore), 
                         predmax = max(plot_allom_csize$RegScore))

# Deformation of the meanshape obtained on the minima and maxima of allometry
Min_allom_sub1 <- tps3d(mesh_meanshape_rhinos_sub1, meanshape_rhinos_sub1, preds_csize$predmin)
shade3d(Min_allom_sub1, col = "ivory")
Max_allom_sub1 <- tps3d(mesh_meanshape_rhinos_sub1, meanshape_rhinos_sub1, preds_csize$predmax)
shade3d(Max_allom_sub1, col = "ivory")

# Heatmap
heatmap_allom_sub1 <- meshDist(Min_allom_sub1, mesh2 = Max_allom_sub1, steps = 200, save = TRUE, file = "Heatmap_allom_sub1")


#################################
########## SUBSAMPLE 2 ##########
#################################

# ACP ----------------------
#morpho

pca2<-procSym(slided2$dataslide, scale = TRUE)

labels_patella<-substr(lab, start=1, stop=80) 
labels_patella<-gsub("_Left", "", labels_patella)
labels_patella<-gsub("_Right", "", labels_patella)
labels_patella<-gsub("_miroir", "", labels_patella)
labels_patella<-gsub("_Patella", "", labels_patella)

couleurs <-c("grey","grey","grey","grey","grey","grey","grey","grey","grey","grey","red","red","red","red","red","red","black","black","black","deeppink3","deeppink3","deeppink3","deeppink3","brown3","slateblue4","slateblue4","magenta4","deeppink4","royalblue4","royalblue4","brown4","salmon3","salmon4","tomato3","blue","blue","green","green","green","green","green","green","darkgoldenrod4","darkgoldenrod3","darkgoldenrod3","darkgoldenrod3","darkgoldenrod3","darkgoldenrod3","wheat3","wheat3","wheat4","wheat4","wheat4","wheat4")
sex <- c("M","M","M","F","F","M","F","M",NA,NA,"M",NA,NA,NA,NA,"M",NA,"M",NA,NA,"F","F","M",NA,"F","M","M","M","M","F","M","F",NA,"M","M",NA,"F","M","F",NA,"F","M",NA,"M","F","M","M","M","M",NA,"M","M",NA,NA)
circ_femur <- c(240.504,227.902,228.179,233.056,218.029,236.607,227.007,223.436,231.777,209.875,168.993,174.063,179.173,163.46,170.566,159.646,NA,185.431,211.147,210.661,106.188,96.433,124.555,128.884,147.547,165.066,141.761,126.331,132.742,138.189,140.94,131.587,139.107,123.258,206.699,209.406,230.095,225.837,224.071,231.211,210.418,208.353,97.3488,107.199,81.848,107.6,113.334,120.607,97.623,NA,111.034,111.125,103.946,80.7567)

# Plot with convex Hulls
PCscores <- as.data.frame(pca2$PCscores[,1:2])
PCscores <- cbind(PCscores, couleurs)

variances <- pca2$Variance
lab1 <- paste0('PC1 = ', round((variances[1,2]), digits = 2), '%')
lab2 <- paste0('PC2 = ', round((variances[2,2]), digits = 2), '%')

hull_data <-   
  PCscores %>%
  drop_na() %>%
  group_by(couleurs) %>% 
  slice(chull(PC1, PC2))

ggplot(PCscores, aes(x = PC1, y = PC2)) +
  geom_point(aes(color=couleurs, size=log(pca2$size)), shape = symbols) +
  geom_polygon(data = hull_data, color = NA, aes(fill = couleurs, color = couleurs), alpha = 0.3, show.legend = FALSE) +
  geom_text(label = labels_patella, size = 2, nudge_x = 0.004) +
  labs(x = lab1, y = lab2) +
  scale_color_manual(values=c("black","blue","brown3","brown4","darkgoldenrod3","darkgoldenrod4","deeppink3","deeppink4","green","grey","magenta4","red","royalblue4","salmon3","salmon4","slateblue4","tomato3","wheat3","wheat4")) +
  scale_fill_manual(values=c("black","blue","brown3","brown4","darkgoldenrod3","darkgoldenrod4","deeppink3","deeppink4","green","grey","magenta4","red","royalblue4","salmon3","salmon4","slateblue4","tomato3","wheat3","wheat4")) +
  coord_fixed() +
  scale_size_identity()

# NEIGHBOUR JOINING ----------------------------------

distances = dist(pca2$PCscores)
arbre_nj = nj(distances)
ggtree(arbre_nj, layout = "ape") +
  geom_tippoint(color=couleurs, size=log(pca2$size), shape = symbols) +
  geom_tiplab(size=2, hjust = 0.5)

# CORRELATION TESTS --------------
#PC1
log_taille_centroide <- as.vector(log(pca2$size))
PC1 <- as.vector(pca2$PCscores[,1])
summary(lm(PC1~log_taille_centroide))
correlation_PC1 <- cor.test(PC1,log_taille_centroide, method="pearson")
correlation_PC1

# PC2
PC2 <- as.vector(pca2$PCscores[,2])
summary(lm(PC2~log_taille_centroide))
correlation_PC2 <- cor.test(PC2,log_taille_centroide, method="pearson")
correlation_PC2

# PERMANOVA
adonis2(pca2$PCscores ~ as.factor(couleurs) + as.factor(sex), method = "euclidean", permutations = 999)

############
# Subsequent code for theoretical shapes, vector vizualisation, allometric tests and heatmaps similar to previous code for subsamples 1 & 2
############

# PHYLOGENETIC SIGNAL ----------
arbre <- read.nexus(file = "./arbre_patellas_perisso_actuels.nex")
plot.phylo(arbre, cex = 0.7)

# Mean per specimen
{
Ceratotherium_simum <- arrMean3(pca2$rotated[,,1:10])
Dicerorhinus_sumatrensis <- arrMean3(pca2$rotated[,,11:16])
Diceros_bicornis <- arrMean3(pca2$rotated[,,17:19])
Equus_africanus_asinus <- arrMean3(pca2$rotated[,,20:23])
Equus_quagga_burchellii <- pca2$rotated[,,24]
Equus_ferus_caballus <- arrMean3(pca2$rotated[,,25:26])
Equus_grevyi <- pca2$rotated[,,27]
Equus_hemionus <- pca2$rotated[,,28]
Equus_ferus_przewalskii <- arrMean3(pca2$rotated[,,29:30])
Equus_quagga_boehmi <- pca2$rotated[,,31]
Equus_quagga_chapmani <- pca2$rotated[,,32]
Equus_quagga_quagga <- pca2$rotated[,,33]
Equus_zebra_hartmannae <- pca2$rotated[,,34]
Rhinoceros_sondaicus <- arrMean3(pca2$rotated[,,35:36])
Rhinoceros_unicornis <- arrMean3(pca2$rotated[,,37:42])
Tapirus_bairdii <- pca2$rotated[,,43]
Tapirus_indicus <- arrMean3(pca2$rotated[,,44:48])
Tapirus_pinchaque <- arrMean3(pca2$rotated[,,49:50])
Tapirus_terrestris <- arrMean3(pca2$rotated[,,51:54])
}

array_specimens_moyens <- bindArr(as.matrix(Ceratotherium_simum),as.matrix(Dicerorhinus_sumatrensis),as.matrix(Diceros_bicornis),as.matrix(Equus_africanus_asinus),as.matrix(Equus_quagga_burchellii),as.matrix(Equus_ferus_caballus),as.matrix(Equus_grevyi),as.matrix(Equus_hemionus),as.matrix(Equus_ferus_przewalskii),as.matrix(Equus_quagga_boehmi),as.matrix(Equus_quagga_chapmani),as.matrix(Equus_quagga_quagga),as.matrix(Equus_zebra_hartmannae ),as.matrix(Rhinoceros_sondaicus),as.matrix(Rhinoceros_unicornis),as.matrix(Tapirus_bairdii),as.matrix(Tapirus_indicus),as.matrix(Tapirus_pinchaque),as.matrix(Tapirus_terrestris), along=3)

# Mean Csize
{
  Ceratotherium_simum_csize <- mean(pca2$size[1:10])
  Dicerorhinus_sumatrensis_csize <- mean(pca2$size[11:16])
  Diceros_bicornis_csize <- mean(pca2$size[17:19])
  Equus_africanus_asinus_csize <- mean(pca2$size[20:23])
  Equus_quagga_burchellii_csize <- mean(pca2$size[24])
  Equus_ferus_caballus_csize <- mean(pca2$size[25:26])
  Equus_grevyi_csize <- mean(pca2$size[27])
  Equus_hemionus_csize <- mean(pca2$size[28])
  Equus_ferus_przewalskii_csize <- mean(pca2$size[29:30])
  Equus_quagga_boehmi_csize <- mean(pca2$size[31])
  Equus_quagga_chapmani_csize <- mean(pca2$size[32])
  Equus_quagga_quagga_csize <- mean(pca2$size[33])
  Equus_zebra_hartmannae_csize <- mean(pca2$size[34])
  Rhinoceros_sondaicus_csize <- mean(pca2$size[35:36])
  Rhinoceros_unicornis_csize <- mean(pca2$size[37:42])
  Tapirus_bairdii_csize <- mean(pca2$size[43])
  Tapirus_indicus_csize <- mean(pca2$size[44:48])
  Tapirus_pinchaque_csize <- mean(pca2$size[49:50])
  Tapirus_terrestris_csize <- mean(pca2$size[51:54])
} 

points_csize_moyenne_csize <- c(Ceratotherium_simum_csize,Dicerorhinus_sumatrensis_csize,Diceros_bicornis_csize,Equus_africanus_asinus_csize,Equus_quagga_burchellii_csize,Equus_ferus_caballus_csize,Equus_grevyi_csize,Equus_hemionus_csize,Equus_ferus_przewalskii_csize,Equus_quagga_boehmi_csize,Equus_quagga_chapmani_csize,Equus_quagga_quagga_csize,Equus_zebra_hartmannae_csize,Rhinoceros_sondaicus_csize,Rhinoceros_unicornis_csize,Tapirus_bairdii_csize,Tapirus_indicus_csize,Tapirus_pinchaque_csize,Tapirus_terrestris_csize)

names(points_csize_moyenne_csize) <- c("Ceratotherium_simum","Dicerorhinus_sumatrensis","Diceros_bicornis","Equus_africanus_asinus","Equus_quagga_burchellii","Equus_ferus_caballus","Equus_grevyi","Equus_hemionus","Equus_ferus_przewalskii","Equus_quagga_boehmi","Equus_quagga_chapmani","Equus_quagga_quagga","Equus_zebra_hartmannae","Rhinoceros_sondaicus","Rhinoceros_unicornis","Tapirus_bairdii","Tapirus_indicus","Tapirus_pinchaque","Tapirus_terrestris")

dimnames(array_specimens_moyens)[[3]] <- c("Ceratotherium_simum","Dicerorhinus_sumatrensis","Diceros_bicornis","Equus_africanus_asinus","Equus_quagga_burchellii","Equus_ferus_caballus","Equus_grevyi","Equus_hemionus","Equus_ferus_przewalskii","Equus_quagga_boehmi","Equus_quagga_chapmani","Equus_quagga_quagga","Equus_zebra_hartmannae","Rhinoceros_sondaicus","Rhinoceros_unicornis","Tapirus_bairdii","Tapirus_indicus","Tapirus_pinchaque","Tapirus_terrestris")

labels_moyenne <- c("Ceratotherium_simum","Dicerorhinus_sumatrensis","Diceros_bicornis","Equus_africanus_asinus","Equus_quagga_burchellii","Equus_ferus_caballus","Equus_grevyi","Equus_hemionus","Equus_ferus_przewalskii","Equus_quagga_boehmi","Equus_quagga_chapmani","Equus_quagga_quagga","Equus_zebra_hartmannae","Rhinoceros_sondaicus","Rhinoceros_unicornis","Tapirus_bairdii","Tapirus_indicus","Tapirus_pinchaque","Tapirus_terrestris")

couleurs_moyenne <- c("grey","red","black","deeppink3","brown3","slateblue4","magenta4","deeppink4","royalblue4","brown4","salmon3","salmon4","tomato3","blue","green","darkgoldenrod4","darkgoldenrod3","wheat3","wheat4")

pca2_moyenne <- procSym(array_specimens_moyens, scale = TRUE)

# Test for phylogenetic signal
gpa <- gpagen(array_specimens_moyens)
signalphy <- physignal(gpa$coords, phy = arbre, iter = 999)
summary(signalphy)

# Blomberg K on Csize
indice_k <- phylosig(arbre, points_csize_moyenne_csize, test=TRUE, method="K", nsim=1000)
indice_k


# EQUIDS ONLY ------------

# ACP 
pca2_equids<-procSym(slided2$dataslide[,,20:34], scale = TRUE)

labels_patella_equids<-substr(lab[20:34], start=1, stop=80) 
labels_patella_equids<-gsub("_Left", "", labels_patella_equids)
labels_patella_equids<-gsub("_Right", "", labels_patella_equids)
labels_patella_equids<-gsub("_miroir", "", labels_patella_equids)
labels_patella_equids<-gsub("_Patella", "", labels_patella_equids) 

# couleurs
couleurs_equids <-c("deeppink3","deeppink3","deeppink3","deeppink3","brown3","slateblue4","slateblue4","magenta4","deeppink4","royalblue4","royalblue4","brown4","salmon3","salmon4","tomato3")
sex_equids <- c(NA,"F","F","M",NA,"F","M","M","M","M","F","M","F",NA,"M")
symbols_equids <- c(19,19,19,19,19,19,19,19,19,19,19,19,19,19,17)

# Plot with convex Hulls
PCscores_equids <- as.data.frame(pca2_equids$PCscores[,1:2])
PCscores_equids <- cbind(PCscores_equids, couleurs_equids)

variances_equids <- pca2_equids$Variance
lab1_equids <- paste0('PC1 = ', round((variances_equids[1,2]), digits = 2), '%')
lab2_equids <- paste0('PC2 = ', round((variances_equids[2,2]), digits = 2), '%')

hull_data_equids <-   
  PCscores_equids %>%
  drop_na() %>%
  group_by(couleurs_equids) %>% 
  slice(chull(PC1, PC2))

ggplot(PCscores_equids, aes(x = PC1, y = PC2)) +
  geom_point(aes(color=couleurs_equids, size=log(pca2_equids$size)), shape = symbols_equids) +
  geom_polygon(data = hull_data_equids, color = NA, aes(fill = couleurs_equids, color = couleurs_equids), alpha = 0.3, show.legend = FALSE) +
  geom_text(label = labels_patella_equids, size = 2, nudge_x = 0.004) +
  labs(x = lab1_equids, y = lab2_equids) +
  scale_color_manual(values=c("brown3","brown4","deeppink3","deeppink4","magenta4","royalblue4","salmon3","salmon4","slateblue4","tomato3")) +
  scale_fill_manual(values=c("brown3","brown4","deeppink3","deeppink4","magenta4","royalblue4","salmon3","salmon4","slateblue4","tomato3")) +
  coord_fixed() +
  scale_size_identity()

# NEIGHBOUR JOINING ----------------------------------

distances_equids = dist(pca2_equids$PCscores)
arbre_nj_equids = nj(distances_equids)
ggtree(arbre_nj_equids, layout = "ape") +
  geom_tippoint(color=couleurs_equids, size=log(pca2_equids$size), shape = symbols_equids) +
  geom_tiplab(size=2, hjust = 0.5)

# THEORETICAL SHAPES ---------------------
meanshape_equids <- as.matrix(pca2_equids$mshape)
findMeanSpec(pca2_equids$rotated)
closest_lm_equids <- as.matrix(slided3$dataslide[,,38])

closest_mesh_equids <- ply2mesh(file="Equus_quagga_chapmani_RBINS_1218_Patella_Left_miroir.ply")
mesh_meanshape_equids <- tps3d(closest_mesh_equids, closest_lm_equids, meanshape_equids)
shade3d(mesh_meanshape_equids, col = "ivory")

# PC1
max_PC1_equids <- restoreShapes(max(pca2_equids$PCscores), pca2_equids$PCs[,1],pca2_equids$mshape)
min_PC1_equids <- restoreShapes(min(pca2_equids$PCscores), pca2_equids$PCs[,1],pca2_equids$mshape)

Patella_equids_Min_PC1 <- tps3d(mesh_meanshape_equids, meanshape_equids, min_PC1_equids)
shade3d(Patella_equids_Min_PC1, col = "ivory")
Patella_equids_Max_PC1 <- tps3d(mesh_meanshape_equids, meanshape_equids, max_PC1_equids)
shade3d(Patella_equids_Max_PC1, col = "ivory")
spheres3d(max_PC1_equids, radius = 0.001)

# PC2
max_PC2_equids <- restoreShapes(max(pca2_equids$PCscores), pca2_equids$PCs[,2],pca2_equids$mshape)
min_PC2_equids <- restoreShapes(min(pca2_equids$PCscores), pca2_equids$PCs[,2],pca2_equids$mshape)

Patella_equids_Min_PC2 <- tps3d(mesh_meanshape_equids, meanshape_equids, min_PC2_equids)
shade3d(Patella_equids_Min_PC2, col = "ivory")
Patella_equids_Max_PC2 <- tps3d(mesh_meanshape_equids, meanshape_equids, max_PC2_equids)
shade3d(Patella_equids_Max_PC2, col = "ivory")

# TAPIRS ONLY ------------

# ACP
pca2_tapirs<-procSym(slided3$dataslide[,,43:54], scale = TRUE)

labels_patella_tapirs<-substr(lab[43:54], start=1, stop=80) 
labels_patella_tapirs<-gsub("_Left", "", labels_patella_tapirs)
labels_patella_tapirs<-gsub("_Right", "", labels_patella_tapirs)
labels_patella_tapirs<-gsub("_miroir", "", labels_patella_tapirs)
labels_patella_tapirs<-gsub("_Patella", "", labels_patella_tapirs) 

# couleurs
couleurs_tapirs <-c("darkgoldenrod4","darkgoldenrod3","darkgoldenrod3","darkgoldenrod3","darkgoldenrod3","darkgoldenrod3","wheat3","wheat3","wheat4","wheat4","wheat4","wheat4")
sex_tapirs <- c(NA,"M","F","M","M","M","M",NA,"M","M",NA,NA)
symbols_tapirs <- c(19,19,19,19,19,19,19,19,19,19,19,19)

# Plot with convex Hulls
PCscores_tapirs <- as.data.frame(pca2_tapirs$PCscores[,1:2])
PCscores_tapirs <- cbind(PCscores_tapirs, couleurs_tapirs)

variances_tapirs <- pca2_tapirs$Variance
lab1_tapirs <- paste0('PC1 = ', round((variances_tapirs[1,2]), digits = 2), '%')
lab2_tapirs <- paste0('PC2 = ', round((variances_tapirs[2,2]), digits = 2), '%')

hull_data_tapirs <-   
  PCscores_tapirs %>%
  drop_na() %>%
  group_by(couleurs_tapirs) %>% 
  slice(chull(PC1, PC2))

ggplot(PCscores_tapirs, aes(x = PC1, y = PC2)) +
  geom_point(aes(color=couleurs_tapirs, size=log(pca2_tapirs$size)), shape = symbols_tapirs) +
  geom_polygon(data = hull_data_tapirs, color = NA, aes(fill = couleurs_tapirs, color = couleurs_tapirs), alpha = 0.3, show.legend = FALSE) +
  geom_text(label = labels_patella_tapirs, size = 2, nudge_x = 0.004) +
  labs(x = lab1_tapirs, y = lab2_tapirs) +
  scale_color_manual(values=c("darkgoldenrod3","darkgoldenrod4","wheat3","wheat4")) +
  scale_fill_manual(values=c("darkgoldenrod3","darkgoldenrod4","wheat3","wheat4")) +
  coord_fixed() +
  scale_size_identity()

# NEIGHBOUR JOINING ----------------------------------

distances_tapirs = dist(pca2_tapirs$PCscores)
arbre_nj_tapirs = nj(distances_tapirs)
ggtree(arbre_nj_tapirs, layout = "ape") +
  geom_tippoint(color=couleurs_tapirs, size=log(pca2_tapirs$size), shape = symbols_tapirs) +
  geom_tiplab(size=2, hjust = 0.5)

# THEORETICAL SHAPES ---------------------
meanshape_tapirs <- as.matrix(pca2_tapirs$mshape)
findMeanSpec(pca2_tapirs$rotated)
closest_lm_tapirs <- as.matrix(slided3$dataslide[,,63])

closest_mesh_tapirs <- ply2mesh(file="Tapirus_terrestris_RBINS_1185D_Patella_Right.ply")
mesh_meanshape_tapirs <- tps3d(closest_mesh_tapirs, closest_lm_tapirs, meanshape_tapirs)
shade3d(mesh_meanshape_tapirs, col = "ivory")

# PC1
max_PC1_tapirs <- restoreShapes(max(pca2_tapirs$PCscores), pca2_tapirs$PCs[,1],pca2_tapirs$mshape)
min_PC1_tapirs <- restoreShapes(min(pca2_tapirs$PCscores), pca2_tapirs$PCs[,1],pca2_tapirs$mshape)

Patella_tapirs_Min_PC1 <- tps3d(mesh_meanshape_tapirs, meanshape_tapirs, min_PC1_tapirs)
shade3d(Patella_tapirs_Min_PC1, col = "ivory")
Patella_tapirs_Max_PC1 <- tps3d(mesh_meanshape_tapirs, meanshape_tapirs, max_PC1_tapirs)
shade3d(Patella_tapirs_Max_PC1, col = "ivory")
spheres3d(max_PC1_tapirs, radius = 0.001)

# PC2
max_PC2_tapirs <- restoreShapes(max(pca2_tapirs$PCscores), pca2_tapirs$PCs[,2],pca2_tapirs$mshape)
min_PC2_tapirs <- restoreShapes(min(pca2_tapirs$PCscores), pca2_tapirs$PCs[,2],pca2_tapirs$mshape)

Patella_tapirs_Min_PC2 <- tps3d(mesh_meanshape_tapirs, meanshape_tapirs, min_PC2_tapirs)
shade3d(Patella_tapirs_Min_PC2, col = "ivory")
Patella_tapirs_Max_PC2 <- tps3d(mesh_meanshape_tapirs, meanshape_tapirs, max_PC2_tapirs)
shade3d(Patella_tapirs_Max_PC2, col = "ivory")