Load required packages

#load packages
library("dplyr")
library("lattice")
library("ggplot2")
library("tidyr")
library("viridis")
library("broom")

1. Load data with info for all samples and the methylation probe data

#load data on shrew
library(readr)
df <- read_csv("SampleSheetAgeN90final.csv")

#remove fetus from both data sets
df_nofet=subset(df, df$Age>0.1)

#add column with location info for island-mainland (Island = BPI and Long Island, mainland = other populations )
df_nofet<- mutate(df_nofet, Location_1 = ifelse(df_nofet$Location == "Bon Portage Island Nova Scotia" | df_nofet$Location == "Long Island Nova Scotia", "yes", "no"))

#remove rows that are doubles of same individual
df_nofet = filter(df_nofet, !(OriginalOrderInBatch %in% c("30", "34", "36", "37", "38", "39", "40", "41", "42", "43", "44")))

2. Linear models for morphology between Island and Mainland groups

#lm for weight in relation age, sex and location
model_weight <- lm(Weight ~ Age + Sex + Location_1, data=df_nofet)

#get model statistics
#summary(model_weight)
#tidy(model_weight)
#glance(model_weight)
library(sjPlot)
tab_model(model_weight, show.ci = 0.95, show.r2 = TRUE, show.re.var=TRUE)
  Weight
Predictors Estimates CI p
(Intercept) 5.07 4.44 – 5.70 <0.001
Age 1.15 0.21 – 2.08 0.018
Sex [M] 0.01 -0.61 – 0.64 0.962
Location 1 [yes] -2.09 -2.79 – -1.39 <0.001
Observations 33
R2 / R2 adjusted 0.571 / 0.526 
#lm for body length in relation age, sex and location
model_length <- lm(Length ~ Age + Sex + Location_1, data=df_nofet)

#get model statistics
#summary(model_length)
tab_model(model_length, show.ci = 0.95, show.r2 = TRUE, show.re.var=TRUE)
  Length
Predictors Estimates CI p
(Intercept) 55.23 51.79 – 58.66 <0.001
Age 8.43 3.36 – 13.50 0.002
Sex [M] -0.20 -3.61 – 3.20 0.903
Location 1 [yes] -10.35 -14.15 – -6.54 <0.001
Observations 32
R2 / R2 adjusted 0.536 / 0.487 
#lm for skull length in relation age, sex and location
model_skull_length <- lm(df_nofet$`Skull length` ~ Age + Sex + Location_1, data=df_nofet)

#get model statistics
#summary(model_skull_length)
tab_model(model_skull_length, show.ci = 0.95, show.r2 = TRUE, show.re.var=TRUE)
## Warning: Using `$` in model formulas can produce unexpected results. Specify your model
##   using the `data` argument instead.
##   Try: `Skull length` ~ Age + Sex +
##   Location_1, data =
  df_nofet$Skull length
Predictors Estimates CI p
(Intercept) 19.16 18.73 – 19.58 <0.001
Age -0.31 -0.93 – 0.32 0.321
Sex [M] -0.21 -0.63 – 0.21 0.319
Location 1 [yes] -1.92 -2.39 – -1.45 <0.001
Observations 31
R2 / R2 adjusted 0.810 / 0.788

3. Boxplots for morphology between Island and Mainland groups

library(ggpubr)

#boxplot for weight
#split boxplots by location and make y axis weight, change color of dots based on age
boxplot_weight <- ggplot(df_nofet, aes(x=Location_1, y=Weight, color=Age)) + 
  #remove outlier extra shape
  geom_boxplot(outlier.shape = NA) +
  #add dots, add jitter to make them spread a bit and make shape change by sex
  geom_point(position=position_jitterdodge(dodge.width=0.5), aes(shape=Sex)) +
  #change color of dots based on age where young is red and blue old
  scale_color_gradient(low = "red", high = "blue") +
  theme_classic()+
  #title, legend and axis labels
  labs(title="Weight",x ="Location", y = "Weight (g)") +
  labs(color='Age (years)') +
  #change y axis scale
  scale_y_continuous(limits=c(2.3,7.75)) +
  #change x axis labels from yes and no to island and mainland
  scale_x_discrete(labels=c("no" = "Mainland", "yes" = "Island")) +
  #put title in middle and include t-test stats on figure
  theme(plot.title = element_text(hjust = 0.5)) + stat_compare_means(method = "t.test", label.x=1.25)
#view plot
boxplot_weight

boxplot_length <- ggplot(df_nofet, aes(x=Location_1, y=Length, color=Age)) + 
  geom_boxplot(outlier.shape = NA) +
  geom_point(position=position_jitterdodge(dodge.width=0.5), aes(shape=Sex)) +
  scale_color_gradient(low = "red", high = "blue") +
  theme_classic()+
  labs(title="Body Length",x ="Location", y = "Length (mm)") +
  labs(color='Age (years)') +
  scale_y_continuous(limits=c(43,70), breaks=seq(45,70,5)) +
  scale_x_discrete(labels=c("no" = "Mainland", "yes" = "Island"))+
  theme(plot.title = element_text(hjust = 0.5)) + stat_compare_means(method = "t.test", label.x=1.25)
boxplot_length

boxplot_skull_length <- ggplot(df_nofet, aes(x=Location_1, y=df_nofet$`Skull length`, color=Age)) + 
  geom_boxplot(outlier.shape = NA) +
  geom_point(position=position_jitterdodge(dodge.width=0.5), aes(shape=Sex)) +
  scale_color_gradient(low = "red", high = "blue") +
  theme_classic()+
  labs(title="Skull Length",x ="Location", y = "Length (mm)") +
  labs(color='Age (years)') +
  scale_y_continuous(limits=c(15.75,19.65)) +
  scale_x_discrete(labels=c("no" = "Mainland", "yes" = "Island"))+
  theme(plot.title = element_text(hjust = 0.5)) + stat_compare_means(method = "t.test", label.x=1.25)
boxplot_skull_length

#take legend form one of the plots and make it an object
legend <- get_legend(boxplot_length) 

#remove legends from all plots to be able to nicely merge them together
boxplot_length_1 <- boxplot_length + theme(legend.position = "none")
boxplot_weight_1 <- boxplot_weight + theme(legend.position = "none")
boxplot_skull_length_1 <- boxplot_skull_length + theme(legend.position = "none")

library(cowplot)
#merge all boxplots to make one figure
#add labels (A, B, C) for each panel, include legend object, use ncol = 4 to put all plots side by side and select widths for each panel (make legend smaller to the side)
final_plot <- plot_grid(boxplot_weight_1, boxplot_length_1, boxplot_skull_length_1, labels = c('A', 'B', 'C'), legend, ncol = 4,  rel_widths = c(0.75, 0.75, 0.75, 0.3))
#view final plot
final_plot

#save manually as pdf and make sure everything fits properly and nothing is cut off in final plot