#leptopoma sexual dimorphism in size, shape and colour patterns


all_morphometric_data=read.csv(choose.files()) # Select and open Open Additional File 4.csv
"Additional File 1 Leptopoma perlucidum shell morphometric data.csv"

# R.4.0 working
library(ggplot2)
library(reshape2)
measurements_data_with_gender_locality <- melt(all_morphometric_data, id.vars = c("Sexes", "Locations"), measure.vars = c("Shell.Height",	"Shell.Width",	"Aperture.Height",	"Aperture.Width",	"Spire.Height"))
levels(measurements_data_with_gender_locality$variable)<-c("Shell.Height",	"Shell.Width",	"Aperture.Height",	"Aperture.Width",	"Spire.Height") 

### Two-Way ANOVA differences in size between male and female and among locations ----
## Boxplot ----
ggplot(measurements_data_with_gender_locality, aes(x = Locations, y = value, fill = Sexes)) + 
  geom_boxplot()+facet_wrap(~variable, scale="free", ncol=2)+
  labs(x="Locations", y= "Size (mm)", fill="Sexes") # 7 X 7.5





## Descriptive statistics ----

combined_factors <- paste(all_morphometric_data$Sexes,all_morphometric_data$Locations)

all_morphometric_data_combined_factors <- cbind(all_morphometric_data[2:8],combined_factors)

all_morphometric_data_combined_factors$Sexes <- factor(all_morphometric_data_combined_factors$Sexes)
all_morphometric_data_combined_factors$Locations <- factor(all_morphometric_data_combined_factors$Locations)
all_morphometric_data_combined_factors$combined_factors <- factor(all_morphometric_data_combined_factors$combined_factors)



str(all_morphometric_data_combined_factors)

summary_SH <- aggregate(all_morphometric_data_combined_factors$Shell.Height, by=list(all_morphometric_data_combined_factors$Sexes, all_morphometric_data_combined_factors$Locations), function(x) c(mean = mean(x), sd = sd(x)))
summary_SW <- aggregate(all_morphometric_data_combined_factors$Shell.Width, by=list(all_morphometric_data_combined_factors$Sexes, all_morphometric_data_combined_factors$Locations), function(x) c(mean = mean(x), sd = sd(x)))
summary_AH <- aggregate(all_morphometric_data_combined_factors$Aperture.Height, by=list(all_morphometric_data_combined_factors$Sexes, all_morphometric_data_combined_factors$Locations), function(x) c(mean = mean(x), sd = sd(x)))
summary_AW <- aggregate(all_morphometric_data_combined_factors$Aperture.Width, by=list(all_morphometric_data_combined_factors$Sexes, all_morphometric_data_combined_factors$Locations), function(x) c(mean = mean(x), sd = sd(x)))
summary_SpH <- aggregate(all_morphometric_data_combined_factors$Spire.Height, by=list(all_morphometric_data_combined_factors$Sexes, all_morphometric_data_combined_factors$Locations), function(x) c(mean = mean(x), sd = sd(x)))

summary_SH
summary_SW
summary_AH
summary_AW
summary_SpH



## Power analysis for two sample test
# sample size of each group
#"Because sampling was opportunistic and effect size (Cohen's d) 
#cannot be determined prior to sampling, we calculated effect size 
#based on the differences between male and female snails in the five 
#shell measurements for each location. The very large effect sizes 
#found in the measurements for both sexes on Tiga Island ranged from 
#0.75 to 1.20, while on Gaya Island, the values ranged from 1.15 to 1.92. 
#Therefore, we checked the power of the tests after the experiment had 
#been conducted by using the calculated effect size (d), significant level 
#of 0.05, and the smaller n for each group,  n = 30 and 9 for Tiga and 
#Gaya Island respectively."


library(pwr)

# two sample test (shell height - between sexes)

# Gaya Island
Gaya_SH <- pwr.t.test(n=9, d=1.92,sig.level=0.05,type="two.sample",alternative="two.sided")
Gaya_SW <- pwr.t.test(n=9, d=1.78,sig.level=0.05,type="two.sample",alternative="two.sided")
Gaya_AH <- pwr.t.test(n=9, d=1.68,sig.level=0.05,type="two.sample",alternative="two.sided")
Gaya_AW <- pwr.t.test(n=9, d=1.15,sig.level=0.05,type="two.sample",alternative="two.sided")
Gaya_SpH <- pwr.t.test(n=9, d=1.73,sig.level=0.05,type="two.sample",alternative="two.sided")
Gaya_SH
Gaya_SW
Gaya_AH
Gaya_AW
Gaya_SpH

# Tiga Island
Tiga_SH <- pwr.t.test(n=30, d=0.99,sig.level=0.05,type="two.sample",alternative="two.sided")
Tiga_SW <- pwr.t.test(n=30, d=1.20,sig.level=0.05,type="two.sample",alternative="two.sided")
Tiga_AH <- pwr.t.test(n=30, d=0.94,sig.level=0.05,type="two.sample",alternative="two.sided")
Tiga_AW <- pwr.t.test(n=30, d=1.20,sig.level=0.05,type="two.sample",alternative="two.sided")
Tiga_SpH <- pwr.t.test(n=30, d=0.75,sig.level=0.05,type="two.sample",alternative="two.sided")
Tiga_SH
Tiga_SW
Tiga_AH
Tiga_AW
Tiga_SpH


## Assumption test (Normality and homogeneity of variances tests)
morphometric_data_female_Gaya <- subset(all_morphometric_data_combined_factors, combined_factors == "female P. Gaya")
morphometric_data_female_Tiga <- subset(all_morphometric_data_combined_factors, combined_factors == "female P. Tiga")
morphometric_data_male_Gaya <- subset(all_morphometric_data_combined_factors, combined_factors == "male P. Gaya") 
morphometric_data_male_Tiga <- subset(all_morphometric_data_combined_factors, combined_factors == "male P. Tiga") 

# Aperture height ----
# normality test
shapiro.test(morphometric_data_female_Gaya$Aperture.Height)
shapiro.test(morphometric_data_female_Tiga$Aperture.Height)
shapiro.test(morphometric_data_male_Gaya$Aperture.Height)
shapiro.test(morphometric_data_male_Tiga$Aperture.Height)
# homogeneity of variances test
library(car)
leveneTest(all_morphometric_data_combined_factors$Aperture.Height~
              all_morphometric_data_combined_factors$combined_factors)


two_way_anova_Aperture.Height <- aov(Aperture.Height ~ Sexes + Locations + Sexes * Locations, data=all_morphometric_data_combined_factors)

summary(two_way_anova_Aperture.Height)


# Aperture Width ----
# normality test
shapiro.test(morphometric_data_female_Gaya$Aperture.Width)
shapiro.test(morphometric_data_female_Tiga$Aperture.Width)
shapiro.test(morphometric_data_male_Gaya$Aperture.Width)
shapiro.test(morphometric_data_male_Tiga$Aperture.Width)
# honogeneity of variances test
library(car)
leveneTest(all_morphometric_data_combined_factors$Aperture.Width~
             all_morphometric_data_combined_factors$combined_factors)

two_way_anova_Aperture.Width <- aov(Aperture.Width~Sexes *Locations, data=all_morphometric_data_combined_factors)

summary(two_way_anova_Aperture.Width)

# Shell Width----
# normality test
shapiro.test(morphometric_data_female_Gaya$Shell.Width)
shapiro.test(morphometric_data_female_Tiga$Shell.Width)
shapiro.test(morphometric_data_male_Gaya$Shell.Width)
shapiro.test(morphometric_data_male_Tiga$Shell.Width)
# honogeneity of variances test
library(car)
leveneTest(all_morphometric_data_combined_factors$Shell.Width~
             all_morphometric_data_combined_factors$combined_factors)

two_way_anova_Shell.Width <- aov(Shell.Width~Sexes *  Locations, data=all_morphometric_data_combined_factors)

summary(two_way_anova_Shell.Width)


# Spire Height----
# normality test
shapiro.test(morphometric_data_female_Gaya$Spire.Height)
shapiro.test(morphometric_data_female_Tiga$Spire.Height)
shapiro.test(morphometric_data_male_Gaya$Spire.Height)
shapiro.test(morphometric_data_male_Tiga$Spire.Height)
# honogeneity of variances test
library(car)
leveneTest(all_morphometric_data_combined_factors$Spire.Height~
             all_morphometric_data_combined_factors$combined_factors)

two_way_anova_Spire.Height<- aov(Spire.Height~Sexes * Locations, data=all_morphometric_data_combined_factors)

summary(two_way_anova_Spire.Height)


# Shell Height----
shapiro.test(morphometric_data_female_Gaya$Shell.Height)
shapiro.test(morphometric_data_female_Tiga$Shell.Height)
shapiro.test(morphometric_data_male_Gaya$Shell.Height)
shapiro.test(morphometric_data_male_Tiga$Shell.Height)

library(car)
leveneTest(all_morphometric_data_combined_factors$Shell.Height~
             all_morphometric_data_combined_factors$combined_factors)

two_way_anova_Shell.Height<- aov(Shell.Height~Sexes * Locations, data=all_morphometric_data_combined_factors)

summary(two_way_anova_Shell.Height)



#geomorph analysis -----
#ONly R-4.0.3 works
library(geomorph)
library(ggplot2)


all_morphometric_data=read.csv(choose.files()) # Select and open Open Additional File 4.csv
"Additional File 1 Leptopoma perlucidum shell morphometric data.csv"

shell_character_for_GM <- all_morphometric_data[,9:26]
names(shell_character_for_GM)

LM_matrix <- as.matrix(shell_character_for_GM)
is.numeric(LM_matrix)
LM_coords <- arrayspecs(LM_matrix, 9, 2)
class(LM_coords)

classifiers <- all_morphometric_data[,2:3]
str(classifiers)

GM_data <- list(LM_coords,classifiers$Sexes,classifiers$Locations)
names(GM_data) <- c("landmarks","Sexes","Locations")

leptopoma.gpa <- gpagen(GM_data$landmarks, print.progress = FALSE)
summary(leptopoma.gpa)
leptopoma.gpa$coords
leptopoma.gpa$Csize

plotAllSpecimens(leptopoma.gpa$coords)

PCA_Leptopoma <- gm.prcomp(leptopoma.gpa$coords)
summary(PCA_Leptopoma)

plot(PCA_Leptopoma)


PCA_score <- cbind(classifiers,data.frame(PCA_Leptopoma$x))

PCA_plot_Leptopoma_all <- ggplot(PCA_score, aes(x=Comp1, y=Comp2)) + 
  geom_point(aes(color = factor(Sexes), shape = factor(Locations)), size = 3) + 
  labs(x = "Principal Component 1 (31.8%)",y = "Principal Component 2 (23.0%)")
PCA_plot_Leptopoma_all # 7.65 X 6

PCA_plot_Leptopoma_by_location <- ggplot(PCA_score, aes(x=Comp1, y=Comp2)) + 
  geom_point(aes(color = factor(Sexes)), size = 3) + 
  labs(x = "Principal Component 1 (31.8%)",y = "Principal Component 2 (23.0%)") +
  facet_wrap(~Locations)
PCA_plot_Leptopoma_by_location #13 X 6


mesh_Leptopoma <- mshape(leptopoma.gpa$coords)
# 5 X 5
plotRefToTarget(PCA_Leptopoma$shapes$shapes.comp1$min, mesh_Leptopoma)
plotRefToTarget(PCA_Leptopoma$shapes$shapes.comp1$max, mesh_Leptopoma)
plotRefToTarget(PCA_Leptopoma$shapes$shapes.comp2$min, mesh_Leptopoma)
plotRefToTarget(PCA_Leptopoma$shapes$shapes.comp2$max, mesh_Leptopoma)

plotRefToTarget(PCA_Leptopoma$shapes$shapes.comp1$min, PCA_Leptopoma$shapes$shapes.comp1$max, method = "vector", mag = 2)
plotRefToTarget(PCA_Leptopoma$shapes$shapes.comp2$min, PCA_Leptopoma$shapes$shapes.comp2$max, method = "vector", mag = 2)

## ANOVA GM----
geomorph_data_frame <- geomorph.data.frame(leptopoma.gpa, Sexes=classifiers$Sexes, Locations=classifiers$Locations)

fit.size <- procD.lm(coords ~ log(Csize), 
                     data = geomorph_data_frame, print.progress = FALSE)
anova(fit.size)

fit.locations_interact_sexes <- procD.lm(coords ~ Sexes*Locations, 
                                     data = geomorph_data_frame, print.progress = FALSE)
anova(fit.locations_interact_sexes)