extended adaptation of models to growth conditions

- define carbon sources with should/shouldn't enable growth of a model
- example enabling growth with glucose and disabling growth with acetate:
  - ./gapseq adapt -m toy/myb71.RDS -w cpd00027:TRUE,cpd00029:F -c toy/myb71-rxnWeights.RDS -g toy/myb71-rxnXgenes.RDS -b toy/myb71-all-Reactions.tbl
- useful to integrate growth data coming for example from biolog
- fixes bug in metacyc EC number matching (refers #101)

src/adapt.R

@@ -5,19 +5,29 @@ spec <- matrix(c(
   'model', 'm', 1, "character", "GapSeq model to be extended (RDS or SBML)",
   'help' , 'h', 0, "logical", "help",
   'add' , 'a', 2, "character", "reactions or pathways that should be added to the model (comma-separated)",
-  'remove' , 'r', 2, "character", "reactions or pathways that should be removed from the model (comma-separated)"
+  'remove' , 'r', 2, "character", "reactions or pathways that should be removed from the model (comma-separated)",
+  'growth' , 'w', 2, "character", "compounds for which growth or non-growth should be ensured (cpd00027:TRUE,cpd00017:FALSE)",
+ '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."
+
 ), ncol = 5, byrow = T)
 
 opt <- getopt(spec)
 
 # Help Screen
-if ( !is.null(opt$help) | is.null(opt$model)){
+if ( !is.null(opt$help) || is.null(opt$model)){
   cat(getopt(spec, usage=TRUE))
   q(status=1)
 }
 
-if( is.null(opt$add) & is.null(opt$remove) ){
-  print("Please specify reactions or pathway that should be added (-a) or removed (-r)")
+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)")
+  q(status=1)
+}
+
+if(!is.null(opt$growth) && (is.null(opt$rxn.weights.file) || is.null(opt$rxnXgene.table) || is.null(opt$rxn.blast.file))){
+  print("For adapting growth for substances please provide the reaction weights table, gene-reaction table and the blast-results table.")
   q(status=1)
 }
 
@@ -60,6 +70,10 @@ if ( is.null(opt$model) ) { opt$model =  paste0(script.dir,"/../toy/myb71.RDS")}
 mod.file            <- opt$model
 ids.add             <- opt$add
 ids.remove          <- opt$remove
+sub.growth          <- opt$growth
+rxn.weights.file    <- opt$rxn.weights.file
+rxnXgene.table      <- opt$rxnXgene.table
+rxn.blast.file      <- opt$rxn.blast.file
 
 # Parameters:
 sbml.export  <- FALSE 
@@ -78,7 +92,7 @@ meta.rxn <- fread(paste0(script.dir, "/../dat/meta_rea.tbl"))
 seed_x_mets   <- fread(paste0(script.dir,"/../dat/seed_metabolites_edited.tsv"), header=T, stringsAsFactors = F, na.strings = c("null","","NA"))
 
 cat("Loading model files", mod.file, "\n")
-mod        <- construct_full_model(script.dir)
+fullmod        <- construct_full_model(script.dir)
 if ( toupper(file_ext(mod.file)) == "RDS" ){
   mod.orig <- readRDS(mod.file)
 }else{ 
@@ -114,7 +128,7 @@ ids2seed <- function(ids){
       ids2seed.dt <- rbind(ids2seed.dt, data.table(id=ids[i], id.type="EC number", db.rxn="", seed=rxn.str))
     }else if ( !is.na(idx.rxn[i]) ){
       rxn <- gsub("\\|","",meta.rxn[idx.rxn[i], id])
-      ec  <- meta.rxn[idx.rxn[i], ec]
+      ec  <- str_remove(meta.rxn[idx.rxn[i], ec], "^EC-")
       rxn.name <- meta.rxn[idx.rxn[i], name]
       rxn.str <- system(paste0(script.dir, "/getDBhit.sh ", paste(rxn, paste0("'",rxn.name,"'"), ec, "seed")), intern=T)
       ids2seed.dt <- rbind(ids2seed.dt, data.table(id=ids[i], id.type="metacyc rxn", db.rxn=ids[i], seed=rxn.str))
@@ -150,13 +164,13 @@ if( !is.null(ids.add) ){
   # remove reaction which are already in model
   rxn.new <- setdiff(rxn.add, str_remove(mod.orig@react_id, "_[a-z]0$"))
   print(paste("Reactions already in model: ", paste0(setdiff(rxn.add, rxn.new), collapse = ",")))
-  if( length(rxn.new)==0 ) stop("No reactions left.")
+  if( length(rxn.new)==0 ) stop("No reactions can be added.")
   rxn.add <- rxn.new
 
   # remove reacitons which are not in full model
-  rxn.avail <- intersect(rxn.add, str_remove(mod@react_id, "_[a-z]0$"))
+  rxn.avail <- intersect(rxn.add, str_remove(fullmod@react_id, "_[a-z]0$"))
   print(paste("Reactions not in gapseq reaction database: ", paste0(setdiff(rxn.add, rxn.avail), collapse = ",")))
-  if( length(rxn.avail)==0 ) stop("No reactions left.")
+  if( length(rxn.avail)==0 ) stop("No reactions can be added.")
   rxn.add <- rxn.avail
 
   mod.out <- add_reaction_from_db(mod.orig, react = rxn.add, gs.origin = 10)
@@ -188,6 +202,44 @@ if( !is.null(ids.remove) ){
   print(paste("Removed reactions: ", paste0(rxn.remove, collapse = ",")))
 }
 
+
+if( !is.null(sub.growth) ){
+    print("adapt model to growth table")
+    growth.dt <- data.table()
+    for(check_input in unlist(str_split(sub.growth, ","), recursive=F)){
+        test <- unlist(str_split(check_input, ":"))
+        if(length(test)==2 && all(!is.na(str_match(test[1], "cpd[0-9]{5}"))) && !is.na(as.logical(test[2]))){
+            growth.dt <- rbind(growth.dt, data.table(sub=test[1], growth=as.logical(test[2])))
+        }else{
+            print(paste("Error in growth state of compound:", check_input))
+        }
+    }
+    if(nrow(growth.dt)==0) quit()
+
+    #source(paste0(script.dir, "/construct_full_model.R"))
+    #fullmod <- construct_full_model(script.dir)
+    source(paste0(script.dir, "/gf.adapt.R"))
+
+    print(growth.dt)
+    mod.out <- mod.orig
+    for(i in 1:nrow(growth.dt)){
+        sub.id <- growth.dt[i, sub]
+        ex.id <- paste0("EX_", sub.id, "_e0")
+        growth <- growth.dt[i, growth]
+        if(growth){
+            if( !ex.id %in% mod.out@react_id) mod.out <- addReact(mod.out, id=ex.id, met=paste0(sub.id,"[e0]"),Scoef=-1,lb=0,ub=1000)
+            mod.out <- add_growth(mod.out, add.met.id=sub.id, weights=rxn.weights.file, genes=rxnXgene.table)
+        }else{
+            if( !ex.id %in% mod.out@react_id){
+                print(paste("Model is already not growing with", sub.id))
+                next
+            }
+            mod.out <- rm_growth(mod.out, del.met.id=sub.id, rxn.blast.file=rxn.blast.file)
+        }
+        #print(paste(sub.id, growth))
+    }
+}
+
 out.id <- gsub(".xml|.RDS|.rds","",gsub("-draft","",basename(mod.file)))
 out.rds <- paste0("./",out.id,"-adapt",".RDS")
 saveRDS(mod.out, file = out.rds)