diff --git a/DESCRIPTION b/DESCRIPTION index 1053921..dcba558 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: BacArena Title: Modeling Framework for Cellular Communities in their Environments -Version: 1.6 +Version: 1.7 Authors@R: c( person("Eugen", "Bauer", , "eugen.bauer@uni.lu", role = c("aut")), person("Johannes", "Zimmermann", , "j.zimmermann@iem.uni-kiel.de", role = c("aut", "cre"))) Author: Eugen Bauer [aut], @@ -11,7 +11,7 @@ Description: Can be used for simulation of organisms living in metabolic models determine the uptake and release of compounds. Biological processes such as movement, diffusion, chemotaxis and kinetics are available along with data analysis techniques. -URL: https://github.com/euba/BacArena +URL: https://BacArena.github.io/ BugReports: https://github.com/euba/BacArena/issues Depends: R (>= 3.0.0), @@ -31,7 +31,6 @@ Imports: plyr, Rcpp Suggests: - sybilSBML, parallel, knitr, rmarkdown diff --git a/NEWS b/NEWS index a8c637e..7ac9d19 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,17 @@ +BacArena v.1.7 (Release date: 2018-03-21) +============== +Changes: +- fine-tuning on growth model + * several duplications per time step + * fixed no growth bug for very high accumulated biomass +- improved handling of non-biomass objectives +- better status printing +- more then one compound in chemotaxis +- wrong unit of diffusion constant +- several small updated in plotting methods + + + BacArena v.1.6 (Release date: 2017-03-21) ============== Changes: diff --git a/R/Arena.R b/R/Arena.R index 0055d0e..ed1fa86 100644 --- a/R/Arena.R +++ b/R/Arena.R @@ -337,13 +337,13 @@ setMethod("addSubs", "Arena", function(object, smax=0, mediac=object@mediac, dif if( !template ) diffmat <- conv * diffmat else - smax <- smax * (arena@n * arena@m) / sum(diffmat!=0) + smax <- smax * (object@n * object@m) / sum(diffmat!=0) } if( is(diffmat, "list") ){ if( !template ) diffmat <- lapply(diffmat, function(x){conv*x}) else - smax <- sapply(seq_along(smax), function(i){smax[i] * (arena@n * arena@m) / sum(diffmat[[i]]) }) + smax <- sapply(seq_along(smax), function(i){smax[i] * (object@n * object@m) / sum(diffmat[[i]]) }) } @@ -497,7 +497,7 @@ setMethod("addEssentialMed", "Arena", function(object, org, only_return=FALSE, l min_val <- -ex@lowbnd[which(ex_max<0)] # get lower bounds also for -INF cases min_val[min_val==Inf] <- 1000 - predefined <- sapply(arena@media, function(x){sum(x@diffmat)}) # only use compounds which are not already in present in the media + predefined <- sapply(object@media, function(x){sum(x@diffmat)}) # only use compounds which are not already in present in the media idx <- which(min_id %in% names(predefined)[which(predefined==0)]) if(length(idx)==0) return(object) @@ -1473,6 +1473,7 @@ setMethod(show, "Arena", function(object){ #' @slot mfluxlist A list of containing highly used metabolic reactions per time step. #' @slot shadowlist A list of containing shadow prices per time step. #' @slot subchange A vector of all substrates with numbers indicating the degree of change in the overall simulation. +#' @slot exchangeslist A list of containing exchanges per time step. setClass("Eval", contains="Arena", representation( @@ -2850,6 +2851,7 @@ setMethod("checkCorr", "Eval", function(object, corr=NULL, tocheck=list()){ #' @param spec_nr Number of the specie #' @param sub_nr Maximal number of substances to be show #' @param cutoff Shadow costs should be smaller than cutoff +#' @param noplot Do not plot #' @details Returns ggplot objects setGeneric("plotShadowCost", function(object, spec_nr=1, sub_nr=10, cutoff=-1, noplot=FALSE){standardGeneric("plotShadowCost")}) #' @export @@ -3089,46 +3091,42 @@ setMethod("plotSubDist2", "Eval", function(object, sub, times=NULL){ #' @return Heatmap (ggplot2) #' #' @examples -#' Do Not Run -#' sim : simulation object -#' speciesA = 1 # E. Coli (wild): sequence number in names(sim@specs) => 1 -#' speciesB = 3 # E. Coli (mutant): sequence number in names(sim@specs) => 3 -#' var_nr = 30 # (manually selected) -#' HeatMapFeeding(object = sim, speciesA = 1, speciesB = 3, var_nr = 30) +#' sim <- sihumi_test +#' HeatMapFeeding(object = sim, speciesA = 1, speciesB = 2, var_nr = 78) #' setGeneric("HeatMapFeeding", function(object, speciesA, speciesB, var_nr){standardGeneric("HeatMapFeeding")}) #' @export #' @rdname HeatMapFeeding setMethod("HeatMapFeeding", "Eval", function(object, speciesA, speciesB, var_nr){ - A<-BacArena::plotSpecActivity(object,spec_list = speciesA, var_nr = var_nr, ret_data = T) - z<-data.frame() - chronos <- (seq_along(object@simlist)-1) - for (t in chronos) - { a<-data.frame() - a <-BacArena::findFeeding3(object, time = t, mets = A$sub[1:var_nr], plot = F) - if (nrow(a)!=0) - {z<-rbind.data.frame(z,a) - } + A <- BacArena::plotSpecActivity(object,spec_list = speciesA, var_nr = var_nr, ret_data = T) + z <- data.frame() + chronos <- 1:(length(object@simlist)-1) + for (t in chronos) { + a<-data.frame() + a <-BacArena::findFeeding3(object, time = t, mets = A$sub[1:var_nr], plot = F) + if (nrow(a)!=0) + z <- rbind.data.frame(z,a) } d <- data.frame() # prod:speciesA cons:speciesB -> 1 # prod:speciesB cons:speciesA -> -1 # NA -> 0 - for (h in 1:nrow(z)) - {if (z[h,1] == names(object@specs)[speciesA]) (d[h,1]=1) - if (z[h,1] == names(object@specs)[speciesB]) (d[h,1]=-1) + for (h in 1:nrow(z)){ + if (z[h,1] == names(object@specs)[speciesA]) + d[h,1] <- 1 + else if (z[h,1] == names(object@specs)[speciesB]) + d[h,1] <- -1 + else + d[h,1] <- NA } colnames(d)[1] <- "status" q <- cbind(z,d) - p <- reshape::cast(q, met ~ sim_step, value = "status") - p[is.na(p)] <- 0 - # GGPLOT METHOD based on https://stackoverflow.com/questions/8406394/how-to-produce-a-heatmap-with-ggplot2 - p2 <- p - p.m <- reshape::melt(p2) - p3 <- ggplot2::ggplot(p.m, ggplot2::aes(sim_step,met)) + - ggplot2::geom_tile(ggplot2::aes(fill = value), colour = "white") + + q$status[is.na(q$status)] <- 0 + p3 <- ggplot2::ggplot(q, ggplot2::aes_string("sim_step","met")) + + ggplot2::geom_tile(ggplot2::aes_string(fill = "status"), colour = "white") + ggplot2::scale_fill_gradient2(name = "Producer", low = "red", mid = "green", high = "blue", breaks=seq(-1,1,by=1), labels = c(names(object@specs)[speciesB], "no feeding", names(object@specs)[speciesA])) + ggplot2::xlab("Simulation Step") + ggplot2::ylab("Exchange Reactions") - return((p3)) }) + return((p3)) +}) diff --git a/R/Organism.R b/R/Organism.R index f904003..5c31a60 100644 --- a/R/Organism.R +++ b/R/Organism.R @@ -120,7 +120,7 @@ Organism <- function(model, algo="fba", ex="EX_", ex_comp=NA, csuffix="\\[c\\]", if(length(pot_biomass)==0) stop("No objection function set in model") print("No objective function set, try set automatically...") print(paste("found:", pot_biomass)) - print(paste("set new biomass function:", pot_biomas[1])) + print(paste("set new biomass function:", pot_biomass[1])) sybil::obj_coef(model)[which(sybil::react_id(model)==pot_biomass[1])] <- 1 rbiomass <- pot_biomass[1] }else{ @@ -392,9 +392,10 @@ setMethod("optimizeLP", "Organism", function(object, lpob=object@lpobj, lb=objec shadow <- glpkAPI::getRowsDualGLPK(lpob@problem@oobj)[ex@met_pos] names(shadow) <- ex@met_id }else if(lpob@problem@solver=="cplexAPI" & solve_ok){ # TODO: cplex shadow costs needs update of optimzeProb call in optimizeLP() which is having side effects ... - ex <- findExchReact(object@model) - shadow <- cplexAPI::getPiCPLEX(lpob@problem@oobj@env, lpob@problem@oobj@lp, 0, lpob@nr-1)[ex@met_pos] - names(shadow) <- ex@met_id + warning("cplex shadow costs are not supported yet") + #ex <- findExchReact(object@model) + #shadow <- cplexAPI::getPiCPLEX(lpob@problem@oobj@env, lpob@problem@oobj@lp, 0, lpob@nr-1)[ex@met_pos] + #names(shadow) <- ex@met_id } names(fbasl$fluxes) <- names(object@lbnd) @@ -822,7 +823,7 @@ setMethod("growth_par", "Bac", function(object, population, j, fbasol, tstep){ #' arena <- Arena(n=20,m=20) #initialize the environment #' arena <- addOrg(arena,bac,amount=10) #add 10 organisms #' arena <- addSubs(arena,40) #add all possible substances -#' chemotaxis(bac,arena,1) +#' chemotaxis(bac,arena,1, "EX_glc(e)") setGeneric("chemotaxis", function(object, population, j, chemo){standardGeneric("chemotaxis")}) #' @export #' @rdname chemotaxis @@ -924,13 +925,14 @@ setMethod("simBac", "Bac", function(object, arena, j, sublb, bacnum, sec_obj="no #' @param lpobject linar programming object (copy of organism@lpobj) that have to be a deep copy in parallel due to pointer use in sybil. #' @param sec_obj character giving the secondary objective for a bi-level LP if wanted. #' @param cutoff value used to define numeric accuracy +#' @param with_shadow True if shadow cost should be stores (default off). #' @return Returns the updated enivironment of the \code{population} parameter with all new positions of individuals on the grid and all new substrate concentrations. #' -setGeneric("simBac_par", function(object, arena, j, sublb, bacnum, lpobject, sec_obj="none", cutoff=1e-6){standardGeneric("simBac_par")}) +setGeneric("simBac_par", function(object, arena, j, sublb, bacnum, lpobject, sec_obj="none", cutoff=1e-6, with_shadow=FALSE){standardGeneric("simBac_par")}) #' @export #' @rdname simBac_par -setMethod("simBac_par", "Bac", function(object, arena, j, sublb, bacnum, lpobject, sec_obj="none", cutoff=1e-6){ - const <- constrain(object, object@medium, lb=-sublb[j,object@medium]/bacnum, ub=object@ubnd*bacnum, #scale to population size +setMethod("simBac_par", "Bac", function(object, arena, j, sublb, bacnum, lpobject, sec_obj="none", cutoff=1e-6, with_shadow=FALSE){ + const <- constrain(object, object@medium, lb=-sublb[j,object@medium]/bacnum, #scale to population size dryweight=arena@orgdat[j,"biomass"], tstep=arena@tstep, scale=arena@scale) lobnd <- const[[1]]; upbnd <- const[[2]] fbasol <- optimizeLP(object, lb=lobnd, ub=upbnd, j=j, sec_obj=sec_obj, cutoff=cutoff, with_shadow=with_shadow) diff --git a/R/Stuff.R b/R/Stuff.R index 8b09391..070fd77 100644 --- a/R/Stuff.R +++ b/R/Stuff.R @@ -172,6 +172,7 @@ lsd <- function(y){lb=mean(y)-stats::sd(y); ifelse(lb<0,0,lb)} #' @param ret_data Set true if data should be returned #' @param num_var Number of varying substances to be shown (if mediac is not specified) #' @param unit Unit for the substances which should be used for plotting (default: mmol) +#' @param useNames Use substance names instead of ids #' #' @return list of three ggplot object for further formating #' diff --git a/man/Eval-class.Rd b/man/Eval-class.Rd index cfda467..17cd955 100644 --- a/man/Eval-class.Rd +++ b/man/Eval-class.Rd @@ -19,5 +19,7 @@ Structure of the S4 class \code{Eval} inheriting from class \code{\link{Arena-cl \item{\code{shadowlist}}{A list of containing shadow prices per time step.} \item{\code{subchange}}{A vector of all substrates with numbers indicating the degree of change in the overall simulation.} + +\item{\code{exchangeslist}}{A list of containing exchanges per time step.} }} diff --git a/man/HeatMapFeeding.Rd b/man/HeatMapFeeding.Rd index 091820b..8e6f449 100644 --- a/man/HeatMapFeeding.Rd +++ b/man/HeatMapFeeding.Rd @@ -27,11 +27,7 @@ The generic function \code{HeatMapFeeding} returns a heatmap between two selecte the producer. The compounds are originated from the most varying substances of the first cell model called "speciesA". } \examples{ - Do Not Run - sim : simulation object - speciesA = 1 # E. Coli (wild): sequence number in names(sim@specs) => 1 - speciesB = 3 # E. Coli (mutant): sequence number in names(sim@specs) => 3 - var_nr = 30 # (manually selected) - HeatMapFeeding(object = sim, speciesA = 1, speciesB = 3, var_nr = 30) +sim <- sihumi_test +HeatMapFeeding(object = sim, speciesA = 1, speciesB = 2, var_nr = 78) } diff --git a/man/chemotaxis.Rd b/man/chemotaxis.Rd index 5565e4a..d60d91e 100644 --- a/man/chemotaxis.Rd +++ b/man/chemotaxis.Rd @@ -32,7 +32,7 @@ bac <- Bac(Ec_core,deathrate=0.05, chem = "EX_o2(e)", arena <- Arena(n=20,m=20) #initialize the environment arena <- addOrg(arena,bac,amount=10) #add 10 organisms arena <- addSubs(arena,40) #add all possible substances -chemotaxis(bac,arena,1) +chemotaxis(bac,arena,1, "EX_glc(e)") } \seealso{ \code{\link{Bac-class}} and \code{\link{emptyHood}} diff --git a/man/plotShadowCost.Rd b/man/plotShadowCost.Rd index c8f9c91..b92a4e4 100644 --- a/man/plotShadowCost.Rd +++ b/man/plotShadowCost.Rd @@ -20,6 +20,8 @@ plotShadowCost(object, spec_nr = 1, sub_nr = 10, cutoff = -1, \item{sub_nr}{Maximal number of substances to be show} \item{cutoff}{Shadow costs should be smaller than cutoff} + +\item{noplot}{Do not plot} } \description{ The generic function \code{plotShadowCost} plots substances have the highest impact on further growth (shadow cost < 0) diff --git a/man/plotSubCurve.Rd b/man/plotSubCurve.Rd index 2a2e6ce..f3b34c9 100644 --- a/man/plotSubCurve.Rd +++ b/man/plotSubCurve.Rd @@ -21,6 +21,8 @@ plotSubCurve(simlist, mediac = NULL, time = c(NULL, NULL), scol = NULL, \item{ret_data}{Set true if data should be returned} \item{num_var}{Number of varying substances to be shown (if mediac is not specified)} + +\item{useNames}{Use substance names instead of ids} } \value{ list of three ggplot object for further formating diff --git a/man/simBac_par.Rd b/man/simBac_par.Rd index aa77165..db2bccc 100644 --- a/man/simBac_par.Rd +++ b/man/simBac_par.Rd @@ -7,10 +7,10 @@ \title{Function for one simulation iteration for objects of Bac class} \usage{ simBac_par(object, arena, j, sublb, bacnum, lpobject, sec_obj = "none", - cutoff = 1e-06) + cutoff = 1e-06, with_shadow = FALSE) \S4method{simBac_par}{Bac}(object, arena, j, sublb, bacnum, lpobject, - sec_obj = "none", cutoff = 1e-06) + sec_obj = "none", cutoff = 1e-06, with_shadow = FALSE) } \arguments{ \item{object}{An object of class Bac.} @@ -28,6 +28,8 @@ simBac_par(object, arena, j, sublb, bacnum, lpobject, sec_obj = "none", \item{sec_obj}{character giving the secondary objective for a bi-level LP if wanted.} \item{cutoff}{value used to define numeric accuracy} + +\item{with_shadow}{True if shadow cost should be stores (default off).} } \value{ Returns the updated enivironment of the \code{population} parameter with all new positions of individuals on the grid and all new substrate concentrations.