ajout script

2-Nutrient_Flux_Bruch.R

+rm(list=ls())
+setwd("C:/users/camille.poulet/Work Folders/Documents/Thse/Script/Nutrients_Fluxes")
+
+###### Attention, perte de masse = 2* perte de gonade du  energie. 
+###### Donc on ne peut pas approximer le poids avant et aprs comme tant la rduction de la gonade, mais comme tant la rduction de la gonade + perte enrgtique du soma ! 
+
+############## Packages ###############
+
+library(tidyverse)
+library(devtools)
+library(ggpmisc)
+
+############# Prparation & Nettoyage du jeu de donnes ########
+
+BDalosesBruch<-read.csv("BDalosesBruch.csv", header= T, sep=";", dec=",") %>%
+  rename(Season = Anne,
+         River = Rivire.pigeage,
+         Sex = Sexe,
+         LT = Lt..cm.,
+         LF = Lf..cm. ,
+         WT = M.tot..g.,
+         WG = M.gonades..g.,
+         Fish_ID = N..poisson) %>% 
+  na.omit(.)  %>% 
+  filter (River != "Dordogne/Garonne")  # On enlve les sites d'chantillonage indtermins 
+
+BDalosesBruch$Sex = str_to_upper(BDalosesBruch$Sex) #Homognistaion des caractres sur la variable "sex": tout est not en MAJ 
+
+#On isole les donnes pr-reproduction
+
+BDalosesBruchPreSpawn<- BDalosesBruch %>% 
+  filter (LOT == "Golfech" | LOT=="Tuilires")  
+
+#On les enlves de la BDD initial pour n'avoir que les donnes pots-reproduction
+BDalosesBruch<- BDalosesBruch %>%   
+  filter (LOT != "Golfech" & LOT != "Tuilires")
+
+  
+######## Calcul des moyennes WT, WG, LT, LF  #########
+
+#============ Donnes Post reproduction ===========
+MeanFeatures <- BDalosesBruch %>%
+  select("Sex","Season","WT","WG","LF","LT") %>% 
+  group_by(Sex) %>%
+  #mutate(LF_Quantile =quantile(.$LF, prob = 0.5)) %>% 
+  summarise_each(.,funs(mean)) 
+
+#============ Donnes Pr-reproduction ============
+MeanFeaturesPreSpawn <- BDalosesBruchPreSpawn %>%
+  select("Sex","Season","WT","WG","LF","LT") %>% 
+  group_by(Sex) %>%
+  summarise_each(.,funs(mean)) 
+
+######## LF #########
+
+tapply(BDalosesBruch$LF, BDalosesBruch[,c("Sex")],min, na.rm = TRUE)
+tapply(BDalosesBruch$LF, BDalosesBruch[,c("Sex")],max, na.rm = TRUE)
+
+
+tapply(BDalosesBruch$LF, BDalosesBruch[,c("Season","Sex")],quantile, na.rm = TRUE, probs= 0.05)
+tapply(BDalosesBruch$LF, BDalosesBruch[,c("Season","Sex")],quantile, na.rm = TRUE, probs= 0.5)
+tapply(BDalosesBruch$LF, BDalosesBruch[,c("Season","Sex")],quantile, na.rm = TRUE, probs= 0.95)
+
+######## WG  #########
+
+#=========== Post reproduction ================
+
+tapply(BDalosesBruch$WG, BDalosesBruch[,c("Season","Sex")],min, na.rm = TRUE)
+tapply(BDalosesBruch$WG, BDalosesBruch[,c("Season","Sex")],max, na.rm = TRUE)
+
+#=========== Pr reproduction ================
+
+tapply(BDalosesBruchPreSpawn$WG,BDalosesBruchPreSpawn [,c("Season","Sex")],min, na.rm = TRUE)
+tapply(BDalosesBruchPreSpawn$WG, BDalosesBruchPreSpawn[,c("Season","Sex")],max, na.rm = TRUE)
+
+########### Relation taille-poids ###############
+
+#========== Vrification des Wpre avec les donnes de Taverny 
+
+#Taverny (p44 livre Baglinire & Elie), 1991
+
+#WT = aLT^b
+
+#Female 
+#a=1.2102e-6
+#b=3.3429
+
+#Male 
+#a= 2.4386e-6
+#b= 3.2252
+
+WeigtSizeRelationshipTaverny = function (a,Lt,b){
+  
+  a * Lt^b 
+}
+
+BDalosesBruchPreSpawn<-BDalosesBruchPreSpawn %>%
+  select("Sex","Season","WT","WG","LT") %>% 
+  group_by(Sex,Season) %>%
+  mutate(LTmm= LT*10) %>% 
+  mutate(WTTaverny= case_when(
+    Sex=="F" ~ WeigtSizeRelationshipTaverny(1.2102e-6,LTmm,3.3429),
+    Sex=="M" ~ WeigtSizeRelationshipTaverny(2.4386e-6,LTmm,3.2252))) %>% 
+  mutate(DiffWeight= abs(WT-WTTaverny))
+
+#======= Graphique relation taille-poids ========
+
+#TabBdMix<-gather(data = BDalosesBruchPreSpawn,key = Weight,value= Value,WTTaverny,WT,factor_key = TRUE)
+
+WT.Taverny<-BDalosesBruchPreSpawn %>%
+  select("WTTaverny")
+ggplot(TabBdMix, aes (LT,Value))+
+  geom_point(stat="identity")+
+  geom_smooth(method="lm")+
+  stat_poly_eq(parse=TRUE, aes(label=paste(stat(eq.label),stat(adj.rr.label), sep="--")), formula = y~x, coef.digits = 5, rr.digits =2)+
+  xlab("Size (cm)")+
+  ylab("Weight (g)")+
+  theme_bw()+
+  facet_grid(Sex~.)+
+  ggtitle("Pre spawning")
+
+
+#Pre reproduction 
+ggplot(BDalosesBruchPreSpawn, aes(log(LTmm),log(WTTaverny)))+
+  geom_point(stat="identity")+
+  geom_smooth(method="lm")+
+  stat_poly_eq(parse=TRUE, aes(label=paste(stat(eq.label),stat(adj.rr.label), sep="--")), formula = y~x, coef.digits = 5, rr.digits =2)+
+  xlab("Size (cm)")+
+  ylab("Weight (g)")+
+  theme_bw()+
+  facet_grid(Sex~.)+
+  ggtitle("Pre spawning")
+
+#Post Reproduction
+ggplot(BDalosesBruch, aes(LF,WT))+
+  geom_point(stat="identity")+
+  geom_smooth(method="lm")+
+  #stat_poly_eq(parse=TRUE, aes(label=..eq.label..), formula = y~x, coef.digits = 5)+
+  xlab("Size (cm)")+
+  ylab("Weight (g)")+
+  theme_bw()+
+  facet_grid(Sex~.)+
+  ggtitle("Post spawning")
+
+
+#========= Vrification de la relation entre LT et LF ========#
+
+# LT = a+bFL (from Gaygusuz et al 2006)
+
+
+ggplot(BDalosesBruchPreSpawn, aes(LF,LT))+
+  geom_point(stat="identity")+
+  geom_smooth(method="lm")+
+  stat_poly_eq(parse=TRUE, aes(label=..eq.label..), formula = y~x, coef.digits = 5)+
+  xlab("LF (cm)")+
+  ylab("LT(cm)")+
+  theme_bw()
+
+
+
+SSRelationship = function (a,b,Lf){
+  
+  a + b*Lf
+}
+
+BDalosesBruch<-BDalosesBruch %>%
+  select("Sex","Season","WT","WG","LF","LT") %>% 
+  group_by(Sex,Season) %>%
+  mutate(LT.SSRelationship= SSRelationship(6.9324,0.9611, LF)) %>% 
+  mutate(DiffWeight=(LT-LT.SSRelationship))
+
+BDalosesBruch<-BDalosesBruch %>%
+  select("Sex","Season","WT","WG","LF","LT") %>% 
+  group_by(Sex,Season) %>%
+  mutate(LT.SSRelationship= case_when(
+    Sex=="F" ~ SSRelationship(exp(0.64858),0.85885, LF),
+    Sex=="M" ~ SSRelationship(exp(0.59854), 0.87097,LF))) %>% 
+  mutate(DiffWeight=(LT-LT.SSRelationship))
+
+
+#Vrification sur les donnes de Bruch 
+
+
+
+
+
+
+
+########## Calcul des GSI ################
+
+# GSI = (WG / (WT-WG))*100
+
+#========== Post reproduction ==============
+
+DataGSI <- BDalosesBruch %>%
+  select("Sex","WT","WG") %>% 
+  group_by("Sex") %>% 
+  #summarise_each(.,funs(mean)) %>% 
+  mutate(GSI = WG/(WT))
+
+######## Calcul des GSI based on Taverny WSR #####
+
+BDalosesBruchPreSpawn %>%
+  select("Sex","WG","DiffWeight") %>% 
+  group_by(Sex) %>%
+  mutate(GSI = WG/(DiffWeight)) %>% 
+  summarise(., avg=mean(GSI))
+
+
+#########" TEST trouver l'erreur #########
+
+hist(DataGSI$GSI[DataGSI$Sex=="F"], breaks=200, col="red")
+abline(v=0.15)
+abline(v=0.12, col="green")
+
+DataGSI <- BDalosesBruch %>%
+  select("Sex","WT","WG", "LT", "LF") %>% 
+  group_by("Sex") %>% 
+  filter(Sex=="F") %>% 
+  #summarise_each(.,funs(mean)) %>% 
+  mutate(GSI = WG/(WT)) %>% 
+  filter(GSI<0.12)
+
+ggplot(DataGSI, aes(log(LF),log(WT)))+
+  geom_point(stat="identity")+
+  geom_smooth(method="lm")+
+  stat_poly_eq(parse=TRUE, aes(label=paste(stat(eq.label),stat(adj.rr.label), sep="--")), formula = y~x, coef.digits = 5, rr.digits =2)+
+  xlab("Size (cm)")+
+  ylab("Weight (g)")+
+  theme_bw()+
+  ggtitle("Post spawning")
+
+
+#========== Pr reproduction ==============
+
+DataGSIPreSpawn <- BDalosesBruchPreSpawn %>%
+  select("Sex","WT","WG") %>% 
+  group_by(Sex) %>%
+  summarise_each(.,funs(mean)) %>% 
+  mutate(GSI = WG/(WT))
+
+#=========== Calcul d'une perte de poids total ==========
+
+WTpre<-tapply(BDalosesBruchPreSpawn$WT, BDalosesBruchPreSpawn[,"Sex"],mean, na.rm = TRUE)
+WTpost<-tapply(BDalosesBruch$WT, BDalosesBruch[,"Sex"],mean, na.rm = TRUE)
+WTloss <- (WTpre - WTpost)
+
+#=========== Calcul d'une perte de poids gonade ==========
+
+#Barber et al, 2018: 
+
+#Wgpr-Wgpost = mass of spawned egg/sperm 
+
+WGpre<-tapply(BDalosesBruchPreSpawn$WG, BDalosesBruchPreSpawn[,"Sex"],mean, na.rm = TRUE)
+WGpost<-tapply(BDalosesBruch$WG, BDalosesBruch[,"Sex"],mean, na.rm = TRUE)
+WGloss <- (WGpre - WGpost) 
+
+#========== Calcul de WT post-reproduction par la relation allomtrie ======== 
+
+# WT = a * Lt ^b /( GSI +1)
+# WT = Wpre/(GSI +1)
+
+#==== Based on Taverny weight-size relationship 
+
+ComputeWeigtPostSpawning = function (a, Lt, b, GSI){
+
+  Wpre <- WeigtSizeRelationshipTaverny (a,Lt,b)
+  
+    return (Wpost = Wpre/(GSI+1))
+  
+}
+
+BDalosesBruch<-BDalosesBruch%>%
+  select("Sex","Season","WT","WG","LF","LT") %>% 
+  group_by(Sex,Season) %>%
+  mutate(LT.mm= LT*10) %>% 
+  mutate(WT.CPL= case_when(
+    Sex=="F" ~ ComputeWeigtPostSpawning(1.2102e-6,LT.mm,3.3429,0.10389710),
+    Sex=="M" ~ ComputeWeigtPostSpawning(2.4386e-6,LT.mm,3.2252,0.05026625))) %>% 
+  mutate(DiffWeight=(WT-WT.CPL))
+
+ggplot(BDalosesBruch, aes(WT.CPL,WT))+
+  geom_point(stat="identity")+
+  xlab("Weight size relationship (Taverny)")+
+  ylab("Weight (g) ")+
+  theme_bw()+
+  facet_grid(Sex~.)
+
+
+BDalosesBruch %>% 
+  filter (LT > 65)
+
+########### Calcul d'une relation d'allomtrie WG et LT #####
+
+#Pre spawning 
+ggplot(BDalosesBruchPreSpawn, aes(log(WT),log(WG)))+
+  geom_point(stat="identity")+
+  geom_smooth(method="lm")+
+  stat_poly_eq(parse=TRUE, aes(label=paste(stat(eq.label),stat(adj.rr.label), sep="--")), formula = y~x, coef.digits = 5, rr.digits =2)+
+  xlab("Size (cm)")+
+  ylab("Weight (g)")+
+  theme_bw()+
+  facet_grid(Sex~.)+
+  ggtitle("Pre spawning")
+
+#Post spawning 
+ggplot(BDalosesBruch, aes(log(WT),log(WG)))+
+  geom_point(stat="identity")+
+  geom_smooth(method="lm")+
+  stat_poly_eq(parse=TRUE, aes(label=paste(stat(eq.label),stat(adj.rr.label), sep="--")), formula = y~x, coef.digits = 5, rr.digits =2)+
+  xlab("Size (cm)")+
+  ylab("Weight (g)")+
+  theme_bw()+
+  facet_grid(Sex~.)+
+  ggtitle("Post spawning")
+
+BDalosesBruch<-BDalosesBruch %>% 
+  filter (log(LF)<4.2)
+
+
+
+#========= Calcul des Wpost et comparaison avec la BDD ======
+
+#totalWeightPost = totalWeightPre * (1-GSIfemalePost)+ totalWeightPost * GSIfemalePost * CoeffLossWeight
+
+GSI <- DataGSI$GSI[,DataGSI$Sex=="F"]
+WTpost = mean (BDalosesBruchPreSpawn$WT)*(1-DataGSI$GSI[,sex=="F"])+ mean(BD)
+
+############ Boxplot exploratoires ##########
+
+par(mfrow=c(1,3))
+boxplot(BDalosesBruch$WG ~ BDalosesBruch$Sex, names=c("F","M"), xlab = "Sex", ylab= "Mean Gonad Mass (g)",main="Gonad Mass")
+boxplot(BDalosesBruch$WT ~ BDalosesBruch$Sex, names=c("F","M"), xlab = "Sex", ylab= "Mean Total Mass (g)", main = "Total Mass")
+boxplot(BDalosesBruch$LF ~ BDalosesBruch$Sex, names=c("F","M"), xlab = "Sex", ylab= "Mean LF (cm)", main = "LF")
+
+
+par(mfrow=c(1,3))
+boxplot(MeanFeatures$WG ~ MeanFeatures$Sex, names=c("F","M"), xlab = "Sex", ylab= "Mean Gonad Mass (g)",main="Gonad Mass")
+boxplot(MeanFeatures$WT ~ MeanFeatures$Sex, names=c("F","M"), xlab = "Sex", ylab= "Mean Total Mass (g)", main = "Total Mass")
+boxplot(MeanFeatures$LF ~ MeanFeatures$Sex, names=c("F","M"), xlab = "Sex", ylab= "Mean LF (cm)", main = "LF")
+
+
+############## Mean Gonad Mass #############
+ggplot(MeanFeatures, aes(Season,(WG)))+
+  geom_bar(stat="identity")+
+  xlab("Season")+
+  ylab("Gonad Mass (g)")+
+  theme_bw()+
+  facet_grid(Sex~.)
+
+############## Mean Total Mass #############
+ggplot(MeanFeatures, aes(Season,(WT)))+
+  geom_bar(stat="identity")+
+  xlab("Season")+
+  ylab("Total Mass (g)")+
+  theme_bw()+
+  facet_grid(Sex~.)
+
+############## Mean Total Lenght #############
+ggplot(MeanFeatures, aes(Season,(LF)))+
+  geom_bar(stat="identity")+
+  xlab("Season")+
+  ylab("Total length (cm)")+
+  theme_bw()+
+  facet_grid(Sex~.)
+
+############## GSI #############
+
+ggplot(BDalosesBruch, aes(Season, (WT-WG)/(WT)))
+geom_bar(stat="identity")+
+  xlab("Season")+
+  ylab("Total length (cm)")+
+  theme_bw()+
+  facet_grid(Sex~.)
+
+
+#TODO 
+
+#Faire un graphique avec les GSI pour regarder si on obtient des donnes similaires  celles trouves dans la littrature 
+#Regarder les relation taille-poids pour dfinir des coefficient a et b 
+#Faire des quantiles 5% sur la taille et le poids pour enlever les donnes abrentes 
+#Comparer les valeurs de poids moyen des gonades avec ceux de la littrature 
+
+################ StringR essay #################
+
+BDalosesBruch$Fish_ID<-as.character(BDalosesBruch$Fish_ID)
+
+## Detect function 
+Pre_Spawn<-str_detect(BDalosesBruch$Fish_ID, regex("Golfech.*$", ignore_case = TRUE))
+Pre_Spawn<-str_detect(BDalosesBruch$Fish_ID, regex("Tuilires.*$", ignore_case = TRUE))
+
+str_count(Pre_Spawn, "TRUE")
+
+## Subset function 
+Pre_Spawn<-str_subset(BDalosesBruch$Fish_ID, "Golfech.*$")
+Pre_Spawn<-str_subset(BDalosesBruch$Fish_ID, regex ("Tuilires.*$", ignore_case = TRUE))
+
+str_count(Pre_Spawn, "TRUE")
+
+