######################################################################### ### Morphometric Analysis from Hartke, Sprenger et al. for Crematogaster ######################################################################### rm(list=ls()); ls() # tyding up workspace # setwd("PATH") ######################################################################### ### Read file with functions source("_Morphometrics_functions.R", chdir = TRUE) ######################################################################### ### Read dataset data <- "Crematogaster_morphometry_measurements.csv" dat0 <- read.csv(data) dat <- dat0 ######################################################################### ### Check data reliability R specimen <- dat$specimen X <- dat[, c(6:28)] reliability <- icc.R(g=specimen, x=X) Rtab <- reliability$Rtable # Supplementary table of reliability R, sort by increasing R (Rtab <- Rtab[with(Rtab, order(R)),]) # Supplementary Figure of reliability R, sort by increasing R R <- Rtab$R; names <- Rtab$varname barplot(Rtab$R, main="Crematogaster: Reliability of characters (ordered by increasing R)", ylab="% reliability (R)", ylim=c(1, 100), names.arg=Rtab$varname, cex.names=1) ######################################################################### ### Aggregate duplicate values by calculating their mean dat <- aggregate(dat[,c(c(6:28))], by=dat[c("specimen", "chemotype", "location")], FUN=mean) ######################################################################### ### Select numeric variables for the analysis X <- dat[, c(4:26)]; X <- X[, order(names(X))] X <- X[, -c(11,18,21)] # removal of spo.d, ppt.l and omc.d with R<85% str(X) ######################################################################### ### Select shapes and colors for graphs shapes <- c(A = 17, B = 16) Cremblue <- rgb(red=19/255,green=82/255,blue=139/255) Crempurple <- rgb(red=164/255,green=92/255,blue=208/255) colours <- c("A"=Cremblue, "B"=Crempurple) ######################################################################### ### Select variables for calculations chemotype <- factor(dat$chemotype); specimen <- dat$specimen; location <- dat$location ######################################################################### ### Calculate isosize and a shape PCA spca <- shapePCA(X, rpc=3) loadings <- round(spca$loadings,2); spca$components <- round(spca$components, 3); pc1_s <- spca$components[2,1]*100; pc2_s <- spca$components[2,2]*100; pc3_s <- spca$components[2,3]*100 isosize <- as.numeric(spca$isosize) shapePC1 <- spca$PCmatrix[,1] shapePC2 <- spca$PCmatrix[,2] df <- data.frame(isosize, shapePC1, shapePC2, chemotype, specimen) ######################################################################### ### Create plots for Figure 4 # Figure 4A plot <- ggplot(df, aes(isosize, shapePC1, shape = chemotype, color = chemotype)) plot <- plot + geom_point(size=3) + theme(aspect.ratio=1) + labs(title="") plot <- plot + geom_smooth(method=lm, se=FALSE, size=1.5) plot <- plot + scale_shape_manual(values=shapes) plot <- plot + scale_colour_manual(values=colours) plot <- plot + xlab("Isometric size") + ylab(paste("Shape PC1 (",pc1_s,"%)", sep="")) plot <- plot + labs(colour = "Chemotype", shape= "Chemotype") a <- plot + theme(axis.title.x = element_text(size=16), axis.title.y = element_text(size=16), plot.title = element_text(size=16), legend.text = element_text(face="italic", size=13)) a # Figure 4B plot <- ggplot(df, aes(isosize, shapePC2, shape = chemotype, color = chemotype)) plot <- plot + geom_point(size=3) + theme(aspect.ratio=1) + labs(title="") plot <- plot + geom_smooth(method=lm, se=FALSE, size=1.5) plot <- plot + scale_shape_manual(values=shapes) plot <- plot + scale_colour_manual(values=colours) plot <- plot + xlab("Isometric size") + ylab(paste("Shape PC2 (",pc2_s,"%)", sep="")) plot <- plot + labs(colour = "Chemotype", shape= "Chemotype") b <- plot + theme(axis.title.x = element_text(size=16), axis.title.y = element_text(size=16), plot.title = element_text(size=16), legend.text = element_text(face="italic", size=13)) b ######################################################################### ### PCA ratio spectrum pcaRS(X, pc=1, sequence=1, labelsize=1.5, barlwd=5, nosize=0.8, bootrep=1000, barcol="#1F78B4") pcaRS(X, pc=2, sequence=1, labelsize=1.5, barlwd=5, nosize=0.8, bootrep=1000, barcol="#1F78B4") ######################################################################### ### Statistical tests # Test for differences in shape SPCAdata=cbind(shapePC1, shapePC2) M=manova(SPCAdata~chemotype+location) summary(M) summary.aov(M) # Test for differences in size t.test(isosize~chemotype)