OD and Luminescence modeling of the INSA_Lyon1 project: FIAT LUX
INSA_Lyon1
2022-10-08
This R markdown document has been done by our INSA_Lyon1 team for the iGEM 2022 competition. The aim of our project FIAT LUX is to create luminescent bacteria to study their behavior inside a living organism. At the lab, we transformed bacteria by adding one of the 4 following plasmids: pSEVA521, pSEVA521-fiatlux, pSEVA531, pSEVA531-fiatlux. Then we measured the optical density (OD) and the luminescence of these cultures. In this modeling page, we will present the work we realized to model and analyze the OD and the luminescence.
If you want more information on our project, feel free to take a look at our wiki.
OD and luminescence modeling and analysis of our project
Set your working directory with the setwd() command.
OD MODELING
Libraries
library(nlstools)
library(dplyr)
library(reshape2)
library(ggplot2)
library(nlme)Data
data_OD_Lyon1 <- read.table("Tecan_23_09_OD.csv", head = T,sep=",")
data_OD_INSA <- read.table("Tecan_20_09_OD.csv", head = T,sep=",")The file Tecan_23_09_OD.csv (variable data_OD_Lyon1) contains the data collected for the September 23 Tecan, taken with the Lyon1 machine whose plate layout is shown below:
Plate layout of the Tecan of September 23
Thus, the resulting files contain:
a column “Time_seconds” containing the time (in seconds).
a column containing the OD value given by the machine for each of the tecan cell. For example, the column pSEVA521_0_A corresponds to the OD values collected by the machine for the first replicate of pSEVA521 without tetracycline, which corresponds to the cell A1 on the plate layout. For the OD values of the plasmid containing fiatlux, we can see on the layout that we have 2 replicates per line for each antibiotic concentration. As a result, we have for example the column pSEVA521-fiatlux_0_E_1 corresponding to the values given by the machine for the plasmid pSEVA521-fitalux with 0 µg/mL in the cell E1, and pSEVA521-fiatlux_0_E_2 corresponding to the same conditions but in the cell E2.
5 columns containing the mean of the OD values given by the machine, for LB medium with a concentration of tetracycline of 0, 3, 5, 7 and 10 µg/mL respectively (mean_LB_0, mean_LB_3, etc.).
a column containing the OD value corrected by the corresponding mean LB value for each cell of the tecan. For example, the column OD_pSEVA521_0_A has been obtained by retrieving the column mean_LB_0 from the column pSEVA_521_0_A. These are the columns we will use for this OD modeling part.
To compute the mean of the OD values of the LB medium for each Tetracycline concentration, we first plotted the graph of the OD against time for each replicate, to see if there was any aberrant value.
plot(LB_0_A~Time_seconds, data=data_OD_INSA,main="OD as a function of time on the LB medium without antibiotics",cex.main=0.8,xlab="Time (s)", ylab = "OD")
We were not supposed to see any growth for the LB replicates, so we did not use the values of the growing curves to compute the mean.
The file Tecan_20_09_OD.csv (variable data_OD_INSA) contains the data collected for the September 20 Tecan, taken with the INSA machine whose plate layout is shown below:
Plate layout of the Tecan of September 20
Implementation
First, we plotted the OD as a function of time for the 4 different plasmids (pSEVA521, pSEVA531, pSEVA521-fiatlux, pSEVA531-fiatlux) to get an insight of the differences. As the values were taken on 2 different machines and the values for pSEVA531-fiatlux were taken only on the INSA machine, we had to use both tables data_OD_Lyon1 and data_OD_INSA. However, data on the Lyon1 machine were collected for a longer time than the ones on the INSA machine. We decided to select only part of the data of the Lyon1 machine to obtain data on the same period of time.
cropped_data_OD_Lyon1=data_OD_Lyon1[0:145,]plot(OD_pSEVA521_3_A~Time_seconds,data=cropped_data_OD_Lyon1,col = "#85CFE8", pch =0,cex=0.9, ylim=c(0,0.7),main="OD as a function of time, antibiotic concentration = 3 µg/mL", xlab= "Time (s)",ylab="OD", cex.main=0.8)
points(OD_pSEVA521_fiatlux_3_F_2~Time_seconds,data=cropped_data_OD_Lyon1,col = "#3000C0",pch =1,cex=0.9)
points(OD_pSEVA531_3_C~Time_seconds,data=cropped_data_OD_Lyon1,col = "#E07575",pch=2,cex=0.9)
points(OD_pSEVA531_fiatlux_3_H~Time_seconds,data=data_OD_INSA,col = "#B30000",pch =3,cex=0.9)
legend(x=90000,y=0.38,legend=c("pSEVA521","pSEVA521-fiatlux","pSEVA531","pSEVA531-fiatlux"),
col=c("#85CFE8","#3000C0","#E07575","#B30000"),pch=c(0,1,2,3),cex=0.8, text.font=c(1,3,1,3))
Below is the code to adjust the 3 different models (Baranyi, Gompertz and Verhulst) on our curves, compare their AIC, and get the parameters:
On the example below, we used the values of the 4th replicate of the plasmid pSEVA521 with 5 µg/mL. First, we plotted the OD as a function of time to see the shape of the curve:
plot(OD_pSEVA521_5_D~Time_seconds,data = data_OD_Lyon1,main = "OD of plasmid pSEVA521 as a function of Time, with 5µg/mL of Tetracycline", xlab="Time (s)", ylab = "OD", cex.main=0.8)
Looking at this graph, we decided to crop the table to avoid the first value that seems to be an outlier.
cropped_data_OD_Lyon1 =data_OD_Lyon1[data_OD_Lyon1$Time_seconds>63,]Baranyi model
baranyi =as.formula(OD_pSEVA521_5_D ~ ymax+log10(-1+exp(mu*lag)+exp(mu*Time_seconds))-log10(exp(mu*Time_seconds)-1+exp(mu*lag)*10^(ymax-y0)))
valinit<-list(y0=0.005, ymax=0.4, mu=0.0001, lag=0.5)
preview(formula=baranyi, data=cropped_data_OD_Lyon1, start=valinit)
##
## RSS: 0.191ajusbar <- nls(formula=baranyi, data=cropped_data_OD_Lyon1, start=valinit, control = nls.control(maxiter = 1000, minFactor = 3.05176e-06, warnOnly = TRUE))
plotfit(ajusbar, smooth=TRUE, main="OD as a function of Time", xlab="Time (s)", ylab = "OD", sub="Plasmid pSEVA521, Antibiotic Concentration = 5µg/mL", cex.main=0.8)
legend(x=70000,y=0.3,legend=c("Experimental","Adjusted Baranyi model"),
col=c("black","red"),pch=c("o"," "),lty=c(0,1),cex=0.8)
overview(ajusbar)##
## ------
## Formula: OD_pSEVA521_5_D ~ ymax + log10(-1 + exp(mu * lag) + exp(mu *
## Time_seconds)) - log10(exp(mu * Time_seconds) - 1 + exp(mu *
## lag) * 10^(ymax - y0))
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## y0 5.648e-02 1.342e-03 42.099 < 2e-16 ***
## ymax 4.207e-01 2.317e-04 1815.582 < 2e-16 ***
## mu 4.222e-05 4.935e-07 85.563 < 2e-16 ***
## lag -3.495e+03 6.578e+02 -5.314 3.24e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.002041 on 175 degrees of freedom
##
## Number of iterations to convergence: 7
## Achieved convergence tolerance: 3.456e-06
##
## ------
## Residual sum of squares: 0.000729
##
## ------
## t-based confidence interval:
## 2.5% 97.5%
## y0 5.383347e-02 5.912913e-02
## ymax 4.202560e-01 4.211707e-01
## mu 4.125056e-05 4.319847e-05
## lag -4.793638e+03 -2.197157e+03
##
## ------
## Correlation matrix:
## y0 ymax mu lag
## y0 1.0000000 -0.1923413 0.5484226 0.7588357
## ymax -0.1923413 1.0000000 -0.5697858 -0.4440230
## mu 0.5484226 -0.5697858 1.0000000 0.9433294
## lag 0.7588357 -0.4440230 0.9433294 1.0000000AIC(ajusbar)## [1] -1703.647Gompertz model
gomp = as.formula(OD_pSEVA521_5_D~alpha*exp(-exp(-k*(Time_seconds-beta))))
valinit<-list(alpha = 0.70, beta = 1.698353, k = 0.00007 )
preview(formula=gomp, data=cropped_data_OD_Lyon1, start=valinit)
##
## RSS: 15ajusgomp <- nls(formula=gomp, data=cropped_data_OD_Lyon1, start=valinit)
plotfit(ajusgomp, smooth=TRUE, main="OD as a function of Time", xlab="Time (s)", ylab = "OD", sub="Plasmid pSEVA521, Antibiotic Concentration = 5µg/mL", cex.main=0.8)
legend(x=70000,y=0.3,legend=c("Experimental","Adjusted Gompertz model"),
col=c("black","red"),pch=c("o"," "),lty=c(0,1),cex=0.8)
overview(ajusgomp)##
## ------
## Formula: OD_pSEVA521_5_D ~ alpha * exp(-exp(-k * (Time_seconds - beta)))
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## alpha 4.191e-01 4.140e-04 1012.42 <2e-16 ***
## beta 9.363e+03 1.346e+02 69.57 <2e-16 ***
## k 4.984e-05 4.371e-07 114.01 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.004208 on 176 degrees of freedom
##
## Number of iterations to convergence: 6
## Achieved convergence tolerance: 2.129e-06
##
## ------
## Residual sum of squares: 0.00312
##
## ------
## t-based confidence interval:
## 2.5% 97.5%
## alpha 4.182871e-01 4.199211e-01
## beta 9.097390e+03 9.628570e+03
## k 4.897434e-05 5.069976e-05
##
## ------
## Correlation matrix:
## alpha beta k
## alpha 1.00000000 0.04093045 -0.4634946
## beta 0.04093045 1.00000000 0.4971135
## k -0.46349458 0.49711354 1.0000000AIC(ajusgomp)## [1] -1445.594Verhulst model
ver = as.formula(OD_pSEVA521_5_D~k/(1+(k/n0-1)*exp(-r0*Time_seconds)))
valinit<-list(n0 = 0.08, r0 = 0.00012, k = 0.7)
preview(formula=ver, data=cropped_data_OD_Lyon1, start=valinit) 
##
## RSS: 13.6ajusver <- nls(formula=ver, data=cropped_data_OD_Lyon1, start=valinit)
plotfit(ajusver, smooth=TRUE, main="OD as a function of Time", xlab="Time (s)", ylab = "OD", sub="Plasmid pSEVA521, Antibiotic Concentration = 5µg/mL", cex.main=0.8)
legend(x=70000,y=0.3,legend=c("Experimental","Adjusted Verhulst model"),
col=c("black","red"),pch=c("o"," "),lty=c(0,1),cex=0.8)
overview(ajusver)##
## ------
## Formula: OD_pSEVA521_5_D ~ k/(1 + (k/n0 - 1) * exp(-r0 * Time_seconds))
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## n0 1.049e-01 1.949e-03 53.80 <2e-16 ***
## r0 6.318e-05 9.851e-07 64.14 <2e-16 ***
## k 4.175e-01 6.830e-04 611.32 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.007218 on 176 degrees of freedom
##
## Number of iterations to convergence: 6
## Achieved convergence tolerance: 9.407e-06
##
## ------
## Residual sum of squares: 0.00917
##
## ------
## t-based confidence interval:
## 2.5% 97.5%
## n0 1.010287e-01 1.087225e-01
## r0 6.124086e-05 6.512908e-05
## k 4.161749e-01 4.188707e-01
##
## ------
## Correlation matrix:
## n0 r0 k
## n0 1.0000000 -0.8485189 0.2101277
## r0 -0.8485189 1.0000000 -0.3818363
## k 0.2101277 -0.3818363 1.0000000AIC(ajusver)## [1] -1252.391The Baranyi model results to be the one with the lowest AIC. We collected the values of the parameters obtained for each condition with the Baranyi model in a table.
STATISTICAL ANALYSIS OF THE PARAMETERS
We compared the parameters of the Baranyi models for each condition thanks to statistical models and tests in order to determine significant differences between our conditions.
Libraries
library(nlme)
library(lme4)
library(lmerTest)
library(WWGbook)
library(emmeans)
library(forcats)
library(dplyr)Global analysis of the parameter
Data
data=read.table('mu.csv',head=T,sep=",")
head(data)## Antibiotic Machine Plasmid Value
## 1 TET0 INSA pSEVA521 5.216e-06
## 2 TET0 INSA pSEVA521 1.776e-06
## 3 TET3 INSA pSEVA521 1.657e-06
## 4 TET3 INSA pSEVA521 1.759e-06
## 5 TET3 INSA pSEVA521 1.897e-06
## 6 TET3 INSA pSEVA521 1.941e-06Variables
machine=as.factor(data$Machine)#2 different measurement machines: INSA and LYON1
plasmid=as.factor(data$Plasmid)#4 different plasmids: pSEVA521, pSEVA531, pSEVA521-fiatlux and pSEVA531-fiatlux
mu=data$Value #growth rate of the bacteria
antibio=as.factor(data$Antibiotic)#antibiotic concentration (µg/mL of medium): 5 different concentrations
antibio=fct_relevel(antibio,"TET0","TET3","TET5","TET7","TET10")TET0= medium without antibiotics antibiotic=selective medium (our plasmid contains resistance for tetracycline)
Objective : we want to study the bacterial behavior according to which plasmid it contains, and in which antibiotic concentration it grows We did measurement during 48h on two different machines. We will take in account that there can be variations in our data depending the machine we used
Here is our experiment plan:
table(antibio, plasmid)## plasmid
## antibio pSEVA521 pSEVA521_fiatlux pSEVA531 pSEVA531_fiatlux
## TET0 5 4 7 4
## TET3 8 8 8 4
## TET5 8 7 8 4
## TET7 4 8 7 4
## TET10 8 8 8 4Graphics
Boxplot representing the growtn rate µ according to the plasmid
par(mfrow=c(2,3))
mu_unit <- expression(µ ~ (s^-1))
boxplot(mu[antibio=="TET0"]~plasmid[antibio=="TET0"], main="0 µg/mL of Tetracycline", cex.main=0.8, xlab="Plasmid",ylab=mu_unit,names=c("pSEVA521","pSEVA521-fiatlux","pSEVA531","pSEVA531-fiatlux"))
boxplot(mu[antibio=="TET3"]~plasmid[antibio=="TET3"], main="3 µg/mL of Tetracycline", cex.main=0.8,xlab="Plasmid",ylab=mu_unit,names=c("pSEVA521","pSEVA521-fiatlux","pSEVA531","pSEVA531-fiatlux"))
boxplot(mu[antibio=="TET5"]~plasmid[antibio=="TET5"], main="5 µg/mL of Tetracycline", cex.main=0.8,xlab="Plasmid",ylab=mu_unit,names=c("pSEVA521","pSEVA521-fiatlux","pSEVA531","pSEVA531-fiatlux"))
boxplot(mu[antibio=="TET7"]~plasmid[antibio=="TET7"], main="7 µg/mL of Tetracycline", cex.main=0.8,xlab="Plasmid",ylab=mu_unit,names=c("pSEVA521","pSEVA521-fiatlux","pSEVA531","pSEVA531-fiatlux"))
boxplot(mu[antibio=="TET10"]~plasmid[antibio=="TET10"], main="10 µg/mL of Tetracycline", cex.main=0.8,xlab="Plasmid",ylab=mu_unit,names=c("pSEVA521","pSEVA521-fiatlux","pSEVA531","pSEVA531-fiatlux"))
par(mfrow=c(1,1))
The apparent heteroscedasticity of the variance on the boxplots is mainly due to the fact that we don’t have the same number of replicates for each plasmid.
mean_of_mu_unit = expression("Mean of µ " ~ (s^-1))
interaction.plot(antibio,plasmid, mu, las=1, main = "Interaction plot between the amount of antibiotics and the plasmid contained in the bacteria",cex.main=0.8, xlab="Tetracycline concentration (µg/mL)",ylab=mean_of_mu_unit,,trace.label = "Plasmid")
The bacteria containing a plasmid with fiatlux seem to growth differently depending on the presence of antibiotic: on medium without Tetracycline (TET0), they have a lower growth rate than with antibiotics A concentration of 10 µg/mL of Tetracycline seems to lower the growth rate of the empty plasmid. Plasmids with fiatlux seem to react similarly to the antibiotic concentration.
coplot(mu~machine|antibio:plasmid, rows=1, ylab = mu_unit, xlab = "Machine (INSA or Lyon1)")
This in an insight of the effect of the machine There is no big appearent effect here, but we cannot conclude Indeed our experimental plan lacks some datas: for instance pSEVA521-fiatlux on INSA machine.
Implementation of the model and results
We use an lmer model of the library lme4, used to fit linear mixed-effects models.
model1 <- lmer(mu~plasmid*antibio+(1|machine))We include the effect of the machine as a random effect. Then, we test the effect of the machine, with ranova, used to test the effect of the random factors:
ranova(model1)## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## mu ~ plasmid + antibio + (1 | machine) + plasmid:antibio
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 22 991.11 -1938.2
## (1 | machine) 21 982.61 -1923.2 17.002 1 3.734e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The p-value (3.734e-05) indicates a significant effect. During the rest of our study, we will take this into account.
Now, let’s test the effect of the fixed factors with anova:
anova(model1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## plasmid 1.5456e-08 5.1519e-09 3 681.42 16.8150 1.511e-10 ***
## antibio 9.2517e-09 2.3129e-09 4 685.06 7.5489 5.925e-06 ***
## plasmid:antibio 1.9163e-08 1.5969e-09 12 684.53 5.2120 2.004e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(model1)#test effect of plasmid and antibio As we supposed by observing the data on graph, the plasmid and the antibiotic concentration affect the growth rate in a significant way (taking a risk of 0.05) We will perform further test to have a better insight of the situation
emmeans.out<-emmeans(model1,~plasmid*antibio)## Cannot use mode = "kenward-roger" because *pbkrtest* package is not installedplot(emmeans.out,main="Estimated marginal means for each condition",xlab=mu_unit)
According to this plot, we can say that: - on medium without antibiotics (TET0), plasmids have around the same growth rate(except pSEVA531-fiatlux that is a little lower) On medium with 3, 5, and 7 µg/mL of Tetracycline, we observe a lower growth rate for the empty plasmids (pSEVA531 and pSEVA521) than for the plasmids with fiatlux Overall, pSEVA531-fiatlux and pSEVA521-fiatlux seem to have the same effect on the growth rate.
Now, we perform a MCP test:
pairs(emmeans.out,adjust="none",by=c("antibio"))## antibio = TET0:
## contrast estimate SE df t.ratio p.value
## pSEVA521 - pSEVA521_fiatlux 8.38e-06 1.19e-05 691 0.706 0.4801
## pSEVA521 - pSEVA531 1.09e-05 1.03e-05 685 1.065 0.2874
## pSEVA521 - pSEVA531_fiatlux 2.59e-05 1.20e-05 698 2.159 0.0312
## pSEVA521_fiatlux - pSEVA531 2.56e-06 1.12e-05 699 0.228 0.8199
## pSEVA521_fiatlux - pSEVA531_fiatlux 1.75e-05 1.31e-05 713 1.343 0.1799
## pSEVA531 - pSEVA531_fiatlux 1.50e-05 1.11e-05 693 1.348 0.1781
##
## antibio = TET3:
## contrast estimate SE df t.ratio p.value
## pSEVA521 - pSEVA521_fiatlux -3.88e-05 9.00e-06 702 -4.315 <.0001
## pSEVA521 - pSEVA531 4.21e-06 8.75e-06 683 0.481 0.6309
## pSEVA521 - pSEVA531_fiatlux -3.51e-05 1.09e-05 696 -3.213 0.0014
## pSEVA521_fiatlux - pSEVA531 4.30e-05 9.00e-06 702 4.783 <.0001
## pSEVA521_fiatlux - pSEVA531_fiatlux 3.75e-06 1.15e-05 718 0.326 0.7448
## pSEVA531 - pSEVA531_fiatlux -3.93e-05 1.09e-05 696 -3.598 0.0003
##
## antibio = TET5:
## contrast estimate SE df t.ratio p.value
## pSEVA521 - pSEVA521_fiatlux -2.28e-05 9.30e-06 700 -2.447 0.0147
## pSEVA521 - pSEVA531 9.92e-06 8.75e-06 683 1.133 0.2576
## pSEVA521 - pSEVA531_fiatlux -3.53e-05 1.09e-05 696 -3.228 0.0013
## pSEVA521_fiatlux - pSEVA531 3.27e-05 9.30e-06 700 3.513 0.0005
## pSEVA521_fiatlux - pSEVA531_fiatlux -1.25e-05 1.17e-05 717 -1.065 0.2874
## pSEVA531 - pSEVA531_fiatlux -4.52e-05 1.09e-05 696 -4.136 <.0001
##
## antibio = TET7:
## contrast estimate SE df t.ratio p.value
## pSEVA521 - pSEVA521_fiatlux -2.64e-05 1.15e-05 718 -2.289 0.0224
## pSEVA521 - pSEVA531 -3.04e-06 1.11e-05 693 -0.273 0.7847
## pSEVA521 - pSEVA531_fiatlux -1.95e-05 1.24e-05 683 -1.576 0.1155
## pSEVA521_fiatlux - pSEVA531 2.33e-05 9.37e-06 705 2.488 0.0131
## pSEVA521_fiatlux - pSEVA531_fiatlux 6.85e-06 1.15e-05 718 0.595 0.5522
## pSEVA531 - pSEVA531_fiatlux -1.65e-05 1.11e-05 693 -1.481 0.1390
##
## antibio = TET10:
## contrast estimate SE df t.ratio p.value
## pSEVA521 - pSEVA521_fiatlux -6.78e-05 9.00e-06 702 -7.534 <.0001
## pSEVA521 - pSEVA531 -4.43e-05 8.75e-06 683 -5.057 <.0001
## pSEVA521 - pSEVA531_fiatlux -3.20e-05 1.09e-05 696 -2.934 0.0035
## pSEVA521_fiatlux - pSEVA531 2.35e-05 9.00e-06 702 2.616 0.0091
## pSEVA521_fiatlux - pSEVA531_fiatlux 3.58e-05 1.15e-05 718 3.107 0.0020
## pSEVA531 - pSEVA531_fiatlux 1.22e-05 1.09e-05 696 1.119 0.2636
##
## Degrees-of-freedom method: satterthwaiteOn medium without TET: we observe no significant effect of the plasmid. On TET3 and TET5: the presence of fiatlux in the plasmid has an effect. On TET7 : the presence of fiatlux in the plasmid has a significant effect on the growth rate for pSEVA521 but not for pSEVA531. On TET10, pSEVA531 and pSEVA531-fiatlux have the same effect. We must be careful when interpreting a result as non-significative because the test lost power because there is a lot of cases to compare.
We check the homogeneity of our variance thanks to this plot:
plot(model1,main="Check the homogeneity of the variance",cex.main=0.6)
Analysis of the parameter for a growth on a medium without antibiotics
As we intend to use our tool in situ, to infect plants, we will focus on the differences of our plasmid on medium without antibiotics. The objective is to find the best vector for our tool between pSEVA531 or pSEVA521.
Data
data1=filter(data,Antibiotic=="TET0")Variables
machine1=as.factor(data1$Machine)#2 differents measurment machines: INSA and LYON1
plasmid1=as.factor(data1$Plasmid)#4 different plasmid: pSEVA521, pSEVA531, pSEVA521-fiatlux and pSEVA531-fiatlux
mu1=data1$Value#growth rate of the bacteriaGraphics
boxplot(mu1~plasmid1, main="Boxplot representing the bacterial growth rate (µ) according to each plasmid, without antibiotics",cex.main=0.8,xlab="Plasmid",ylab=mu_unit)
The growth rate seems lower for both pSEVA531 plasmids than for the pSEVA521 plamsids.
Implementation of the model and results
model2 <- lmer(mu1~plasmid1+(1|machine1))anova(model2)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## plasmid1 1.3618e-09 4.5394e-10 3 9.1638 5.6776 0.01787 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The plasmid has a significant effect on the growth rate.
ranova(model2)## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## mu1 ~ plasmid1 + (1 | machine1)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 6 157.65 -303.31
## (1 | machine1) 5 141.87 -273.74 31.566 1 1.928e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The machine has a significant effect on the growth rate.
MCP test:
emmeans.out1<-emmeans(model2,~plasmid1)## Cannot use mode = "kenward-roger" because *pbkrtest* package is not installedplot(emmeans.out1,ylab="Plasmid",xlab=mu_unit)
pairs(emmeans.out1,adjust="none")## contrast estimate SE df t.ratio p.value
## pSEVA521 - pSEVA521_fiatlux 2.52e-05 6.35e-06 9.16 3.975 0.0031
## pSEVA521 - pSEVA531 3.71e-06 5.31e-06 9.14 0.699 0.5022
## pSEVA521 - pSEVA531_fiatlux 6.31e-07 6.77e-06 9.18 0.093 0.9277
## pSEVA521_fiatlux - pSEVA531 -2.15e-05 6.35e-06 9.18 -3.392 0.0078
## pSEVA521_fiatlux - pSEVA531_fiatlux -2.46e-05 8.20e-06 9.22 -3.002 0.0145
## pSEVA531 - pSEVA531_fiatlux -3.08e-06 6.03e-06 9.16 -0.510 0.6218
##
## Degrees-of-freedom method: satterthwaiteAccording to the plot and the results, we can say that the growth rate of the bacteria containing pSEVA521-fiatlux is significantly higher than those of the bacteria containing pSEVA531-fiatlux. On medium without antibiotics, pSEVA521-fiatlux has a higher growth rate.
Analysis of the parameter for data collected on a single machine
As we work on data measured on the same machine, we can get rid of the influence of the machine.
Data
data2_1=filter(data,Machine=="INSA")
data2 = filter(data2_1,Plasmid=="pSEVA531" | Plasmid=="pSEVA531_fiatlux")Variables
antibio2=as.factor(data2$Antibiotic)#5 different concentration of antibiotics
antibio2=fct_relevel(antibio2,"TET0","TET3","TET5","TET7","TET10")
plasmid2=as.factor(data2$Plasmid)#2 different plasmids: pSEVA531 and pSEVA531-fiatlux
mu2=data2$Value #growth rate of the bacteriaPlan of the experiment:
table(antibio2,plasmid2)#4 replicates of pSEVA531 and pSEVA531-fiatlux## plasmid2
## antibio2 pSEVA531 pSEVA531_fiatlux
## TET0 4 4
## TET3 4 4
## TET5 4 4
## TET7 4 4
## TET10 4 4Graphics
plot(plasmid2,mu2,main = "Boxplot of the growth rate according to the plasmid, regardless of the antibiotic concentration", xlab="Plasmid",ylab=mu_unit,cex.main=0.8,names=c("pSEVA531","pSEVA531-fiatlux"))
The growth rate for the bacteria containing pSEVA531-fiatlux seems higher.
interaction.plot(antibio2,plasmid2, mu2, main = "Interaction plot between the amount of antibiotics and the plasmid contained in the bacteria",cex.main=0.8, xlab="Tetracycline concentration (µg/mL)",ylab=mean_of_mu_unit,trace.label = "Plasmid")
We can suppose that there is an interaction between the presence of fiatlux in the plasmid, and the amount of antibiotics. On medium without antibiotics, the growth rate seems lower.
Implementation of the model and results
model3=lm(mu2~antibio2*plasmid2)
anova(model3)## Analysis of Variance Table
##
## Response: mu2
## Df Sum Sq Mean Sq F value Pr(>F)
## antibio2 4 1.1851e-08 2.9628e-09 9.1440 5.926e-05 ***
## plasmid2 1 7.2770e-10 7.2772e-10 2.2460 0.1444
## antibio2:plasmid2 4 3.3188e-09 8.2971e-10 2.5607 0.0588 .
## Residuals 30 9.7204e-09 3.2401e-10
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There is a significative effect of the plasmids and significative effect of the antibiotics. There is also an interaction between the concentration of antibiotic and the plasmid
Let’s perform a Tukey test:
model3a=aov(mu2~plasmid2-1)
TukeyHSD(model3a)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mu2 ~ plasmid2 - 1)
##
## $plasmid2
## diff lwr upr p adj
## pSEVA531_fiatlux-pSEVA531 8.530655e-06 -7.853337e-06 2.491465e-05 0.29852model3b=aov(mu2[antibio2=="TET0"]~plasmid2[antibio2=="TET0"]-1)
TukeyHSD(model3b)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mu2[antibio2 == "TET0"] ~ plasmid2[antibio2 == "TET0"] - 1)
##
## $`plasmid2[antibio2 == "TET0"]`
## diff lwr upr p adj
## pSEVA531_fiatlux-pSEVA531 2.315525e-06 -5.193094e-06 9.824144e-06 0.4790582The difference in the growth rate between pSEVA531 and pSEVA531-fiatlux is not significant.
Check the linear hypotheses:
par(mfrow=c(2,2))
plot(model3)
par(mfrow=c(1,1))All the hypotheses are checked : - no funnel effect - normality of the residuals - no outliers (Cook’s distances < 3) - data are independant - the model is linear
We realized the same analysis on the other machine (Lyon1).
We can perform exactly the same study for the Nmax parameter, by taking the file Nmax.csv instead of mu.csv.
LUMINESCENCE MODELING
Apart from bacterial growth, we also wanted to have information on the luminescence produced by our fiatlux bacteria. For that we plotted several representations of luminescence as a function of OD and time.
Data
data_Lumi_OD_521 <- read.table("Tecan_23_09_Lumi_OD.csv", head = T,sep=",")
data_Lumi_OD_531 <- read.table("Tecan_20_09_Lumi_OD.csv", head = T,sep=",")These files contain:
a column “Time_seconds” containing the time (in seconds).
a column containing the OD value corrected by the corresponding mean LB value for each cell of the tecan corresponding to a plasmid with fiatlux. These are the columns we used for the OD modeling part. For example, the column OD_pSEVA521_fiatlux_0_A contains the OD values of the first replicate of the bacteria that contain pSEVA521-fiatlux and grew without antibiotics.
a column containing the luminescence value for each cell of the tecan corresponding to a plasmid with fiatlux. For example, the column Lumi_pSEVA521_fiatlux_0_A contains the luminescence values of the first replicate of the bacteria that contain pSEVA521-fiatlux and grew without antibiotics.
Graphics
Below, we plot the graph of the luminescence/OD ratio of Dickeya with pSEVA521-fiatlux against the time, for different Tetracycline concentrations.
plot(Lumi_pSEVA521_fiatlux_10_E_1/OD_pSEVA521_fiatlux_10_E_1~Time_seconds,data=data_Lumi_OD_521,col="#006600",cex=0.8,main="Luminescence/OD as a function of time for Dickeya with pSEVA521-fiatlux, for different antibiotic concentrations",cex.main=0.8,xlab="Time (s)",ylab="Luminescence/OD")
points(Lumi_pSEVA521_fiatlux_7_E_1/OD_pSEVA521_fiatlux_7_E_1~Time_seconds,data=data_Lumi_OD_521,col="#339900",cex=0.8)
points(Lumi_pSEVA521_fiatlux_5_E_1/OD_pSEVA521_fiatlux_5_E_1~Time_seconds,data=data_Lumi_OD_521,col="#66CC33",cex=0.8)
points(Lumi_pSEVA521_fiatlux_3_E_1/OD_pSEVA521_fiatlux_3_E_1~Time_seconds,data=data_Lumi_OD_521,col="#99FF00",cex=0.8)
points(Lumi_pSEVA521_fiatlux_0_E_1/OD_pSEVA521_fiatlux_0_E_1~Time_seconds,data=data_Lumi_OD_521,col="#CCFF00",cex=0.8)
legend(x=100000,y=47000,title="Antibiotic concentrations",legend = c("0 µg/mL","3 µg/mL","5 µg/mL","7 µg/mL","10 µg/mL"),
col = c("#CCFF00","#99FF00","#66CC33","#339900","#006600"),pch=1,cex=0.6)
abline(v=35000,col="#006600")
abline(v=28000,col="#339900")
abline(v=21500,col="#66CC33")
abline(v=18300,col="#99FF00")
abline(v=13000,col="#CCFF00")
As we observe that the moment of the peak of intensity vary according to the antibiotic concentration, we can suppose that the shape of the bacterial growth also varies and that the lag time increases with the concentration of antibiotics. This would explain why the peak of intensity is reached later with 10 µg/mL of Tetracycline than with lower concentrations. Below, we plot the OD as a function of time for exactly the same plasmid and concentrations to check this hypothesis. On the above graph, we added vertical lines where the peak was reached, and we will add them on the graph with the OD as well to check that the peak of luminescence is reached during the exponential phase of the bacterial growth.
plot(OD_pSEVA521_fiatlux_0_E_1~Time_seconds,data=data_Lumi_OD_521,col="#CCFF00",cex=0.8,main="OD as a function of time for Dickeya with pSEVA521-fiatlux, for different antibiotic concentrations",cex.main=0.8,xlab="Time (s)",ylab="OD")
points(OD_pSEVA521_fiatlux_10_E_1~Time_seconds,data=data_Lumi_OD_521,col="#006600",cex=0.8)
points(OD_pSEVA521_fiatlux_7_E_1~Time_seconds,data=data_Lumi_OD_521,col="#339900",cex=0.8)
points(OD_pSEVA521_fiatlux_5_E_1~Time_seconds,data=data_Lumi_OD_521,col="#66CC33",cex=0.8)
points(OD_pSEVA521_fiatlux_3_E_1~Time_seconds,data=data_Lumi_OD_521,col="#99FF00",cex=0.8)
legend('topright',title="Antibiotic concentrations",legend = c("0 µg/mL","3 µg/mL","5 µg/mL","7 µg/mL","10 µg/mL"),
col = c("#CCFF00","#99FF00","#66CC33","#339900","#006600"),pch=1,cex=0.65)
abline(v=35000,col="#006600")
abline(v=28000,col="#339900")
abline(v=21500,col="#66CC33")
abline(v=18300,col="#99FF00")
abline(v=13000,col="#CCFF00")
As we expected, the higher is the antibiotic concentration, the longer is the lag time, so the later is the exponential phase. Thanks to the vertical lines, we confirm that the peak of intensity is reached during the exponential phase of the bacterial growth.
Now, we plot the graph of the luminescence/OD as a function of time for pSEVA531-fiatlux.
plot(Lumi_pSEVA531_fiatlux_10_E/OD_pSEVA531_fiatlux_10_E~Time_seconds,data=data_Lumi_OD_531,col="#006600",cex=0.8,main="Luminescence/OD as a function of time for Dickeya with pSEVA531-fiatlux, for different antibiotic concentrations",cex.main=0.8,xlab="Time (s)",ylab="Luminescence/OD",ylim=c(0,80000))
points(Lumi_pSEVA531_fiatlux_7_E/OD_pSEVA531_fiatlux_7_E~Time_seconds,data=data_Lumi_OD_531,col="#339900",cex=0.8)
points(Lumi_pSEVA531_fiatlux_5_E/OD_pSEVA531_fiatlux_5_E~Time_seconds,data=data_Lumi_OD_531,col="#66CC33",cex=0.8)
points(Lumi_pSEVA531_fiatlux_3_E/OD_pSEVA531_fiatlux_3_E~Time_seconds,data=data_Lumi_OD_531,col="#99FF00",cex=0.8)
points(Lumi_pSEVA531_fiatlux_0_E/OD_pSEVA531_fiatlux_0_E~Time_seconds,data=data_Lumi_OD_531,col="#CCFF00",cex=0.8)
abline(v=58000,col="#006600")
abline(v=2000,col="#006600")
legend('topright',title="Antibiotic concentrations",legend = c("0 µg/mL","3 µg/mL","5 µg/mL","7 µg/mL","10 µg/mL"),
col = c("#CCFF00","#99FF00","#66CC33","#339900","#006600"),pch=1,cex=0.65)
The curves on this graph have a peak at the beginning like the other graph, but these peaks are less visible since the lower values are higher than the lower values of the other graph. However, we see 2 peaks of intensity for an antibiotic concentration of 10 µg/mL. To make sure that the luminescence work properly for this replicate, we plot its OD as a function of time and add the 2 vertical lines corresponding to the luminescence peak to check that one of the 2 is during the exponential phase of the growth.
plot(OD_pSEVA531_fiatlux_0_E~Time_seconds,data=data_Lumi_OD_531,col="#CCFF00",cex=0.8,main="OD as a function of time for Dickeya with pSEVA531-fiatlux, for different antibiotic concentrations",cex.main=0.8,xlab="Time (s)",ylab="OD")
points(OD_pSEVA531_fiatlux_10_E~Time_seconds,data=data_Lumi_OD_531,col="#006600",cex=0.8)
points(OD_pSEVA531_fiatlux_7_E~Time_seconds,data=data_Lumi_OD_531,col="#339900",cex=0.8)
points(OD_pSEVA531_fiatlux_5_E~Time_seconds,data=data_Lumi_OD_531,col="#66CC33",cex=0.8)
points(OD_pSEVA531_fiatlux_3_E~Time_seconds,data=data_Lumi_OD_531,col="#99FF00",cex=0.8)
abline(v=58000,col="#006600")
abline(v=2000,col="#006600")
legend('topright',title="Antibiotic concentrations",legend = c("0 µg/mL","3 µg/mL","5 µg/mL","7 µg/mL","10 µg/mL"),
col = c("#CCFF00","#99FF00","#66CC33","#339900","#006600"),pch=1,cex=0.65)
The second luminescence peak is during the exponential phase as expected.
Finaly, we plot the luminescence against the OD for Dickeya with pSEVA531-fiatlux, for different Tetracycline concentrations.
plot(Lumi_pSEVA531_fiatlux_10_E~OD_pSEVA531_fiatlux_10_E,data=data_Lumi_OD_531,col="#006600",cex=0.8,main="Luminescence as a function of OD for Dickeya with pSEVA531-fiatlux, for different antibiotic concentrations",cex.main=0.8,xlab="OD",ylab="Luminescence")
points(Lumi_pSEVA531_fiatlux_7_E~OD_pSEVA531_fiatlux_7_E,data=data_Lumi_OD_531,col="#339900",cex=0.8)
points(Lumi_pSEVA531_fiatlux_5_E~OD_pSEVA531_fiatlux_5_E,data=data_Lumi_OD_531,col="#66CC33",cex=0.8)
points(Lumi_pSEVA531_fiatlux_3_E~OD_pSEVA531_fiatlux_3_E,data=data_Lumi_OD_531,col="#99FF00",cex=0.8)
points(Lumi_pSEVA531_fiatlux_0_E~OD_pSEVA531_fiatlux_0_E,data=data_Lumi_OD_531,col="#CCFF00",cex=0.8)
legend(x=0.02,y=26000,title="Antibiotic concentrations",legend = c("0 µg/mL","3 µg/mL","5 µg/mL","7 µg/mL","10 µg/mL"),
col = c("#CCFF00","#99FF00","#66CC33","#339900","#006600"),pch=1,cex=0.6)
We notice that the first part of the curves are linear, starting at 0, which suggests that there is a relation between the luminescence and the OD. One of our next goals is to characterize.