initial release of v2 (based on Rmd notebook)

Author: pierremillard

README.md

@@ -2,27 +2,31 @@
 
 ## Overview
 
-These R scripts perform all analyses detailed in the following publication:
+This R code perform all analyses detailed in the following publication:
 
 > From toxic waste to beneficial nutrient: acetate boosts *Escherichia coli* growth at low glycolytic flux.
 >
-> P. Millard, S. Uttenweiler-Joseph, B. Enjalbert. bioRxiv, 2022, doi: [10.1101/2022.09.20.506926](https://doi.org/10.1101/2022.09.20.506926)
+> P. Millard, T. Gosselin-Monplaisir, S. Uttenweiler-Joseph, B. Enjalbert. bioRxiv, 2022, doi: [10.1101/2022.09.20.506926](https://doi.org/10.1101/2022.09.20.506926)
 
 The model is available in COPASI format in the `./model/` directory, and in SBML format from the Biomodels database (http://www.ebi.ac.uk/biomodels/) under identifier MODEL2005050001. Details on the model can be found in the original publication (Millard et al., 2021, eLife, doi: [10.7554/eLife.63661](https://doi.org/10.7554/eLife.63661)) and at https://github.com/MetaSys-LISBP/acetate_regulation. Experimental data are provided in the `./data/` directory. Calculation results and figures are saved in the `./results/` directory.
 
-Details on the calculations can be found in the [original publication](https://doi.org/10.1101/2022.09.20.506926) and in the R scripts.
+Details on the calculations can be found in the [original publication](https://doi.org/10.1101/2022.09.20.506926) and in the Rmd notebook.
+
+We distribute an HTML file showing the [notebook’s output after execution](https://htmlpreview.github.io/?https://github.com/MetaSys-LISBP/glucose_acetate_interplay/blob/main/notebook.html)
 
 The code is open-source and available under [GPLv3 license](https://www.gnu.org/licenses/gpl-3.0.txt).
 
 ## Dependencies
 
-Some R packages are required.
+If not yet done, download and install [Rstudio](https://posit.co/downloads/).
+
+Some R packages are also required.
 
-`RColorBrewer`, `gplots`, `pso` and `optimx` can be installed
-by running the following command in an R console:
+`RColorBrewer`, `gplots`, `pso`, `colorspace` and `optimx` can be installed
+by running the following command in the Rstudio console:
 
 ```bash
-install.packages(c("RColorBrewer", "gplots", "pso", "optimx"))
+install.packages(c("RColorBrewer", "gplots", "pso", "optimx", "colorspace"))
 ```
 
 `CoRC` can be installed
@@ -40,38 +44,22 @@ Additional information on CoRC installation and usage are available from the CoR
 
 ## Usage
 
-To run all analyses detailed in the publication and reproduce Figures 1-4, 6 and 7:
+To run all analyses detailed in the publication and reproduce Figures 1-4, 6-7 and EV1-2:
 
 - download the [GitHub repository tarball](https://github.com/pierremillard/glucose_acetate_interplay/archive/refs/heads/master.zip) and unpack it somewhere on your disk
 
-- open a terminal
-
-- go to the code directory, e.g.:
-
-```bash
-cd /home/usr/data/glucose_acetate_interplay/
-```
-
-- open an R session:
+- start Rstudio
 
-```bash
-R
-```
-
-- run simulations, generate all figures and perform statistical analyses with:
-
-```bash
-source("run_calculations.R")
-```
+- open the file `./notebook.Rmd`
 
-All results will be saved in the `./results/` directory.
+The notebook will open in Rstudio, in each cell you can read/modify/execute a proposed code as well as read accompanying comments.
 
 ## Questions
 If you have any question, please open a new *issue*
 to [our GitHub issue tracker](https://github.com/MetaSys-LISBP/glucose_acetate_interplay/issues) so we could discuss together. 
 
 ## How to cite
-Millard P., Uttenweiler-Joseph S., Enjalbert B. From toxic waste to beneficial nutrient: acetate boosts *Escherichia coli* growth at low glycolytic flux. bioRxiv, 2022, doi: [10.1101/2022.09.20.506926](https://doi.org/10.1101/2022.09.20.506926)
+Millard P., Gosselin-Monplaisir T. Uttenweiler-Joseph S., Enjalbert B. From toxic waste to beneficial nutrient: acetate boosts *Escherichia coli* growth at low glycolytic flux. bioRxiv, 2022, doi: [10.1101/2022.09.20.506926](https://doi.org/10.1101/2022.09.20.506926)
 
 ## Authors
 Pierre Millard

data/data_E_coli_K12_BW25113_Figure_6-7.txt

@@ -71,3 +71,15 @@ E_coli_K12_BW25113	delta_pta	galactose	13	0	0	0.2643
 E_coli_K12_BW25113	delta_pta	galactose	13	0	10	0.2369
 E_coli_K12_BW25113	delta_pta	galactose	14	0	0	0.2556
 E_coli_K12_BW25113	delta_pta	galactose	14	0	10	0.2389
+E_coli_K12_BW25113	delta_pgi	glucose	19	0	0	0.2117
+E_coli_K12_BW25113	delta_pgi	glucose	19	0	10	0.2925
+E_coli_K12_BW25113	delta_pfkA	glucose	19	0	0	0.0739
+E_coli_K12_BW25113	delta_pfkA	glucose	19	0	10	0.143
+E_coli_K12_BW25113	delta_pgi	glucose	20	0	0	0.1793
+E_coli_K12_BW25113	delta_pgi	glucose	20	0	10	0.2991
+E_coli_K12_BW25113	delta_pfkA	glucose	20	0	0	0.0706
+E_coli_K12_BW25113	delta_pfkA	glucose	20	0	10	0.1658
+E_coli_K12_BW25113	delta_pgi	glucose	21	0	0	0.1816
+E_coli_K12_BW25113	delta_pgi	glucose	21	0	10	0.3156
+E_coli_K12_BW25113	delta_pfkA	glucose	21	0	0	0.0758
+E_coli_K12_BW25113	delta_pfkA	glucose	21	0	10	0.1691

data/data_E_coli_K12_BW25113_Figure_EV1.txt

@@ -0,0 +1,25 @@
+strain	genotype	carbon_source	experiment_number	aMG	ACE	delta_pH
+E_coli_K12_BW25113	wt	glucose	1	0	0	-0.23
+E_coli_K12_BW25113	wt	glucose	1	100	0	-0.31
+E_coli_K12_BW25113	wt	glucose	1	100	10	-0.21
+E_coli_K12_BW25113	wt	glucose	1	100	100	-0.14
+E_coli_K12_BW25113	wt	glucose	2	0	0	-0.35
+E_coli_K12_BW25113	wt	glucose	2	100	0	-0.25
+E_coli_K12_BW25113	wt	glucose	2	100	10	-0.19
+E_coli_K12_BW25113	wt	glucose	2	100	100	-0.14
+E_coli_K12_BW25113	wt	glucose	3	0	0	-0.27
+E_coli_K12_BW25113	wt	glucose	3	100	0	-0.24
+E_coli_K12_BW25113	wt	glucose	3	100	10	-0.15
+E_coli_K12_BW25113	wt	glucose	3	100	100	-0.12
+E_coli_K12_BW25113	wt	glycerol	4	0	0	-0.3
+E_coli_K12_BW25113	wt	glycerol	4	0	10	-0.22
+E_coli_K12_BW25113	wt	galactose	5	0	0	-0.25
+E_coli_K12_BW25113	wt	galactose	5	0	10	-0.17
+E_coli_K12_BW25113	wt	glycerol	6	0	0	-0.21
+E_coli_K12_BW25113	wt	glycerol	6	0	10	-0.15
+E_coli_K12_BW25113	wt	galactose	7	0	0	-0.17
+E_coli_K12_BW25113	wt	galactose	7	0	10	-0.14
+E_coli_K12_BW25113	wt	glycerol	8	0	0	-0.16
+E_coli_K12_BW25113	wt	glycerol	8	0	10	-0.14
+E_coli_K12_BW25113	wt	galactose	9	0	0	-0.15
+E_coli_K12_BW25113	wt	galactose	9	0	10	-0.12

data/data_E_coli_K12_BW25113_Figure_EV2.txt

@@ -0,0 +1,30 @@
+genotype	carbon_source	experiment_number	aMG	ACETIC_ACID	mu	delta_pH
+E_coli_K12_BW25113	glucose	1	0	0	0.6283	-0.25
+E_coli_K12_BW25113	glucose	1	100	0	0.3905	-0.22
+E_coli_K12_BW25113	glucose	1	100	10	0.439	-0.18
+E_coli_K12_BW25113	glycerol	1	0	0	0.3807	-0.21
+E_coli_K12_BW25113	glycerol	1	0	10	0.4168	-0.1
+E_coli_K12_BW25113	galactose	1	0	0	0.1509	-0.11
+E_coli_K12_BW25113	galactose	1	0	10	0.3063	0
+E_coli_K12_BW25113	glucose	2	0	0	0.59	-0.26
+E_coli_K12_BW25113	glucose	2	0	10	0.5678	-0.19
+E_coli_K12_BW25113	glucose	3	0	0	0.6029	-0.27
+E_coli_K12_BW25113	glucose	3	0	10	0.5793	-0.15
+E_coli_K12_BW25113	glucose	4	0	0	0.6038	-0.29
+E_coli_K12_BW25113	glucose	4	0	10	0.5786	-0.18
+E_coli_K12_BW25113	glucose	5	0	0	0.6234	-0.31
+E_coli_K12_BW25113	glucose	5	100	0	0.4272	-0.38
+E_coli_K12_BW25113	glucose	5	100	10	0.4555	-0.23
+E_coli_K12_BW25113	glucose	6	0	0	0.6238	-0.27
+E_coli_K12_BW25113	glucose	6	100	0	0.3949	-0.21
+E_coli_K12_BW25113	glucose	6	100	10	0.4264	-0.16
+E_coli_K12_BW25113	glycerol	7	0	0	0.3856	-0.19
+E_coli_K12_BW25113	glycerol	7	0	10	0.4065	-0.15
+E_coli_K12_BW25113	galactose	8	0	0	0.1897	-0.08
+E_coli_K12_BW25113	galactose	8	0	10	0.3186	-0.04
+E_coli_K12_BW25113	glucose	9	0	0	0.632250927	-0.33
+E_coli_K12_BW25113	glucose	9	100	0	0.423267115	-0.23
+E_coli_K12_BW25113	glycerol	10	0	0	0.389361548	-0.29
+E_coli_K12_BW25113	glycerol	10	0	10	0.402106384	-0.23
+E_coli_K12_BW25113	galactose	11	0	0	0.183135751	-0.16
+E_coli_K12_BW25113	galactose	11	0	10	0.318044338	0.03

misc.R

@@ -3,10 +3,10 @@
 # load dependencies
 
 library(RColorBrewer)
-library(gplots)
 library(CoRC)
 library(pso)
 library(optimx)
+library(colorspace)
 
 
 # polynomial fit
@@ -28,7 +28,7 @@ cost <- function(par, x, ace, func){
   if (anyNA(xpred)){
     xpred <- rep(1e10, length(x))
   }
-  res <- (xpred - x)**2
+  res <- ((xpred - x)/0.1)**2
   return(sum(res))
 }
 
@@ -39,13 +39,17 @@ sim_sigmoid <- function(par, ace){
   return(xpred)
 }
 
-fit_sigmoid <- function(x_val, y_val, par_ini_weibull, lower_sigmoid, upper_sigmoid){
+fit_sigmoid <- function(x_val, y_val, par_ini_sigmoid, lower_sigmoid, upper_sigmoid){
   ord <- order(x_val)
   x <- y_val[ord]
   ace <- x_val[ord]
   ace <- ace[!is.na(x)]
   x <- x[!is.na(x)]
-  res_fit <- psoptim(par=par_ini_sigmoid, fn=cost, x=x, ace=ace, func=sim_sigmoid, lower=lower_sigmoid, upper=upper_sigmoid, control = list(maxit=100000, abstol=1e-6, trace.stats=TRUE))
+  res_fit <- psoptim(par=par_ini_sigmoid, fn=cost, x=x, ace=ace, func=sim_sigmoid, lower=lower_sigmoid, upper=upper_sigmoid, control = list(maxit=20000, abstol=1e-6, trace.stats=TRUE))
+  
+  res_stats <- optim(par=res_fit$par, fn=cost, x=x, ace=ace, func=sim_sigmoid, method="L-BFGS-B", hessian = TRUE, lower=lower_sigmoid, upper=upper_sigmoid, control = list(maxit=100))
+  res_stats$sd <- sqrt(diag(solve(res_stats$hessian)))
+  res_fit$res_stats <- res_stats
   return(res_fit)
 }
 
@@ -199,13 +203,56 @@ plot_with_ci_2 <- function(x1, y1, y2, x2, y3, sd_y3, col, h=NULL, ...){
   plot_points(x2, y3, sd_y3, offset=0.002, col=col)
 }
 
-plot_with_ci_3 <- function(sim_results, x, specie, col, ...){
-  plot(sim_results[1,,x], apply(sim_results[,,specie], 2, mean), col=col, type="l", ...)
+plot_with_ci_3 <- function(sim_results, x, specie, col, add_to_plot=FALSE, ...){
+  if (add_to_plot){
+    lines(sim_results[1,,x], apply(sim_results[,,specie], 2, mean), col=col, ...)
+  }else{
+    plot(sim_results[1,,x], apply(sim_results[,,specie], 2, mean), col=col, type="l", ...)
+  }
   polygon(x=c(sim_results[1,,x], rev(sim_results[1,,x])),
           y=c(apply(sim_results[,,specie], 2, max), rev(apply(sim_results[,,specie], 2, min))),
           col=paste(col, "55", sep=""), border=NA)
 }
 
+plot_with_ci_4 <- function(sim_results, x, specie, col, add_to_plot=FALSE, ...){
+  if (add_to_plot){
+    lines(sim_results[1,,x], apply(sim_results[,,specie], 2, mean), col=col, ...)
+  }else{
+    plot(sim_results[1,,x], apply(sim_results[,,specie], 2, mean), col=col, type="l", ...)
+  }
+  polygon(x=c(sim_results[1,,x], rev(sim_results[1,,x])),
+          y=c(apply(sim_results[,,specie], 2, mean)+apply(sim_results[,,specie], 2, sd), rev(apply(sim_results[,,specie], 2, mean)-apply(sim_results[,,specie], 2, sd))),
+          col=paste(col, "55", sep=""), border=NA)
+}
+
+plot_with_ci_5 <- function(simulation_results, n, x, specie, col, add_to_plot=FALSE, ...){
+  if (add_to_plot){
+    lines(x, rev(apply(simulation_results[,,n,specie], 2, mean)), col=col, ...)
+  }else{
+    plot(x, rev(apply(simulation_results[,,n,specie], 2, mean)), col=col, type="l", axes=FALSE, ...)
+    axis(1, at = seq(20,100,by=20), labels = paste(seq(100,20, by=-20)), tick = TRUE)
+    axis(2, tick = TRUE, las=1)
+    box()
+  }
+  polygon(x=c(x, rev(x)),
+          y=c(rev(apply(simulation_results[,,n,specie], 2, mean))+rev(apply(simulation_results[,,n,specie], 2, sd)), rev(rev(apply(simulation_results[,,n,specie], 2, mean))-rev(apply(simulation_results[,,n,specie], 2, sd)))),
+          col=paste(col, "55", sep=""), border=NA)
+}
+
+#plot_with_ci_5 <- function(simulation_results, n, x, specie, col, add_to_plot=FALSE, ...){
+#  if (add_to_plot){
+#    lines(x, rev(apply(simulation_results[,,n,specie], 2, mean)), col=col, ...)
+#  }else{
+#    plot(x, rev(apply(simulation_results[,,n,specie], 2, mean)), col=col, type="l", axes=FALSE, ...)
+#    axis(1, at = seq(20,100,by=20), labels = paste(seq(100,20, by=-20)), tick = TRUE)
+#    axis(2, tick = TRUE, las=1)
+#    box()
+#  }
+#  polygon(x=c(x, rev(x)),
+#          y=c(rev(apply(simulation_results[,,n,specie], 2, min)), rev(rev(apply(simulation_results[,,n,specie], 2, max)))),
+#          col=paste(col, "55", sep=""), border=NA)
+#}
+
 plot_no_ci <- function(fit_results, cond, specie, col, ...){
   if (specie %in% dimnames(fit_results[[cond]]$simulations)$specie){
     specie_id <- specie
@@ -231,3 +278,8 @@ error.bar <- function(x, y, upper, lower=upper, length=0.05,...){
   arrows(x,y+upper, x, y-lower, angle=90, code=3, length=length, ...)
 }
 
+get_ace_threshold <- function(ace_flux, ace_conc){
+  idx <- which.min(abs(ace_flux))
+  return(ace_conc[idx])
+}
+