#Analysing flower shape variation on the VIRTUAL3D virtual3D dataset

# Script written by Sara van de Keker (saravandekerke@outlook.com)
# Written December 2018, updated October 2019
# Accompanying van de Kerke et al. 2019 Capturing variation in floral shape; a virtual3D based morphospace for Pelargonium, PeerJ preprint

library(geomorph) # using version 3.1.2
library(ggplot2) # using version 3.0.0
library(wesanderson) # using version 0.3.6

### set working directory 
setwd("C:/Users/kerke001/OneDrive - WageningenUR/morphometrics_chapter_2-3/Rscripten/")



###########################################
############### Data import ###############
# Import the datalist created in the earlier steps (S3 Supplementary R script Rescale&Combine) 
load("pelargonium_VIRTUAL3D_initial.Rda")


# Import variablesfile containing only individuals for those species that show overlap between the datasets and variablesfile containing
# VIRTUAL3D variable info, and rename columns
classifier.variables.PETAL.matching <- read.csv("Variables_PETAL_overlap.csv")
classifier.individuals <- read.csv("individuals68species3.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 = "_")



#####################################################################
############### Generalised Procrustes Analysis (GPA) ############### 
# Final GPA analysis on the VIRTUAL3D virtual3D dataset
pelargonium.VIRTUAL3D.GPA <- gpagen(pelargonium.VIRTUAL3D.initial$Coords, 
                                   surfaces=NULL, PrinAxes=TRUE, max.iter=NULL, 
                                   ProcD=TRUE, Proj=TRUE, print.progress = TRUE)


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



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


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

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

save(pelargonium.VIRTUAL3D.anova, file = "Pelargonium_VIRTUAL3D_ANOVA.Rda")


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


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





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


plotAllSpecimens(pelargonium.VIRTUAL3D.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.VIRTUAL3D.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 VIRTUAL3D in a list to facilitate Geomorph. We combine the coordinate data of the VIRTUAL3D dataset with numerous 
# variable data 

# VIRTUAL3D dataset
pelargonium.VIRTUAL3D <- list(pelargonium.VIRTUAL3D.GPA$coords, 
                             pelargonium.VIRTUAL3D.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.VIRTUAL3D) <- c("Coords", "Csize", "Clade", "Pollinator", 
                                 "PetalNumber", "number", "ID", 
                                 "Species","Combination", "speciesnumber")


save(pelargonium.VIRTUAL3D, file = "Pelargonium_VIRTUAL3D.Rda")
load("Pelargonium_VIRTUAL3D.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. 


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

save(PCA_VIRTUAL3D, file = "PCA_VIRTUAL3D.Rda")
load("PCA_VIRTUAL3D.Rda")


########################################
############### FIGURE 5 ###############
########################################

##### PCA plots #####

# Convert geomorph dataframe containing PCA results to a standard dataframe suitable for ggplot
data.C2F4.initial <- as.data.frame(PCA_VIRTUAL3D$pc.scores[,1:6])

# Add Species name to the dataframe
data.C2F4 <- cbind(data.C2F4.initial, pelargonium.VIRTUAL3D$Species)

# We add an extra variable to be able to highligh the individuals of specific species
data.C2F4$colour.VIRTUAL3D <- "a"
data.C2F4[data.C2F4$`pelargonium.VIRTUAL3D$Species` %in% "P. mutans",]$colour.VIRTUAL3D <- "P. mutans"
data.C2F4[data.C2F4$`pelargonium.VIRTUAL3D$Species` %in% "P. multibracteatum",]$colour.VIRTUAL3D <- "P. multibracteatum"
data.C2F4[data.C2F4$`pelargonium.VIRTUAL3D$Species` %in% "P. myrrhifolium",]$colour.VIRTUAL3D <- "P. myrrhifolium"
data.C2F4[data.C2F4$`pelargonium.VIRTUAL3D$Species` %in% "P. crithmifolium",]$colour.VIRTUAL3D <- "P. crithmifolium"
data.C2F4[data.C2F4$`pelargonium.VIRTUAL3D$Species` %in% "P. pseudoglutinosum",]$colour.VIRTUAL3D <- "P. pseudoglutinosum"
data.C2F4[data.C2F4$`pelargonium.VIRTUAL3D$Species` %in% "P. mutans",]$colour.VIRTUAL3D <- "P. mutans"
data.C2F4[data.C2F4$`pelargonium.VIRTUAL3D$Species` %in% "P. patulum",]$colour.VIRTUAL3D <- "P. patulum"
data.C2F4[data.C2F4$`pelargonium.VIRTUAL3D$Species` %in% "P. crispum",]$colour.VIRTUAL3D <- "P. crispum"

# Create a colour palette based on the number of species to be highlighted in the new variable above
colourCount.VIRTUAL3D = length(unique(data.C2F4$colour.VIRTUAL3D))
getPalette.VIRTUAL3D = colorRampPalette(wes_palette(name="Zissou1"))

# Create PCA plot for VIRTUAL3D dataset (Change x and y values to show PC axes)
VIRTUAL3D <-
  ggplot(data = data.C2F4,
         mapping = aes(x = PC1,
                       y = PC2,
                       group = colour.VIRTUAL3D,
                       colour = colour.VIRTUAL3D)) +
  geom_point(size = 3, colour = "black", alpha = 0.4) +
  geom_point(data = subset(data.C2F4, colour.VIRTUAL3D != "a"),
             size = 8, alpha = 0.5) +
  scale_color_manual(values = getPalette.VIRTUAL3D(colourCount.VIRTUAL3D)) +
  geom_text(data=subset(data.C2F4, colour.VIRTUAL3D != "a"),
            aes(label=colour.VIRTUAL3D), 
            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%)")
VIRTUAL3D


##### Outlineframes #####
# Create the outlineframes illustrating the shapes on the extremes of PC axes
ref <- mshape(pelargonium.VIRTUAL3D.GPA$coords)

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

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

plotRefToTarget(ref.VIRTUAL,PCA_VIRTUAL3D$pc.shapes$PC1min,
                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")