diff --git a/RAnalysis/Output/DryWeights/F1/DryWeights_bFactorLengthpCO2.pdf b/RAnalysis/Output/DryWeights/F1/DryWeights_bFactorLengthpCO2.pdf index 687f639..52e7b00 100644 Binary files a/RAnalysis/Output/DryWeights/F1/DryWeights_bFactorLengthpCO2.pdf and b/RAnalysis/Output/DryWeights/F1/DryWeights_bFactorLengthpCO2.pdf differ diff --git a/RAnalysis/Scripts/F1_Phys/F1_DryWeights_analysis.Rmd b/RAnalysis/Scripts/F1_Phys/F1_DryWeights_analysis.Rmd index 9d3fd9e..1a7f03f 100644 --- a/RAnalysis/Scripts/F1_Phys/F1_DryWeights_analysis.Rmd +++ b/RAnalysis/Scripts/F1_Phys/F1_DryWeights_analysis.Rmd @@ -268,7 +268,7 @@ dev.off() * here we do no require a bfactor becuase we are already scalling to the individual! ### Build poroprtion dataset -```{r Poroportion calculate} +```{r Proportion calculate} Proportion_master <- F1_dryweights_calc.long %>% @@ -2109,4 +2109,12 @@ print(ggarrange(Gonad_final_plot, Adductor_final_plot, dev.off() +pdf("Output/DryWeights/F1/DryWeights_bFactorLengthpCO2.pdf", height=7.5, width =10) +print(ggarrange(Gonad_final_plot, Adductor_final_plot, + Somatic_final_plot, Tissue_final_plot_SUBSET, + Shell_final_plot_SUBSET, + ncol = 3, nrow = 2 + )) +dev.off() + ``` diff --git a/RAnalysis/Scripts/F1_Phys/F1_Lengths_analysis.Rmd b/RAnalysis/Scripts/F1_Phys/F1_Lengths_analysis.Rmd index 68cd0b9..3c23c39 100644 --- a/RAnalysis/Scripts/F1_Phys/F1_Lengths_analysis.Rmd +++ b/RAnalysis/Scripts/F1_Phys/F1_Lengths_analysis.Rmd @@ -21,13 +21,14 @@ library(dplyr) library(ggplot2) library(forcats) library(lme4) -library(lmerTest) -library(performance) +# library(lmerTest) +# library(performance) +# install.packages('') library(car) library(ggpubr) -library(SciViews) +# library(SciViews) library(Rmisc) -library(hrbrthemes) +# library(hrbrthemes)lmerTest # SET WORKING DIRECTORY ::::::::::::::::::::::::::::::::::::::::::::::: knitr::opts_knit$set(root.dir = "C:/Users/samjg/Documents/Github_repositories/Airradians_multigen_OA/RAnalysis") @@ -226,8 +227,47 @@ for (i in 1:nrow(ANOVA_Age_DPF)) { print(AOVdf_total) # print to monitor progress } View(AOVdf_total) # view all the anova tests within data -View(AOVdf_total %>% dplyr::filter(P_val < 0.05)) +View(AOVdf_total %>% dplyr::filter(P_val < 0.05)) # 243, 275, 318 are significnatly different via one0-way anova write.csv(AOVdf_total, "C:/Users/samjg/Documents/Github_repositories/Airradians_multigen_OA/RAnalysis/Output/Length/F1/F1_Length_ANOVA_table.csv") + +# (2) You can easitly do #1 by piping in the formula here BUT this does not have contingency based on assumptions +F1_length_MEANS_growout %>% + group_by(Age_DPF) %>% + t_test(data =., Length_mm ~ pCO2) %>% #243, 275, 318 are significnatly different via one0-way anova + adjust_pvalue(method = "bonferroni") %>% + add_significance("p.adj") # with bonferroni JUST 273 is significantly different! + + +# (3) GLMM +F1_length_MEANS_growout$ID <- paste0(F1_length_MEANS_growout$Replicate, + '_', + F1_length_MEANS_growout$pH) # ID variable to run the GLMM withput averaging! + +F1_GLMM <- glmer( + Length_mm ~ + pCO2 * Age_DPF + (1 | ID), + data = F1_length_MEANS_growout, + # gaussian(link = "identity") + Gamma(link="identity") + # Gamma(link="log") + ) + +summary(F1_GLMM) + +# check the distribution of residuals +simulationOutput <- simulateResiduals(fittedModel = F1_GLMM, use.u = T) +residuals <- residuals(F1_GLMM, type = "response", retype="normalized") + +# pdf(file="Output/1_Survival_Growth/size_growth/Experiment4/Experiment4_Length_GLMM_QQ_Resid.pdf", height = 4, width = 7) +plot(simulationOutput) +# dev.off() + +# pdf(file="Output/1_Survival_Growth/size_growth/Experiment4/Experiment4_Length_GLMM_Resid.pdf", height = 5, width = 5) +plot(residuals) +# dev.off() + +# output the anova table +car::Anova(F1_GLMM, type=3) # JUST AN EFFECT OF TIME ``` \ No newline at end of file diff --git a/RAnalysis/Scripts/F2_Phys/F2_Lengths_analysis.Rmd b/RAnalysis/Scripts/F2_Phys/F2_Lengths_analysis.Rmd index 4a8b953..5684276 100644 --- a/RAnalysis/Scripts/F2_Phys/F2_Lengths_analysis.Rmd +++ b/RAnalysis/Scripts/F2_Phys/F2_Lengths_analysis.Rmd @@ -222,4 +222,46 @@ View(AOVdf_total %>% dplyr::filter(P_val < 0.05)) write.csv(AOVdf_total, "C:/Users/samjg/Documents/Github_repositories/Airradians_multigen_OA/RAnalysis/Output/Length/F2/F2_Length_ANOVA_table.csv") + + +# (2) You can easitly do #1 by piping in the formula here BUT this does not have contingency based on assumptions +F2_length_MEANS_bytank %>% + group_by(Age_DPF) %>% + t_test(data =., Length_mm ~ pCO2) %>% #243, 275, 318 are significnatly different via one0-way anova + adjust_pvalue(method = "bonferroni") %>% + add_significance("p.adj") # with bonferroni JUST 273 is significantly different! +# NO SIGNIFICANT DIFFERENCES WHETHER WITH CONSERVATIVE ADJUSTMENT OR NOT! + +# (3) GLMM +F2_length_MEANS_bytank$ID <- paste0(F2_length_MEANS_bytank$Replicate, + '_', + F2_length_MEANS_bytank$pH) # ID variable to run the GLMM withput averaging! + +F2_GLMM <- glmer( + Length_mm ~ + pCO2 * Age_DPF + (1 | ID), + data = F2_length_MEANS_bytank, + # gaussian(link = "identity") + Gamma(link="identity") + # Gamma(link="log") + ) + +summary(F2_GLMM) + +# check the distribution of residuals +simulationOutput <- simulateResiduals(fittedModel = F2_GLMM, use.u = T) +residuals <- residuals(F2_GLMM, type = "response", retype="normalized") + +# pdf(file="Output/1_Survival_Growth/size_growth/Experiment4/Experiment4_Length_GLMM_QQ_Resid.pdf", height = 4, width = 7) +plot(simulationOutput) +# dev.off() + +# pdf(file="Output/1_Survival_Growth/size_growth/Experiment4/Experiment4_Length_GLMM_Resid.pdf", height = 5, width = 5) +plot(residuals) +# dev.off() + +# output the anova table +car::Anova(F2_GLMM, type=3) # JUST AN EFFECT OF TIME + + ``` \ No newline at end of file diff --git a/RAnalysis/Scripts/Popgen/F1Juveniles_F3Juvneiles_LOWvMOD_workbook.Rmd b/RAnalysis/Scripts/Popgen/F1Juveniles_F3Juvneiles_LOWvMOD_workbook.Rmd index 773518f..f6d902c 100644 --- a/RAnalysis/Scripts/Popgen/F1Juveniles_F3Juvneiles_LOWvMOD_workbook.Rmd +++ b/RAnalysis/Scripts/Popgen/F1Juveniles_F3Juvneiles_LOWvMOD_workbook.Rmd @@ -23,7 +23,9 @@ library(pcadapt) library(adegenet) library(qqman) library(ggpubr) -install.packages("AFLPsim") +library(radiator) # https://thierrygosselin.github.io/radiator/articles/get_started.html +# run bayscan in R using radiator https://rdrr.io/github/thierrygosselin/radiator/man/run_bayescan.html +?run_bayescan # library(pegas) # library(LDlinkR) # library(ldsep) # need a Seqinfo object Airradians not registered in package @@ -37,6 +39,8 @@ install.packages("AFLPsim") * check out thise use of bayscan for neutral blaancing and outlier loci https://github.com/laurabenestan/Bayescan * tutorial to contrast the outliers detected from bayscan and pcadapt https://rpubs.com/lbenestan/outlier +* Narum and Hess 2011 " Comparisons included simulation procedures (FDIST2, ARLEQUIN v.3.5 and BAYESCAN) as well as more conventional tools such as global FST histograms. Of the three simulation methods, FDIST2 and BAYESCAN typically had the lowest type II error, BAYESCAN had the least type I error and Arlequin had highest type I and II error. High error rates in Arlequin with a hierarchical approach were partially because of confounding scenarios where patterns of adaptive variation were contrary to neutral structure; however, Arlequin consistently had highest type I and type II error in all four simulation scenarios tested in this study. Given the results provided here, it is important that outlier loci are interpreted cautiously and error rates of various methods are taken into consideration in studies of adaptive molecular variation, especially when hierarchical structure is included." + ## load files * vcf.gz files