#Load packages library(readr) library(lmerTest) library(lme4) library(lmerTest) library(emmeans) library(ggplot2) library(gridExtra) library(grid) library(dplyr) #Set factors in the dataset Linear_mixed_effect_model_input_data$Strain<-factor(Linear_mixed_effect_model_input_data$Strain) Linear_mixed_effect_model_input_data$Family<-factor(Linear_mixed_effect_model_input_data$Family) Linear_mixed_effect_model_input_data$Dam<-factor(Linear_mixed_effect_model_input_data$Dam) Linear_mixed_effect_model_input_data$Sire<-factor(Linear_mixed_effect_model_input_data$Sire) Linear_mixed_effect_model_input_data$SW<-factor(Linear_mixed_effect_model_input_data$SW) Linear_mixed_effect_model_input_data$Background<-factor(Linear_mixed_effect_model_input_data$Background) #Remove NAs from data Linear_mixed_effect_model_input_data_no_NA<-na.omit(Linear_mixed_effect_model_input_data) #Set factors in the dataset Linear_mixed_effect_model_input_data_no_NA$Strain<-factor(Linear_mixed_effect_model_input_data_no_NA$Strain) Linear_mixed_effect_model_input_data_no_NA$Family<-factor(Linear_mixed_effect_model_input_data_no_NA$Family) Linear_mixed_effect_model_input_data_no_NA$Dam<-factor(Linear_mixed_effect_model_input_data_no_NA$Dam) Linear_mixed_effect_model_input_data_no_NA$Sire<-factor(Linear_mixed_effect_model_input_data_no_NA$Sire) Linear_mixed_effect_model_input_data_no_NA$SW<-factor(Linear_mixed_effect_model_input_data_no_NA$SW) Linear_mixed_effect_model_input_data_no_NA$Background<-factor(Linear_mixed_effect_model_input_data_no_NA$Background) #Regression between fork length and kype height kype_height_resid_no_log<- lm(kype_height_cm ~ Body_length_mature_cm, data=Linear_mixed_effect_model_input_data) summary(kype_height_resid_no_log) #Regression between fork length and kype length kype_length_resid_no_log<- lm(kype_length_5_4 ~ Body_length_mature_cm, data=Linear_mixed_effect_model_input_data) summary(kype_length_resid_no_log) #Regression between log10 transformed fork length and log10 transformed kype height kype_height_resid<- lm(log_kype_height ~ log_body_length_mature_cm, data=Linear_mixed_effect_model_input_data) summary(kype_height_resid) kype_height_resid<-resid(kype_height_resid) head(kype_height_resid) #Regression between log10 transformed fork length and log10 transformed kype length kype_length_resid<- lm(log_5_4 ~ log_body_length_mature_cm, data=Linear_mixed_effect_model_input_data) summary(kype_length_resid) kype_length_resid<-resid(kype_length_resid) head(kype_length_resid) #Full AKL LME Full_AKL<-lmer(AKL~ Strain*SW*GSI+ (1|Strain:Family)+ (1|Dam) +(1|Sire), data = Linear_mixed_effect_model_input_data_no_NA) #Reducing model using step function s1_AKL<-step(Full_AKL) s1_AKL #Reduced model after step function Final_model_AKL <- lmer(formula = AKL ~ Strain + SW + (1 |Sire), data = Linear_mixed_effect_model_input_data) #Effect of factors within LME anova(Final_model_AKL) summary(Final_model_AKL) #Caluclating means from LME emmeans (Final_model_AKL, pairwise ~Strain, adjust="tukey") emmeans (Final_model_AKL, pairwise ~SW, adjust="tukey") #Full AKH LME Full_AKH_model<-lmer(AKH~ SW*Strain*GSR+ (1|Strain:Family)+ (1|Dam) +(1|Sire), data = Linear_mixed_effect_model_input_data_no_NA) #Reducing model using step function s1_AKH<-step(Full_AKH_model) s1_AKH #Reduced model after step function Final_model_AKH <- lmer(formula = AKH ~ SW + Strain + (1 | Sire), data = Linear_mixed_effect_model_input_data) #Effect of factors within LME anova(Final_model_AKH) summary(Final_model_AKH) #Caluclating means from LME emmeans (Final_model_AKH, pairwise ~Strain, adjust="tukey") emmeans (Final_model_AKH, pairwise ~SW, adjust="tukey") #Jitter boxplots showing AKH/AKL means for SW/strains emmeans_AKL_SW_ms$SW<-factor(emmeans_AKL_SW_ms$SW) ggplot(data=emmeans_AKL_SW_ms,aes(x=SW,y=emmean,color=SW))+ geom_point()+ geom_jitter(data = Linear_mixed_effect_model_input_data, aes(x=SW,y=AKL),size=1)+ geom_pointrange(data=emmeans_AKL_SW_ms, mapping=aes(x=SW, y=emmean, ymin=upper.CL, ymax=lower.CL),shape=22,color='red',fill='white')+ theme_bw(base_size=20)+ scale_color_manual(values=c("#CC79A7","#009E73","#000000"))+theme(legend.position="none")+ coord_cartesian(ylim = c(-0.4,0.5)) emmeans_AKL_strains_ms$Strain<-factor(emmeans_AKL_strains_ms$Strain,levels= c('Arna','Vosso','Figgjo','Hybrid_FM','Hybrid_MF','Farmed'),ordered=TRUE) ggplot(data=emmeans_AKL_strains_ms,aes(x=Strain,y=emmean,color=Strain))+ geom_point()+ geom_jitter(data = Linear_mixed_effect_model_input_data, aes(x=Strain,y=AKL),size=1)+ geom_pointrange(data=emmeans_AKL_strains_ms, mapping=aes(x=Strain, y=emmean, ymin=upper.CL, ymax=lower.CL),shape=22,color='red',fill='white')+ theme_bw(base_size=20)+ scale_color_manual(values=c("#999999","#0072B2","#56B4E9","#F0E442","#E69F00","#D55E00"))+theme(legend.position="none")+ coord_cartesian(ylim = c(-0.4,0.5)) emmeans_AKH_SW_ms$SW<-factor(emmeans_AKH_SW_ms$SW) ggplot(data=emmeans_AKH_SW_ms,aes(x=SW,y=emmean,color=SW))+ geom_point()+ geom_jitter(data = Linear_mixed_effect_model_input_data, aes(x=SW,y=AKH),size=1)+ geom_pointrange(data=emmeans_AKH_SW_ms, mapping=aes(x=SW, y=emmean, ymin=upper.CL, ymax=lower.CL),shape=22,color='red',fill='white')+ theme_bw(base_size=20)+ scale_color_manual(values=c("#CC79A7","#009E73","#000000"))+theme(legend.position="none")+ coord_cartesian(ylim = c(-0.4,0.5)) emmeans_AKH_strains_ms$Strain<-factor(emmeans_AKH_strains_ms$Strain,levels= c('Arna','Vosso','Figgjo','Hybrid_FM','Hybrid_MF','Farmed'),ordered=TRUE) ggplot(data=emmeans_AKH_strains_ms,aes(x=Strain,y=emmean,color=Strain))+ geom_point()+ geom_jitter(data = Linear_mixed_effect_model_input_data, aes(x=Strain,y=AKH),size=1)+ geom_pointrange(data=emmeans_AKH_strains_ms, mapping=aes(x=Strain, y=emmean, ymin=upper.CL, ymax=lower.CL),shape=22,color='red',fill='white')+ theme_bw(base_size=20)+ scale_color_manual(values=c("#999999","#0072B2","#56B4E9","#F0E442","#E69F00","#D55E00"))+theme(legend.position="none")+ coord_cartesian(ylim = c(-0.4,0.5)) library("geomorph") # Geometric morphometrics 6 landmarks(lm) landmarks6lm <- readland.tps("Mature_female_and_males_6lm.TPS",specID=c("ID")) corrected6lm<-gpagen(landmarks6lm, PrinAxes=TRUE) plotAllSpecimens(corrected6lm$coords, mean=TRUE, label=TRUE) geometric_PCA6lm<-plotTangentSpace(corrected6lm$coords, axis1=1, axis2=2, label = TRUE) PCA1_6lm<-plotTangentSpace(landmarks6lm, axis1=1, axis2=2, label = TRUE) write.csv(geometric_PCA6lm$pc.summary$importance, "PCAxes_geometric_all_6lm.csv") write.csv(geometric_PCA6lm$pc.scores, "PCScores_geometric_all_6lm.csv") # extract values needed from the PCA tables PCScores_geometric <- read_csv("PCScores_geometric_all_6lm.csv") PCAxes_geometric <- read_csv("PCAxes_geometric_all_6lm.csv") PCScores_geometric_sub <- PCScores_geometric[,1:3] PCAxes_geometric_sub <- PCAxes_geometric[2,2:3] write.csv(PCScores_geometric_sub, "PCScores_geometric6lm.csv") write.csv(PCAxes_geometric_sub, "PCAxes_geometric6lm.csv") #Add SW column manually. ID 0-444 = SW1, ID 0-97 = SW2, ID 1-62 = SW3, ID 0-95 = Female PCScores_geometric6lm <- read_csv("PCScores_geometric6lm.csv") PCScores_geometric6lm$SW<-factor(PCScores_geometric6lm$SW) #PCA plot for 6lm and including females with frequency plots on x and y axis pdf("Males_female_6lm_PCA_plot.pdf") PCA_geometric6lm<-ggplot(PCScores_geometric6lm,aes(PC1,PC2, color=SW))+ theme_bw()+ geom_point()+ scale_color_manual(values = c("#CC79A7","#009E73","#000000","red")) + coord_fixed()+ stat_ellipse(type = "t",size=1)+ theme(legend.position = "none",axis.title.x = element_blank(),axis.title.y = element_blank()) xdensity <- ggplot(PCScores_geometric6lm, aes(PC1, fill=SW)) + theme_bw()+ geom_density(alpha=.5) + theme(axis.title.y = element_blank(),axis.text.y = element_blank(),axis.ticks = element_blank(),axis.title.x = element_blank())+ scale_fill_manual(values = c("#CC79A7","#009E73","#000000","red")) + theme(legend.position = "none") ydensity <- ggplot(PCScores_geometric6lm, aes(PC2, fill=SW)) + theme_bw()+ geom_density(alpha=.5) + theme(axis.text.x = element_blank(),axis.ticks = element_blank(),axis.title.x = element_blank())+ scale_fill_manual(values = c("#CC79A7","#009E73","#000000","red")) + theme(legend.position = "none")+ coord_flip() blankPlot <- ggplot()+geom_blank(aes(1,1))+ theme(plot.background = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank()) grid.arrange(xdensity, blankPlot, PCA_geometric6lm, ydensity, ncol=2, nrow=2, widths=c(4, 1.4), heights=c(1.8, 4)) dev.off()