#################################################################################################################

#install packages

if("geiger" %in% rownames(installed.packages()) == FALSE){install.packages("geiger")
} else {print (paste0("'geiger' has already been installed in library"))}
if("maps" %in% rownames(installed.packages()) == FALSE){install.packages("maps")
} else {print (paste0("'maps' has already been installed in library"))}
if("phytools" %in% rownames(installed.packages()) == FALSE){install.packages("phytools")
} else {print (paste0("'phytools' has already been installed in library"))}
if("pacman" %in% rownames(installed.packages()) == FALSE){install.packages("pacman")
} else {print (paste0("'pacman' has already been installed in library"))}
if("geomorph" %in% rownames(installed.packages()) == FALSE){install.packages("geomorph")
} else {print (paste0("'geomorph' has already been installed in library"))}
if("ape" %in% rownames(installed.packages()) == FALSE){install.packages("ape")
} else {print (paste0("'ape' has already been installed in library"))}


#load packages
pacman::p_load("geomorph", "geiger", "ape", "phytools", "maps")



#################################################################################################################
############################################# INPUT FILES #######################################################

# 1- tps file
# 2- metadata file
# 3- tree file

data(pupfish)
str(pupfish)


#################################################################################################################
######################### Phylogenetic Generalized Least Squares (PLGS): Linear Models  #########################

plethtree <- read.tree('plethtree.tre')
plethland <- readland.tps('PlethodonLand.tps',specID = "ID",
                          warnmsg = FALSE)

gps <- read.csv('PlethGps.csv', header=TRUE, row.names=1)
Y.gpa <- gpagen(plethland, print.progress = FALSE)
M <- mshape(Y.gpa$coords)
size <- Y.gpa$Csize

shape <- Y.gpa$coords
shape.test <- treedata(phy = plethtree, data = two.d.array(shape), warnings = TRUE)


data.matched <- treedata(phy = plethtree, data = gps, warnings=FALSE)
elev <- as.factor(data.matched$data); names(elev) <- row.names(data.matched$data)

gdf <- geomorph.data.frame(shape=shape, size=size,elev = elev, plethtree=plethtree)

links <- matrix(c(4,3,2,1,1,6,7,8,9,10,1,1,11,5,5,4,2,3,7,8,9,10,11,9,10,1),
                ncol=2,byrow=FALSE)
plot(ladderize(plethtree),edge.width=3)
axisPhylo(1)


pgls.reg <- procD.pgls(f1 = shape ~ size, phy = plethtree, data = gdf, print.progress = FALSE)
summary(pgls.reg)

allom.plot <- plot(pgls.reg, type = "regression", predictor = gdf$size,
                   reg.type ="RegScore", pch=19, cex=1.5, xlab = "size") # make sure to have a predictor 
fit.line <- lm(allom.plot$RegScore ~ gdf$size)
abline(fit.line, col = "red")

### Phylogenetic Signal

PS.shape <- physignal(A=shape,phy=plethtree,iter=999, print.progress = FALSE)
summary(PS.shape)
plot(PS.shape)

### Phylogenetic Ordination

#### Phylomorphospace
PCA <- gm.prcomp(shape, phy=plethtree)
plot(PCA, phylo = TRUE, pch=21, bg=gdf$elev, cex=2, phylo.par = list(tip.labels = FALSE, node.labels = FALSE) )
legend("topleft", pch=21, pt.bg = unique(gdf$elev), legend = levels(gdf$elev))

#### Phylogenetic PCA (pPCA): With GLS-centered residuals
pPCA <- gm.prcomp(shape, phy=plethtree, GLS = TRUE, transform = TRUE)
plot(pPCA, phylo = FALSE, pch=21, bg=gdf$elev, cex=2, phylo.par = list(tip.labels = FALSE, node.labels = FALSE) )
legend("topleft", pch=21, pt.bg = unique(gdf$elev), legend = levels(gdf$elev))

#### Phylogenetic PCA (pPCA): With GLS-transformed residuals
pPCA2 <- gm.prcomp(shape, phy=plethtree, GLS = TRUE, transform = TRUE)
plot(pPCA2, phylo = TRUE, pch=21, bg=gdf$elev, cex=2, phylo.par = list(tip.labels = FALSE, node.labels = FALSE) )
legend("topleft", pch=21, pt.bg = unique(gdf$elev), legend = levels(gdf$elev))

#### Phylogenetically-Aligned Components Analysis (PACA)
PACA <- gm.prcomp(shape, phy=plethtree, align.to.phy = TRUE)
plot(PACA, phylo = TRUE, pch=21, bg=gdf$elev, cex=2, phylo.par = list(tip.labels = FALSE, node.labels = FALSE) )
legend("topleft", pch=21, pt.bg = unique(gdf$elev), legend = levels(gdf$elev))