# Analysing flower shape variation on TUBE and PETAL datasets separate

# 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(scales) # using version 0.5.0
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 and check ###############

# We can import both landmark coordinate data in a .tsp file, as well as additional 
# variable data from a .csv  file (Excel based).
# We can also import a .csv file indicating which landmarks are semi/sliding landmarks.



### Import TPS file with landmark data of the front dataset, 
### which has only 'normal' landmarks.

pelargonium.PETAL.raw <- readland.tps("werkbestand_PETAL_new.TPS", 
                                      specID = "image", warnmsg = TRUE)



### import TPS file with landmark data of the side dataset, 
### which has 'normal' and 'sliding' landmarks

pelargonium.TUBE.raw <- readland.tps("werkbestand_TUBE_new_update4.TPS", 
                                     specID = "ID", readcurves = TRUE, warnmsg = TRUE)



### Import .csv file containing sliding landmark info. 
### Make sure the file indicates for each sliding landmark in between 
### which landmarks it is placed.

curvepts <- read.csv("semilandmarks.csv", header = TRUE)
show(curvepts)



### Check the dimensions of the data. Our data is formatted in a (3D) array, because we have 
### coordinates for each individual separate.  

dim(pelargonium.TUBE.raw)
is.matrix(pelargonium.PETAL.raw)
is.array(pelargonium.PETAL.raw)
View(pelargonium.PETAL.raw)




###########################################################################
############### Import a .csv file containing variable data ###############
# For each of the datasets (TUBE and PETAL) we upload a separate file containing information 
# for the full dataset (_complete) and containing information for those species that show overlap 
# between the datasets (_overlap)

### Import .csv file for the front dataset
classifier.variables.PETAL.complete <- read.csv("Variables_PETAL_complete.csv", 
                                                header=TRUE)
classifier.variables.PETAL.matching <- read.csv("Variables_PETAL_overlap.csv")

### Import .csv file for the side dataset
classifier.variables.TUBE.complete <- read.csv("variables_TUBE_complete_2.csv", 
                                               header=TRUE)
classifier.variables.TUBE.matching <- read.csv("variables_TUBE_overlap.csv", header=TRUE)





######################################################
############### Combine data in a list ###############
######################################################
# We combine the raw coordinate data and the species names from the variables file in a list to 
# perform the first GPA and the PCA analyses for the separate datasets


pelargonium.PETAL.initial <- list(pelargonium.PETAL.raw, 
                                  classifier.variables.PETAL.complete$Species)
names(pelargonium.PETAL.initial) <- c("Coords.raw", "Species")



# side dataset
pelargonium.TUBE.initial <- list(pelargonium.TUBE.raw, 
                                 classifier.variables.TUBE.complete$Species) 
names(pelargonium.TUBE.initial) <- c("Coords.raw", "Species")




save(pelargonium.TUBE.initial, file = "Pelargonium_TUBE_initial.Rda")
save(pelargonium.PETAL.initial, file = "Pelargonium_PETAL_initial.Rda")


#####################################################################
############### Generalised Procrustes Analysis (GPA) ############### 
#####################################################################
# GPA analysis on the separate TUBE and PETAL datasets

pelargonium.TUBE.initial.GPA <- gpagen(pelargonium.TUBE.raw, curves=NULL, 
                                       surfaces=NULL, PrinAxes=TRUE, max.iter=10, 
                                       ProcD=TRUE, Proj=TRUE, print.progress = TRUE)



pelargonium.PETAL.initial.GPA <- gpagen(pelargonium.PETAL.raw, curves=NULL, 
                                        surfaces=NULL, PrinAxes=TRUE, max.iter=10, 
                                        ProcD=TRUE, Proj=TRUE, print.progress = TRUE)



save(pelargonium.TUBE.initial.GPA, file = "Pelargonium_TUBE_initial_GPA.Rda")
save(pelargonium.PETAL.initial.GPA, file = "Pelargonium_PETAL_initial_GPA.Rda")


load("Pelargonium_TUBE_initial_GPA.Rda")
load("Pelargonium_PETAL_initial_GPA.Rda")




##################################################
############### Check for outliers ###############



plotOutliers(pelargonium.TUBE.initial.GPA$coords, inspect.outliers = FALSE)

outliers <- plotOutliers(pelargonium.TUBE.initial.GPA$coords)

plotRefToTarget(mshape(pelargonium.TUBE.initial.GPA$coords),
                pelargonium.TUBE.initial.GPA$coords[,,outliers[128]],
                method="vector", label = T)





#######################################################################
############### PCA on separate TUBE and PETAL datasets ###############


### TUBE

PCA_TUBE <- plotTangentSpace(pelargonium.TUBE.initial.GPA$coords, axis1 = 1, axis2 = 2, 
                             warpgrids = FALSE, mesh = NULL, legend = FALSE)
#,label = pelargonium.TUBE.initial$Species, ) 


### PETAL


PCA_PETAL <- plotTangentSpace(pelargonium.PETAL.initial.GPA$coords, axis1 = 1, axis2 = 2, 
                              warpgrids = FALSE, mesh = NULL, legend = FALSE)
#label = pelargonium.PETAL.initial$Species,)

### Save

save(PCA_TUBE, file = "PCA_TUBE.Rda")
save(PCA_PETAL, file = "PCA_PETAL.Rda")

load("PCA_PETAL.Rda")





#################################################################
############### S6 Supplementary Figure PC scores ###############
# Making the barplots of for the TUBE and PETAL datasets 

### TUBE

pvar_TUBE <- (PCA_TUBE$sdev^2)/(sum(PCA_TUBE$sdev^2))
names(pvar_TUBE) <- seq(1:length(pvar_TUBE))
barplot(pvar_TUBE, xlab= "Principal Components", ylab = "% Variance")

### PETAL

pvar_PETAL <- (PCA_PETAL$sdev^2)/(sum(PCA_PETAL$sdev^2))
names(pvar_PETAL) <- seq(1:length(pvar_PETAL))
barplot(pvar_PETAL, xlab= "Principal Components", ylab = "% Variance")




######################################################
############### Combine data in a list ###############
# We combine the raw coordinate data, the species names from the variables file, as well as 
# the new GPA coordinates


### PETAL
pelargonium.PETAL.rotated <- list(pelargonium.PETAL.raw, 
                                  classifier.variables.PETAL.complete$Species,
                                  pelargonium.PETAL.initial.GPA)

names(pelargonium.PETAL.rotated) <- c("Coords.raw", "Species", "Coords.GPA")



### TUBE
pelargonium.TUBE.rotated <- list(pelargonium.TUBE.raw, 
                                 classifier.variables.TUBE.complete$Species,
                                 pelargonium.TUBE.initial.GPA)

names(pelargonium.TUBE.rotated) <- c("Coords.raw", "Species", "Coords.GPA")


### Save

save(pelargonium.PETAL.rotated, file = "pelargonium_PETAL_rotated.Rda")
save(pelargonium.TUBE.rotated, file = "pelargonium_TUBE_rotated.Rda")

load("pelargonium_PETAL_rotated.Rda")
load("pelargonium_TUBE_rotated.Rda")


#############################################################################################################
############### S5 Supplementary Figure mean consensus configuration and Procrustes residuals ###############
#############################################################################################################
# Load files indicating which landmarks have to be connected and plot all coordinates for all specimens

links <- read.csv("Links_TUBE.csv")
links.PETAL <- read.csv("Links_PETAL.csv")

p <- plotAllSpecimens(pelargonium.TUBE.initial.GPA$coords, label = F, links = links, mean = T,
                 plot.param = list(pt.bg = "grey", pt.cex = 1, pt.col = "#DEA109", pch = 5, 
                                   mean.cex=0.5, 
                                   link.col="#3182bd", link.lwd = 3,
                                   txt.pos=3, txt.cex=1.5))

p <- plotAllSpecimens(pelargonium.PETAL.initial.GPA$coords, label = F, links = links.PETAL,
                      plot.param = list(pt.bg = "grey", pt.cex = 1, pt.col = "#DEA109", pch = 5, 
                                        mean.cex=0.5, 
                                        link.col="#3182bd", link.lwd = 3,
                                        txt.pos=3, txt.cex=1.5))




##################################################
############### Pedicel/Tube ratio ###############
##################################################

### Calculate interlandmark distance

# Make a matrix defining interlandmark distances for the tube and pedicel
lmks <- data.frame(TUBE = c(1,2), PEDICEL = c(6,10), 
                   row.names = c("start", "end")) 
A <- pelargonium.TUBE.raw
lineardists <- interlmkdist(A, lmks)

# If measured value for tube is larger than measure value for the pedicel, we replace the tube value with the
# pedicel value because we know it is biologically not possible to have a larger tube.
lineardists[,1][lineardists[,1]>lineardists[,2]] <- lineardists[,2][lineardists[,1]>lineardists[,2]]

# Calculate the ratio between the pedicel and tube and add value to the matrix
lineardists.ratio.initial <- lineardists[,1]/lineardists[,2]
lineardists<- cbind(lineardists, lineardists.ratio.initial)



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

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

###TUBE

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

# Add pedicel/tube ratio value to the dataframe
data.C2F3.initial <- cbind(data.C2F3.raw, lineardists[,3])

# Order data and variables in alphabetical order and combine
data.C2F3.sorted <- data.C2F3.initial[order(dimnames(data.C2F3.initial)[[1]]),]
classifier.TUBE.sorted <- classifier.variables.TUBE.complete[order(classifier.variables.TUBE.complete$File_name),]
data.C2F3 <- cbind(data.C2F3.sorted, classifier.TUBE.sorted$Species)

# We add an extra variable to be able to highligh the individuals of specific species
data.C2F3$colour.TUBE <- "a"
data.C2F3[data.C2F3$`classifier.TUBE.sorted$Species` %in% "P. triste",]$colour.TUBE <- "P. triste"
data.C2F3[data.C2F3$`classifier.TUBE.sorted$Species` %in% "P. crithmifolium",]$colour.TUBE <- "P. crithmifolium"
data.C2F3[data.C2F3$`classifier.TUBE.sorted$Species` %in% "P. pseudoglutinosum",]$colour.TUBE <- "P. pseudoglutinosum"
data.C2F3[data.C2F3$`classifier.TUBE.sorted$Species` %in% "P. mutans",]$colour.TUBE <- "P. mutans"
data.C2F3[data.C2F3$`classifier.TUBE.sorted$Species` %in% "P. patulum",]$colour.TUBE <- "P. patulum"

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

# Create PCA plot for TUBE dataset
TUBE <-
  ggplot(data = data.C2F3,
         mapping = aes(x = PC1,
                       y = PC2,
                       group = colour.TUBE,
                       colour = colour.TUBE)) +
  geom_point(size = 3, colour = "black", alpha = data.C2F3$`lineardists[, 3]`) +
  geom_point(data = subset(data.C2F3, colour.TUBE != "a"),
             size = 8) +
  scale_color_manual(values = getPalette.TUBE(colourCount.TUBE)) +
  geom_text(data=subset(data.C2F3, colour.TUBE != "a"),
            aes(label=colour.TUBE), 
           size=3,
            nudge_x = 0.02,
            nudge_y = 0.02) +
  theme_bw(base_size = 20) +
  theme(legend.position = "none") +
  labs(x = "PC1 (47%)", y = "PC2 (28%)")
TUBE


###PETAL

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

# Order data and variables in alphabetical order and combine
data.C2F3.sorted <- data.C2F3.initial[order(dimnames(data.C2F3.initial)[[1]]),]
classifier.PETAL.sorted <- classifier.variables.PETAL.complete[order(classifier.variables.PETAL.complete$ID),]
data.C2F3.PETAL <- cbind(data.C2F3.sorted, classifier.PETAL.sorted$Species)

# We add an extra variable to be able to highligh the individuals of specific species
data.C2F3.PETAL$colour.PETAL <- "a"
data.C2F3.PETAL[data.C2F3.PETAL$`classifier.PETAL.sorted$Species` %in% "P. mutans",]$colour.PETAL <- "P. mutans"
data.C2F3.PETAL[data.C2F3.PETAL$`classifier.PETAL.sorted$Species` %in% "P. multibracteatum",]$colour.PETAL <- "P. multibracteatum"
data.C2F3.PETAL[data.C2F3.PETAL$`classifier.PETAL.sorted$Species` %in% "P. myrrhifolium",]$colour.PETAL <- "P. myrrhifolium"

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

# Create PCA plot for PETAL dataset
PETAL <-
  ggplot(data = data.C2F3.PETAL,
         mapping = aes(x = PC1,
                       y = PC2,
                       group = colour.PETAL,
                       colour = colour.PETAL)) +
  geom_point(size = 3, colour = "black") +
  geom_point(data = subset(data.C2F3.PETAL, colour.PETAL != "a"),
             size = 8) +
  scale_color_manual(values = getPalette.PETAL(colourCount.PETAL)) +
  geom_text(data=subset(data.C2F3.PETAL, colour.PETAL != "a"),
            aes(label=colour.PETAL), 
            size=3,
            nudge_x = 0.02,
            nudge_y = 0.02) +
  theme_bw(base_size = 20) +
  theme(legend.position = "none") +
  labs(x = "PC1 (40.5%)", y = "PC2 (13.1%)")
PETAL


##### Outlineframes #####
# Create the outlineframes illustrating the shapes on the extremes of PC axes

### TUBE
ref.TUBE <- mshape(pelargonium.TUBE.initial.GPA$coords)
pellinks.TUBE <- read.csv("Links_TUBE.csv")

plotRefToTarget(ref,PCA_TUBE$pc.shapes$PC1max,
                gridPars=gridPar(pt.bg = "grey",
                                 link.col = "grey",
                                 tar.link.col="black",
                                 tar.link.lwd=2,
                                 pt.size = 0.3, 
                                 tar.pt.bg = "green",
                                 out.col = "grey",
                                 tar.pt.size = 0.6),
                method = "points",
                links = pellinks)


### PETAL
ref.PETAL <- mshape(pelargonium.PETAL.initial.GPA$coords)
pellinks.PETAL <- read.csv("Links_PETAL.csv")

dev.new(width=50, height=11)
plotRefToTarget(ref.PETAL,PCA_PETAL$pc.shapes$PC1max,
                gridPars=gridPar(pt.bg = "grey",
                                 link.col = "grey",
                                 tar.link.col="black",
                                 tar.link.lwd=2,
                                 pt.size = 0.3, 
                                 tar.pt.bg = "green",
                                 out.col = "grey",
                                 tar.pt.size = 0.6), 
                method = "points",
                links = pellinks.PETAL,
                label = F)
save(PETAL_1max, file = "PETAL_1max.Rda")

pelargonium.PETAL.initial.GPA$coords[,,128]
                               


##################################################################
############### Extract symmetric component PETAL ################ 
##################################################################
# Check whether symmetry in floral shape is of any influence in our analyses by extracting the symmetric component and
# creating an additional PCA plot. 

landpairs <- read.csv("Pairs.csv")
classifier.variables.PETAL.complete <- read.csv("Variables_PETAL_complete.csv", 
                                                header=TRUE)

gdf.sym <- geomorph.data.frame(shape=pelargonium.PETAL.rotated$Coords.GPA$coords, ind = classifier.variables.PETAL.complete$ID)

PETAL.symm <- bilat.symmetry(A = shape, data = gdf.sym, object.sym = TRUE, land.pairs = landpairs)

save(PETAL.symm, file = "PETAL_symm.Rda")
load("PETAL_symm.Rda")

PCA_SYM <- plotTangentSpace(PETAL.symm$symm.shape, axis1 = 1, axis2 = 2, 
                            warpgrids = T, mesh = NULL, legend = FALSE)