#Morphological variation of long-legged velvet mite larvae in deep time

#Sofía I. Arce, Florian Braig, Simon J. Linhart, Jehan LeCadre, Olympia Salvamoser, Patrick Müller, Carolin Haug, Joachim T. Haug.

#Supplementary materials, File 2: R-Code used for the outline analysis and quantification.

#Required packages
library(Momocs)
library(dispRity)
library(ggplot2)
library(RColorBrewer)
library(readxl)

set.seed(1234)

### Import data ----------------------------------------------------------------

lf <- list.files("FILE_PATH", 
                 full.names = TRUE)
coo <- import_jpg(lf, auto.notcentered = TRUE, fun.notcentered = NULL, 
                  threshold = 0.8) 

Momocs_mites_1 <- read.csv("DATA_SHEET.csv", header=T, sep=",")
tail(Momocs_mites_1)
Momocs_mites <- subset(Momocs_mites_1, included.in.analysis=="yes")

df1 <- data.frame(fauna = as.factor(Momocs_mites$fauna),
                  drawn_by = as.factor(Momocs_mites$drawn_by),
                  parasite_free = as.factor(Momocs_mites$parasite_free),
                  INGROUP_II = as.factor(Momocs_mites$INGROUP_II),
                  id = as.factor(Momocs_mites$ID))


comp <- Out(coo, fac = df1)
panel(comp, names = TRUE, cex.names = 0.5)


# ELLIPTIC FOURIER ANALYSIS AND PCA --------------------------------------------

stack(comp, first.point = TRUE)
comp <- coo_center(comp)
stack(comp, first.point = TRUE)
comp <- coo_scale(comp)
stack(comp, first.point = TRUE)
comp.l <- coo_slidedirection(comp, direction = "down")
stack(comp.l, first.point = TRUE)
inspect(comp.l)

calibrate_harmonicpower_efourier(comp.l, id = 1:99, 100, drop = 1, 
                                 thresh = c(95, 99, 99.9))
calibrate_deviations_efourier(comp.l, range = c(5, 9, 16, 37), 
                              norm.centsize = TRUE, dist.method = edm_nearest, 
                              interpolate.factor = 1, dist.nbpts = 60)

comp.e <- efourier(comp.l, nb.h = 16, norm = FALSE, start = TRUE) 
boxplot(comp.e, drop = 1)
comp.pca <- PCA(comp.e)
scree(comp.pca, nax = c(1:20))
comp.scores <- comp.pca %>% as_df(13)
PCcontrib(comp.pca, nax = c(1:5))

# PLOTS OF THE MORPHOSPACE ------------------------------------------------------

#Simple plot of morphospace with according shapes in background

plot(comp.pca, fac = "fauna", col = c("orange", "cyan"), 
     fill = c("orange", "cyan"))

plot(comp.pca, fac = "parasite_free")


# Plot with every data point having its respective shape, custom colours, custom 
# polygons and custom axis

# extant Vs. fossils - PC 1 and PC 2
plot_PCA(comp.pca, ~ fauna, points = TRUE, chull = TRUE, chullfilled = FALSE,
         axes = c(1, 2), zoom = 1, morphospace_position = "xy",
         palette = pal_manual(c('red', 'orange')))

plot_PCA(comp.pca, ~ fauna, points = TRUE, chull = TRUE, chullfilled = FALSE,
         axes = c(3, 4), zoom = 1, 
         palette = pal_manual(c('red', 'orange')))


####### using ggplot

###

# extant Vs. fossils - PC 1 and PC 2
ggplot(comp.scores, aes(x = PC1, y = PC2, color = fauna, shape= fauna)) + 
  geom_point(size = 3) +
  #scale_fill_hue() +
  #scale_fill_hue(c=45, l=80) +
  scale_color_manual(legend_title, values = c("#d7191c", "#ffd700")) +
  scale_shape_manual(legend_title, values = c(19, 17)) +
  coord_fixed (ratio = 1, xlim = NULL, ylim = NULL, expand = TRUE) +
  xlab("PC1 (40.2%)") + ylab("PC2 (26.6%)") + 
  theme(lab = element_text(size = 20)) + 
  theme_light(base_size = 10)

# extant Vs. fossils - PC 3 and PC 4

ggplot(comp.scores, aes(x = PC3, y = PC4, color = fauna, shape= fauna)) + 
  geom_point(size = 3) +
  #scale_fill_hue() +
  #scale_fill_hue(c=45, l=80) +
  scale_color_manual(legend_title, values = c("#d7191c", "#ffd700")) +
  scale_shape_manual(legend_title, values = c(19, 17)) +
  coord_fixed (ratio = 1, xlim = NULL, ylim = NULL, expand = TRUE) +
  xlab("PC3") + ylab("PC4") + 
  theme(lab = element_text(size = 20)) + 
  theme_light(base_size = 10)



### to check ID in the plot

plot_label <- ggplot(comp.scores, aes(x = PC1, y = PC2, colour=fauna,  label=id)) + 
  geom_text(aes(label = id), hjust = 0.5, vjust = 0.5, size = 3) +
  geom_point(size = 2) 


# DISPRITY ANALYSIS ------------------------------------

data <- data.frame(comp.scores[, 6:15])
group.v <- data.frame(df1$fauna)
rownames(data) <- 1:nrow(data)
comp_subsets <- custom.subsets(data, group = group.v)
minimum_size <- min(size.subsets(comp_subsets))
comp_bootstrapped <- boot.matrix(data = comp_subsets, rarefaction = minimum_size,
                                   bootstraps = 10000)

disp.metrics <- test.metric(comp_bootstrapped, metric = c(mean, displacements), 
                            shifts = c("random", "size"))
summary(disp.metrics)
plot(disp.metrics)

comp_pos <- dispRity(data = comp_bootstrapped, 
                        metric = c(mean, displacements))
summary(comp_pos)
plot(comp_pos)

test.dispRity(comp_pos, test = t.test,  rarefaction = minimum_size)
test.dispRity(comp_pos, test = wilcox.test,  rarefaction = minimum_size)


disp.metrics <- test.metric(comp_bootstrapped, metric = c(sum, variances), 
                            shifts = c("random", "size"))
summary(disp.metrics)
plot(disp.metrics)

comp_size <- dispRity(data = comp_bootstrapped, 
                    metric = c(sum, variances))
summary(comp_size)
plot(comp_size)

test.dispRity(comp_size, test = t.test,  rarefaction = minimum_size)
test.dispRity(comp_size, test = wilcox.test,  rarefaction = minimum_size)