# growth rate analysis for Duyi's exp - loss of biomass AND biodiversity library(nlsMicrobio) library(data.table) library(plyr) library(lattice) library(ggplot2) library(devtools) library(cowplot) require(nlme) require(minpack.lm) library(reshape2) library(devtools) library(MuMIn) library(emmeans) library(AICcmodavg) library("papeR") library(janitor) #need to run the NLS Loop from the previous script! #set wd to ("/Users/elisaschaum/Dropbox/Duyi data") setwd("~/Dropbox/Duyi data and dilution paper") #read file Biod<-read.csv("final_round_grand.csv") #inspect head(Biod) str(Biod) #make cells.mL and day numeric Biod$cells.ml<-as.numeric(Biod$cells.ml) Biod$day<-as.numeric(Biod$day) #make an identifier containing the dilution, temp, and region Biod <- within(Biod, id.01 <- as.character(factor(Region):factor(Temp):factor(dilution_real):factor(biorep))) # if not needed for biorep, use this one Biod <- within(Biod, id.01 <- as.character(factor(Region):factor(Temp):factor(dilution_real))) #we need to make a column with the counts as log values Biod$logcounts<-log10(Biod$cells.ml) #we also get rid of all NAs Biod<-Biod[complete.cases(Biod), ] #check if it did the one job it had head(Biod) summary(Biod) str(Biod) #now, exclude the FULL treatment from the analysis, because it did not follow a real growth curve but remained largely stable. can also run a quick plot to double check. qplot(day, logcounts, data=Biod, facets=Region~Temp, colour=dilution_real)+theme_classic() Biod<-subset(Biod, dilution_real!=1) # we only keep the column for day, log counts, and the identifier. You can look at effects on size etc in a different script if you want to! Biod_sub<-Biod[,c(11,17,18)] colnames(Biod_sub) <-c('t','id.01','LOG10N') # rename columns so they fit the model parameters as R likes them id.01 <- unique(Biod_sub$id.01) # find all region/ident combinations that are unique (so we pool by bioreplicate) head(Biod_sub) #this is the function for the most simple logistic model - we will probably use a slightly more sophisticated one gm <- function(r, K, t, T0.biom){ gm <- K /(1 + ((K - T0.biom)/min(T0.biom))*exp(-r*t)) return(gm) } #get time at which there is mumax - note that you may need to add lag phase to that later on. atmu<-function (r, K, T0.biom = 2) { return(t<-((T0.biom/(1-T0.biom/K)))/r) } #if there is no lag or carrying capacity, don't use the gompertz model. you can check ?nlsMicrobio for recommendations ### find parameters that fit the curves best - maybe try one that does not have a lag phase res_logis_lag <- nlsLoop::nlsLoop(Biod_sub, #this will throw hissy fits until it finds a good fit! Error messages are normal as not all of the curves it tries out are going to fit model = baranyi_without_lag, tries = 1000, id_col = 'id.01', param_bds = c(4,8,1,2,0.3,1.8), #parameter space na.action = na.omit, lower = c(LOG10Nmax=4,LOG10N0=1,mumax=0.3), upper = c(LOG10Nmax=8,LOG10N0=2,mumax=1.8)) head(res_logis_lag$params) growth_params_lag<-res_logis_lag$params growth_params_lag$atmum<-atmu(r=growth_params_lag$mumax, K=growth_params_lag$LOG10Nmax,T0.biom=3) head(res_logis_lag$predictions) #these are the predictions that we will use to fit the curve growth_preds_lag<-res_logis_lag$predictions #now we want to get our temp and region and isolate back into dataframe growth_params_lag$Region<- as.vector(t(sapply(as.character(growth_params_lag$id.01), function(y) strsplit(y,split=":")[[1]][1]))) growth_params_lag$Temp<- as.vector(t(sapply(as.character(growth_params_lag$id.01), function(y) strsplit(y,split=":")[[1]][2]))) growth_params_lag$Dilution<- as.vector(t(sapply(as.character(growth_params_lag$id.01), function(y) strsplit(y,split=":")[[1]][3]))) growth_params_lag$Biorep<- as.vector(t(sapply(as.character(growth_params_lag$id.01), function(y) strsplit(y,split=":")[[1]][4]))) #check that this worked head(growth_params_lag) #YAY #now we also get the details back for the predictions head(growth_preds_lag) growth_preds_lag$Region<- as.vector(t(sapply(as.character(growth_preds_lag$id.01), function(y) strsplit(y,split=":")[[1]][1]))) growth_preds_lag$Temp<- as.vector(t(sapply(as.character(growth_preds_lag$id.01), function(y) strsplit(y,split=":")[[1]][2]))) growth_preds_lag$Dilution<- as.vector(t(sapply(as.character(growth_preds_lag$id.01), function(y) strsplit(y,split=":")[[1]][3]))) growth_preds_lag$Biorep<- as.vector(t(sapply(as.character(growth_preds_lag$id.01), function(y) strsplit(y,split=":")[[1]][4]))) #check that this worked head(growth_preds_lag) # YAY!! colnames(growth_preds_lag)[3] <- 'pred_rate' #we rename "Born" back to Bornholm growth_params_lag$Region[growth_params_lag$Region == "Born"] <- "Bornholm" growth_preds_lag$Region[growth_preds_lag$Region == "Born"] <- "Bornholm" Biod$Region<-as.character(Biod$Region) Biod$Region[Biod$Region == "Born"] <- "Bornholm" #bring into correct order growth_params_lag$Region <- factor(growth_params_lag$Region, levels = c('Kiel', 'Arkona', 'Bornholm')) growth_preds_lag$Region <- factor(growth_preds_lag$Region, levels = c('Kiel', 'Arkona', 'Bornholm')) Biod$Region<-factor(Biod$Region, levels = c('Kiel', 'Arkona', 'Bornholm')) #make plots for inspection, then stats! head(growth_preds_lag) #you might want to save the preds and params data frames. As well as the model output so you don't have to start from scratch each time. #make plots with predicted fits - is for SI only subReg<-subset(growth_params_lag, Region!='Arkona') subRegpred<-subset(growth_preds_lag, Region!='Arkona') subBiod<-subset(Biod, Region!='Arkona') subReg$Region <- factor(subReg$Region, levels = c('Kiel', 'Bornholm')) subRegpred$Region <- factor(subRegpred$Region, levels = c('Kiel', 'Bornholm')) subBiod$Region <- factor(subBiod$Region, levels = c('Kiel', 'Bornholm')) colnames(subRegpred)<-c("id.01", "t", "pred_rate", "Region", "Temp" ,"dilution_real","Biorep") #this plot has been arranged so as to work when id.01 does NOT include the bioreps, as it is working on the averages! #which we still have to make subBiod_av<-ddply(subBiod,c("Region","Temp","dilution_real","day"), function(df) return(c(logcounts.avg=mean(df$logcounts),logcounts.se=sd(df$logcounts)/sqrt(length(subBiod))))) #need to include day 1 though subBiod_av<-read.csv("growth_out_Duyi.csv") subBiod_av$Region <- factor(subBiod_av$Region, levels = c('Kiel', 'Bornholm')) qplot(as.numeric(day), logcounts.avg, data=subBiod_av, colour=Region,alpha=as.factor((dilution_real)))+ theme_classic(base_size = 14, base_family = "Helvetica")+ facet_wrap(Region~Temp)+ geom_errorbar(aes(ymin=logcounts.avg-logcounts.se, ymax=logcounts.avg+logcounts.se))+ scale_colour_manual("Region",values=c("orange","cornflowerblue"))+ geom_line(data = subRegpred, aes(x = t, y = pred_rate,colour = Region, alpha=as.factor(dilution_real)))+ theme(legend.position='top') # YES MUCH BETTER qplot(as.numeric(day), logcounts.avg, data=subBiod_av, colour=Region,alpha=as.factor((dilution_real)))+ theme_classic(base_size = 14, base_family = "Helvetica")+ facet_wrap(Region~Temp)+ geom_errorbar(aes(ymin=logcounts.avg-logcounts.se, ymax=logcounts.avg+logcounts.se))+ scale_colour_manual("Region",values=c("orange","cornflowerblue"))+ scale_fill_manual("Region",values=c("orange","cornflowerblue"))+ geom_smooth (aes(group=1, colour=Region, fill=Region),method="gam",formula= y ~ s(x, k = 3))+ theme(legend.position='top') # YES MUCH BETTER #plot r(same as µmax) - colouring etc is up to you...by temp/region/dilution? r_plot01<-ggplot(data=growth_params_lag,aes(Dilution, mumax))+ geom_boxplot(aes(fill=as.factor(Temp)))+ geom_jitter(aes(shape=Biorep, colour=Temp))+ scale_colour_manual("Region",values=c("red","blue","darkgreen"))+ theme_classic(base_size = 14, base_family = "Helvetica")+ facet_wrap(Temp~Region)+ theme(legend.position='top') + ylab("cell doubling (day^-1)") r_plot01 #need to make a new identifier growth_params_lag<-within(growth_params_lag, id.plot<-as.character(factor(Region):factor(Temp))) r_plot01.1<-ggplot(data=growth_params_lag,aes(log(as.numeric(Dilution)), (mumax)))+ geom_jitter(aes(colour=Region))+ scale_colour_manual("Region",values=c("cornflowerblue","purple","lightgreen"))+ scale_fill_manual("Region",values=c("cornflowerblue","purple","lightgreen"))+ theme_classic(base_size = 14, base_family = "Helvetica")+ facet_wrap(Region~Temp)+ #geom_smooth ( aes(group=id.plot, colour=Region, fill=Region),method="gam",formula= y ~ s(x, k = 3))+ theme(legend.position='top')+ geom_smooth(aes(group=id.plot, colour=Region, fill=Region), method="lm")+ theme(legend.position='top') + ylab("Growth rate µ (day-1)")+ xlab("Dilution (ln)") r_plot01.1 r_plot01.1K<-ggplot(data=growth_params_lag,aes(log(as.numeric(Dilution)), (LOG10Nmax)))+ geom_jitter(aes(colour=Region))+ scale_colour_manual("Region",values=c("cornflowerblue","purple","lightgreen"))+ scale_fill_manual("Region",values=c("cornflowerblue","purple","lightgreen"))+ theme_classic(base_size = 14, base_family = "Helvetica")+ facet_wrap(Region~Temp)+ #geom_smooth ( aes(group=id.plot, colour=Region, fill=Region),method="gam",formula= y ~ s(x, k = 3))+ theme(legend.position='top')+ geom_smooth(aes(group=id.plot, colour=Region, fill=Region), method="lm")+ theme(legend.position='top') + ylab("carrying capacity")+ xlab("Dilution (ln)") r_plot01.1K r_plot01.1.1<-ggplot(data=growth_params_lag,aes(((Dilution)), (mumax)))+ geom_point(aes(colour=as.factor(Region)))+ geom_jitter(aes(colour=Region))+ scale_colour_manual("Region",values=c("cornflowerblue","purple","lightgreen"))+ scale_fill_manual("Region",values=c("cornflowerblue","purple","lightgreen"))+ theme_classic(base_size = 14, base_family = "Helvetica")+ facet_wrap(Region~Temp, scales="free_y")+ geom_smooth(aes(group=id.plot, colour=Region, fill=Region), method="lm")+ theme(legend.position='top') + ylab("Growth rate µ")+ xlab(" Dilution") r_plot01.1.1 r_plot01.2<-ggplot(data=growth_params_lag,aes(log(as.numeric(Dilution)), (LOG10Nmax)))+ geom_point(aes(colour=as.factor(Region)))+ geom_jitter(aes(colour=Region))+ scale_colour_manual("Region",values=c("cornflowerblue","purple","lightgreen"))+ scale_fill_manual("Region",values=c("cornflowerblue","purple","lightgreen"))+ theme_classic(base_size = 14, base_family = "Helvetica")+ facet_wrap(Region~Temp)+ geom_smooth(aes(group=id.plot, colour=Region, fill=Region), method="lm")+ theme(legend.position='top') + ylab("log carrying capacity")+ xlab("ln Dilution") r_plot01.2 r_plot01.2.1<-ggplot(data=growth_params_lag,aes(((Dilution)), (LOG10Nmax)))+ geom_point(aes(colour=as.factor(Region)))+ geom_jitter(aes(colour=Region))+ scale_colour_manual("Region",values=c("cornflowerblue","purple","lightgreen"))+ scale_fill_manual("Region",values=c("cornflowerblue","purple","lightgreen"))+ theme_classic(base_size = 14, base_family = "Helvetica")+ facet_wrap(Region~Temp, scales="free_y")+ geom_smooth ( aes(group=id.plot, colour=Region, fill=Region),method="gam",formula= y ~ s(x, k = 3))+ theme(legend.position='top') + ylab("log carrying capacity")+ xlab(" Dilution") r_plot01.2.1 #also want this with biomass - will have to take a separate script that retains the biomass column? Or merge with size data somehow using the id.01 column #first make a size column in Biod Biod$sizeum<-Biod$size/10000+1.3 Biod$volume<- 1.333333*3.14159265359*((Biod$sizeum/2)^3) #assuming spheres #now merge - need to keep id.01 in Biod and volume column Biod_for_merge<-Biod[,c(17,19,20)] Merged<-merge(Biod_for_merge,growth_params_lag) head(Merged) Merged$biomass<-exp(Merged$LOG10Nmax)*Merged$volume #should convert into carbon # As a rule of thumb, we assume a sphere, and use y =ax^b. With, x your cell volume in um^3, a in this case = 0.165 and b=1.046. Not the most accurate science, I am afraid. This is roughly after Montagnes, D.J.S., J.A. Berges, P. J. Harrison, and F. J. R. Taylor. 1994, Table 3 Merged$carbon_bio<-0.165*(Merged$sizeum^1.046) Merged$carbon_bioALL<-Merged$carbon_bio*exp(Merged$LOG10Nmax) r_plot01.2.2<-ggplot(data=Merged,aes(log(as.numeric(Dilution)),log(carbon_bioALL)))+ geom_point(aes(colour=as.factor(Region)))+ scale_colour_manual("Region",values=c("cornflowerblue","purple","lightgreen"))+ scale_fill_manual("Region",values=c("cornflowerblue","purple","lightgreen"))+ theme_classic(base_size = 14, base_family = "Helvetica")+ facet_wrap(Region~Temp, scales="free_y")+ geom_smooth ( aes(group=id.plot, colour=Region, fill=Region),method="gam",formula= y ~ s(x, k = 3))+ theme(legend.position='top') + ylab(" ln biomass in at carrying capacity[pg C]")+ xlab(" ln Dilution") r_plot01.2.2 r_plot01.1.3<-ggplot(data=Merged,aes(log(as.numeric(Dilution)), (mumax)))+ geom_point(aes(colour=as.factor(Region)),position=position_dodge(width=0.5))+ scale_colour_manual("Region",values=c("cornflowerblue","purple","lightgreen"))+ scale_fill_manual("Region",values=c("cornflowerblue","purple","lightgreen"))+ theme_classic(base_size = 14, base_family = "Helvetica")+ facet_wrap(Region~Temp, scales="free_y")+ geom_smooth ( aes(group=id.plot, colour=Region, fill=Region),method="lm")+ theme(legend.position='top') + ylab(" umax")+ xlab(" ln Dilution") r_plot01.1.3