## Evolution of Syndromes
rm(list = ls()) 
Sys.setenv(LANG = "en")

#### Mode of Evolution and test for adaptation to different phenotypic optima

# Mode of Evo - Pagel's tree modification

library(ape)
library(geiger)
library(phytools)

setwd("~/Documents/Uni/PHD/Statistics/Evo/Review")
tree <- read.nexus("Reginato_9-gene_Dated_fixed_topology.tre") # load tree

tree=tree$ALL_FOSSILS 

dat <- read.csv("MelasTraitsMean_Rev.csv", header=T, sep=";", row.names = 1) # load data

# remove species with NAs
dat=dat[c(1:60,62:393,395:413),]

# Name syndromes
syndrome <- dat$Pred.Syn.Disp 
names(syndrome) <- dat$Phylo.Name 

#bring them into correct order -> remove tips not in data and data not in tips
tree1<-treedata(tree,syndrome,sort=T,warnings=T)$phy
data<-treedata(tree,syndrome,sort=T,warnings=T)$data


### check for evo mode

# no tree trans
no_transARD <- fitDiscrete(tree1, data, model = "ARD", transform="none") # no transoftmation - rate is constant
no_transARD # AIC = 283.185858, 

no_transER <- fitDiscrete(tree1, data, model = "ER", transform="none")
no_transER # AIC = 317.998184

no_transSYM <- fitDiscrete(tree1, data, model = "SYM", transform="none")
no_transSYM # AIC = 304.743996

aicsNOTrans <- c(no_transARD=no_transARD$opt$aicc, no_transER=no_transER$opt$aicc, no_transSYM=no_transSYM$opt$aicc)
aic.w(aicsNOTrans) 
# ARD lowest AICc - 283.185858 no_transARD  no_transER no_transSYM 
#                               0.99997024  0.00000004  0.00002972 

# early burst transformation
EB_ARD <- fitDiscrete(tree1, data, model = "ARD", transform="EB") # model where the rate of evolution increases or decreases exponentially through time
EB_ARD # frequency of best fit = NA

EB_ER <- fitDiscrete(tree1, data, model = "ER", transform="EB")
EB_ER # frequency of best fit = NA

EB_SYM <- fitDiscrete(tree1, data, model = "SYM", transform="EB")
EB_SYM # frequency of best fit = NA

aicsEB <- c(EB_ARD=EB_ARD$opt$aicc, EB_ER=EB_ER$opt$aicc, EB_SYM=EB_SYM$opt$aicc)
aic.w(aicsEB) 

# kappa
kappa_ARD <- fitDiscrete(tree1, data, transform="kappa",  model="ARD", bounds=list(kappa=c(0,1))) # model of trait evolution where character divergence is related to the number of speciation events between two species
kappa_ARD # kappa = 0, AIC = 285.548453

kappa_ER <- fitDiscrete(tree1, data, transform="kappa",  model="ER", bounds=list(kappa=c(0,1)))
kappa_ER # kappa = 0.037505, AIC = 291.061805

kappa_SYM <- fitDiscrete(tree1, data, transform="kappa",  model="SYM", bounds=list(kappa=c(0,1)))
kappa_SYM # kappa = 0.092256, AIC = 283.817545

aicsKapp <- c(kappa_ARD=kappa_ARD$opt$aicc, kappa_ER=kappa_ER$opt$aicc, kappa_SYM=kappa_SYM$opt$aicc)
aic.w(aicsKapp) # SYM  = 0.73265131 => trait change occurs at events of speciation
                #  kappa_ARD   kappa_ER  kappa_SYM 
#                   0.21575122 0.02368044 0.76056834 
# delta
delta_ARD <- fitDiscrete(tree1, data, transform="delta", model="ARD") # relative contributions of early versus late evolution in the tree to the covariance of species trait value
delta_ARD # Where delta is greater than 1, recent evolution has been relatively fast; if delta is less than 1, recent evolution has been comparatively slow
# delta = 2.999999, AIC = 279.336399

delta_ER <- fitDiscrete(tree1, data, transform="delta", model="ER")
delta_ER # delta = 2.999999, AIC = 308.860349

delta_SYM <- fitDiscrete(tree1, data, transform="delta", model="SYM")
delta_SYM # delta = 2.999999, AIC = 297.496185

aicsDelta <- c(delta_SYM=delta_SYM$opt$aicc, delta_ER=delta_ER$opt$aicc, delta_ARD=delta_ARD$opt$aicc)# transformed model have higher AIC than the original - so no transformation to explain the diet of of species better than the gradual BM model
aic.w(aicsDelta) #c ARD delta = 0.99983032=> recent evolution has been fast
                  # delta_SYM   delta_ER  delta_ARD 
                  # 0.00016901 0.00000067 0.99983032 
# lambda
lambda_ARD <- fitDiscrete(tree1, data,transform="lambda",  model="ARD", bounds=list(lambda=c(0,1))) # models that fits the extent to which the phylogeny predicts covariance among trait values for species.
lambda_ARD # # lambda = 1.000000, AIC = 301.780507, frequency of best fit = NA

lambda_ER <- fitDiscrete(tree1, data,transform="lambda",  model="ER", bounds=list(lambda=c(0,1)))
lambda_ER # lambda = 0.944656, AIC = 315.850471, frequency of best fit = NA

lambda_SYM <- fitDiscrete(tree1, data,transform="lambda",  model="SYM", bounds=list(lambda=c(0,1)))
lambda_SYM # lambda = 0.974223, AIC = 305.703791, frequency of best fit = NA

aicsLambda <- c(lambda_ER=lambda_ER$opt$aicc, lambda_ARD=lambda_ARD$opt$aicc, lambda_SYM=lambda_SYM$opt$aicc)# transformed model have higher AIC than the original - so no transformation to explain the diet of of species better than the gradual BM model
aic.w(aicsLambda) # ARD lambda = 0.99986673 => strong phylo signal 
                  # lambda_ER lambda_ARD lambda_SYM 
                  # 0.00000097 0.99986673 0.00013231 
# total
aicsTotal <- c(EB_SYM=EB_SYM$opt$aicc, kappa_SYM=kappa_SYM$opt$aicc, delta_ARD=delta_ARD$opt$aicc, lambda_ARD=lambda_ARD$opt$aicc, no_transARD=no_transARD$opt$aicc)# transformed model have higher AIC than the original - so no transformation to explain the diet of of species better than the gradual BM model
aic.w(aicsTotal) # delta_ARD  => recent evolution has been fast
                # EB_SYM   kappa_SYM   delta_ARD  lambda_ARD no_transARD 
                #0.00000000  0.11803317  0.74773818  0.01576946  0.11845920 

#################
# OU Estimation #
#################
# Adaptation to different phenotypic optima

library(devtools)
install_github("glmgen/genlasso")
install_github("khabbazian/l1ou")
library(genlasso)
library(l1ou)
library(phytools)
library(ape)
library(geiger)
library(parallel)

Data <- read.csv("MelasTraitsMean_Rev.csv", sep=";",  dec = ",")
tree <- read.nexus("Reginato_9-gene_Dated_fixed_topology.tre")

tree <- tree$ALL_FOSSILS

## remove species only NAs!
tabSpace=Data[c(1:60,62:393,395:413),]

#subset the traits from M1.3
SpaceCHAR=tabSpace[,c(5:43,45:47,52,53,55)]

numeric_columns <- which(sapply(SpaceCHAR, function(x) class(x) == "numeric")) # get all numeric variables

Spaceochar=numeric_columns 

source("~/Documents/Uni/PHD/Statistics/PCoA/F_DISPARITY_functions vs4.R") 
library(vegan)

dis=MD(as.matrix(SpaceCHAR),Spaceochar,removinvar=F) # produces large list that contains dissimilarity matrix we need dis$FOOTE
mapco=PCO(dis$FOOTE, correction=" ")

PC1_3 <- mapco$pco[,1:3] # returnes the PC values 1,2,3 for all species
row.names(PC1_3)<-tabSpace$Phylo.Name # Use Phylo names to lable PC
tree<-treedata(tree,PC1_3)$phy
PC1_3<-treedata(tree,PC1_3,sort=T,warnings=T)$data

#02 - PREPARE THE DATA
data<-adjust_data(tree,PC1_3)

#03 - estimate and plot shifts in phenotypic optima - OU Model
fM0<-estimate_shift_configuration(data$tree,data$Y,criterion = "pBIC",nCores=1) # all 3 axis
fM1<-estimate_shift_configuration(data$tree,data$Y[,1],criterion = "pBIC",nCores=1) #(adaptation 'strength' alpha and drift variance sigma^2) # 1st axis
fM2<-estimate_shift_configuration(data$tree,data$Y[,2],criterion = "pBIC",nCores=1) #(adaptation 'strength' alpha and drift variance sigma^2)
fM3<-estimate_shift_configuration(data$tree,data$Y[,3],criterion = "pBIC",nCores=1) #(adaptation 'strength' alpha and drift variance sigma^2)

save(fM0, file="Output/fM0InclMedi.Rdata")
save(fM1, file="Output/fM1InclMedi.Rdata")
save(fM2, file="Output/fM2InclMedi.Rdata")
save(fM3, file="Output/fM3InclMedi.Rdata")

######
quartz()
pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Evo_PC1-3InclMedi.pdf")
plot(fM0,show.tip.label=FALSE)
dev.size(units = c("cm"))
dev.off()

pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Evo_PC1InclMedi.pdf")
plot(fM1,show.tip.label=FALSE)
dev.size(units = c("cm"))
dev.off()

pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Evo_PC2InclMedi.pdf")
plot(fM2,show.tip.label=FALSE)
dev.size(units = c("cm"))
dev.off()

pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Evo_PC3InclMedi.pdf")
plot(fM3,show.tip.label=FALSE)
dev.size(units = c("cm"))
dev.off()

summary(fM0)


#04 - estimate convergent shifts 
fM0conv<-estimate_convergent_regimes(fM0,criterion="pBIC",nCores=1)
fM0conv

quartz()
pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Evo_PC1-3_Conv_Shifts_PhenoOpt..pdf")
plot(fM0conv, show.tip.label=FALSE)
dev.size(units = c("cm"))
dev.off()

save(fM0conv, file="Output/fM0conv.Rdata")

sink("Output/Sum_Shifts_PhenoOpt_Pc_1-3Fin.txt")
summary(fM0)
summary(fM0conv)
sink(NULL)

pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/TEST_tree.pdf", height=20)
plot(tree, show.tip.label=FALSE)
edgelabels()
dev.size(units = c("cm"))
dev.off()

#
fM1conv<-estimate_convergent_regimes(fM1,criterion="pBIC",nCores=1)
fM1conv

quartz()
pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Evo_PC1_Conv_Shifts_PhenoOpt..pdf")
plot(fM1conv, show.tip.label=FALSE)
dev.size(units = c("cm"))
dev.off()

save(fM1conv, file="Output/fM1conv.Rdata")

sink("Output/Sum_Shifts_PhenoOpt_Pc_1Fin.txt")
summary(fM1)
summary(fM1conv)
sink(NULL)

#
fM2conv<-estimate_convergent_regimes(fM2,criterion="pBIC",nCores=1) 
fM2conv

quartz()
pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Evo_PC2_Conv_Shifts_PhenoOpt..pdf")
plot(fM2conv, show.tip.label=FALSE)
dev.size(units = c("cm"))
dev.off()

save(fM2conv, file="Output/fM2conv.Rdata")

sink("Output/Sum_Shifts_PhenoOpt_Pc_2Fin.txt")
summary(fM2)
summary(fM2conv)
sink(NULL)

#
fM3conv<-estimate_convergent_regimes(fM3,criterion="pBIC",nCores=1)
fM3conv

quartz()
pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Evo_PC3_Conv_Shifts_PhenoOpt..pdf")
plot(fM3conv, show.tip.label=FALSE)
dev.size(units = c("cm"))
dev.off()

save(fM3conv, file="Output/fM3conv.Rdata")

sink("Output/Sum_Shifts_PhenoOpt_Pc_3Fin.txt")
summary(fM3)
summary(fM3conv)
sink(NULL)

### the lower the BIC (pBIC) the better the model fits

# fit_OU - Fits an OU model based on a given configuration - configuration obtained from estimate_shift_configuration
shift.config<-fM0[["shift.configuration"]] # use best Model 290,339,400,105,654,32,298,669
fM_0000<-fit_OU(data$tree,data$Y,shift.config, criterion="pBIC") # = fm0 - only if one wants to compare this model with models based on knowledge of taxons

# plot best model
quartz()
pdf("~/Documents/Uni/PHD/Statistics/Evo/Evo_PC1-3InclMedi.pdf")
plot(fM0,show.tip.label=FALSE)
dev.size(units = c("cm"))
dev.off()

### plot ACE with OU model PC1-3 all traits
setwd("~/Documents/Uni/PHD/Statistics/Evo/Review/Output")

load("fM0InclMedi.RData")
head(fM0)
fM0$tree

edgelength <- tree1$edge.length
fM0$tree$edge.length <- edgelength

edge <- tree1$edge
fM0$tree$edge <- edge

shift.config<-fM0$shift.configuration 
fM_0000<-fit_OU(fM0$tree,fM0$Y,shift.config, criterion="pBIC")

fM_0000$tree$tip.label <- rep(c(""), times = 336)
##########


shiftcol <- c("khaki4", "darkorange", "forestgreen", "orangered", "navy", "olivedrab", "slateblue", "lightcyan3", "mediumseagreen", "darkslateblue","firebrick", "gray")
shiftcol1 <- c("gray0", "gray24", "gray29", "gray50", "gray39", "gray49", "gray54", "gray79", "gray64", "gray69", "gray34", "gray84", "gray89")

quartz()
pdf("~/Documents/Uni/PHD/Statistics/MCMC/MuSSE_CharacterMapping_SYM_and_OUPC1-3allTraits_with_lables.pdf", height = 15, width = 12)
plot(fM0, palette = shiftcol, type="phylogram",show.tip.label = T, asterisk = F, edge.shift.ann = F, label.offset = 3)
nodelabels(pie=st,piecol=piecol, cex=0.5) # st from MCMC
tiplabels(pch = 22, bg = tip_col, cex = 1, adj = 1)
axisPhylo()
dev.size(units = c("cm"))
dev.off()

#### Only 7 shifted tribes together (Astronieae,Blakeeae / Pyxidantheae,Melastomateae,Merianieae,Miconieae,Olisbeoideae,Sonerileae)

# select tribest that shifted
tabSpaceShift <- subset(tabSpace,Tribe...Subfamily=="Astronieae"|Tribe...Subfamily=="Blakeeae / Pyxidantheae"|Tribe...Subfamily=="Melastomateae"|Tribe...Subfamily=="Merianieae"|Tribe...Subfamily=="Miconieae"|Tribe...Subfamily=="Olisbeoideae"|Tribe...Subfamily=="Sonerileae")

#subset the traits from M1.3
SpaceCHAR=tabSpaceShift[,c(5:43,45:47,52,53,55)]

numeric_columns <- which(sapply(SpaceCHAR, function(x) class(x) == "numeric")) # get all numeric variables

Spaceochar=numeric_columns 

source("~/Documents/Uni/PHD/Statistics/PCoA/F_DISPARITY_functions vs4.R") #script has to be in same folder - runs in background
library(vegan)

dis=MD(as.matrix(SpaceCHAR),Spaceochar,removinvar=F) # produces large list that contains dissimilarity matrix we need dis$FOOTE
mapco=PCO(dis$FOOTE, correction=" ")

PC1_3 <- mapco$pco[,1:3] # returnes the PC values 1,2,3 for all species
row.names(PC1_3)<-tabSpaceShift$Phylo.Name # Use Phylo names to lable PC
tree<-treedata(tree,PC1_3)$phy
PC1_3<-treedata(tree,PC1_3,sort=T,warnings=T)$data

####
##### pcoa
polllist3<-as.vector(tabSpaceShift$Pred.Syn.PcoA)

colcons3<-polllist3
colcons3[which(colcons3=="Buzz")]<-"blue"
colcons3[which(colcons3=="Buzz Pred")]<-"cornflowerblue"
colcons3[which(colcons3=="Passerine")]<-"darkgoldenrod"
colcons3[which(colcons3=="Nectar foraging vertebrates")]<-"deeppink"
colcons3[which(colcons3=="Nectar foraging vertebrates Pred")]<-"lightpink"
colcons3[which(colcons3=="Generalist")]<-"chartreuse"
colcons3[which(colcons3=="Generalist Pred")]<-"darkseagreen4"
              
names(colcons3)<-tabSpaceShift$Taxon
            
Eig=mapco$PCAEval[mapco$PCAEval>0]
var1.2=Eig[1]*100/sum(Eig)
var2.2=Eig[2]*100/sum(Eig)
var3.2=Eig[3]*100/sum(Eig)
            
c(var1.2, var2.2, var3.2) #var1 to 3 give the percentage of variance explained by the axes
            
sum(var1.2,var2.2,var3.2)
            
pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/PCoA_Only_Shifted_Tribes_M1_Rev.pdf", width=10, height=10)
quartz()
#par(mfrow=c(1,2)) #to split window
plot(mapco$pco, col = colcons3, pch=16, cex=4, main="Only Shifted Tribes M1 Expl. Var. 35.49%", xlab = "PC1_20.16%", ylab = "PC2_8.76%") #colour for presumed pollinator
legend("bottomleft", legend = c("Buzz-Bee","Buzz Pred","Nectar","Nectar Pred","Passerine", "Generalis", "Generalist Pred"),
col =  c("blue","cornflowerblue", "deeppink","lightpink", "darkgoldenrod", "chartreuse", "darkseagreen4"),
pch = 16, cex = 1 )
            
dev.size(units = c("cm"))
dev.off()
            
#02 - PREPARE THE DATA
data<-adjust_data(tree,PC1_3)
            
#03 - estimate and plot shifts in phenotypic optima - OU Model
fM0Shif<-estimate_shift_configuration(data$tree,data$Y,criterion = "pBIC",nCores=1) # all 3 axis
fM1Shif<-estimate_shift_configuration(data$tree,data$Y[,1],criterion = "pBIC",nCores=1) #(adaptation 'strength' alpha and drift variance sigma^2) # 1st axis
fM2Shif<-estimate_shift_configuration(data$tree,data$Y[,2],criterion = "pBIC",nCores=1) #(adaptation 'strength' alpha and drift variance sigma^2)
fM3Shif<-estimate_shift_configuration(data$tree,data$Y[,3],criterion = "pBIC",nCores=1) #(adaptation 'strength' alpha and drift variance sigma^2)
            
save(fM0Shif,file="~/Documents/Uni/PHD/Statistics/Evo/Review/Output/fM0_onlyShift.Rdata")
save(fM1Shif, file="~/Documents/Uni/PHD/Statistics/Evo/Review/Output/fM1_onlyShift.Rdata")
save(fM2Shif, file="~/Documents/Uni/PHD/Statistics/Evo/Review/Output/fM2_onlyShift.Rdata")
save(fM3Shif, file="~/Documents/Uni/PHD/Statistics/Evo/Review/Output/fM3_onlyShift.Rdata")
            
            
quartz()
pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Evo_PC1-3_onlyShift.pdf")
plot(fM0Shif,show.tip.label=FALSE)
dev.size(units = c("cm"))
dev.off()
            
pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Evo_PC1_onlyShift.pdf")
plot(fM1Shif,show.tip.label=FALSE)
dev.size(units = c("cm"))
dev.off()
            
pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Evo_PC2_onlyShift.pdf")
plot(fM2Shif,show.tip.label=FALSE)
dev.size(units = c("cm"))
dev.off()
            
pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Evo_PC3_onlyShift.pdf")
plot(fM3Shif,show.tip.label=FALSE)
dev.size(units = c("cm"))
dev.off()
            
summary(fM0Shif)
            
#04 - estimate convergent shifts - same output as estimate_shift_configuration - OU Model
fM0convShif<-estimate_convergent_regimes(fM0Shif,criterion="pBIC",nCores=1)
fM0convShif
            
quartz()
pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Evo_PC1-3_Conv_Shifts_PhenoOpt_onlyShift.pdf")
plot(fM0convShif, show.tip.label=FALSE)
dev.size(units = c("cm"))
dev.off()
            
save(fM0convShif, file="~/Documents/Uni/PHD/Statistics/Evo/Review/Output/fM0conv_onlyShift.Rdata")
            
sink("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Sum_Shifts_PhenoOpt_Pc_1-3_onlyShift.txt")
summary(fM0Shif)
summary(fM0convShif)
sink(NULL)

####
fM1convShif<-estimate_convergent_regimes(fM1Shif,criterion="pBIC",nCores=1)
fM1convShif

quartz()
pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Evo_PC1_Conv_Shifts_PhenoOpt_onlyShift.pdf")
plot(fM1convShif, show.tip.label=FALSE)
dev.size(units = c("cm"))
dev.off()

save(fM1convShif, file="~/Documents/Uni/PHD/Statistics/Evo/Review/Output/fM1conv_onlyShift.Rdata")

sink("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Sum_Shifts_PhenoOpt_Pc_1_onlyShift.txt")
summary(fM1Shif)
summary(fM1convShif)
sink(NULL)

####
fM2convShif<-estimate_convergent_regimes(fM2Shif,criterion="pBIC",nCores=1) # no shifts
fM2convShif

quartz()
pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Evo_PC2_Conv_Shifts_PhenoOpt_onlyShift.pdf")
plot(fM2convShif, show.tip.label=FALSE)
dev.size(units = c("cm"))
dev.off()

save(fM2convShif, file="~/Documents/Uni/PHD/Statistics/Evo/Review/Output/fM2conv_onlyShift.Rdata")

sink("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Sum_Shifts_PhenoOpt_Pc_2_onlyShift.txt")
summary(fM2Shif)
summary(fM2convShif)
sink(NULL)

####
fM3convShif<-estimate_convergent_regimes(fM3Shif,criterion="pBIC",nCores=1)
fM3convShif

quartz()
pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Evo_PC3_Conv_Shifts_PhenoOpt_onlyShift.pdf")
plot(fM3convShif, show.tip.label=FALSE)
dev.size(units = c("cm"))
dev.off()

save(fM3convShif, file="~/Documents/Uni/PHD/Statistics/Evo/Review/Output/fMconv_onlyShift.Rdata")

sink("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Sum_Shifts_PhenoOpt_Pc_3_onlyShift.txt")
summary(fM3Shif)
summary(fM3convShif)
sink(NULL)

#################### now OU Models for all tribes that shifted pollinators separately
### Astronieae,Blakeeae / Pyxidantheae,Melastomateae,Merianieae,Miconieae,Olisbeoideae,Sonerileae

tabSpaceShift <- subset(tabSpace,Tribe...Subfamily=="Sonerileae") #|Tribe...Subfamily=="Blakeeae / Pyxidantheae"|Tribe...Subfamily=="Melastomateae"|Tribe...Subfamily=="Merianieae"|Tribe...Subfamily=="Miconieae"|Tribe...Subfamily=="Olisbeoideae"|Tribe...Subfamily=="Sonerileae")

#subset the traits from M1.3
SpaceCHAR=tabSpaceShift[,c(5:43,45:47,52,53,55)]

numeric_columns <- which(sapply(SpaceCHAR, function(x) class(x) == "numeric")) # get all numeric variables

Spaceochar=numeric_columns 

source("~/Documents/Uni/PHD/Statistics/PCoA/F_DISPARITY_functions vs4.R") #script has to be in same folder - runs in background
library(vegan)

dis=MD(as.matrix(SpaceCHAR),Spaceochar,removinvar=F) # produces large list that contains dissimilarity matrix we need dis$FOOTE
mapco=PCO(dis$FOOTE, correction=" ")

PC1_3 <- mapco$pco[,1:3] # returnes the PC values 1,2,3 for all species
row.names(PC1_3)<-tabSpaceShift$Phylo.Name # Use Phylo names to lable PC
tree<-treedata(tree,PC1_3)$phy
PC1_3<-treedata(tree,PC1_3,sort=T,warnings=T)$data

####
##### pcoa
polllist3<-as.vector(tabSpaceShift$Pred.Syn.PcoA)

colcons3<-polllist3
colcons3[which(colcons3=="Buzz")]<-"blue"
  colcons3[which(colcons3=="Buzz Pred")]<-"cornflowerblue"
    colcons3[which(colcons3=="Passerine")]<-"darkgoldenrod"
      colcons3[which(colcons3=="Nectar foraging vertebrates")]<-"deeppink"
        colcons3[which(colcons3=="Nectar foraging vertebrates Pred")]<-"lightpink"
          colcons3[which(colcons3=="Generalist")]<-"chartreuse"
            colcons3[which(colcons3=="Generalist Pred")]<-"darkseagreen4"
              
            names(colcons3)<-tabSpaceShift$Taxon
            
            Eig=mapco$PCAEval[mapco$PCAEval>0]
            var1.2=Eig[1]*100/sum(Eig)
            var2.2=Eig[2]*100/sum(Eig)
            var3.2=Eig[3]*100/sum(Eig)
            
            c(var1.2, var2.2, var3.2) #var1 to 3 give the percentage of variance explained by the axes
            
            sum(var1.2,var2.2,var3.2)
            
            pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/PCoA_Only_Shifted_Tribes_M1_Son_Rev_Fin.pdf", width=10, height=10)
            quartz()
            #par(mfrow=c(1,2)) #to split window
            plot(mapco$pco, col = colcons3, pch=16, cex=4, main="Only Shifted Tribes M1 Sonerileae Expl. Var. 36.75%", xlab = "PC1_16.02%", ylab = "PC2_10.64%") #colour for presumed pollinator
            #legend("bottomleft", legend = c("Buzz-Bee","Buzz Pred","Nectar","Nectar Pred","Passerine", "Generalis", "Generalist Pred"),
                  # col =  c("blue","cornflowerblue", "deeppink","lightpink", "darkgoldenrod", "chartreuse", "darkseagreen4"),
                  # pch = 16, cex = 1 )
            
            dev.size(units = c("cm"))
            dev.off()
            
#############02 - PREPARE THE DATA
            data<-adjust_data(tree,PC1_3)
            
            #03 - estimate and plot shifts in phenotypic optima - OU Model
            fM0Shif<-estimate_shift_configuration(data$tree,data$Y,criterion = "pBIC",nCores=1) # all 3 axis
            fM1Shif<-estimate_shift_configuration(data$tree,data$Y[,1],criterion = "pBIC",nCores=1) #(adaptation 'strength' alpha and drift variance sigma^2) # 1st axis
            fM2Shif<-estimate_shift_configuration(data$tree,data$Y[,2],criterion = "pBIC",nCores=1) #(adaptation 'strength' alpha and drift variance sigma^2)
            fM3Shif<-estimate_shift_configuration(data$tree,data$Y[,3],criterion = "pBIC",nCores=1) #(adaptation 'strength' alpha and drift variance sigma^2)
            
            save(fM0Shif,file="~/Documents/Uni/PHD/Statistics/Evo/Review/Output/fM0_onlyShift_Son.Rdata")
            save(fM1Shif, file="~/Documents/Uni/PHD/Statistics/Evo/Review/Output/fM1_onlyShift_Son.Rdata")
            save(fM2Shif, file="~/Documents/Uni/PHD/Statistics/Evo/Review/Output/fM2_onlyShift_Son.Rdata")
            save(fM3Shif, file="~/Documents/Uni/PHD/Statistics/Evo/Review/Output/fM3_onlyShift_Son.Rdata")
          
            quartz()
            pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Evo_PC1-3_onlyShift_Son.pdf")
            plot(fM0Shif,show.tip.label=FALSE)
            dev.size(units = c("cm"))
            dev.off()
            
            pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Evo_PC1_onlyShift_Son.pdf")
            plot(fM1Shif,show.tip.label=FALSE)
            dev.size(units = c("cm"))
            dev.off()
            
            pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Evo_PC2_onlyShift_Son.pdf")
            plot(fM2Shif,show.tip.label=FALSE)
            dev.size(units = c("cm"))
            dev.off()
            
            pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Evo_PC3_onlyShift_Son.pdf")
            plot(fM3Shif,show.tip.label=FALSE)
            dev.size(units = c("cm"))
            dev.off()
            
            summary(fM0Shif)
            
            #04 - estimate convergent shifts - same output as estimate_shift_configuration - OU Model
            fM0convShif<-estimate_convergent_regimes(fM0Shif,criterion="pBIC",nCores=1)
            fM0convShif
            
            quartz()
            pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Evo_PC1-3_Conv_Shifts_PhenoOpt_onlyShift_Son.pdf")
            plot(fM0convShif, show.tip.label=FALSE)
            dev.size(units = c("cm"))
            dev.off()
            
            save(fM0convShif, file="~/Documents/Uni/PHD/Statistics/Evo/Review/Output/fM0conv_onlyShift_Son.Rdata")
            
            sink("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Sum_Shifts_PhenoOpt_Pc_1-3_onlyShift_Son.txt")
            summary(fM0Shif)
            summary(fM0convShif)
            sink(NULL)
            
            ####
            fM1convShif<-estimate_convergent_regimes(fM1Shif,criterion="pBIC",nCores=1)
            fM1convShif
            
            quartz()
            pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Evo_PC1_Conv_Shifts_PhenoOpt_onlyShift_Son.pdf")
            plot(fM1convShif, show.tip.label=FALSE)
            dev.size(units = c("cm"))
            dev.off()
            
            save(fM1convShif, file="~/Documents/Uni/PHD/Statistics/Evo/Review/Output/fM1conv_onlyShift_Son.Rdata")
            
            sink("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Sum_Shifts_PhenoOpt_Pc_1_onlyShift_Son.txt")
            summary(fM1Shif)
            summary(fM1convShif)
            sink(NULL)
            
            ####
            fM2convShif<-estimate_convergent_regimes(fM2Shif,criterion="pBIC",nCores=1) # no shifts
            fM2convShif
            
            quartz()
            pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Evo_PC2_Conv_Shifts_PhenoOpt_onlyShift_Son.pdf")
            plot(fM2convShif, show.tip.label=FALSE)
            dev.size(units = c("cm"))
            dev.off()
            
            save(fM2convShif, file="~/Documents/Uni/PHD/Statistics/Evo/Review/Output/fM2conv_onlyShift_Son.Rdata")
            
            sink("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Sum_Shifts_PhenoOpt_Pc_2_onlyShift_Son.txt")
            summary(fM2Shif)
            summary(fM2convShif)
            sink(NULL)
            
            ####
            fM3convShif<-estimate_convergent_regimes(fM3Shif,criterion="pBIC",nCores=1)
            fM3convShif
            
            quartz()
            pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Evo_PC3_Conv_Shifts_PhenoOpt_onlyShift_Son.pdf")
            plot(fM3convShif, show.tip.label=FALSE)
            dev.size(units = c("cm"))
            dev.off()
            
            save(fM3convShif, file="~/Documents/Uni/PHD/Statistics/Evo/Review/Output/fMconv_onlyShift_Son.Rdata")
            
            sink("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Sum_Shifts_PhenoOpt_Pc_3_onlyShift_Son.txt")
            summary(fM3Shif)
            summary(fM3convShif)
            sink(NULL)
           
####################
# Plot tree with lables for Mel PC1-3 and PC1 and Pyx PC1-3 + Conv
synd <- tabSpaceShift$Pred.Syn.Disp 
names(synd) <- tabSpaceShift$Phylo.Name 

datcol1<-treedata(tree,synd,sort=T,warnings=T)$data
datcol <- as.data.frame(datcol1)

# Color syndromes
synd1 <- as.data.frame(datcol$V1)
tip_col <- synd1$`datcol$V1`

tip_col[which(tip_col=="Buzz")]<-"cornflowerblue" # cornflowerblue
tip_col[which(tip_col=="Passerine")]<-"darkgoldenrod" # darkgoldenrod
tip_col[which(tip_col=="Nectar foraging vertebrates")]<-"deeppink" # deeppink
tip_col[which(tip_col=="Generalist")]<-"darkorchid3" # chartreuse
        
      
pdf("~/Documents/Uni/PHD/Statistics/Evo/Review/Output/Shifts_WherePC1_conv_onlyShift_Mel.pdf", width = 20, height=80)
quartz()
plot(fM1convShif, type="phylogram")
tiplabels(pch = 22, bg = tip_col, cex = 1, adj = 1)
dev.size(units = c("cm"))
dev.off()



