# Analysing Flower shape variation on COMBINED dataset

# Script written by Sara van de Keker (sara.vandekerke@wur.nl)
# Written December 2018, updated December 2018
# Accompanying van de Kerke et al. 

library(geomorph) # using version ...
library(geoR)

### set working directory 
setwd("D:/PhD/OneDrive - WageningenUR/morphometrics_chapter_2-3/Rscripten")



### Load Pelargonium.combined.initial dataset

load("pelargonium_combined_initial.Rda")




#####################################################################
############### Generalised Procrustes Analysis (GPA) ############### 

pelargonium.combined.GPA <- gpagen(pelargonium.combined.initial$Coords, 
                                   surfaces=NULL, PrinAxes=TRUE, max.iter=NULL, 
                                   ProcD=TRUE, Proj=TRUE, print.progress = TRUE)


save(pelargonium.combined.GPA, file = "Pelargonium_combined_GPA.Rda")
load("Pelargonium_combined_GPA.Rda")


classifier.variables.PETAL.matching <- read.csv("Variables_PETAL_overlap.csv")
classifier.individuals <- read.csv("individuals68species.csv", header=FALSE)
classifier.individuals$numbers <- paste(classifier.individuals$V1, 
                                        classifier.individuals$V2, sep = "_")
classifier.individuals$species <- classifier.individuals$V6

classifier.individuals$combination <- 
  paste(classifier.individuals$V4,classifier.individuals$V5, sep = "_")


###############################################################################
############### Check whether data is normally distributed ####################


gdf3 <- geomorph.data.frame(pelargonium.combined.GPA) # make geomorph data frame

pelargonium.combined.anova <- procD.lm(coords ~ Csize, data = gdf3, iter = 999,
                                       RRPP = TRUE, print.progress = TRUE)
plot(pelargonium.combined.anova) # diagnostic plots

save(pelargonium.combined.anova, file = "Pelargonium_combined_ANOVA.Rda")


######################################################################
############### Check for outliers and their direction ###############


plotOutliers(pelargonium.combined.GPA$coords)
outliers <- plotOutliers(pelargonium.combined.GPA$coords)
plotRefToTarget(mshape(pelargonium.combined.GPA$coords),pelargonium.combined.GPA$coords[,,outliers[61]],
                method="vector", label = T)





########################################################################################
############### Plot all landmarks making use of the new GPA coordinates ############### 
########################################################################################


plotAllSpecimens(pelargonium.combined.GPA$coords, mean=TRUE, label=TRUE, 
                 plot.param=list(pt.bg="gray", txt.col=1, txt.cex=1.4, 
                                 pt.cex=1, mean.cex=2, mean.bg=1))






###############################################################################################
############### Plot individual specimens making use of the new GPA coordinates ############### 
###############################################################################################


plot3d(pelargonium.combined.GPA$coords[,,2][,c(1,2,3)], size = 8) 

plot(pelargonium.SPUR.raw[,,72])
plot(pelargonium.front.raw[,,1])




######################################################
############### Combine data in a list ############### 
######################################################


# All data has to be combined in a list to facilitate Geomorph. We combine the coordinate data (pelargonium.raw) with the 
# variable data (classifier.variables.front.extended)

# Combined dataset
pelargonium.combined <- list(pelargonium.combined.GPA$coords, 
                             pelargonium.combined.GPA$Csize, 
                             classifier.variables.PETAL.matching$Clade, 
                             classifier.variables.PETAL.matching$Pollinator, 
                             classifier.variables.PETAL.matching$PetalNumber,
                             classifier.variables.PETAL.matching$number,
                             classifier.variables.PETAL.matching$ID, 
                             classifier.individuals$species,
                             classifier.individuals$combination,
                             classifier.individuals$numbers)
names(pelargonium.combined) <- c("Coords", "Csize", "Clade", "Pollinator", 
                                 "PetalNumber", "number", "ID", 
                                 "Species","Combination", "speciesnumber")


save(pelargonium.combined, file = "Pelargonium_combined.Rda")
load("Pelargonium_combined.Rda")


#########################################################################
############### Plot TangentSpace to view results of GPA ################ 
#########################################################################



### The TangentSpace uses a PCA analysis to visualise the shape variation among the individuals in the
### GPA aligned dataset along their principal axes. 
### We can colour individuals using grouping variables as defined in the imported .csv file. 


# combined dataset


PCA_COMBINED <- plotTangentSpace(pelargonium.combined$Coords, axis1 = 1, axis2 = 2, 
                             warpgrids = FALSE, mesh = NULL, legend = FALSE,
                             label = pelargonium.combined$Species)









########################################
############### FIGURE 4 ###############
########################################

### PCA plot of all individuals coloured by main clade

data.C2F4.initial <- as.data.frame(PCA_COMBINED$pc.scores[,1:6])
data.C2F4 <- cbind(data.C2F4.initial, pelargonium.combined$Species)

data.C2F4$colour.COMBINED <- "a"
data.C2F4[data.C2F4$`pelargonium.combined$Species` %in% "P. mutans",]$colour.COMBINED <- "P. mutans"
data.C2F4[data.C2F4$`pelargonium.combined$Species` %in% "P. multibracteatum",]$colour.COMBINED <- "P. multibracteatum"
data.C2F4[data.C2F4$`pelargonium.combined$Species` %in% "P. myrrhifolium",]$colour.COMBINED <- "P. myrrhifolium"
data.C2F4[data.C2F4$`pelargonium.combined$Species` %in% "P. crithmifolium",]$colour.COMBINED <- "P. crithmifolium"
data.C2F4[data.C2F4$`pelargonium.combined$Species` %in% "P. pseudoglutinosum",]$colour.COMBINED <- "P. pseudoglutinosum"
data.C2F4[data.C2F4$`pelargonium.combined$Species` %in% "P. mutans",]$colour.COMBINED <- "P. mutans"
data.C2F4[data.C2F4$`pelargonium.combined$Species` %in% "P. patulum",]$colour.COMBINED <- "P. patulum"
data.C2F4[data.C2F4$`pelargonium.combined$Species` %in% "P. crispum",]$colour.COMBINED <- "P. crispum"

colourCount.COMBINED = length(unique(data.C2F4$colour.COMBINED))
getPalette.COMBINED = colorRampPalette(wes_palette(name="Zissou1"))

Combined <-
  ggplot(data = data.C2F4,
         mapping = aes(x = PC3,
                       y = PC2,
                       group = colour.COMBINED,
                       colour = colour.COMBINED)) +
  geom_point(size = 3, colour = "black", alpha = 0.4) +
  geom_point(data = subset(data.C2F4, colour.COMBINED != "a"),
             size = 8, alpha = 0.5) +
  scale_color_manual(values = getPalette.COMBINED(colourCount.COMBINED)) +
  geom_text(data=subset(data.C2F4, colour.COMBINED != "a"),
            aes(label=colour.COMBINED), 
            size=3,
            nudge_x = 0.02,
            nudge_y = 0.02) +
  #geom_text(aes(label=pelargonium.combined$Species), 
  #          size=3,
  #          nudge_x = 0.02,
  #          nudge_y = 0.02) +
  theme_bw(base_size = 20) +
  theme(legend.position = "none") +
  labs(x = "PC3 (14%)", y = "PC2 (19%)")
Combined


#############################################
############### Create shapes ###############
#############################################
ref <- mshape(pelargonium.combined.GPA$coords)

pellinks <- read.csv("Links_VIRTUAL3D.csv")



plotRefToTarget(ref,pelargonium.co,
                               gridPars=gridPar(tar.pt.bg = "#465670", tar.link.col="#465670",tar.link.lwd=2,
                                                pt.size = 0.1, tar.pt.size = 0.3), 
                               method="points", 
                               links = pellinks)

 

ref.VIRTUAL <- mshape(pelargonium.combined.GPA$coords)
pellinks.VIRTUAL <- read.csv("Links_VIRTUAL3D.csv")

plotRefToTarget(ref.VIRTUAL,PCA_COMBINED$pc.shapes$PC2min,
                gridPars=gridPar(pt.bg = "#696969",
                                 link.col = "#696969",
                                 tar.link.col="black",
                                 tar.link.lwd=5,
                                 pt.size = 0.1, 
                                 tar.pt.bg = "black",
                                 out.col = "#696969",
                                 tar.pt.size = 0.1,
                                 link.lwd = 3), 
                method = "points",
                links = pellinks.VIRTUAL,
                label = F)
rgl.postscript("PC3max.pdf",fmt="pdf")