Actinobacteria_growthcurves

Packages

if(require(growthrates)){library(growthrates)}else{install.packages("growthrates");library(growthrates)}

## Loading required package: growthrates

## Loading required package: lattice

## Loading required package: deSolve

if(require(data.table)){library(data.table)}else{install.packages("data.table");library(data.table)}

## Loading required package: data.table

if(require(ggplot2)){library(ggplot2)}else{install.packages("ggplot2");library(ggplot2)}

## Loading required package: ggplot2

RHODOCOCCUS RHODNII

BHI

Data

Import

data <- read.csv("~/Dropbox/Conjugation_in_other_species/growthrates/R.rodnii/R.rodni_BHI.csv")

Visualization

xyplot(OD~time | strain, data= data, groups=replicate)

Growth curve

Construction

means.s <- c()
sd.s <- c()
times <- data$time[data$strain=="S"]

for (i in 1:length(times)){
  means.s[i] <- mean(data$OD[data$strain == "S" & data$time==times[i]])
  sd.s[i] <- sd(data$OD[data$strain == "S" & data$time==times[i]]) / sqrt(length((data$OD[data$strain == "S" & data$time==times[i]])))
}

susceptibles <- data.frame(time=times, strain=c("S"), mean.OD=means.s, se.OD=sd.s)



means.r <- c()
sd.r <- c()

for (i in 1:length(times)){
  means.r[i] <- mean(data$OD[data$strain == "R" & data$time==times[i]])
  sd.r[i] <- sd(data$OD[data$strain == "R" & data$time==times[i]]) / sqrt(length((data$OD[data$strain == "R" & data$time==times[i]])))
}

resistants <- data.frame(time=times, strain=c("R"), mean.OD=means.r, se.OD=sd.r)

data.plot <- rbind(susceptibles, resistants)

Visualization

ggplot(data.plot, aes(x=time, y=mean.OD, group=strain, col=strain))+
  geom_point()+
  geom_line()+
  geom_errorbar(aes(ymin=mean.OD-se.OD, ymax=mean.OD+se.OD))+
  labs(x="Time", y="OD 600", title="R.rhodnii growth in BHI media")+
  scale_color_discrete(name="Strain", labels=c("Susceptible", "Resistant"))

Estimating the model

Non-parametric estimation of max growth rate (to estimate maxgrowth values that we will use in next model)

fitall <- all_splines(OD ~ time | strain + replicate, data=data, spar = 0.5)
maxgrowth <- data.frame(coef(fitall))
maxgrowth # non-parametric model parameter estimates

##               y0     mumax
## R:1 0.0043955046 0.2958175
## S:1 0.0069859075 0.2991656
## R:2 0.0002007227 1.2012965
## S:2 0.0001691623 1.2892352
## R:3 0.0017902722 0.3454103
## S:3 0.0023940110 0.4001799

Plot production

par(mfrow = c(3,2))
par(mar = c(2.5,4,2,1))
plot(fitall)

Parametric fit of growth using non-parametric data

mean.mumax <- mean(unlist(maxgrowth$mumax))
mean.y_0 <- mean(unlist(maxgrowth$y0))

##initial parameters and box constraints## 
p <- c(y0 = mean.y_0, mumax = mean.mumax, K = 0.5, h0 = 2)
lower <- c(y0 = 0.00001, mumax = 1e-3, K = 0.05, h0 = 0)
upper <- c(y0 = 0.6, mumax = 2, K = 2, h0 = 15)

Fitting Barayani-Roberts growth model

many_baranyi <- all_growthmodels(
                                  OD ~ grow_baranyi(time, parms) | strain + replicate, 
                                  data = data, 
                                  p = p, 
                                  lower = lower, 
                                  upper = upper, 
                                  log = "y")
 
new.params <- as.data.frame(coef(many_baranyi))
new.params

##             y0     mumax         K           h0
## R:1 0.02839954 0.1274673 0.6787982 3.653051e-08
## S:1 0.05032289 0.1075876 1.0640609 3.622217e-08
## R:2 0.02982710 0.1355132 0.7425861 2.652790e-08
## S:2 0.05258908 0.1203697 1.1245876 6.847910e-08
## R:3 0.02117541 0.1308403 0.6541136 3.652134e-08
## S:3 0.04303626 0.1138580 1.0263450 2.311012e-08

Adjusting parameter h0

new.params$h0 <- mean(new.params$h0)
new.params

##             y0     mumax         K           h0
## R:1 0.02839954 0.1274673 0.6787982 3.789852e-08
## S:1 0.05032289 0.1075876 1.0640609 3.789852e-08
## R:2 0.02982710 0.1355132 0.7425861 3.789852e-08
## S:2 0.05258908 0.1203697 1.1245876 3.789852e-08
## R:3 0.02117541 0.1308403 0.6541136 3.789852e-08
## S:3 0.04303626 0.1138580 1.0263450 3.789852e-08

Fitting Barayani-Roberts growth model with adjusted parameters

many_baranyi.adjusted <- all_growthmodels(
                                  OD ~ grow_baranyi(time, parms) | strain + replicate, 
                                  data = data, 
                                  p = new.params, 
                                  lower = lower, 
                                  upper = upper, 
                                  log = "y")
 
param.estimates <- as.data.frame(coef(many_baranyi.adjusted))
param.estimates # parametric model parameter estimates

##             y0     mumax         K           h0
## R:1 0.02839954 0.1274673 0.6787982 3.789853e-08
## S:1 0.05032291 0.1075876 1.0640611 1.816131e-11
## R:2 0.02982710 0.1355132 0.7425861 2.670281e-11
## S:2 0.05258911 0.1203697 1.1245878 1.291550e-11
## R:3 0.02117541 0.1308402 0.6541136 1.889322e-11
## S:3 0.04303632 0.1138579 1.0263454 1.276146e-11

Plot production

par(mfrow = c(3,2))
par(mar = c(2.5,4,2,1))
plot(many_baranyi.adjusted)

** Statistical Analysis**

t.test(x=param.estimates[c(2,4,6),]$mumax, y=param.estimates[c(1,3,5),]$mumax, paired=TRUE)

## 
##  Paired t-test
## 
## data:  param.estimates[c(2, 4, 6), ]$mumax and param.estimates[c(1, 3, 5), ]$mumax
## t = -12.575, df = 2, p-value = 0.006265
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.02326671 -0.01140371
## sample estimates:
## mean of the differences 
##             -0.01733521

t.test(x=param.estimates[c(2,4,6),]$K, y=param.estimates[c(1,3,5),]$K, paired=TRUE)

## 
##  Paired t-test
## 
## data:  param.estimates[c(2, 4, 6), ]$K and param.estimates[c(1, 3, 5), ]$K
## t = 97.018, df = 2, p-value = 0.0001062
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.3629871 0.3966772
## sample estimates:
## mean of the differences 
##               0.3798322

** Box/dot plot**

param.estimates$strain <- substr(rownames(param.estimates),1,1)

boxplot(param.estimates$mumax~param.estimates$strain)
stripchart(param.estimates$mumax~param.estimates$strain, vertical = TRUE,
    method = "jitter", add = TRUE, pch = 20, col = 'blue')

boxplot(param.estimates$K~param.estimates$strain)
stripchart(param.estimates$K~param.estimates$strain, vertical = TRUE,
    method = "jitter", add = TRUE, pch = 20, col = 'blue')

>#####Minimal media

Data

Import

data <- read.csv("~/Dropbox/Conjugation_in_other_species/growthrates/R.rodnii/R.rodni_MM.csv")

Visualization

xyplot(OD~time | strain, data= data, groups=replicate)

Growth curve

Construction

means.s <- c()
sd.s <- c()
times <- data$time[data$strain=="S"]

for (i in 1:length(times)){
  means.s[i] <- mean(data$OD[data$strain == "S" & data$time==times[i]])
  sd.s[i] <- sd(data$OD[data$strain == "S" & data$time==times[i]]) / sqrt(length((data$OD[data$strain == "S" & data$time==times[i]])))
}

susceptibles <- data.frame(time=times, strain=c("S"), mean.OD=means.s, se.OD=sd.s)



means.r <- c()
sd.r <- c()

for (i in 1:length(times)){
  means.r[i] <- mean(data$OD[data$strain == "R" & data$time==times[i]])
  sd.r[i] <- sd(data$OD[data$strain == "R" & data$time==times[i]]) / sqrt(length((data$OD[data$strain == "R" & data$time==times[i]])))
}

resistants <- data.frame(time=times, strain=c("R"), mean.OD=means.r, se.OD=sd.r)

data.plot <- rbind(susceptibles, resistants)

Visualization

ggplot(data.plot, aes(x=time, y=mean.OD, group=strain, col=strain))+
  geom_point()+
  geom_line()+
  geom_errorbar(aes(ymin=mean.OD-se.OD, ymax=mean.OD+se.OD))+
  labs(x="Time", y="OD 600", title="R.rhodnii growth in minimal media")+
  scale_color_discrete(name="Strain", labels=c("Susceptible", "Resistant"))

Estimating the model

Non-parametric estimation of max growth rate (to estimate maxgrowth values that we will use in next model)

fitall <- all_splines(OD ~ time | strain + replicate, data=data, spar = 0.5)
maxgrowth <- data.frame(coef(fitall))
maxgrowth # non-parametric model parameter estimates

##               y0     mumax
## R:1 0.0000415686 0.9308069
## S:1 0.0002029309 0.8577552
## R:2 0.0001424462 0.8239395
## S:2 0.0042707954 0.4392541
## R:3 0.0038733321 0.3874977
## S:3 0.0001663033 0.4059563

Plot production

par(mfrow = c(3,2))
par(mar = c(2.5,4,2,1))
plot(fitall)

Parametric fit of growth using non-parametric data

mean.mumax <- mean(unlist(maxgrowth$mumax))
mean.y_0 <- mean(unlist(maxgrowth$y0))

##initial parameters and box constraints## 
p <- c(y0 = mean.y_0, mumax = mean.mumax, K = 0.5, h0 = 2)
lower <- c(y0 = 0.00001, mumax = 1e-2, K = 0.05, h0 = 0)
upper <- c(y0 = 0.6, mumax = 3, K = 1, h0 = 15)

Fitting Barayani-Roberts growth model

many_baranyi <- all_growthmodels(
                                  OD ~ grow_baranyi(time, parms) | strain + replicate, 
                                  data = data, 
                                  p = p, 
                                  lower = lower, 
                                  upper = upper, 
                                  log = "y")
 
new.params <- as.data.frame(coef(many_baranyi))
new.params

##              y0     mumax         K        h0
## R:1 0.011970518 0.8608436 0.4245219 12.040819
## S:1 0.015485414 0.9447999 0.4509619 13.786392
## R:2 0.007834689 0.3845183 0.3360463  3.932289
## S:2 0.010709828 0.3469516 0.3835736  3.315672
## R:3 0.021721023 0.8609228 0.4295527 12.692632
## S:3 0.022040978 0.9572198 0.4820828 14.008130

Adjusting parameter h0

new.params$h0 <- mean(new.params$h0)
new.params

##              y0     mumax         K       h0
## R:1 0.011970518 0.8608436 0.4245219 9.962656
## S:1 0.015485414 0.9447999 0.4509619 9.962656
## R:2 0.007834689 0.3845183 0.3360463 9.962656
## S:2 0.010709828 0.3469516 0.3835736 9.962656
## R:3 0.021721023 0.8609228 0.4295527 9.962656
## S:3 0.022040978 0.9572198 0.4820828 9.962656

Fitting Barayani-Roberts growth model with adjusted parameters

many_baranyi.adjusted <- all_growthmodels(
                                  OD ~ grow_baranyi(time, parms) | strain + replicate, 
                                  data = data, 
                                  p = new.params, 
                                  lower = lower, 
                                  upper = upper, 
                                  log = "y")
 
param.estimates <- as.data.frame(coef(many_baranyi.adjusted))
param.estimates # parametric model parameter estimates

##              y0     mumax         K        h0
## R:1 0.011970582 0.8608474 0.4245218 12.040890
## S:1 0.015485500 0.9448064 0.4509618 13.786516
## R:2 0.007834579 0.3845156 0.3360465  3.932224
## S:2 0.010709860 0.3469518 0.3835736  3.315680
## R:3 0.021721177 0.8609165 0.4295530 12.692544
## S:3 0.022041234 0.9572376 0.4820826 14.008458

Plot production

par(mfrow = c(3,2))
par(mar = c(2.5,4,2,1))
plot(many_baranyi.adjusted)

** Statistical Analysis**

t.test(x=param.estimates[c(2,4,6),]$mumax, y=param.estimates[c(1,3,5),]$mumax, paired=TRUE)

## 
##  Paired t-test
## 
## data:  param.estimates[c(2, 4, 6), ]$mumax and param.estimates[c(1, 3, 5), ]$mumax
## t = 1.1137, df = 2, p-value = 0.3813
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1362254  0.2313696
## sample estimates:
## mean of the differences 
##              0.04757211

t.test(x=param.estimates[c(2,4,6),]$K, y=param.estimates[c(1,3,5),]$K, paired=TRUE)

## 
##  Paired t-test
## 
## data:  param.estimates[c(2, 4, 6), ]$K and param.estimates[c(1, 3, 5), ]$K
## t = 5.2745, df = 2, p-value = 0.03412
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.007768857 0.076562353
## sample estimates:
## mean of the differences 
##               0.0421656

** Box/dot plot**

param.estimates$strain <- substr(rownames(param.estimates),1,1)

boxplot(param.estimates$mumax~param.estimates$strain)
stripchart(param.estimates$mumax~param.estimates$strain, vertical = TRUE,
    method = "jitter", add = TRUE, pch = 20, col = 'blue')

boxplot(param.estimates$K~param.estimates$strain)
stripchart(param.estimates$K~param.estimates$strain, vertical = TRUE,
    method = "jitter", add = TRUE, pch = 20, col = 'blue')

Actinobacteria_growthcurves_analysis

Dr. Sonsiray Alvarez-Narvaez

November 18, 2019

Packages

RHODOCOCCUS RHODNII

BHI