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

#install packages

if("vegan" %in% rownames(installed.packages()) == FALSE){install.packages("vegan")
} else {print (paste0("'vegan' has already been installed in library"))}
if("fossil" %in% rownames(installed.packages()) == FALSE){install.packages("fossil")
} else {print (paste0("'fossil' has already been installed in library"))}
if("ecodist" %in% rownames(installed.packages()) == FALSE){install.packages("ecodist")
} else {print (paste0("'ecodist' has already been installed in library"))}
if("RStoolbox" %in% rownames(installed.packages()) == FALSE){install.packages("RStoolbox")
} else {print (paste0("'RStoolbox' has already been installed in library"))}
if("raster" %in% rownames(installed.packages()) == FALSE){install.packages("raster")
} else {print (paste0("'raster' has already been installed in library"))}
if("pacman" %in% rownames(installed.packages()) == FALSE){install.packages("pacman")
} else {print (paste0("'pacman' has already been installed in library"))}
if("geomorph" %in% rownames(installed.packages()) == FALSE){install.packages("geomorph")
} else {print (paste0("'geomorph' has already been installed in library"))}
if("ape" %in% rownames(installed.packages()) == FALSE){install.packages("ape")
} else {print (paste0("'ape' has already been installed in library"))}



library(geomorph)
library(ape)
library(ecodist)
library(geodata)
library(raster)
library(RStoolbox)
library(fossil)
library(vegan)

#################################################################################################################
############################################# INPUT FILES #######################################################
# 1- cvs file with the geographic coordinates of each individual
# 2- tps file
# 3- Environmental layers used in the ENMs 



#################################################################################################################
########################################## Morphometric Distance ################################################

setwd(" ") # Set working directory with the tps


# 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), ]

coord<- metadata_final[,7:8] #select only the Long/Lat column

mDIST <- earth.dist(coord, dist = TRUE)
mDIST <- as.dist(mDIST)
mDIST_vec <- as.vector(mDIST)
hist(mDIST_vec)

pcnm<-pcnm(mDIST, dist.ret = FALSE)

mDIST2<-as.matrix(dist(mDIST))
#write.csv(mDIST2, "geo_dist.csv") # save the geographic distance of the genetic points


PCA_mor <- gm.prcomp(gpa$coords)


PCs_mor<- PCA_mor$x[,1:5] ### ~70% of the variation

PCA_mor_dist<-ecodist::distance(PCs_mor, method = "mahalanobis", sprange=NULL, spweight=NULL) 

PCA_mor_dist2<-as.matrix(dist(PCA_mor_dist))

#write.csv(PCA_mor_dist2, "PCA_mor_dist.csv") # save the PCA genetic distance


#################################################################################################################
######################################### Environmental Variable ################################################

dir.create("env_variables") #creating a new directory
setwd("./env_variables") #set working directory with the environmental variables for the entire area


# download variables 
bio_data <- worldclim_global(var = "bio", res = 10 , path = tempdir())

# Visualize
plot(bio_data[[1]])  # First bioclim variable

# Save each variable as separate ASCII file
dir.create("./bioclim_ascii")
output_dir <- "bioclim_ascii"

for(i in 1:nlyr(bio_data)) {
  output_file <- file.path(output_dir, paste0(names(bio_data)[i], ".asc"))
  
  # Correct way to write ASCII format:
  writeRaster(
    bio_data[[i]],
    filename = output_file,
    datatype = "FLT4S",  # 4-byte floating point
    NAflag = -9999,
    overwrite = TRUE
  )
  
  message("Saved: ", output_file)
}


setwd("./bioclim_ascii/")
file.remove(list.files(pattern ='asc.aux.xml'))
file.remove(list.files(pattern ='.prj'))

#reading current bioclim data
files <- list.files(".", pattern='asc', full.names=TRUE )

predictors<- raster::stack(files)
names(predictors)

#Perform the PCA:
pca_env <- rasterPCA(predictors, nComp=1, norm=T, spca=F)
summary(pca_env$model) # PC1 78% of the total variance
pca_env$model$loadings


PC1_mor<-extract(pca_env$map$PC1, coord)

#write.csv(PC1_mor, "PC1_env_mor.csv") # save the PCA env distance

#################################################################################################################
########################################### Redundancy Analysis  ################################################


### Morphometric

PC1_mor[is.na(PC1_mor)] <- 0

dbrda4<-capscale(PCA_mor_dist2 ~ PC1_mor[,2], dist = "man") 
summary(dbrda4)
sig4<- anova(dbrda4)
sig4

dbrda5<-capscale(PCA_mor_dist2 ~ pcnm$vectors, dist = "man") 
dbrda5
summary(dbrda5)
sig5<- anova(dbrda5)
sig5

dbrda6<-capscale(PCA_mor_dist2 ~ PC1_mor[,2] + Condition(pcnm$vectors), dist = "man") 
summary(dbrda6)
sig6<- anova(dbrda6)
sig6