###############################################################################################################################
### R script for ordination methods, convergence statistics, phylogenetic signal, macroevolutionary landscapes, disparity   ###
### Written by V Fischer, 2018-2020                                                 		          						###
### Related paper: Fischer et al. 2020. The macroevolutionary landscape of short-necked plesiosaurians. Scientific Reports. ###
###############################################################################################################################

#=====================================================#
#===================== Libraries =====================#
#=====================================================#

#Data treatment
library(ape)
library(Claddis)
library(vegan)
library(paleotree)
library(pvclust)
library(dispRity)
library(MASS)
library(psych)
library(phytools)
library(adephylo)
library(convevol)
library(raster)
library(RRphylo)

#Data viz
library(ggplot2)
library(ggdendro)
library(ggrepel)
library(dendextend)
library(plotly)
library(gridExtra)
library(htmlwidgets)
library(orca) #used to extract .pdf from plotly; might requires installation of command-line software

#The red and blue colors are taken from the "Life Aquatic" palette in Wes Anderson palettes: https://wesandersonpalettes.tumblr.com 

#====================================================================#
#===================== Set options and get data =====================#
#====================================================================#

#Set "region" to choose the body region to be used. Use "#" in front of unselected region.
region <- "whole"
region <- "craniodental"

#Select the ordination method that will be used to compute the morphospaces. Use "#" in front of unselected method
ordin_method <- "PCoA"
ordin_method <- "NMDS"

#Define the completeness threshold
threshold <- 0.45

#Get data and apply z-transform
raw_data <- as.data.frame(read.csv("Supp_file_bauplan_data.csv",na.strings="NA"))
#for testing:
#raw_data <- as.data.frame(read.csv("~/Dropbox (EDDy Lab)/A_Research/A_Projects/A_Cretaceous/Plesio Landscape/bauplan_data2.csv",na.strings="NA"))
#

bauplan_data_full <- data.frame(raw_data[,-1],row.names=raw_data[,1])

#taxa <- bauplan_data_full$Taxon; bauplan_data_full <- bauplan_data_full[,-1]; rownames(bauplan_data_full) <- taxa

if (region== "craniodental"){
  bauplan_data <- bauplan_data_full[,3:10] #restricting the dataset to craniodental features
  } else {
    bauplan_data <- bauplan_data_full[,-c(1:2)] #if not, just remove the column about clades
  } 
mandible_data <- bauplan_data_full[,2]; names(mandible_data) <- rownames(bauplan_data_full)

#Lists of taxa for each major clade
Pliosauridae <- rownames(bauplan_data_full[bauplan_data_full$Clade!="Polycotylidae",])
Thalassophonea <- rownames(bauplan_data_full[bauplan_data_full$Clade=="Thalassophonea",])
Polycotylidae <- rownames(bauplan_data_full[bauplan_data_full$Clade=="Polycotylidae",])

#Lists of taxa within each temporal bin (Early Jurassic, Mid Jurassic, Late Jurassic, Early Creatceous, Late Cretaceous)
time_windows <- list(EJ = c("Thalassiodracon_hawkinsii","Hauffiosaurus_longirostris","Hauffiosaurus_tomistomimus","Hauffiosaurus_zanoni"), MJ = c("Marmornectes_candrewi","Peloneustes_philarchus","Simolestes_vorax","Pliosaurus_andrewsi","Liopleurodon_ferox","Anguanax_zignoi"), LJ = c("Pliosaurus_funkei","Pliosaurus_westburyensis","Pliosaurus_carpenteri","Pliosaurus_brachydeirus","Pliosaurus_macromerus","Pliosaurus_brachyspondylus","Pliosaurus_cf_kevani","Pliosaurus_kevani","Gallardosaurus_iturraldei","Pliosaurus_rossicus","Pliosaurus_irgisensis","Pliosaurus_almanzaensis","Pliosaurus_patagonicus"), EK = c("Kronosaurus_queenslandicus","QM_F51291","Makhaira_rossica","Luskhan_itilensis","Acostasaurus_pavachoquensis","Sachicasaurus_vitae","Kronosaurus_boyacensis","Edgarosaurus_muddi","Richmond_pliosaur"), LK = c("Polyptychodon_sp","Megacephalosaurus_eulerti","Brachauchenius_lucasi","Brachauchenius_MNA","Plesiopleurodon_wellesi","Palmulasaurus_quadratus","Pahasapasaurus_haasi","Polycotylus_latipinnis","Thililua_longicollis","Trinacromerum_bentonianum","Dolichorhynchops_osborni","Dolichorhynchops_bonneri","Sulcusuchus_erraini","Mauriciosaurus_fernandezi","Dolichorhynchops_tropicensis","Eopolycotylus_rankini","Georgiasaurus_penzensis","Dolichorhynchops_ROM_29010","Dolichorhynchops_herschelensis","Manemergus_anguirostris"))

#Amount of missing Bauplan data
paste("Amount of missing data:",100*(sum(is.na(bauplan_data))/(dim(bauplan_data)[1]*dim(bauplan_data)[2])),"%")

#Completeness per Bauplan feature
100*(colMeans(!is.na(bauplan_data)))

#Amount of missing Bauplan data after application of completeness threshold per taxa
completeness <- rowMeans(!is.na(bauplan_data))
paste("Amount of missing data:",100*(sum(is.na(bauplan_data[completeness>threshold,]))/(dim(bauplan_data[completeness>threshold,])[1]*dim(bauplan_data[completeness>threshold,])[2])),"%")

#Scale the data (z-transform) foe each taxa passing the completeness threshold and compute the euclidean distance matrix
data.s <- scale(bauplan_data[completeness>threshold,]) 
d <- dist(data.s)

#Get phylogenetic data and drop tip according to the taxa passing the completeness threshold
trees <- read.nexus("Supp_file_trees.nex") 
consensus_tree <- consensus(trees)
ranges <- read.table("Supp_file_taxon_ages.txt",header=TRUE)
random <- sample(1:length(trees),1)
tree <- drop.tip(trees[[random]],trees[[random]]$tip.label[!trees[[random]]$tip.label %in% rownames(data.s)])
tree_eq <- timePaleoPhy(tree,timeData=ranges[rownames(ranges) %in% rownames(data.s),],type= "equal",vartime = 3)
tree_mbl <- timePaleoPhy(tree,timeData=ranges[rownames(ranges) %in% rownames(data.s),],type= "mbl",vartime=3)

#How well resolved is the strict consensus of pruned trees? (takes a long time)
if (region== "craniodental"){
pruned_trees <- c()
for (i in 1:length(trees)){
  pruned_trees[[i]] <- drop.tip(trees[[i]],trees[[i]]$tip.label[!trees[[i]]$tip.label %in% rownames(data.s)])
}
consensus_tree_pruned <- timePaleoPhy(consensus(pruned_trees),timeData=ranges[rownames(ranges) %in% rownames(data.s),],type= "equal",vartime = 3)
pdf("Consensus_pruned.pdf",width=8, height=5)
geoscalePhylo(ladderize(consensus_tree_pruned),ranges,cex.ts=0.5,cex.tip=0.6)
dev.off()
}


#==============================================================#
#===================== Cluster dendrogram =====================#
#==============================================================#

#Cluster dendrogram
clust_result <- pvclust(as.data.frame(t(data.s)), method.dist="eu", method.hclust="ward.D2", use.cor="pairwise.complete.obs", nboot=1000, r=seq(0.5,10,by=.5)) # "euclidean distances" for method distance and "ward" for clustering method r= different size of data subsets 50% to 1000%

pdf(paste("Robustness cluster_",region,".pdf",sep=""), width=12, height=12)
plot(clust_result)
dev.off()

#Number of groups you want from the cluster dendrogram results
cut <- 2

#Add cladistic and cluster appartenance information to the dendrogram data
dendro <- dendro_data(clust_result$hclust,type="rectangle")
cut_result <- cutree(clust_result$hclust,k=cut)
clust_df <- data.frame(label=names(cut_result), cluster=factor(cut_result))
dendro$labels <- merge(dendro$labels,clust_df, by="label")
dendro$labels$clade[dendro$labels$label %in% Polycotylidae] <- "Polycotylidae"
dendro$labels$clade[dendro$labels$label %in% Pliosauridae] <- "Pliosauridae"
dendro$labels$clade[dendro$labels$label %in% Thalassophonea] <- "Thalassophonea"

#Use ggplot2 to visualise the cluster dendrogram
clade_colors <- c('#f54b1a','#01abe9', '#ff952c'); names(clade_colors) <- levels(as.factor(dendro$labels$clade))
col_clades <- scale_colour_manual(name = "clade",values = clade_colors)
cluster_colors <- c('#23CC8F','#FF0D68'); names(cluster_colors) <-  c("Longirostrines","Latirostrines")
col_clusters <- scale_colour_manual(name = "cluster",values = cluster_colors)
fill_clusters <- scale_fill_manual(name = "cluster",values = cluster_colors)

cluster <- ggplot()+
  geom_segment(data=segment(dendro),aes(x=x, y=y, xend=xend, yend=yend))+
  geom_text(data=label(dendro), aes(x, y, label=gsub("[_]"," ",label(dendro)$label),fontface="italic", hjust=0,color=clade), size=2.5) +
  annotate(geom="text", x=c(17,5), y=11.5, label=c("Longirostrines","Latirostrines"),color=cluster_colors,size=3,hjust=0,angle=45)+
  annotate(geom="text", x=c(30,29,28), y=11.5, label=c("Early pliosaurids","Thalassophonea","Polycotylidae"),color=clade_colors[c(1,3,2)],size=3,hjust=0)+
  col_clades+
  labs(y="",x="")+
  expand_limits(y=c(1.1*max(clust_result$hclust$height),-0.4*max(clust_result$hclust$height)))+
  scale_x_continuous(breaks=NULL)+
  coord_flip() +
  scale_y_reverse(breaks=c(max(dendro$segment$yend),0),minor_breaks=NULL)+
  theme(legend.position="none", panel.background = element_blank())
cluster

#explore the distribution of Bauplan variables per major cluster with psych package
if (region== "craniodental"){
pdf(paste("FigS1_correlations_",region,"_groupings.pdf",sep=""),width=10, height=10)
pairs.panels(bauplan_data_full[completeness>threshold,-c(1:2)],method="pearson",bg=cluster_colors[cut_result],pch=21,lm=TRUE,hist.col="grey",density=FALSE,ellipses=FALSE,stars=TRUE)
dev.off()
}


#=====================================================#
#===================== PERMANOVA =====================#
#=====================================================#

#test the statistical significance of the clusters (the groups were defined above using the "cut" variable and cutree function)
permanova <- adonis2(dist(data.s)~cut_result,data=as.data.frame(cut_result),permutations=1000)
permanova
write.table(permanova, file=paste("permanova_clusters_",region,".csv",sep=""))

#=====================================================================#
#===================== Morphospaces: PCoA & NMDS =====================#
#=====================================================================#

#PCoA
pcoa <- pcoa(d,correction="cailliez")

#NMDS
nmds <- metaMDS(d,k=2,trymax=100) #2 dimensions
pdf(paste("NMDS_stressplot_",region,".pdf",sep=""), width=7, height=7)
stressplot(nmds)
dev.off()

#Standardise data and temporarily homogenise axes labels and add phylogenetic information
if (ordin_method=="PCoA"){
  to_plot <- as.data.frame(pcoa$vectors.cor[,1:2])
} else if (ordin_method=="NMDS"){
  to_plot <- as.data.frame(nmds$points)
} else {
  print("error in ordination method used")
}

colnames(to_plot) <- c("Axis1","Axis2")
to_plot$Mandible <- mandible_data[completeness>threshold]
to_plot$clade <- bauplan_data_full[completeness>threshold,"Clade"]

#Draw convex hulls using the Split-Apply-Combine strategy
splitData <- split(to_plot[,c(1,2,4)], to_plot$clade)
cHull <- lapply(splitData, function(df){
  df[chull(df), ]  # chull computes convex hulls for each subgroup.
})
cHull <- do.call(rbind, cHull) #combines the coordinates of the points (i.e. taxa) needed to draw the convex hull (i.e. only the extremes ones)

#Visualise the morphospace using ggplot
eco_fig <- ggplot(data = to_plot,aes(x = Axis1, y = Axis2))+
  geom_polygon(data = cHull,aes(x = Axis1, y = Axis2, fill = clade),alpha = 1/2) +
  scale_fill_manual(values=clade_colors)+
  geom_point(aes(color=clade,size=log(Mandible))) +
  scale_color_manual(values=clade_colors)+
  geom_text(aes(label=rownames(to_plot),color=clade))+
  coord_fixed(ratio=1,xlim=c(min(to_plot[,1])*1.3,max(to_plot[,1])*1.3),ylim=c(min(to_plot[,2])*1.3,max(to_plot[,2])*1.3))
if(ordin_method=="PCoA"){
  eco_fig <- eco_fig+labs(title = paste("Bauplan data –",ordin_method,region), x = paste("Principal Coordinates", colnames(to_plot)[1],"(",100*(round(pcoa$values[1,"Rel_corr_eig"],digits=3)),"%)"), y = paste("Principal Coordinates", colnames(to_plot)[2],"(",100*(round(pcoa$values[2,"Rel_corr_eig"],digits=3)),"%)"))
} else if (ordin_method=="NMDS"){
  eco_fig <- eco_fig+labs(title = paste("Bauplan data –",ordin_method,region), x = "NMDS1", y = "NMDS2")
} else {
  print("wrong ordination method used")
}
eco_fig
ggsave(paste("Morphospace_",ordin_method,"_",region,".pdf",sep=""),device="pdf",units="mm",width=240,height=240)


#===============================================================#
#===================== Phylogenetic signal =====================#
#===============================================================#

#Tanglegram
comp_clust <- as.dendrogram(clust_result$hclust) #extract cluster from pvclust result
comp_tree_eq <- as.dendrogram(force.ultrametric(tree = tree_eq,method="extend")) #force a phylogenetic tree as a dendrogram, which requires an ultrametric tree
if(region=="craniodental"){
  pdf(paste("Tanglegram_",region,".pdf",sep=""), width=12, height=7)
  par(mfrow=c(2,1))
  tanglegram(comp_tree_eq,comp_clust, fast = TRUE, margin_inner = 13,main_left="Phylogenetic tree",main_right="Cluster dendrogram",axes=FALSE,cex_main=1.5,color_lines = c(rep(clade_colors[[1]],2),rep(clade_colors[[3]],16),rep(clade_colors[[2]],11))) #creates tanglegram
  dev.off()
}

#Mantel test: phenotypic distance VS phylogenetic distance
phylo_dist_eq <- distTips(tree_eq,method="patristic")
phylo_dist_mbl <- distTips(tree_mbl,method="patristic")

mantel_phenoVSphylo_eq <- mantel(xdis=d,ydis=phylo_dist_eq,permutations=1000)
mantel_phenoVSphylo_mbl <- mantel(xdis=d,ydis=phylo_dist_mbl,permutations=1000)

mantel_table <- rbind(cbind(statistic=mantel_phenoVSphylo_eq$statistic,pval=mantel_phenoVSphylo_eq$signif),cbind(statistic=mantel_phenoVSphylo_mbl$statistic,pval=mantel_phenoVSphylo_mbl$signif));rownames(mantel_table) <- c("equal","mbl")

write.table(mantel_table,file=paste("mantel_phenoVSphylo_",region,".csv",sep=""))

#Simple mantel test: phenotypic distance VS temporal distance
time_dist <- dist(rowMeans(ranges[rownames(ranges) %in% rownames(data.s),]))
mantel_phenoVStime <- mantel(xdis=d,ydis=time_dist,permutations=1000)

write.table(cbind(statistic=mantel_phenoVStime$statistic,pval=mantel_phenoVStime$signif),file=paste("mantel_phenoVStime_",region,".csv",sep=""))


#===================================================================================================#
#===================== Mesures of convergence (Stayton 2015; convevol package) =====================#
#===================================================================================================#

nsim <- 1000
plot(pcoa$values$Rel_corr_eig,type="b") #plot relative corrected eigenvalues

# Putative interclade convergence within the longirostrine cluster: Luskhan, Marmornectes, Trinacromerum
convtips_longirostrines <- c("Marmornectes_candrewi","Luskhan_itilensis","Trinacromerum_bentonianum")

# Putative interclade convergence within the latirostrine cluster: Plesiopleurodon, Acostasaurus
convtips_latirostrines <- c("Plesiopleurodon_wellesi","Acostasaurus_pavachoquensis")


## "Equal" time-scaled tree (might take a long time)

#C1,C2,C3,C4 across the first two axes of PCoA - longirostrines
S_metrics_longirostrines_2 <- convratsig(tree_eq,pcoa$vectors.cor[,1:2],convtips_longirostrines,nsim=nsim)

#C1,C2,C3,C4 accross the first five axes of the PCoA - longirostrines
S_metrics_longirostrines_5 <- convratsig(tree_eq,pcoa$vectors.cor[,1:5],convtips_longirostrines,nsim=nsim)

#C1,C2,C3,C4 accross all the axes of the PCoA - longirostrines
S_metrics_longirostrines_all <- convratsig(tree_eq,pcoa$vectors.cor,convtips_longirostrines,nsim=nsim)

#C1,C2,C3,C4 across the first two axes of PCoA - latirostrines
S_metrics_latirostrines_2 <- convratsig(tree_eq,pcoa$vectors.cor[,1:2], convtips_latirostrines,nsim=nsim)

#C1,C2,C3,C4 across the first five axes of PCoA - latirostrines
S_metrics_latirostrines_5 <- convratsig(tree_eq,pcoa$vectors.cor[,1:5], convtips_latirostrines,nsim=nsim)

#C1,C2,C3,C4 across all axes of the PCoA - latirostrines
S_metrics_latirostrines_all <- convratsig(tree_eq,pcoa$vectors.cor, convtips_latirostrines,nsim=nsim)

#C stats table
c_table_eq <- rbind(cbind(S_metrics_longirostrines_2$ObservedCs[1],S_metrics_longirostrines_2$Pvals[1],S_metrics_longirostrines_2$ObservedCs[2],S_metrics_longirostrines_2$Pvals[2],S_metrics_longirostrines_2$ObservedCs[3],S_metrics_longirostrines_2$Pvals[3],S_metrics_longirostrines_2$ObservedCs[4],S_metrics_longirostrines_2$Pvals[4]),cbind(S_metrics_longirostrines_5$ObservedCs[1],S_metrics_longirostrines_5$Pvals[1],S_metrics_longirostrines_5$ObservedCs[2],S_metrics_longirostrines_5$Pvals[2],S_metrics_longirostrines_5$ObservedCs[3],S_metrics_longirostrines_5$Pvals[3],S_metrics_longirostrines_5$ObservedCs[4],S_metrics_longirostrines_5$Pvals[4]),cbind(S_metrics_longirostrines_all$ObservedCs[1],S_metrics_longirostrines_all$Pvals[1],S_metrics_longirostrines_all$ObservedCs[2],S_metrics_longirostrines_all$Pvals[2],S_metrics_longirostrines_all$ObservedCs[3],S_metrics_longirostrines_all$Pvals[3],S_metrics_longirostrines_all$ObservedCs[4],S_metrics_longirostrines_all$Pvals[4]),cbind(S_metrics_latirostrines_2$ObservedCs[1],S_metrics_latirostrines_2$Pvals[1],S_metrics_latirostrines_2$ObservedCs[2],S_metrics_latirostrines_2$Pvals[2],S_metrics_latirostrines_2$ObservedCs[3],S_metrics_latirostrines_2$Pvals[3],S_metrics_latirostrines_2$ObservedCs[4],S_metrics_latirostrines_2$Pvals[4]),cbind(S_metrics_latirostrines_5$ObservedCs[1],S_metrics_latirostrines_5$Pvals[1],S_metrics_latirostrines_5$ObservedCs[2],S_metrics_latirostrines_5$Pvals[2],S_metrics_latirostrines_5$ObservedCs[3],S_metrics_latirostrines_5$Pvals[3],S_metrics_latirostrines_5$ObservedCs[4],S_metrics_latirostrines_5$Pvals[4]),cbind(S_metrics_latirostrines_all$ObservedCs[1],S_metrics_latirostrines_all$Pvals[1],S_metrics_latirostrines_all$ObservedCs[2],S_metrics_latirostrines_all$Pvals[2],S_metrics_latirostrines_all$ObservedCs[3],S_metrics_latirostrines_all$Pvals[3],S_metrics_latirostrines_all$ObservedCs[4],S_metrics_latirostrines_all$Pvals[4]))
colnames(c_table_eq) <- c("C1","p-value","C2","p-value","C3","p-value","C4","p-value") ; rownames(c_table_eq) <- c("Longi_PCo1-2_eq","Longi_PCo1-5_eq","Longi_All_eq","Lati_PCo1-2_eq","Lati_PCo1-5_eq","Lati_All_eq")

## "Mbl" time-scaled tree (might take a long time)

#C1,C2,C3,C4 across the first two axes of PCoA - longirostrines
S_metrics_longirostrines_2 <- convratsig(tree_mbl,pcoa$vectors.cor[,1:2],convtips_longirostrines,nsim=nsim)

#C1,C2,C3,C4 accross the first five axes of the PCoA - longirostrines
S_metrics_longirostrines_5 <- convratsig(tree_mbl,pcoa$vectors.cor[,1:5],convtips_longirostrines,nsim=nsim)

#C1,C2,C3,C4 accross all the axes of the PCoA - longirostrines
S_metrics_longirostrines_all <- convratsig(tree_mbl,pcoa$vectors.cor,convtips_longirostrines,nsim=nsim)

#C1,C2,C3,C4 across the first two axes of PCoA - latirostrines
S_metrics_latirostrines_2 <- convratsig(tree_mbl,pcoa$vectors.cor[,1:2], convtips_latirostrines,nsim=nsim)

#C1,C2,C3,C4 across the first five axes of PCoA - latirostrines
S_metrics_latirostrines_5 <- convratsig(tree_mbl,pcoa$vectors.cor[,1:5], convtips_latirostrines,nsim=nsim)

#C1,C2,C3,C4 across all axes of the PCoA - latirostrines
S_metrics_latirostrines_all <- convratsig(tree_mbl,pcoa$vectors.cor, convtips_latirostrines,nsim=nsim)

#C stats table
c_table_mbl <- rbind(cbind(S_metrics_longirostrines_2$ObservedCs[1],S_metrics_longirostrines_2$Pvals[1],S_metrics_longirostrines_2$ObservedCs[2],S_metrics_longirostrines_2$Pvals[2],S_metrics_longirostrines_2$ObservedCs[3],S_metrics_longirostrines_2$Pvals[3],S_metrics_longirostrines_2$ObservedCs[4],S_metrics_longirostrines_2$Pvals[4]),cbind(S_metrics_longirostrines_5$ObservedCs[1],S_metrics_longirostrines_5$Pvals[1],S_metrics_longirostrines_5$ObservedCs[2],S_metrics_longirostrines_5$Pvals[2],S_metrics_longirostrines_5$ObservedCs[3],S_metrics_longirostrines_5$Pvals[3],S_metrics_longirostrines_5$ObservedCs[4],S_metrics_longirostrines_5$Pvals[4]),cbind(S_metrics_longirostrines_all$ObservedCs[1],S_metrics_longirostrines_all$Pvals[1],S_metrics_longirostrines_all$ObservedCs[2],S_metrics_longirostrines_all$Pvals[2],S_metrics_longirostrines_all$ObservedCs[3],S_metrics_longirostrines_all$Pvals[3],S_metrics_longirostrines_all$ObservedCs[4],S_metrics_longirostrines_all$Pvals[4]),cbind(S_metrics_latirostrines_2$ObservedCs[1],S_metrics_latirostrines_2$Pvals[1],S_metrics_latirostrines_2$ObservedCs[2],S_metrics_latirostrines_2$Pvals[2],S_metrics_latirostrines_2$ObservedCs[3],S_metrics_latirostrines_2$Pvals[3],S_metrics_latirostrines_2$ObservedCs[4],S_metrics_latirostrines_2$Pvals[4]),cbind(S_metrics_latirostrines_5$ObservedCs[1],S_metrics_latirostrines_5$Pvals[1],S_metrics_latirostrines_5$ObservedCs[2],S_metrics_latirostrines_5$Pvals[2],S_metrics_latirostrines_5$ObservedCs[3],S_metrics_latirostrines_5$Pvals[3],S_metrics_latirostrines_5$ObservedCs[4],S_metrics_latirostrines_5$Pvals[4]),cbind(S_metrics_latirostrines_all$ObservedCs[1],S_metrics_latirostrines_all$Pvals[1],S_metrics_latirostrines_all$ObservedCs[2],S_metrics_latirostrines_all$Pvals[2],S_metrics_latirostrines_all$ObservedCs[3],S_metrics_latirostrines_all$Pvals[3],S_metrics_latirostrines_all$ObservedCs[4],S_metrics_latirostrines_all$Pvals[4]))
rownames(c_table_mbl) <- c("Longi_PCo1-2_mbl","Longi_PCo1-5_mbl","Longi_All_mbl","Lati_PCo1-2_mbl","Lati_PCo1-5_mbl","Lati_All_mbl")
write.table(rbind(c_table_eq,c_table_mbl),file=paste("Table1_Stayton_Cs_",region,".csv",sep=""))


#============================================================================================================#
#===================== Mesures of convergence (Castiglione et al 2019; RRphylo package) =====================#
#============================================================================================================#

## "Equal" timescaled tree

#Longirostrines (Luskhan itilensis + Marmonectes candrewi + Trinacromerum bentonianum)
states <- rep("nostate",length(tree_eq$tip.label));names(states) <- tree_eq$tip.label #create state vector
states[convtips_longirostrines] <- "longirostrines"
conv_longi_2 <- search.conv(tree=tree_eq, y=pcoa$vectors.cor[,1:2], state=states,declust=TRUE, foldername = "Castiglione_longi_eq/2/")
conv_longi_5 <- search.conv(tree=tree_eq, y=pcoa$vectors.cor[,1:5], state=states,declust=TRUE, foldername = "Castiglione_longi_eq/5/")
conv_longi_all <- search.conv(tree=tree_eq, y=pcoa$vectors.cor, state=states,declust=TRUE, foldername = "Castiglione_longi_eq/all/")


#latirostrines (Plesiopleurodon wellesi + Acostasaurus pavachoquensis)
states <- rep("nostate",length(tree_eq$tip.label));names(states) <- tree_eq$tip.label #create state vector
states[convtips_latirostrines] <- "latirostrines"
conv_lati_2 <- search.conv(tree=tree_eq, y=pcoa$vectors.cor[,1:2], state=states,declust=TRUE, foldername = "Castiglione_lati_eq/2/")
conv_lati_5 <- search.conv(tree=tree_eq, y=pcoa$vectors.cor[,1:5], state=states,declust=TRUE, foldername = "Castiglione_lati_eq/5/")
conv_lati_all <- search.conv(tree=tree_eq, y=pcoa$vectors.cor, state=states,declust=TRUE, foldername = "Castiglione_lati_eq/all/")

Castiglione_eq <- rbind(conv_longi_2,conv_longi_5,conv_longi_all,conv_lati_2,conv_lati_5,conv_lati_all)
rownames(Castiglione_eq) <- rownames(c_table_eq)

## "Mbl" timescaled tree

#Longirostrines (Luskhan itilensis + Marmonectes candrewi + Trinacromerum bentonianum)
states <- rep("nostate",length(tree_eq$tip.label));names(states) <- tree_mbl$tip.label #create state vector
states[convtips_longirostrines] <- "longirostrines"
conv_longi_2 <- search.conv(tree=tree_mbl, y=pcoa$vectors.cor[,1:2], state=states,declust=TRUE, foldername = "Castiglione_longi_mbl/2/")
conv_longi_5 <- search.conv(tree=tree_mbl, y=pcoa$vectors.cor[,1:5], state=states,declust=TRUE, foldername = "Castiglione_longi_mbl/5/")
conv_longi_all <- search.conv(tree=tree_mbl, y=pcoa$vectors.cor, state=states,declust=TRUE, foldername = "Castiglione_longi_mbl/all/")


#latirostrines (Plesiopleurodon wellesi + Acostasaurus pavachoquensis)
states <- rep("nostate",length(tree_eq$tip.label));names(states) <- tree_mbl$tip.label #create state vector
states[convtips_latirostrines] <- "latirostrines"
conv_lati_2 <- search.conv(tree=tree_mbl, y=pcoa$vectors.cor[,1:2], state=states,declust=TRUE, foldername = "Castiglione_lati_mbl/2/")
conv_lati_5 <- search.conv(tree=tree_mbl, y=pcoa$vectors.cor[,1:5], state=states,declust=TRUE, foldername = "Castiglione_lati_mbl/5/")
conv_lati_all <- search.conv(tree=tree_mbl, y=pcoa$vectors.cor, state=states,declust=TRUE, foldername = "Castiglione_lati_mbl/all/")

Castiglione_mbl <- rbind(conv_longi_2,conv_longi_5,conv_longi_all,conv_lati_2,conv_lati_5,conv_lati_all)
rownames(Castiglione_mbl) <- rownames(c_table_mbl)
write.table(rbind(Castiglione_eq,Castiglione_mbl),file=paste("Table2_Castiglione_",region,".csv",sep=""))


#==========================================================================================================#
#===================== Differences of disparity between groups using dispRity package =====================#
#==========================================================================================================#

## 1: Per clade ##

#Sampled taxa per clade (Thalassophonea and Polycotylidae)
clades <- list(as.character(dendro$labels[dendro$labels$clade=="Polycotylidae",1]),as.character(dendro$labels[dendro$labels$clade=="Thalassophonea",1]))
names(clades) <- c("Polycotylidae","Thalassophonea")

#Create subsets (split), bootstrap each subset, compute disparity, and test for differences using Wilcoxon test
subsets <- custom.subsets(data=pcoa$vectors.cor,group=clades)
number_bootstraps <- 1000
bootstraps <- boot.matrix(subsets,bootstraps = number_bootstraps)
disparity_clades <- dispRity(bootstraps,metric=c(sum,variances))
plot(disparity_clades)
test.dispRity(disparity_clades,test=wilcox.test)

#Create a data frame extracting the bootstraped data from dispRity (this will be used for the geom_jitter function in ggplot)
disp_mat <- data.frame(disp=double(),clade=character(),stringsAsFactors =FALSE)
for (i in 1:length(clades)){
  disp_mat[((i*number_bootstraps)-number_bootstraps+1):(i*number_bootstraps),"disp"] <- disparity_clades$disparity[[i]][[2]][1,]
  disp_mat[((i*number_bootstraps)-number_bootstraps+1):(i*number_bootstraps),"clade"] <- names(clades)[i]
}

#Visualise the results using ggplot
disparity_clades <- ggplot(data=disp_mat,aes(clade,disp))+
  geom_jitter(width=0.2,aes(color=clade,alpha=0.5))+
  col_clades+
  geom_boxplot(aes(alpha=0.5),outlier.color = NA)+
  labs(title=paste("Clade disparity,",region,"data,", number_bootstraps,"bootstraps \nWilcoxon test p value=",signif(test.dispRity(disparity_clades,test=wilcox.test)[[2]]$p.value,digits=2)),y="Disparity",x="")+
  theme(legend.position="none",panel.grid.major = element_blank(), panel.grid.minor = element_blank())


## 2: Per major Bauplan cluster ##

#Sampled taxa per cluster identified above via the cutree (k=2)
clusters <- list(as.character(dendro$labels[dendro$labels$cluster==1,1]),as.character(dendro$labels[dendro$labels$cluster==2,1]))
names(clusters) <- c("Longirostrines","Latirostrines")

#Create subsets (split), bootstrap each subset, compute disparity, and test for differences using Wilcoxon test
subsets <- custom.subsets(data=pcoa$vectors.cor,group=clusters)
number_bootstraps <- 1000
bootstraps <- boot.matrix(subsets,bootstraps = number_bootstraps)
disparity_clusters <- dispRity(bootstraps,metric=c(sum,variances))
plot(disparity_clusters)
test.dispRity(disparity_clusters,test=wilcox.test)

#Create a data frame extracting the bootstraped data from dispRity (this will be used for the geom_jitter function in ggplot)
disp_mat <- data.frame(disp=double(),cluster=character(),stringsAsFactors =FALSE)
for (i in 1:length(clusters)){
  disp_mat[((i*number_bootstraps)-number_bootstraps+1):(i*number_bootstraps),"disp"] <- disparity_clusters$disparity[[i]][[2]][1,]
  disp_mat[((i*number_bootstraps)-number_bootstraps+1):(i*number_bootstraps),"cluster"] <- names(clusters)[i]
}

#Visualise the results using ggplot
disparity_clusters <- ggplot(data=disp_mat,aes(cluster,disp))+
  geom_jitter(width=0.2,aes(color=cluster,alpha=0.5))+
  col_clusters+
  scale_x_discrete(limits = (unique(disp_mat$cluster)))+
  geom_boxplot(aes(alpha=0.5),outlier.color = NA)+
  labs(title=paste("Major cluster disparity,",region,"data,", number_bootstraps,"bootstraps \nWilcoxon test p value=",signif(test.dispRity(disparity_clusters,test=wilcox.test)[[2]]$p.value,digits=2)),y="Disparity",x="")+
  theme(legend.position="none",panel.grid.major = element_blank(), panel.grid.minor = element_blank())

## Mandible size differences
mandi_clust <- as.data.frame(cbind(log(mandible_data[completeness>threshold]),cut_result)) ; colnames(mandi_clust) <- c("mandible_length","cluster") ; mandi_clust[mandi_clust$cluster==1,2] <- "Longirostrines" ; mandi_clust[mandi_clust$cluster==2,2] <- "Latirostrines" ; mandi_clust$cluster <- as.factor(mandi_clust$cluster)
mandible_distrib <- ggplot(mandi_clust, aes(x=mandible_length, fill=cluster,color=cluster)) +
  geom_density(alpha=0.4, position = 'identity') +
  fill_clusters+
  col_clusters+
  geom_vline(xintercept=c(log(1000),log(2000)),linetype="dotted",color = "black")+
  annotate(geom="text", x=c(log(1000),log(2000)), y=1.1, label=c("1m","2m"),color="black")+
  labs(title="Mandible size distribution per major cluster",x="log10 mandible length")+
  theme(legend.position="none",panel.grid.major = element_blank(), panel.grid.minor = element_blank())
mandible_distrib

##  Arrange graphs into Fig 2 (cluster, disparity):
titlefontsize <- 9

if (region=="craniodental"){
  grid.arrange(cluster+labs(title=paste("Cluster dendogram \n",region," data",sep=""),tags="A")+
                 theme(plot.title = element_text(size = titlefontsize)),mandible_distrib+labs(tags="B")+
                 theme(plot.title = element_text(size = titlefontsize)),disparity_clades+labs(tags="C")+
                 theme(plot.title = element_text(size = titlefontsize)),disparity_clusters+labs(tags="D")+
                 theme(plot.title = element_text(size = titlefontsize)),nrow=2,widths=c(1.5,1.5))
  save_fig <- arrangeGrob(cluster+labs(title=paste("Cluster dendrogram \n",region," data",sep=""),tags="A")+
                            theme(plot.title = element_text(size = titlefontsize)),mandible_distrib+labs(tags="B")+
                            theme(plot.title = element_text(size = titlefontsize)),disparity_clades+labs(tags="C")+
                            theme(plot.title = element_text(size = titlefontsize)),disparity_clusters+labs(tags="D")+
                            theme(plot.title = element_text(size = titlefontsize)),nrow=2,widths=c(1.5,1.5))
  ggsave(file=paste("Fig2_clust_disp_",region,".pdf",sep=""),save_fig,device="pdf",units="mm",width=270,height=200)
  
} else if (region=="whole"){
  grid.arrange(cluster+labs(title=paste("Cluster dendrogram \n",region," data",sep=""),tags="A")+
                 theme(plot.title = element_text(size = titlefontsize)),disparity_clades+labs(tags="B")+
                 theme(plot.title = element_text(size = titlefontsize)),disparity_clusters+labs(tags="C")+
                 theme(plot.title = element_text(size = titlefontsize)),nrow=1,widths=c(1.5,1,1))
  save_fig <- arrangeGrob(cluster+labs(title=paste("Cluster dendrogram \n",region," data",sep=""),tags="A")+
                            theme(plot.title = element_text(size = titlefontsize)),disparity_clades+labs(tags="B")+
                            theme(plot.title = element_text(size = titlefontsize)),disparity_clusters+labs(tags="C")+
                            theme(plot.title = element_text(size = titlefontsize)),nrow=1,widths=c(1.5,1,1))
  ggsave(file=paste("FigS2_clust_disp_",region,".pdf",sep=""),save_fig,device="pdf",units="mm",width=340,height=130)
  
} else{
  print("wrong option")
}


#==================================================================================#
#===================== Disparity and overdispersion over time =====================#
#==================================================================================#

## Characterisation of overdispersion by comparing sum of range/sum of variance with
## 1000 bootstraped values sampling the whole dataset, using the same number of taxa per bin

time_windows_pruned <- time_windows
overdisp <- matrix(ncol=length(time_windows),nrow=6)
rownames(overdisp) <- c("SoR","Boot_SoR_mean-2sd","Boot_SoR_mean+2sd","SoV","Boot_SoV_mean-2sd","Boot_SoV_mean+2sd")
random <- 1000
SoR_random <- SoV_random <- matrix(ncol=length(time_windows),nrow=random)
colnames(overdisp) <- colnames(SoR_random) <- colnames(SoV_random) <- names(time_windows)
for (i in 1:length(time_windows)){
  time_windows_pruned[[i]] <- time_windows[[i]][time_windows[[i]] %in% rownames(pcoa$vectors.cor)]
  if (length(time_windows_pruned[[i]])>2){
    bin_taxa <- pcoa$vectors.cor[rownames(pcoa$vectors.cor) %in% time_windows[[i]],]
    overdisp[1,i] <- sum(apply(bin_taxa,2,max)-apply(bin_taxa,2,min)) # sum of ranges along all PCo axes
    overdisp[4,i] <- sum(apply(bin_taxa,2,var)) # sum of variance along all PCo axes
    
    for (j in 1:random){
      sample <- pcoa$vectors.cor[sample(nrow(pcoa$vectors.cor),nrow(bin_taxa)),]
      SoR_random[j,i] <- sum(apply(sample,2,max)-apply(sample,2,min)) # sum of ranges along all PCo axes
      SoV_random[j,i] <- sum(apply(sample,2,var)) # sum of variance along all PCo axes
    }
    overdisp[2,i] <- mean(SoR_random[,i])-2*sd(SoR_random[,i])
    overdisp[3,i] <- mean(SoR_random[,i])+2*sd(SoR_random[,i])
    overdisp[5,i] <- mean(SoV_random[,i])-2*sd(SoV_random[,i])
    overdisp[6,i] <- mean(SoV_random[,i])+2*sd(SoV_random[,i])
  }
}
write.csv(overdisp,file="TableS3_overdispersion.csv")


## Disparity over time with bootstrapped data, using dispRity

time_subsets <- custom.subsets(data=as.matrix(pcoa$vectors.cor),group=time_windows_pruned)
time_bootstraps <- boot.matrix(time_subsets,bootstraps = 1000)
disparity_var <- dispRity(time_bootstraps,metric=c(sum,variances),verbose=TRUE)
write.csv(summary(disparity_var),file=paste("temporal_disparity_",region,".csv",sep=""))
summary <- cbind(1:5,summary(disparity_var)); colnames(summary) <- c("bins","time","n","SoV_obs","SoV_median_boot","min2","min25","max75","max97")

temporal_disp <- ggplot(data=summary)+
  geom_ribbon(aes(x=bins,ymin=min2,ymax=max97),fill="#216472",alpha=0.5)+
  geom_ribbon(aes(x=bins,ymin=min25,ymax=max75),fill="#216472",alpha=0.5)+
  geom_line(aes(x=bins,y=SoV_median_boot),color="#216472",size=3,alpha=1)+
  annotate("text",y=80,x=1,label="NA",alpha=1)+
  ylab("Sum of variance")+
  theme(axis.line.x=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank(),axis.title.x=element_blank(),axis.text.x=element_blank(),axis.ticks.x=element_blank())

#=================================================================================================================#
#===================== Approximation of macroevolutionary landscape using phenotypic density =====================#
#=================================================================================================================#

#change bandwidth (if needed) by a fraction/multiplication of the absolute lengths of axes
factor <- 1
bandw <- factor*c(max(to_plot[,"Axis1"]),max(to_plot[,"Axis2"]))

over_time_MS=function (X){
  bin_taxa <- to_plot[rownames(to_plot) %in% X,]
  graph <- ggplot(data = to_plot,aes(x = Axis1, y = Axis2))+
  stat_density_2d(aes(alpha = ..level..),h=bandw, geom = "polygon")+
  geom_point(data = bin_taxa, aes(x = Axis1, y = Axis2,color=clade,size=4)) +
  col_clades+
  scale_x_continuous(limits = c(min(to_plot[,1])*1.5,max(to_plot[,1])*1.5))+
  scale_y_continuous(limits = c(min(to_plot[,2])*1.5,max(to_plot[,2])*1.5))+
  coord_fixed(ratio=1)+
  theme(legend.position = "none",panel.grid.major = element_blank(), panel.grid.minor = element_blank())
  if(ordin_method=="PCoA"){
    graph <- graph+labs(title = paste("Bauplan data –",ordin_method,region), x = "Principal Coordinate 1", y = "Principal Coordinate 2")
  } else if (ordin_method=="NMDS"){
    graph <- graph+labs(title = paste("Bauplan data –",ordin_method,region), x = "NMDS1", y = "NMDS2")
  } else {
    print("wrong ordination method used")
  }
  return(graph)
}
over_time_plots <- lapply(time_windows,over_time_MS)

## 3D plot a Kernel density estimation with MASS and plot_ly
kd <- kde2d(x=to_plot[,"Axis1"],y=to_plot[,"Axis2"],h= bandw,n=50,lims=1.6*c(min(to_plot[,"Axis1"]),max(to_plot[,"Axis1"]),min(to_plot[,"Axis2"]),max(to_plot[,"Axis2"])))
scene = list(camera = list(eye = list(x = 1, y = 0, z = 1)),tickcolor="white",xaxis=list(title = "",zeroline = FALSE,showline = FALSE,showticklabels = FALSE,showgrid = FALSE),yaxis=list(title = "",zeroline = FALSE,showline = FALSE,showticklabels = FALSE,showgrid = FALSE),zaxis=list(title = "",zeroline = FALSE,showline = FALSE,showticklabels = FALSE,showgrid = FALSE))

interactive_plot <- plot_ly(x = -kd$x, y = kd$y, z = kd$z) %>% add_surface(showscale=FALSE) %>% layout(scene = scene)
interactive_plot

#save as interactive .html and .svg
saveWidget(interactive_plot, paste("HTML_landscape_",region,"_",ordin_method,".html",sep=""))
orca(interactive_plot,paste("temp_landscape_",region,"_",ordin_method,".svg",sep=""))

#find coordinates of maxima
raster <- raster(kd)
max_function <- function(X) max(X, na.rm=TRUE) ## Find the maximum value within the 15*15-cell neighborhood of each cell ; smaller neighborhods will result in more local maxima
moving_matrix <- matrix(1, nrow=15, ncol=15) ## Weight matrix for cells in moving window
localmax <- focal(raster, fun=max_function, w=moving_matrix, pad=TRUE, padValue=NA)
raster2 <- raster==localmax ## Does each cell have the maximum value in its neighborhood?
maxXY <- xyFromCell(raster2, Which(raster2==1, cells=TRUE)) ## Get x-y coordinates of those cells that are local maxima


#=======================================================================================================#
#===================== Phylomorphospace by W. Andrew Barr, modified to add density =====================#
#=======================================================================================================#

source("Supp_file_Phylo_density.R", chdir = TRUE)
to_plot$taxon <- rownames(to_plot)
phylo_scape <- ggphylomorpho(tree=tree_eq,tipinfo=to_plot,names=TRUE,repel=TRUE,xvar=Axis1,yvar=Axis2,factorvar=clade,labelvar=taxon,h=bandw,sizedata=to_plot$Mandible,lim=1.5*max(abs(to_plot[,1:2])))
phylo_scape <- phylo_scape + col_clades + theme(legend.position='right',panel.grid.major = element_blank(), panel.grid.minor = element_blank())
if(ordin_method=="PCoA"){
  phylo_scape <- phylo_scape+labs(x = paste("Principal Coordinates", colnames(to_plot)[1],"(",100*(round(pcoa$values[1,"Rel_corr_eig"],digits=3)),"%)"), y = paste("Principal Coordinates", colnames(to_plot)[2],"(",100*(round(pcoa$values[2,"Rel_corr_eig"],digits=3)),"%)"))
  
} else { phylo_scape <- phylo_scape+labs(x = "NMDS1", y = "NMDS2")
}
#interactive plot:
ggplotly(phylo_scape)


#Arrange  graphs for Fig 3/S12-14:
grid.arrange(phylo_scape+labs(tags="A"),
             over_time_plots[[1]]+labs(title="Early Jurassic",tags="C"),
             over_time_plots[[2]]+labs(title="Mid Jurassic",tags="D"),
             over_time_plots[[3]]+labs(title="Late Jurassic",tags="E"),
             over_time_plots[[4]]+labs(title="Early Cretaceous",tags="F"),
             over_time_plots[[5]]+labs(title="Late Cretaceous",tags="G"),
             temporal_disp+labs(tags="H"),
             heights=c(3,1.6,1),layout_matrix = rbind(c(1,1,1,NA,NA),c(2,3,4,5,6),c(7,7,7,7,7)))
save_fig <- arrangeGrob(phylo_scape+labs(tags="A"),
                        over_time_plots[[1]]+labs(title="Early Jurassic",tags="C"),
                        over_time_plots[[2]]+labs(title="Mid Jurassic",tags="D"),
                        over_time_plots[[3]]+labs(title="Late Jurassic",tags="E"),
                        over_time_plots[[4]]+labs(title="Early Cretaceous",tags="F"),
                        over_time_plots[[5]]+labs(title="Late Cretaceous",tags="G"),
                        temporal_disp+labs(tags="H"),
                        heights=c(3,1.6,1),layout_matrix = rbind(c(1,1,1,NA,NA),c(2,3,4,5,6),c(7,7,7,7,7)))
ggsave(file=paste("MEL_",ordin_method,"_",region,".pdf",sep=""),save_fig,device="pdf",units="mm",width=340,height=340)
