##########################################################################################
####################################### Joyce Prado ######################################

# Script to analyze geometric morphometric data

##########################################################################################

#--------------------------------------#
#         Pacotes e leitura de dados   #
#--------------------------------------#

# Instala o pacote 'geomorph' se necessário
if (!requireNamespace("geomorph", quietly = TRUE)) {
  install.packages("geomorph")
} else {
  message("'geomorph' já está instalado.")
}

if (!requireNamespace("ggplot2", quietly = TRUE)) {
  install.packages("ggplot2")
} else {
  message("'ggplot2' já está instalado.")
}

library(geomorph)
library(ggplot2)

# Leitura do arquivo TPS com coordenadas dos marcos anatômicos
landmarks_raw <- readland.tps("Alouatta_frontal_landmarks_calculados.TPS", specID = "ID")
dim(landmarks_raw)
landmarks_raw[, , 1]

# Leitura da tabela com informações dos indivíduos
metadata <- read.csv("Variables_alouatta_pepe_final.csv", sep = ";")
head(metadata)
str(metadata)

# Definição das conexões      

# Define os pares de marcos conectados
landmark_links <- read.csv("linksfrontalAlouatta.csv", sep = ",")
landmark_links <- landmark_links[,2:3]

#--------------------------------------#
#   Seleção de adultos e subadultos    #
#--------------------------------------#

# Seleciona apenas adultos e subadultos
adult_sub_id <- metadata[metadata$Age %in% c("Adult", "Subadult"), ]


# Converte para matriz 2D, filtra e reconstrói array 3D
landmarks_2d <- two.d.array(landmarks_raw)
landmarks_adult_sub <- landmarks_2d[adult_sub_id$ID, ]
landmarks_array <- arrayspecs(landmarks_adult_sub, p = dim(landmarks_raw)[1], k = 2)

# Confirma estrutura dos dados
dim(landmarks_array)
str(landmarks_array)

#--------------------------------------#
#           Alinhamento GPA            #
#--------------------------------------#

gpa <- gpagen(landmarks_array)
summary(gpa)

# If you have semilandmarks you can add the argument curves

# Plota todos os espécimes com links
plotAllSpecimens(gpa$coords, links = landmark_links)

#--------------------------------------#
#           Detecção de outliers       #
#--------------------------------------#

plotOutliers(gpa$coords, inspect.outliers = TRUE)

# IDs a remover
outliers <- c("ALMAMA19166AM", "ALCLFS11121PM")

# Remove outliers do GPA
to_keep <- !dimnames(gpa$coords)[[3]] %in% outliers
gpa$coords <- gpa$coords[, , to_keep]
gpa$Csize <- gpa$Csize[to_keep]
gpa$ProcD <- gpa$ProcD[to_keep]

# Atualiza metadados após remoção dos outliers

metadata_final <- adult_sub_id[match(dimnames(gpa$coords)[[3]], adult_sub_id$ID), ]


#--------------------------------------#
#      Criação do data.frame final     #
#--------------------------------------#

df_geomorph <- geomorph.data.frame(
  coord  = gpa$coords,
  size   = gpa$Csize,
  species = metadata_final$Species,
  sex     = metadata_final$Sex,
  age     = metadata_final$Age,
  biome   = metadata_final$biome
)

str(df_geomorph)


#### Principal Components Analysis (PCA)
PCA <- gm.prcomp(gpa$coords)
plot(PCA)

gps <- as.factor(paste(df_geomorph$species))  #define some groups for plotting
plot.new()
plot(PCA, pch=22, cex = 1.5, bg = gps) 
legend("topleft", pch=22, pt.bg = unique(gps), legend = unique(gps), ncol=2)

gps <- as.factor(paste(df_geomorph$sex))  #define some groups for plotting
plot.new()
plot(PCA, pch=22, cex = 1.5, bg = gps) 
legend("topleft", pch=22, pt.bg = unique(gps), legend = unique(gps), ncol=2)

gps <- as.factor(paste(df_geomorph$biome))  #define some groups for plotting
plot.new()
plot(PCA, pch=22, cex = 1.5, bg = gps) 
legend("topleft", pch=22, pt.bg = unique(gps), legend = unique(gps), ncol=2)

ggplot(metadata_final, aes(x=PCA$x[,1], y=PCA$x[,2], colour=Species))+geom_point(colour= 'transparent', size=8)+
  geom_jitter(aes(col=Species))+ 
  stat_ellipse(geom = "polygon",alpha = 0.1, aes(fill = Species))+
  scale_fill_manual(values=c("red3","black","darkslategray3","tan2","chocolate","violet","greenyellow","navyblue","springgreen4","coral4"))+
  labs(title = "PCA_Alouatta_species", x= "PCA 1 (29.59%)", y= "PCA 2 (18.48%)")+ theme_bw() 

#ggsave("PCA_frontal_Alouatta_species.png", width=20, height=12, units= "cm", dpi=300)


ggplot(metadata_final, aes(x=PCA$x[,1], y=PCA$x[,2], colour=biome))+geom_point(colour= 'transparent')+
  geom_jitter(aes(shape=biome, col=biome),size=5)+
  scale_shape_manual(values=c(3,16,17,18))+
  #stat_ellipse(geom = "polygon",alpha = 0, aes(fill = Species))+
  scale_color_manual(values=c("red3","black","darkslategray3","green4"))+
  labs(title = "PCA_Skull_frontal_Alouatta_biomes", x= "PCA 1 (29.59%)", y= "PCA 2 (18.48%)")+ theme_bw()


##### Shape Visualization

M <- mshape(gpa$coords)

par(mfrow=c(1,2))
plotRefToTarget(M,gpa$coords[,,39], 
                links = landmark_links, 
                method="vector", mag=3)
mtext("Vector Displacements")
plotRefToTarget(M, gpa$coords[,,39], 
                links = landmark_links, 
                gridPars = gridPar(pt.bg="red", 
                                   link.col="green", 
                                   pt.size = 1), 
                method="vector", mag=3)
mtext("Vector Displacements: Other Options")

#Predict a shape and plot this using shape.predictor (very useful!). This time along PC1:

PC <- PCA$x[,1]
preds <- shape.predictor(gpa$coords, x= PC, Intercept = FALSE, 
                         pred1 = min(PC),  
                         pred2 = max(PC)) # PC 1 extremes, more technically
par(mfrow=c(1, 2))
plotRefToTarget(M, preds$pred1, links = landmark_links)
mtext("PC1 - Min.")
plotRefToTarget(M, preds$pred2, links = landmark_links)
mtext("PC1 - Max.")

#### Average shape per species
shape.2d.means<-rowsum(two.d.array(df_geomorph$coord),metadata_final$Species)/as.vector(table(metadata_final$Species))
shape.means<-arrayspecs(shape.2d.means,dim(landmarks_raw)[1],dim(landmarks_raw)[2])
shape.means

#### Average size per species
size.means<-rowsum(df_geomorph$size,metadata_final$Species)/as.vector(table(metadata_final$Species))
size.means
#write.csv(size.means,file = "size_means_frontal_Alouatta_adultsubadult.csv") 


##### GLM: Regression

# Procrustes ANOVA with permutation
fit.null <- procD.lm(coords ~ 1, data = gpa, 
                     iter = 9999)
fit.alt <- procD.lm(coords ~ log(Csize), data = gpa,
                    turbo = FALSE, iter = 9999)
anova(fit.alt)

plot.new()
par(mfcol = c(1, 1))
plot(fit.alt, type = "regression", 
     reg.type = "RegScore", 
     predictor = gpa$Csize)


#### Centroid size analyses
size <- gpa$Csize
metadata_final<- cbind(metadata_final,size)
table_size_adu_sub_males<-subset(metadata_final, Sex=="male")
table_size_adu_sub_female<-subset(metadata_final, Sex=="female")


p <- ggplot(table_size_adu_sub_males, aes(x=Species, y=size)) + 
  geom_boxplot() + labs(title="Box_Plot_size_males") 
p

p <- ggplot(table_size_adu_sub_female, aes(x=Species, y=size)) + 
  geom_boxplot() + labs(title="Box_Plot_size_frontal_females") 
p

ggplot(metadata_final, aes(x = Species, y = size, fill = Sex)) +
  geom_boxplot(position = position_dodge(0.8)) +
  labs(title = "Boxplots por espécie e sexo",
       x = "Espécie",
       y = "Valor") +
  scale_fill_manual(values = c("male" = "skyblue", "female" = "pink")) +
  theme_minimal()


# Changing background
ggplot(table_size_adu_sub_males, aes(x=Species, y=size)) +
  geom_boxplot(fill='#A4A4A4', color="black")+
  theme_bw()+ labs(title="Box_Plot_size_frontal_males") 

# Change box plot colors by groups
p<-ggplot(metadata_final, aes(x=Species, y=size, fill=Species)) +
  geom_boxplot() + 
  labs(title="Box_Plot_size_frontal_Alouatta") + 
  theme(axis.text.x = element_text(face="italic",size=9, angle=0, vjust=0.8))#+
#geom_point(size=1,color="red")  
#geom_text(aes(label=Order, size=4),hjust=0, vjust=0) + theme_bw() 
p

#ggsave(file="Frontal_Alouatta.svg", plot=p, width=15, height=15)
#ggsave("Frontal_Alouatta.png", width=20, height=12, units= "cm", dpi=300)



### ANOVA

fit <- procD.lm(coord ~ biome * log(size), 
                data = df_geomorph, iter = 999)
anova(fit)

fitm <- manova.update(fit, print.progress = FALSE)
summary(fitm, test = "Pillai")


### Pairwise Comparisons: SIMPLE EXAMPLE

group <- interaction(df_geomorph$species, df_geomorph$sex)
PW <- pairwise(fit, groups = group)
summary(PW)



###### Allometry

plotAllSpecimens(gpa$coord)

df_geomorph$logSize <- log(gpa$Csize)
fit <- procD.lm(coord ~ logSize, 
                data = df_geomorph, 
                print.progress = FALSE) 

plot(fit, type = "regression", 
     reg.type = "RegScore", 
     predictor = df_geomorph$logSize, 
     pch = 19)

anova(fit)

fit$residuals

PCA_corrected <- gm.prcomp(fit$residuals)
plot(PCA_corrected)

gps <- as.factor(paste(df_geomorph$sex))  #define some groups for plotting
plot.new()
plot(PCA_corrected, pch=22, cex = 1.5, bg = gps) 
legend("topleft", pch=22, pt.bg = unique(gps), legend = unique(gps), ncol=2)

#### Comparing Allometric Trajectories

df_geomorph$Group <- interaction(metadata_final$Sex, metadata_final$Age)

fit.common <- procD.lm(coord ~ logSize + Group, 
                       data = df_geomorph, print.progress = FALSE) 
fit.unique <- procD.lm(coord ~ logSize * Group, 
                       data = df_geomorph, print.progress = FALSE)  
anova(fit.unique)



### Plot
par(mfcol = c(1, 2))

plot(fit.common, type = "regression", 
     predictor = df_geomorph$logSize, 
     reg.type = "PredLine", 
     pch=19, 
     col = df_geomorph$Group)

legend("topleft", 
       legend = unique(df_geomorph$Group), 
       pch = 21, 
       pt.bg = unique(df_geomorph$Group))

mtext("Common Slopes")

plot(fit.unique, 
     type = "regression", 
     predictor = df_geomorph$logSize, 
     reg.type = "PredLine", 
     pch=19, 
     col = df_geomorph$Group)

legend("topleft", 
       legend = unique(df_geomorph$Group), 
       pch = 21, 
       pt.bg = unique(df_geomorph$Group))

mtext("Unique Slopes")


save.image("meu_projeto_morpho_geo.RData")