library(geomorph) library(Morpho) library(Rmisc) library(ggplot2) library(dplyr) library(lme4) library(lsmeans) library(multcompView) data_path = "./data/" plot_path = './plots/procrustes_shape/' dir.create(plot_path, recursive = TRUE) #read in TPS file landmarks <- readland.tps(paste0(data_path, "PCA.txt"), specID = c("None", "ID", "imageID"), readcurves = TRUE, warnmsg = TRUE) #set semilandmarks in plots window #define.sliders(landmarks, 36, surfsliders = NULL, write.file = TRUE) #read in semilandmarks generated by define.sliders curves <- read.csv(paste0(data_path, "curveslide.csv")) #read identifier column identifier <- read.csv(paste0(data_path, "GMM_metadata.csv"), header=T) colours <- c('DMSO'='#f2f2f2', 'CHIR-99021'='#a6a6a6', 'IWR-1-endo'='black') identifier$colour = unname(colours[identifier$Treatment]) #add all columns together gmm <- list(landmarks = landmarks, sample = identifier$Sample, tooth.position = identifier$Tooth.position, treatment = identifier$Treatment, colour = identifier$colour) #procrustes superimposition of landmark values procrustes <- gpagen(gmm$landmarks, curves = curves, surfaces = NULL, PrinAxes = TRUE, max.iter = NULL, ProcD = TRUE, Proj = TRUE, print.progress = TRUE) summary(procrustes) # Plot procrustes summary png(paste0(plot_path, 'procrustes.png'), width = 15, height = 15, units = 'cm', res = 400) plot(procrustes) graphics.off() ##################################### #Variance test on PC points - not averaged for sample png(paste0(plot_path, 'tangent_space_all.png'), width = 15, height = 15, units = 'cm', res = 400) PC<-plotTangentSpace(procrustes$coords, groups = gmm$colour, warpgrids = T) graphics.off() var_data <- data.frame(PC1 = PC$pc.scores[,1], treatment = identifier$Treatment) var.test(x = filter(var_data, treatment == "DMSO")$PC1, y = filter(var_data, treatment == "CHIR-99021")$PC1) var.test(x = filter(var_data, treatment == "DMSO")$PC1, y = filter(var_data, treatment == "IWR-1-endo")$PC1) #Plot PCA1 and 2 for all measurements irrespective of sample to show variance xlab <- paste("Principal Component 1 ", "(", round(PC$pc.summary$importance[2,1]*100, 1), "%)", sep="") ylab <- paste("Principal Component 2 ", "(", round(PC$pc.summary$importance[2,2]*100, 1), "%)", sep="") png(paste0(plot_path, 'pca_all_teeth.png'), width = 20, height = 15, units = 'cm', res = 400) plot(PC$pc.scores[,1], PC$pc.scores[,2], pch=c('DMSO'=22, 'IWR-1-endo'=24, 'CHIR-99021'=21)[as.array(identifier$Treatment)], cex=1, bg = colours[as.array((identifier$Treatment))], col="black", xlab=xlab, ylab=ylab, axes=TRUE, cex.lab=1.2, cex.axis=1.2, cex.main=1.2, cex.sub=1.2, asp=T) legend("bottomright", legend= unique(identifier$Treatment), pch=c('DMSO'=22, 'IWR-1-endo'=24, 'CHIR-99021'=21)[as.array(unique(identifier$Treatment))], pt.bg = colours[as.array(unique(identifier$Treatment))], col="black", cex = 1, pt.cex = 1.1, y.intersp = 1, bty = "n") graphics.off() ############################################## ###subset procrustes data set and average based on sample metadata_average <- unique(identifier[,!names(identifier) == 'Tooth.position']) av_list <- lapply(metadata_average$Sample, function(x){ crr_sample = x positions = which(identifier$Sample == crr_sample) crr_average = apply(procrustes$coords[,,positions] , c(1,2), FUN=mean) return(crr_average) }) names(av_list) <- metadata_average$Sample #make array from average proc data array.from.list <- array(unlist(av_list), dim = c(nrow(av_list[[1]]), ncol(av_list[[1]]), length(av_list))) #PCA and plot average proc data png(paste0(plot_path, 'tangent_space_average.png'), width = 15, height = 15, units = 'cm', res = 400) plotTangentSpace(array.from.list, groups = metadata_average$colour, warpgrids = T) graphics.off() ################################### ################################### #Plot final warp figure on sample averages gdf.av <- geomorph.data.frame(shape = array.from.list, treatment = metadata_average$Treatment, sample = metadata_average$Sample, colour = metadata_average$colour) PCA.1 <- plotTangentSpace(gdf.av$shape, groups = gdf.av$colour, warpgrids = T) xlab <- paste("Principal Component 1 ", "(", round(PCA.1$pc.summary$importance[2,1]*100, 1), "%)", sep="") ylab <- paste("Principal Component 2 ", "(", round(PCA.1$pc.summary$importance[2,2]*100, 1), "%)", sep="") plot_data <- data.frame(PCA.1$pc.scores[,1:2], treatment = gdf.av$treatment) plot_data$treatment <- factor(plot_data$treatment, levels = c('DMSO', 'IWR-1-endo', 'CHIR-99021')) png(paste0(plot_path, 'pca_average.png'), width = 20, height = 15, units = 'cm', res = 400) ggplot(plot_data, aes(x = PC1, y = PC2, shape=treatment, fill = treatment)) + geom_point(size=5) + stat_ellipse(aes(colour = treatment), linetype = 2, size = 1, level = .90) + scale_shape_manual(values=c(22, 24, 21)) + theme_classic() + xlab(xlab) + ylab(ylab) + theme(legend.position = c(0.87, 0.15), legend.title = element_blank(), legend.key.size = unit(0.9, 'cm'), legend.text = element_text(size=16), panel.border = element_rect(colour = "black", fill=NA, size=1), axis.text=element_text(size=16), axis.title=element_text(size=18), axis.text.y = element_text(angle = 90, hjust = 0.5)) graphics.off() ref <- mshape(gdf.av$shape)# assign mean shape for use with plotRefToTarget below # Item 2 to plot, the first TPS grid; here we use the outline option to add to the visualisation png(paste0(plot_path, 'shape_x0.png'), width = 12, height = 12, units = 'cm', res = 400) plotRefToTarget(ref,PCA.1$pc.shapes$PC1min,outline=gdf.av$shape) graphics.off() png(paste0(plot_path, 'shape_x1.png'), width = 12, height = 12, units = 'cm', res = 400) plotRefToTarget(ref,PCA.1$pc.shapes$PC1max,outline=gdf.av$shape) graphics.off() png(paste0(plot_path, 'shape_y0.png'), width = 12, height = 12, units = 'cm', res = 400) plotRefToTarget(ref,PCA.1$pc.shapes$PC2min,outline=gdf.av$shape) graphics.off() png(paste0(plot_path, 'shape_y1.png'), width = 12, height = 12, units = 'cm', res = 400) plotRefToTarget(ref,PCA.1$pc.shapes$PC2max,outline=gdf.av$shape) graphics.off() ################################## #variance_test var_data <- data.frame(PC1 = PCA.1$pc.scores[,1], treatment = gdf.av$treatment) var.test(x = filter(var_data, treatment == "DMSO")$PC1, y = filter(var_data, treatment == "IWR-1-endo")$PC1) var.test(x = filter(var_data, treatment == "DMSO")$PC1, y = filter(var_data, treatment == "CHIR-99021")$PC1) #PCA and plot ALL proc data PCA <- plotTangentSpace(procrustes$coords, groups = identifier$colour, warpgrids = T) summary(PCA) ### linear model on procrustes data points reveals significant effect of treatment on overall tooth shape gdf <- geomorph.data.frame(shape = procrustes$coords, treatment = identifier$Treatment, sample = identifier$Sample) out <- procD.lm(shape ~ treatment + sample, data = gdf,iter=999, RRPP = F) summary(out) ###################################################