#Maryam Nouri-Aiin
#PCoA Codes

library(polysat)
library(dplyr)
library(ggplot2)
library(ggfortify)


# Getting started: at tutorial
## Creating a dataset
### Accessor and Replacement Functions for "gendata" Objects
getwd()
agrestis <- read.GeneMapper("/Users/maryam/Desktop/MS Fall 2021/Polysat/A. agrestis for polysat.txt")
summary(agrestis)
Samples(agrestis)
Loci(agrestis)
myloci <- Loci(agrestis)
Description(agrestis) <- "Dataset for the tutorial"
PopNames(agrestis) <- c("PopA", "PopB", "PopC", "PopD", "PopE", "PopF")
PopInfo(agrestis) <- rep(1:6, times= c(26,14,45,25,22,64))
Ploidies(agrestis) <- 4
Usatnts(agrestis) <- c(4, 4, 4, 3, 4, 4, 4)
PopInfo(agrestis)
mysamples <- PopInfo(agrestis)
Samples(agrestis, populations = "PopF")
Usatnts(agrestis)
save(agrestis, file="agrestis.RData")

### Lynch.distance, meandistance.matrix2 16 point
testmat2 <- meandistance.matrix2(agrestis, distmetric = Lynch.distance)
pca5 <- cmdscale(testmat2)
summary(pca5)
pca <- princomp(pca5)

tiff("plot222", width = 5, height = 5, units = 'in', res = 300)

Plot2 <- plot(jitter(pca5[,1], 1500), jitter(pca5[,2], 1000), col=PopInfo(agrestis),
              main = "B                            PCoA with Lynch distance                          ", xlab= "PC1", ylab="PC2", pch=16, cex= 0.75, frame.plot = TRUE)

legend("topleft", col=1:6, legend = c("CRN","AU","CW","MTF","HG","HF"),pch = 20, bty='n', cex=0.75)
dev.off() 

### Bruvo.distance, meandistance.matrix2 19 point same as above
testmat1 <- meandistance.matrix2(agrestis)
pca1 <- cmdscale(testmat1)
pca3 <- princomp(pca1)
tiff("PCoA with Bruvo distance (A. agrestis)", width = 5, height = 5, units = 'in', res = 300)

Mainplot1 <- plot(jitter(pca1[,1],1500), jitter(pca1[,2], 1000), col=PopInfo(agrestis),
                  main = "A                           PCoA with Bruvo distance                          ", xlab= "PC1", ylab="PC2", pch=16, cex= 0.75, frame.plot = TRUE)

legend("topleft", col=1:6, legend = c("CRN","AU","CW","MTF","HG","HF"),pch = 20, bty='n', cex=0.75)

dev.off()