###########################################################################################################################################
### Source code for: Madzia et al. Historical significance and taxonomic status of Ischyrodon meriani... PeerJ                      	###
### Adapted by D. Madzia from a code written by V. Fischer (Université de Liege, Belgium).                                              ###
###########################################################################################################################################

library(ape)

#===================== CLUSTER ANALYSES & PCOA =====================
library(cluster) 

#Get data, apply completeness threshold and scale pruned dataset
data<-read.table("C:/.../data_meriani.txt", header=TRUE,as.is=TRUE) # Modify the path.
threshold <- 0.5
completeness.t <- sort(apply(!is.na(data),1,sum),decreasing=TRUE)/length(data[1,])
taxa.comp <- names(completeness.t)[completeness.t>=threshold] 
data.s<-scale(data[taxa.comp,]) # scale data using std, taking ONLY taxa passing the completeness threshold.

#Compute distance and cluster dendrogram
gower.d <- daisy(data.s, metric="gower") #using Gower (1971 method), which is better for mixed data types (continuous and discrete).
fit <- hclust(gower.d, method="ward.D2")
plot(fit)
warnings()
read(table)
#Export .pdf 
pdf("C:/.../cluster_dendrogram.pdf", width=10, height=7)
plot(fit)
dev.off()

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

pdf("C:/.../PCoA.pdf", width=14, height=14)
biplot(pcoa, Y=NULL, plot.axes = c(2,1), dir.axis1=1, dir.axis2=1)
dev.off()