library(vegan)
library(fossil)
library(fpc)
library(ape)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Read data, convert to presence/absence matrix
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

FaunalList <- read.csv ("El-Sayed_et_al_faunal_data.csv")

PresenceAbsence <- create.matrix  (FaunalList, tax.name = "Clade",  locality = "Site")

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Forbes metric from code from:
#
# Brocklehurst, N., Day, M.O., Rubidge, B.S., and Fröbisch, J., 2017, 
# Olson’s Extinction and the latitudinal biodiversity gradient of tetrapods in the Permian: 
# Proceedings of the Royal Society, Series B, vol. 284, p. 20170231.
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Alroy_Forbes<-function(occs)
{
	output<-matrix(nrow=ncol(occs),ncol=ncol(occs))
	rownames(output)<-colnames(occs)
	colnames(output)<-colnames(occs)
	for(i in 1:nrow(output))
	{
		for (j in 1:ncol(output))
		{
			species.A<-rownames(occs)[which(occs[,rownames(output)[i]]!=0)]
			species.B<-rownames(occs)[which(occs[,colnames(output)[j]]!=0)]

			n<-length(unique(c(species.A,species.B)))

			a<-length(intersect(species.A,species.B))
			b<-length(which(species.A%in%species.B==F))
			c<-length(which(species.B%in%species.A==F))
			
			output[i,j]<-1-((a*(n+sqrt(n)))/((a*(n+sqrt(n)))+(3/2*b*c)))
			
		}
	}
	output<-as.dist(output)
	return(output)
}

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Calculate Forbes index for comparisons
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

distances <- Alroy_Forbes(PresenceAbsence)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Export distances as .csv
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

write.csv(as.matrix(distances), file = "El-Sayed_et_al_distances.csv")

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Test for number of clusters in data (over range from 2-10)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

pamk (distances, krange=2:10, usepam=TRUE)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# NMDS
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

NMDS <- metaMDS (distances, k=2, trymax = 1000)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# PCOA with Lingoes correction
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

PCOA <- pcoa (distances, correction = "lingoes")

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Basic plotting
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

ages <- FaunalList[(match(rownames(NMDS$points), FaunalList [, "Site"])), "Age"]

colors <- vector (length = length (ages))

colors [which(ages == "Priabonian")] <- rgb (255/255, 201/255, 153/255, 1)
colors [which(ages == "Bartonian")] <- rgb (255/255, 201/255, 153/255, 1) 
colors [which(ages == "Lutetian")] <- rgb (255/255, 201/255, 153/255, 1) 
colors [which(ages == "Ypresian")] <- rgb (255/255, 201/255, 153/255, 1) 
colors [which(ages == "PETM")] <- rgb(204/255, 0, 0, 1)
colors [which(ages == "Selandian")] <- rgb(239/255, 144/255, 80/255, 1)
colors [which(ages == "Danian")] <- rgb(239/255, 144/255, 80/255, 1)  
colors [which(ages == "Maastrichtian")] <- rgb(166/255, 216/255, 74/255, 1) 

symbols <- vector (length = length (ages))

symbols [which(ages == "Priabonian")] <- 25 
symbols [which(ages == "Bartonian")] <- 25
symbols [which(ages == "Lutetian")] <- 25
symbols [which(ages == "Ypresian")] <- 25
symbols [which(ages == "PETM")] <- 21
symbols [which(ages == "Selandian")] <- 23 
symbols [which(ages == "Danian")] <- 23
symbols [which(ages == "Maastrichtian")] <- 22


#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Output basic .pdf plots for visualization
# (Templates for Fig. 4, Fig. S7)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 

pdf("El-Sayed_et_al_NMDS.pdf", width=3.5, height=3.5)

plot (NMDS$points[,1], NMDS$points[,2], pch = symbols, col = "transparent", bg = colors, cex.axis = 0.8, cex.lab = 0.8, cex = 1.25, xlab = "NMDS 1", ylab = "NMDS 2", tck=-0.035, asp = 1, xlim = c(-0.56, 0.54), ylim = c(-0.35, 0.35))

dev.off()


pdf("El-Sayed_et_al_PCOA.pdf", width=3.5, height=3.5)

plot (PCOA$vectors.cor[,1], PCOA$vectors.cor[,2], pch = symbols, col = "transparent", bg = colors, cex.axis = 0.8, cex.lab = 0.8, cex = 1.25, xlab = "PCO 1", ylab = "PCO 2", tck=-0.035, asp = 1, xlim = c(-0.4, 0.9), ylim = c(-0.6, 0.6))

dev.off()