# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
                  # # # #                                                                 # # # #
                  # # # #         New Odobenids remains phylogeny + morphometry           # # # #
                  # # # #                 Written by N. chatar 2024                       # # # #
                  # # # #            Related paper: Boisville et al. 2024                 # # # # 
                  # # # #  New species of Ontocetus (Pinnipedia: Odobenidae) from the     # # # # 
                  # # # # Lower Pleistocene of the North Atlantic shows similar feeding   # # # # 
                  # # # # adaptation independent to the extant walrus (Odobenus rosmarus) # # # # 
                  # # # #                                                                 # # # # 
                  # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 


# Set wd to folder containing script
setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) 

# # # Phylogeny, implied weighting # # #
install.packages("strap")
install.packages("paleotree")
install.packages("phytools")
install.packages("Claddis")

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

tree <- read.nexus("./Phylogeny/New_sp_Ontocetus_IW12.nex")

taxa_list <- c("Enaliarctos_emlongi", "Pteronarctos_goedertae", "Callorhinus_ursinus", "Desmatophoca_oregonensis", 
               "Allodesmus_kernensis", "Erignathus_barbatus", "Monachus_monachus", "Prototaria_primigena", 
               "Prototaria_planicephala", "Proneotherium_repenningi", "Neotherium_mirum",  "Kamtschatarctos_sinelnikovae", 
               "Pseudotaria_muramotoi", "Archaeodobenus_akamatsui", "Imagotaria_downsi", "Pelagiarctos", 
               "Nanodobenus_arandai","Pontolis_magnus", "Pontolis_barroni", "Pontolis_kohnoi", "Titanotaria_orangensis", 
               "Osodobenus_eodon","Dusignathus_santacruzensis","Dusignathus_seftoni", "Gomphotaria_pugnax", 
               "Aivukus_cedrosensis", "Protodobenus_japonicus", "Ontocetus_sp.", "Ontocetus_emmonsi", 
               "Odobenus_rosmarus", "Pliopedia_pacifica", "Valenictus_imperialensis", "Valenictus_sheperdi", 
               "Valenictus_chulavistensis", "Valenictus_UCMP_137426", "Ontocetus_posti")

consensus <- consensus(tree)

# Replace number exported in the Nexus file by the correct taxa name
for (i in 1:length(consensus$tip.label)) 
  {
   to_replace <- which(consensus$tip.label == i)
   consensus$tip.label[to_replace] <- taxa_list[i]
}

# Import taxa ages
ages_taxa <- read.csv("./Phylogeny/ages.csv", dec = ".", header = FALSE, row.names = 1)
colnames(ages_taxa) <- c("FAD", "LAD")

# Timescaling
timescaling_cons_tree <- timePaleoPhy(consensus, timeData = ages_taxa, type = "equal", vartime = 1)
geoscalePhylo(ladderize(timescaling_cons_tree,right=FALSE),ages_taxa,cex.ts=1,cex.tip=1, width = 2.5)

# # # Morphometry # # #

install.packages("ggplot2")
install.packages('ggfortify')
install.packages("ggthemes")
install.packages("ggrepel")
install.packages("ggpubr")
install.packages("viridis")

library(ggplot2)
library(ggfortify)
library(ggthemes) 
library(ggrepel)
library(ggpubr)
library(viridis)

setwd("C:/Users/Narimane Chatar/EDDy Lab Dropbox/Narimane Chatar/DECAF data/Collab/New_species_Ontocetus/Morphometry")

Measurments <- read.csv("Measurments_odobenini.csv", header = TRUE)
# Extract spe
species <- Measurments$full_name

shape_to_plot <- c("Ontocetus_posti_male"     = 0, 
                   "Ontocetus_posti_female"   = 15,
                   "Ontocetus_emmonsi_male"   = 1, 
                   "Ontocetus_emmonsi_female" = 16,
                   "Odobenus_rosmarus_male"   = 2,
                   "Odobenus_rosmarus_female" = 17)

color_to_plot <- c("Ontocetus_posti_male"     = inferno(3)[1], 
                   "Ontocetus_posti_female"   = inferno(3)[1],
                   "Ontocetus_emmonsi_male"   = inferno(3)[2], 
                   "Ontocetus_emmonsi_female" = inferno(3)[2],
                   "Odobenus_rosmarus_male"   = viridis(3)[2],
                   "Odobenus_rosmarus_female" = viridis(3)[2])

# keep specimen above the completeness threshold 50%
threshold <- 0.5
Measurments_thresold <- Measurments[apply(!is.na(Measurments[, 7:44]),1,sum)/dim(Measurments[, 7:44])[2]>threshold,]

# sclae numeric part of the matrix 
Measurments_thresold[, 7:44] <- scale(Measurments_thresold[, 7:44])

#Dissimilarity matrix
dissimilarity_measurments <- dist(Measurments_thresold[, 7:44], method = "euclidean")

cluster_measurments <- hclust(dissimilarity_measurments)
cluster_measurments$labels <- Measurments_thresold$to_plot

plot(cluster_measurments)

#Compute PCoA and plot it
pcoa_measurments <- pcoa(dissimilarity_measurments, correction = "cailliez")

df_pcoa_measurments <-cbind(as.data.frame(pcoa_measurments$vectors), 
                            "Species" = Measurments_thresold$full_name, 
                            "To_plot" = Measurments_thresold$to_plot)

plot_pcoa_measurments <-  ggplot(data=df_pcoa_measurments, aes(x=Axis.1,y=Axis.2, colour=To_plot, shape=To_plot)) +
                    theme_minimal() +
                    geom_point(size= 5,alpha=0.8) +
                    scale_color_manual(values = color_to_plot) +
                    scale_fill_manual(values = color_to_plot) +
                    scale_shape_manual(values=shape_to_plot) +
                    geom_rangeframe(color = "black") +
                    labs(x = paste0('Axis 1 = ', round(((pcoa_measurments$values$Rel_corr_eig[1])*100), digits = 2), '%'), 
                         y = paste0('Axis 2 = ', round(((pcoa_measurments$values$Rel_corr_eig[2])*100), digits = 2), '%')) 
plot_pcoa_measurments