library(geomorph)
library(Rmisc)
library(ggplot2)
library(dplyr)
library(lme4)
library(lsmeans)
library(multcompView)

data_path = "./data/"
plot_path = './plots/csize/'
dir.create(plot_path, recursive = TRUE)

#read in TPS file
#writeland.tps(landmarks, "PCA_landmarks_scale.TPS", scale = 0.000407, specID = T)
landmarks <- readland.tps(paste0(data_path, "PCA_landmarks_scale.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)

#add all columns together
gmm <- list(landmarks = landmarks, sample = identifier$Sample, 
            tooth.position = identifier$Tooth.position,
            treatment = identifier$Treatment)

# Generalized procrustes analysis
procrustes <- gpagen(gmm$landmarks, curves = curves, surfaces = NULL, PrinAxes = TRUE,
                     max.iter = NULL, ProcD = TRUE, Proj = TRUE, print.progress = TRUE)


size <- data.frame(Size = procrustes$Csize, Treatment = identifier$Treatment,
                   Sample = identifier$Sample, stringsAsFactors = TRUE)

# Test difference in centroid size between treatments with sample as co-variate
size$Treatment <- relevel(size$Treatment, ref = "DMSO")
size$Sample <- relevel(size$Sample, ref = "DMSO_1")
mod <- lm(Size ~ Treatment+Sample, data = size)

summary(mod)
anova(mod)

# Run post-hoc tukey test to test for pairwise differences between treatments
mod2 <- aov(Size ~ Treatment+Sample, data = size)
TukeyHSD(mod2)


# Plot centroid size values
SE.sum <- summarySE(data = size, measurevar = "Size", groupvars = "Treatment",
                    na.rm = FALSE, conf.interval = 0.95, .drop = TRUE)

plot.dat <- data.frame(SE = SE.sum$se, Size = SE.sum$Size, Treatment = SE.sum$Treatment)

size$Treatment <- factor(x = size$Treatment , levels = c("DMSO", "IWR-1-endo", "CHIR-99021"))

png(paste0(plot_path, 'centroid_size.png'), width = 10, height = 15, units = 'cm', res = 400)
ggplot(plot.dat, aes(x = Treatment, y = Size, colour = Treatment)) +
  geom_point(data = size, aes(x = Treatment, y = Size, colour = Treatment), size = 1.5, alpha = 0.1) +
  geom_point (size = 4) +
  geom_errorbar(aes(ymax = Size + SE,
                    ymin = Size - SE), width = 0.3) +
  labs(x = "Treatment", y = "Centroid Size (mm)") +
  theme_bw() +
  theme(legend.position = "none")+
  theme(axis.line = element_line(colour = "black"),
        panel.grid = element_blank(),
        panel.border = element_blank())
graphics.off()