# Analysing Flower shape variation on SPUR and PETAL datasets separate # 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(scales) # using version ... ### set working directory setwd("D:/PhD/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.SPUR.raw <- readland.tps("werkbestand_SPUR_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 inbetween ### 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.front.raw) is.matrix(pelargonium.front.raw) is.array(pelargonium.front.raw) View(pelargonium.front.raw) ################################################################################## ############### How to import a .csv file containing variable data ############### ### 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.SPUR.complete <- read.csv("variables_SPUR_complete_2.csv", header=TRUE) classifier.variables.SPUR.matching <- read.csv("variables_SPUR_overlap.csv", header=TRUE) ### Check the content of the imported .csv file (#1) or make it pop-up (#2) show() View() ### Check specific parts of the imported .csv file attributes() ### returns a list of components in the .csv file dimnames() ### returns a list of all dimension names (rows and columns) in the file ...$... ### returns a list of all states of that variable for each individual ### Check if a variable is a factor use (factor is a different word for category) ### and if not save as factor: class(...$...) ...$....converted <- as.factor(...$...) ###################################################### ############### Combine data in a list ############### ###################################################### pelargonium.PETAL.initial <- list(pelargonium.PETAL.raw, classifier.variables.PETAL.complete$Species) names(pelargonium.PETAL.initial) <- c("Coords.raw", "Species") # side dataset pelargonium.SPUR.initial <- list(pelargonium.SPUR.raw, classifier.variables.SPUR.complete$Species) names(pelargonium.SPUR.initial) <- c("Coords.raw", "Species") save(pelargonium.SPUR.initial, file = "Pelargonium_SPUR_initial.Rda") save(pelargonium.PETAL.initial, file = "Pelargonium_PETAL_initial.Rda") ##################################################################### ############### Generalised Procrustes Analysis (GPA) ############### ##################################################################### pelargonium.SPUR.initial.GPA <- gpagen(pelargonium.SPUR.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.SPUR.initial.GPA, file = "Pelargonium_SPUR_initial_GPA.Rda") save(pelargonium.PETAL.initial.GPA, file = "Pelargonium_PETAL_initial_GPA.Rda") load("Pelargonium_SPUR_initial_GPA.Rda") load("Pelargonium_PETAL_initial_GPA.Rda") ################################################## ############### Check for outliers ############### plotOutliers(pelargonium.SPUR.initial.GPA$coords, inspect.outliers = FALSE) outliers <- plotOutliers(pelargonium.SPUR.initial.GPA$coords) plotRefToTarget(mshape(pelargonium.SPUR.initial.GPA$coords), pelargonium.SPUR.initial.GPA$coords[,,outliers[128]], method="vector", label = T) ######################################################## ############### PCA on separate datasets ############### ### SPUR PCA_SPUR <- plotTangentSpace(pelargonium.SPUR.initial.GPA$coords, axis1 = 1, axis2 = 2, warpgrids = FALSE, mesh = NULL, legend = FALSE) #,label = pelargonium.SPUR.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_SPUR, file = "PCA_SPUR.Rda") save(PCA_PETAL, file = "PCA_PETAL.Rda") load("PCA_PETAL.Rda") ########################################### ############### Barplot PCA ############### ### SPUR pvar_SPUR <- (PCA_SPUR$sdev^2)/(sum(PCA_SPUR$sdev^2)) names(pvar_SPUR) <- seq(1:length(pvar_SPUR)) barplot(pvar_SPUR, 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 ############### ### 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") ### SPUR pelargonium.SPUR.rotated <- list(pelargonium.SPUR.raw, classifier.variables.SPUR.complete$Species, pelargonium.SPUR.initial.GPA) names(pelargonium.SPUR.rotated) <- c("Coords.raw", "Species", "Coords.GPA") ### Save save(pelargonium.PETAL.rotated, file = "pelargonium_PETAL_rotated.Rda") save(pelargonium.SPUR.rotated, file = "pelargonium_SPUR_rotated.Rda") load("pelargonium_PETAL_rotated.Rda") load("pelargonium_SPUR_rotated.Rda") ################################################## ############### Chapter 2 FIGURE 3 ############### ################################################## links <- read.csv("Links_SPUR.csv") links.PETAL <- read.csv("Links_PETAL.csv") p <- plotAllSpecimens(pelargonium.SPUR.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)) ################################################## ############### Chapter 2 FIGURE 4 ############### #SPUR classifier.variables.SPUR.complete$File_name[classifier.variables.SPUR.complete$File_name %in% dimnames(pelargonium.SPUR.raw)[[3]]==F] data.C2F3.raw <- as.data.frame(PCA_SPUR$pc.scores[,1:6]) data.C2F3.initial <- cbind(data.C2F3.raw, lineardists[,3]) data.C2F3.sorted <- data.C2F3.initial[order(dimnames(data.C2F3.initial)[[1]]),] classifier.SPUR.sorted <- classifier.variables.SPUR.complete[order(classifier.variables.SPUR.complete$File_name),] data.C2F3 <- cbind(data.C2F3.sorted, classifier.SPUR.sorted$Species) data.C2F3$colour.SPUR <- "a" data.C2F3[data.C2F3$`classifier.SPUR.sorted$Species` %in% "P. triste",]$colour.SPUR <- "P. triste" data.C2F3[data.C2F3$`classifier.SPUR.sorted$Species` %in% "P. crithmifolium",]$colour.SPUR <- "P. crithmifolium" data.C2F3[data.C2F3$`classifier.SPUR.sorted$Species` %in% "P. pseudoglutinosum",]$colour.SPUR <- "P. pseudoglutinosum" data.C2F3[data.C2F3$`classifier.SPUR.sorted$Species` %in% "P. mutans",]$colour.SPUR <- "P. mutans" data.C2F3[data.C2F3$`classifier.SPUR.sorted$Species` %in% "P. patulum",]$colour.SPUR <- "P. patulum" colourCount.SPUR = length(unique(data.C2F3$colour.SPUR)) getPalette.SPUR = colorRampPalette(wes_palette(name="Darjeeling1")) SPUR <- ggplot(data = data.C2F3, mapping = aes(x = PC1, y = PC2, group = colour.SPUR, colour = colour.SPUR)) + geom_point(size = 3, colour = "black", alpha = data.C2F3$`lineardists[, 3]`) + geom_point(data = subset(data.C2F3, colour.SPUR != "a"), size = 8) + scale_color_manual(values = getPalette.SPUR(colourCount.SPUR)) + geom_text(data=subset(data.C2F3, colour.SPUR != "a"), aes(label=colour.SPUR), size=3, nudge_x = 0.02, nudge_y = 0.02) + #geom_text(aes(label=classifier.SPUR.sorted$Species), # 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%)") SPUR #PETAL classifier.variables.PETAL.complete$ID[classifier.variables.SPUR.complete$ID %in% dimnames(pelargonium.PETAL.raw)[[3]]==F] dimnames(pelargonium.PETAL.raw)[[3]][dimnames(pelargonium.PETAL.raw)[[3]] %in% classifier.variables.SPUR.complete$ID==F] data.C2F3.initial <- as.data.frame(PCA_PETAL$pc.scores[,1:6]) 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) 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" colourCount.PETAL = length(unique(data.C2F3$colour.PETAL)) getPalette.PETAL = colorRampPalette(wes_palette(name="Cavalcanti1")) 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) + #geom_text(aes(label=classifier.PETAL.sorted$Species), # 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 ref <- mshape(pelargonium.SPUR.initial.GPA$coords) pellinks <- read.csv("Links_SPUR.csv") plotRefToTarget(ref,PCA_SPUR$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) 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] ############################################# ############### Create shapes ############### ############################################# M <- mshape(pelargonium.SPUR.initial.GPA$coords) PC <- PCA_SPUR$pc.scores[,1:2] preds <- shape.predictor(pelargonium.SPUR.initial.GPA$coords, x= PC, Intercept = FALSE, pred1 = c(0.2,0.3)) plot <- plotRefToTarget(M,preds$pred1, 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) ################################################## ############### Pedicel/Spur ratio ############### ################################################## ### Calculate interlandmark distance # Make a matrix defining three interlandmark distances lmks <- data.frame(SPUR = c(1,2), PEDICEL = c(6,10), row.names = c("start", "end")) A <- pelargonium.SPUR.raw lineardists <- interlmkdist(A, lmks) lineardists[,1][lineardists[,1]>lineardists[,2]] <- lineardists[,2][lineardists[,1]>lineardists[,2]] lineardists.ratio.initial <- lineardists[,1]/lineardists[,2] lineardists<- cbind(lineardists, lineardists.ratio.initial)