added way to adapt growth rate

- 'gapseq fill' and 'gapseq adapt' allow to set minimum growth rate now
- use '-k' parameter to set growth rate (default is 0.05)
- gap filling will try to reach at least the given growth rate

src/adapt.R

@@ -10,7 +10,8 @@ spec <- matrix(c(
   'rxn.weights.file', 'c', 2, "character", "Reaction weights table generated by gapseq function \"generate_GSdraft.R\" (RDS format; needed for 'growth' option).",
   'rxnXgene.table','g', 2, "character", "Table with gene-X-reaction associations as generated by the \"generate_GSdraft.R\" (RDS format; needed for 'growth; option)",
   'rxn.blast.file', 'b', 1, "character", "Blast-results table generated by gapseq find.",
-  'output.dir', '-f', 2, "character", "Path to directory, where output files will be saved (default: current directory)",
+  'output.dir', 'f', 2, "character", "Path to directory, where output files will be saved (default: current directory)",
+  'min.growth', 'k', 2, "numeric", "Set minimum growth rate that should be achieved by gap-filling.",
   'verbose', 'v', 2, "logical", "Verbose output and printing of debug messages. Default: FALSE"
 ), ncol = 5, byrow = T)
 
@@ -22,8 +23,8 @@ if ( !is.null(opt$help) || is.null(opt$model)){
   q(status=1)
 }
 
-if( is.null(opt$add) && is.null(opt$remove) && is.null(opt$growth) ){
-  print("Please specify reactions or pathway that should be added (-a) or removed (-r) or list growth compounds (-w)")
+if( is.null(opt$add) && is.null(opt$remove) && is.null(opt$growth) && is.null(opt$min.growth) ){
+  print("Please specify reactions or pathway that should be added (-a) or removed (-r) or list growth compounds (-w) or set new growth rate (-k)")
   q(status=1)
 }
 
@@ -32,6 +33,11 @@ if(!is.null(opt$growth) && (is.null(opt$rxn.weights.file) || is.null(opt$rxnXgen
   q(status=1)
 }
 
+if(!is.null(opt$min.growth) && (is.null(opt$rxn.weights.file) || is.null(opt$rxnXgene.table))){
+  print("For adapting growth rate please provide the reaction weights table, gene-reaction table.")
+  q(status=1)
+}
+
 # get current script path
 if (!is.na(Sys.getenv("RSTUDIO", unset = NA))) {
   # RStudio specific code
@@ -79,6 +85,7 @@ rxnXgene.table      <- opt$rxnXgene.table
 rxn.blast.file      <- opt$rxn.blast.file
 output.dir          <- opt$output.dir
 verbose             <- opt$verbose
+min.growth          <- opt$min.growth
 
 # Parameters:
 sbml.export  <- FALSE 
@@ -250,6 +257,12 @@ if( !is.null(sub.growth) ){
     }
 }
 
+if( !is.null(min.growth) ){
+    source(paste0(script.dir, "/gf.adapt.R"))
+    mod.out <- increase_growth(mod.orig, min.obj.val=min.growth, weights=rxn.weights.file, genes=rxnXgene.table, fullmod=fullmod, verbose=verbose)
+}
+
+
 if(mod.out@react_num == mod.orig@react_num && all(mod.out@react_id == mod.orig@react_id)){
     warning("No changes made")
     quit()

src/gf.adapt.R

@@ -175,7 +175,7 @@ add_growth <- function(model.orig, add.met.id=NA, weights=NA, genes=NA, verbose=
                               script.dir = script.dir,
                               core.only = only.core,
                               verbose=verbose,
-                              gs.origin = 3,
+                              gs.origin = 10,
                               rXg.tab = rXg.tab) 
     #))
     new.reactions <- mod.adapt.lst$rxns.added
@@ -348,3 +348,46 @@ addReactSrc <- function(mod, src, rxn.id){
                        lb=rxn.lb, ub=rxn.ub, obj=rxn.obj)
   return(model.new)
 }
+
+
+increase_growth <- function(model.orig, min.obj.val, weights=NA, genes=NA, verbose=F, only.core=F, fullmod){
+  # This here is needed if another draft than GapSeq's own draft networks are gapfilled
+  if((!"gs.origin" %in% colnames(model.orig@react_attr))) {
+    model.orig@react_attr <- data.frame(seed      = gsub("_.0","",model.orig@react_id),
+                                        gs.origin = 0,
+                                        stringsAsFactors = F)
+  }
+  if( !is.na(weights) ) rxn.weights <- readRDS(weights)
+  if( !is.na(genes) )   rXg.tab     <- readRDS(genes)
+
+  sol <- optimizeProb(model.orig, retOptSol=F)
+  cat("Old growth rate:", round(sol$obj,6), "\n")
+  if(sol$stat == ok & sol$obj >= min.obj.val){
+    warning(paste("Model already has growth rate >=",min.obj.val))
+    return(invisible(model.orig))
+  }else{
+    cat("\nTry to increase growth", "\n")
+    mod.adapt.lst <- gapfill4(mod.orig = model.orig, 
+                              mod.full = fullmod, 
+                              rxn.weights = copy(rxn.weights),  
+                              min.gr = min.obj.val,
+                              bcore = bcore,
+                              dummy.weight = dummy.weight,
+                              script.dir = script.dir,
+                              core.only = only.core,
+                              verbose=verbose,
+                              gs.origin = 10,
+                              rXg.tab = rXg.tab) 
+
+    new.reactions <- mod.adapt.lst$rxns.added
+    if( length(new.reactions) > 0 ){
+      #if( verbose ) cat("Added reactions:", new.reactions, "\n")
+      mod.adapt <- mod.adapt.lst$model
+    }
+  }
+  sol <- optimizeProb(mod.adapt, retOptSol=F)
+  cat("New growth rate:", round(sol$obj,6), "\n")
+  rxn.added <- setdiff(mod.adapt@react_id, model.orig@react_id)
+  printReaction(mod.adapt, react=rxn.added)
+  return(invisible( mod.adapt ))
+}

src/gf.suite.R

@@ -18,7 +18,8 @@ spec <- matrix(c(
   'verbose', 'v', 2, "logical", "Verbose output and printing of debug messages. Default: FALSE",
   'relaxed.constraints', 'r', 2, "logical", "Save final model as unconstraint network (i.e. all exchange reactions are open). Default: FALSE",
   'environment', 'e', 2, "character", "Adjusting reaction directions according to specific environmental conditions. See documentation for details. CAUTION: experimental option!",
-  'write.cs.ferm', 'w', 2, "logical", "Write a list with found carbon sources and fermentation products"
+  'write.cs.ferm', 'w', 2, "logical", "Write a list with found carbon sources and fermentation products",
+  'min.obj.val', 'k', 2, "numeric", "Minimum growth rate that should be achieved by gap-filling. Default: 0.05"
 ), ncol = 5, byrow = T)
 
 opt <- getopt(spec)
@@ -70,6 +71,7 @@ if ( is.null(opt$no.core) ) { opt$no.core = F }
 if ( is.null(opt$relaxed.constraints) ) { opt$relaxed.constraints = F }
 if ( is.null(opt$environment) ) { opt$environment = "" }
 if ( is.null(opt$write.cs.ferm) ) { opt$write.cs.ferm = F }
+if ( is.null(opt$min.obj.val) ) { opt$min.obj.val = 0.05 }
 
 # deprecation notice for flag '-o'
 if(!is.null(opt$depr.output.dir)) {
@@ -93,11 +95,10 @@ no.core             <- opt$no.core
 relaxed.constraints <- opt$relaxed.constraints
 env                 <- opt$environment
 write.cs.ferm       <- opt$write.cs.ferm
-
+min.obj.val         <- opt$min.obj.val
 
 # Parameters:
 dummy.weight <- 100
-min.obj.val  <- 0.05
 sbml.export  <- FALSE 
 
 # create output directory if not already there