#######################################################################################################################################################################
### Script for morphological data treatment and plotting  																											###
### Written by V. Fischer, June 2025          																														###
### Associated paper: Fischer, Weis et al. Vampyromorph coleoid predation by an ichthyosaurian from the Early Jurassic Lagerstätte of Bascharage, Luxembourg. PeerJ	###
#######################################################################################################################################################################

### Libraries

library(readxl)
library(ape)
library(ggplot2)
library(viridis)
library(scales)
library(ggpubr)


###	1. Original data, with z-transform

#Get and tidy data
df <- as.data.frame(read_excel("~/Desktop/Data_Maxwell2012.xlsx",sheet=1,na="NA")); rownames(df) <- df[,1]; df <- df[,-1] ; df$Species <- as.factor(df$Species)

species_colors <- viridis(n=8)[c(6,8,4,2)];names(species_colors) <- levels(df$Species)


#Log measurement data
df[,2:10] <- log(df[,2:10])

#Scale measurement data
df[,2:10] <- scale(df[,2:10])

#PCoA, full dataset
dist <- dist(df[,2:10])
pcoa <- pcoa(dist)
to_plot <- data.frame(Species=df$Species,PCoA1=pcoa$vectors[,1],PCoA2=pcoa$vectors[,2])

split_species <- split(to_plot[,c("PCoA1","PCoA2","Species")], to_plot$Species)
cHull_species <- lapply(split_species, function(df){
  df[chull(df), ]  # chull computes convex hulls for each subgroup.
})
cHull_species <- do.call(rbind, cHull_species) #combines the coordinates of the points (i.e. taxa) needed to draw the convex hull (i.e. only the extremes ones)

#Plot
plot_pcoa_original <- ggplot(data=to_plot,aes(x=PCoA1, y=PCoA2))+
	geom_polygon(data=cHull_species,aes(fill=Species),alpha=1/2)+
	geom_point(aes(color=Species))+
	scale_fill_manual(values= species_colors)+
	scale_colour_manual(values= species_colors)+
	labs(x=paste("PCoA1 ","(",100*(round(pcoa$values$Rel_corr_eig[1],digits=3)),"%)",sep=""), y=paste("PCoA2 ","(",100*(round(pcoa$values$Rel_corr_eig[2],digits=3)),"%)",sep=""))+
	ggtitle("PCoA, original data + z-transform")

	

#PCA, without pelvic girdle data
df <- df[,1:8]

#PCA of reduced dataset
pca <- prcomp(df[,2:8])
summaryPCA <- summary(pca)
to_plot2 <- data.frame(Species=df$Species,PCA1=pca$x[,1],PCA2=pca$x[,2])

split_species <- split(to_plot2[,c("PCA1","PCA2","Species")], to_plot2$Species)
cHull_species <- lapply(split_species, function(df){
  df[chull(df), ]  # chull computes convex hulls for each subgroup.
})
cHull_species <- do.call(rbind, cHull_species) #combines the coordinates of the points (i.e. taxa) needed to draw the convex hull (i.e. only the extremes ones)

#Plot
plot_pca_original <- ggplot(data=to_plot2,aes(x=PCA1, y=PCA2))+
	geom_polygon(data=cHull_species,aes(fill=Species),alpha=1/2)+
	geom_point(aes(color=Species))+
	scale_fill_manual(values= species_colors)+
	scale_colour_manual(values= species_colors)+
	labs(x=paste("PC1 ","(",100*(round(summaryPCA$importance[2,1],digits=3)),"%)",sep=""), y=paste("PC2 ","(",100*(round(summaryPCA$importance[2,2],digits=3)),"%)",sep=""))+
	ggtitle("PCA, original data (no pelvic) + z-transform")


#############################################@

###	2. Data scaled humerus length + z-transform

#Get and tidy data
df_scaled <- as.data.frame(read_excel("~/Desktop/Data_Maxwell2012.xlsx",sheet=2,na="NA")); rownames(df) <- df_scaled[,1]; df_scaled <- df_scaled[,-1] ; df_scaled$Species <- as.factor(df_scaled$Species)

species_colors <- viridis(n=8)[c(6,8,4,2)];names(species_colors) <- levels(df$Species)


#Log measurement data
df_scaled[,2:9] <- log(df_scaled[,2:9])

#Scale measurement data
df_scaled[,2:9] <- scale(df_scaled[,2:9])

#PCoA of full dataset
dist_scaled <- dist(df_scaled[,2:9])
pcoa_scaled <- pcoa(dist_scaled,correction="cailliez")
to_plot_scaled <- data.frame(Species=df_scaled$Species,PCoA1=pcoa_scaled$vectors[,1],PCoA2=pcoa_scaled$vectors[,2])

split_species_scaled <- split(to_plot_scaled[,c("PCoA1","PCoA2","Species")], to_plot_scaled$Species)
cHull_species_scaled <- lapply(split_species_scaled, function(df){
  df[chull(df), ]  # chull computes convex hulls for each subgroup.
})
cHull_species_scaled <- do.call(rbind, cHull_species_scaled) #combines the coordinates of the points (i.e. taxa) needed to draw the convex hull (i.e. only the extremes ones)

#Plot
plot_pcoa_scaled <- ggplot(data=to_plot_scaled,aes(x=PCoA1, y=PCoA2))+
	geom_polygon(data=cHull_species_scaled,aes(fill=Species),alpha=1/2)+
	geom_point(aes(color=Species))+
	scale_fill_manual(values= species_colors)+
	scale_colour_manual(values= species_colors)+
	labs(x=paste("PCoA1 ","(",100*(round(pcoa_scaled$values$Rel_corr_eig[1],digits=3)),"%)",sep=""), y=paste("PCoA2 ","(",100*(round(pcoa_scaled$values$Rel_corr_eig[2],digits=3)),"%)",sep=""))+
	ggtitle("PCoA, scaled data + z-transform")
	

#PCA, without pelvic girdle measurements

pca_scaled <- prcomp(df_scaled[,2:7],center=FALSE)
summaryPCA_scaled <- summary(pca_scaled)
to_plot_scaled2 <- data.frame(Species=df_scaled$Species,PCA1=pca_scaled$x[,1],PCA2=pca_scaled$x[,2])

split_species_scaled2 <- split(to_plot_scaled2[,c("PCA1","PCA2","Species")], to_plot_scaled2$Species)
cHull_species_scaled2 <- lapply(split_species_scaled2, function(df){
  df[chull(df), ]  # chull computes convex hulls for each subgroup.
})
cHull_species_scaled2 <- do.call(rbind, cHull_species_scaled2) #combines the coordinates of the points (i.e. taxa) needed to draw the convex hull (i.e. only the extremes ones)

#Plot
plot_pca_scaled <- ggplot(data=to_plot_scaled2,aes(x=PCA1, y=PCA2))+
	geom_polygon(data=cHull_species_scaled2,aes(fill=Species),alpha=1/2)+
	geom_point(aes(color=Species))+
	scale_fill_manual(values= species_colors)+
	scale_colour_manual(values= species_colors)+
	labs(x=paste("PC1 ","(",100*(round(summaryPCA_scaled$importance[2,1],digits=3)),"%)",sep=""), y=paste("PC2 ","(",100*(round(summaryPCA_scaled$importance[2,2],digits=3)),"%)",sep=""))+
	ggtitle("PCA, scaled data (no pelvic) + z-transform")


ggarrange(plot_pcoa_original, plot_pca_original, plot_pcoa_scaled, plot_pca_scaled,common.legend=TRUE, legend="bottom",nrow=2,ncol=2)
ggsave("FigS1.pdf",device="pdf",units="mm",width=240,height=180)