gf.suite.R

library(getopt)

# get options first
spec <- matrix(c(
  'model', 'm', 1, "character", "gapseq-Draft-Model to be gapfilled (RDS or SBML)",
  'help' , 'h', 0, "logical", "help",
  'media', 'n', 1, "character", "tab- or komma separated table for media components. Requires three named columns: 1 - \"compounds\" (for metab. IDs), 2 - \"name\" (metab. name), 3 - \"maxFlux\" (maximum inflow flux)",
  'target.metabolite', 't', 2, "character", "ID (without compartment suffix) of metabolite that shall be produced. Default: cpd11416 (Biomass)",
  'rxn.weights.file', 'c', 1, "character", "Reaction weights table generated by gapseq function \"generate_GSdraft.R\" (RDS format).",
  'rxnXgene.table','g', 1, "character", "Table with gene-X-reaction associations as generated by the \"generate_GSdraft.R\" (RDS format)",
  'bcore', 'b', 2, "numeric", "Minimum bitscore for reaction associated blast hits to consider reactions as core/candidate reactions. Default: 50",
  'output.dir', '-f', 2, "character", "Path to directory, where output files will be saved (default: current directory)",
  'depr.output.dir', 'o', 2, "character", "deprecated. Use flag\"-f\" instead",
  'sbml.no.output', 's', 2, "logical", "Do not save gapfilled model as sbml file. Default: Save as SBML",
  'quick.gf','q', 2, "logical", "perform only step 1 and 2. Default: FALSE",
  'limit', 'l', 2, "character", "Test metabolite to which search is limitted",
  'no.core', 'x', 2, "logical", "Use always all reactions instead of core reactions, which have sequence evidence. Default: FALSE",
  '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",
  '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)

# Help Screen
if ( !is.null(opt$help) | is.null(opt$model)){
  cat(getopt(spec, usage=TRUE))
  q(status=1)
}

# get current script path
if (!is.na(Sys.getenv("RSTUDIO", unset = NA))) {
    # RStudio specific code
    script.dir    <- dirname(rstudioapi::getSourceEditorContext()$path)
} else{
    initial.options <- commandArgs(trailingOnly = FALSE)
    script.name <- sub("--file=", "", initial.options[grep("--file=", initial.options)])
    script.dir  <- dirname(script.name)
}

if( "cplexAPI" %in% installed.packages() )
  suppressMessages(library(cplexAPI))
if( "sybilSBML" %in% installed.packages() )
  suppressMessages(library(sybilSBML))
suppressMessages(library(sybil))
suppressMessages(library(data.table)); setDTthreads(1)
suppressMessages(library(stringr))
suppressMessages(library(methods))
suppressMessages(library(tools))

# select solver
if( "cplexAPI" %in% rownames(installed.packages()) ){
  sybil::SYBIL_SETTINGS("SOLVER","cplexAPI"); ok <- 1
}else{
  warning("glpkAPI is used but cplexAPI is recommended because it is much faster")
  sybil::SYBIL_SETTINGS("SOLVER","glpkAPI"); ok <- 5
}

# Setting defaults if required
if ( is.null(opt$target.metabolite) ) { opt$target.metabolite = "cpd11416" }
if ( is.null(opt$output.dir) ) { opt$output.dir = "." }
if ( is.null(opt$sbml.no.output) ) { opt$sbml.no.output = F } else { opt$sbml.no.output = T }
#if ( is.null(opt$model) ) { opt$model = "" }
if ( is.null(opt$media) ) { opt$media = paste0(script.dir,"/../dat/media/MM_glu.csv") }
if ( is.null(opt$verbose) ) { opt$verbose = F }
if ( is.null(opt$quick.gf) ) { opt$quick.gf = F }
if ( is.null(opt$bcore) ) { opt$bcore = 50 }
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)) {
  warning("Deprecation notice: Flag '-o' is now replaced with flag '-f'. Please adjust your script(s), as the flag will be removed in a future release.")
  if(is.null(opt$output.dir)) 
    opt$output.dir <- opt$depr.output.dir
}
                        
# Arguments:
mod.file            <- opt$model
media.file          <- opt$media
target.met          <- opt$target.metabolite
rxn.weights.file    <- opt$rxn.weights.file
rxnXgene.table      <- opt$rxnXgene.table
output.dir          <- opt$output.dir
verbose             <- opt$verbose
quick.gf            <- opt$quick.gf
bcore               <- opt$bcore
met.limit           <- opt$limit
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
sbml.export  <- FALSE 

# create output directory if not already there
dir.create(output.dir, recursive = TRUE, showWarnings = FALSE)
if (!dir.exists(output.dir) || file.access(output.dir, mode = 2) == -1)
  stop(paste("Output directory",output.dir,"cannot be created or is not writable."))

# Get list of core-reactions from one or more files (given as comma seperated string by -c)
rxn.weights <- readRDS(rxn.weights.file)
rXg.tab     <- readRDS(rxnXgene.table)

# parsing environment specification string
env <- unlist(str_split(env, ","))
env <- env[env %in% c("","highH2")] # Filter for currently supported environment specifications only

# Little helpers
source(paste0(script.dir,"/add_missing_exRxns.R"))
source(paste0(script.dir,"/constrain.model.R"))
source(paste0(script.dir,"/gapfill4.R"))
source(paste0(script.dir,"/generate_rxn_stoich_hash.R"))
source(paste0(script.dir,"/get_gene_logic_string.R"))
source(paste0(script.dir,"/addMetAttr.R"))
source(paste0(script.dir,"/addReactAttr.R"))
source(paste0(script.dir,"/media_check.R"))
source(paste0(script.dir,"/adjust_model_env.R"))
source(paste0(script.dir,"/construct_full_model.R"))

rm.na <- function(vec){
  idx <- which(is.na(vec))
  if (length(idx) > 0) 
    return(vec[-idx])
  else return(vec)
}

if( no.core ){
  use.core = FALSE
}else{
  use.core = TRUE
}

# database files
carbon.source <- fread(paste0(script.dir, "/../dat/subex.tbl"))
seed_x_mets   <- fread(paste0(script.dir,"/../dat/seed_metabolites_edited.tsv"), header=T, stringsAsFactors = F, na.strings = c("null","","NA"))
seed_x_metCyc <- fread(paste0(script.dir,"/../dat/mnxref_seed-other.tsv"), header = T)

# potentially limit carbon.source
if ( length(met.limit) > 0 ){
  carbon.source <- carbon.source[str_extract(seed, "cpd[0-9]+") == met.limit | seed == met.limit | tolower(name) %like% tolower(met.limit)]
  #print(carbon.source)
  if( nrow(carbon.source)==0 ){
    stop("Limitation of carbon sources failed, nothing found!")
  }
}

# read full model & target model
cat("Loading model files", mod.file, "\n")
mod <- construct_full_model(script.dir)

if ( toupper(file_ext(mod.file)) == "RDS" ){
  mod.orig <- readRDS(mod.file)
}else{ 
  mod.orig <- readSBMLmod(mod.file)}

bu_mod_attr <- mod.orig@mod_attr

# Add annotation column to model attributes if not already there
if(!("annotation" %in% colnames(mod.orig@mod_attr))) {
  bm_ind <- which(mod.orig@react_id == "bio1")
  annostr <- ""
  if(grepl("Bacteria",mod.orig@react_name[bm_ind]))
    annostr <- "tax_domain:Bacteria"
  if(grepl("Archaea",mod.orig@react_name[bm_ind]))
    annostr <- "tax_domain:Archaea"
  
  mod.orig@mod_attr <- cbind(mod.orig@mod_attr,
                             data.frame(annotation = annostr))
}

# adjust environment if needed
if(env[1] != "")
  mod <- adjust_model_env(mod, env, script.dir)

# block reactions in domain, which do not have these reactions
if(any(grepl("tax_domain:", mod.orig@mod_attr[,"annotation"], fixed = T))) {
  domain.rxn.exclusions <- fread(paste0(script.dir, "/../dat/biomass/excluded_reactions.tsv"), header=T, stringsAsFactors = F)
  
  anno_ind   <- which(grepl("tax_domain:", mod.orig@mod_attr[,"annotation"], fixed = T))
  org.domain <- str_match(mod.orig@mod_attr[anno_ind,"annotation"],"tax_domain\\:\\s*(.*)\\s*")[2]
  
  domain.rxn.exclusions <- domain.rxn.exclusions[domain == org.domain, exclude.reaction]
  if(length(domain.rxn.exclusions) > 0) {
    excl_rxns <- paste(domain.rxn.exclusions, collapse = "|")
    
    block_rxns_ids_orig <- mod.orig@react_id[grep(excl_rxns, mod.orig@react_id)]
    if(length(block_rxns_ids_orig) > 0) {
      mod.orig <- changeBounds(mod.orig, react = block_rxns_ids_orig, 
                               lb = rep(0, length(block_rxns_ids_orig)),
                               ub = rep(0, length(block_rxns_ids_orig)))
    }
    
    block_rxns_ids_full <- mod@react_id[grep(excl_rxns, mod@react_id)]
    if(length(block_rxns_ids_full) > 0) {
      mod <- changeBounds(mod, react = block_rxns_ids_full, 
                          lb = rep(0, length(block_rxns_ids_full)),
                          ub = rep(0, length(block_rxns_ids_full)))
    }
  }
}

# This here is needed if another draft than gapseq's own draft networks are gapfilled
if((!"gs.origin" %in% colnames(mod.orig@react_attr))) {
  mod.orig@react_attr <- data.frame(seed      = gsub("_.0","",mod.orig@react_id),
                                    gs.origin = 0,
                                    stringsAsFactors = F)
}

mod.orig   <- add_missing_exchanges(mod.orig)

# add diffusion reactions
mod.orig       <- add_missing_diffusion(mod.orig)
#mod.orig       <- changeBounds(mod.orig, react="EX_cpd11640_e0", lb=0, ub=1) # TODO: Limit the hydrogen evolution rate. Check if it's necessary.

# create complete medium
cat("using media file", media.file, "\n")
if( media.file == "complete" ){
  met.pos <- findExchReact(mod.orig)@met_pos
  met.id  <- gsub("\\[.0\\]","",mod.orig@met_id[met.pos])
  met.name<- mod.orig@met_name[met.pos]
  media <- data.frame(compounds=met.id, name=met.name, maxFlux=100)
  media.file <- paste0(script.dir,"/../dat/media/ALLmed.csv")
  write.csv(media, media.file, quote = F, row.names = F)
}else{
  # perform plausible check of user-defined medium file provide warning if necessary
  media_check(media.file, mod.orig, seed_x_mets)
}


# constrain model
mod.orig <- constrain.model(mod.orig, media.file = media.file)
mod.orig@obj_coef <- rep(0,mod.orig@react_num)

# add metabolite objective + sink
mod.orig <- add_met_sink(mod.orig, target.met, obj = 1) # TODO: add gs.origin

# Perform gapfill
cat("\n\n1. Initial gapfilling: Make model grow on given media using all reactions\n")
mod.fill.lst <- gapfill4(mod.orig = mod.orig, 
                         mod.full = mod, 
                         rxn.weights = copy(rxn.weights), 
                         min.gr = min.obj.val,
                         bcore = bcore,
                         dummy.weight = dummy.weight,
                         script.dir = script.dir,
                         verbose = verbose,
                         gs.origin = 1,
                         rXg.tab = rXg.tab,
                         env = env)

mod.fill1 <- constrain.model(mod.fill.lst$model, media.file = media.file, scaling.fac = 1)
mod.out <- mod.fill1

# get currently present reactions
pres.rxns <- mod.orig@react_id
pres.rxns <- gsub("_.*","",pres.rxns)

# Check if anymore core reactions could be added
mseed.t <- fread(paste0(script.dir, "/../dat/seed_reactions_corrected.tsv"), header=T, stringsAsFactors = F)
mseed   <- copy(mseed.t)
mseed.t <- mseed.t[gapseq.status %in% c("approved","corrected")]
mseed.t <- mseed.t[!(id %in% pres.rxns)]
mseed.t <- mseed.t[id %in% rxn.weights[bitscore > bcore, seed]]

if(nrow(mseed.t)==0)
  cat("No more core reactions in list, that could be added to the model. Skipping gapfilling-steps 2,2b,3, and 4.\n")

core.plus.present.rxns <- c(rxn.weights[bitscore > bcore, seed], pres.rxns)

# A function that tests if the metabolite that is subject for gafillling steps 3 and 4 are part of 
# any reaction in the core reaction list or reactions, that are already included in the model. If not
# gapfilling is not possible and can be skipped.
checkIfGapfillIsPossible <- function(met) {
  # metabolite in [e0] compartment ?
  do.gf <- any(mseed[id %in% core.plus.present.rxns, grepl(gsub("\\[e0","\\[1",met),equation, fixed = T)])
  return(do.gf)
} 

if(nrow(mseed.t)>0) { # Skip steps 2,2b,3, and 4 if core-reaction list does not contain any new reactions.
  
  cat("\n\n2. Biomass gapfilling using core reactions only\n")
  
  mod.orig2 <- mod.out
  
  # load minimal medium and add available carbon sources
  media2 <- fread(paste0(script.dir,"/../dat/media/MM_glu.csv"))
  src.met <- carbon.source[group %in% c("Carbohydrates", "Polymers", "Carboxylic acids", "Amino acids") & seed %in% mod.orig2@react_id, .(seed,name,group)] 
  if( nrow(src.met) == 0){
    warnings("No carbon source exchange reactions found in model, considering all available.")
    src.met <- carbon.source[seed %in% mod.orig2@react_id, .(seed,name,group)]
  }
  # if glucose is not usable then add other carbon source(s)
  if( !any(grepl("alpha-D-Glucose", src.met$name)) ){
    src.carbo <- src.met[group=="Carbohydrates"]
    if( nrow(src.carbo)>0 )
      src.add <- src.carbo # if no glucose is there, then add all other available carbohydrates
    else
      src.add <- src.met # if no carbohydrates is avaiable, then take everything else (probably amino acid biosynthesis is not gapfilled because amino acids are part of the medium)
    media2 <- rbind(media2, data.table(compounds=gsub("^EX_|_e0$","",src.add$seed), name=src.add$name, maxFlux=100))
  }
  
  
  # constrain model  
  mod.orig2 <- constrain.model(mod.orig2, media = media2)
  mod.orig2@obj_coef <- rep(0,mod.orig2@react_num)
  
  bm.ind      <- which(mod.orig2@react_id == "bio1")
  bm.met.inds <- which(mod.orig2@S[,bm.ind]<0)
  bm.met      <- gsub("\\[.0\\]","",mod.orig2@met_id[bm.met.inds])
  bm.met.name <- mod.orig2@met_name[bm.met.inds]
  mod.fill2    <- mod.orig2
  mod.fill2.counter <- 0
  mod.fill2.names <- c()
  
  if( !verbose ) options(warn=-1)
  for( i in seq_along(bm.met) ){
    cat("\r",i,"/",length(bm.met))
    target.new <- bm.met[i]
    
    # add metabolite objective + sink
    rm.sink = TRUE
    if( paste0("EX_",target.new,"_c0") %in% react_id(mod.fill2) )
      rm.sink = FALSE
    mod.fill2  <- add_met_sink(mod.fill2, target.new, obj = 1)
    
    sol <- optimizeProb(mod.fill2, retOptSol=F)
    
    if(sol$stat == ok & sol$obj >= 1e-6){
      mod.fill2@obj_coef <- rep(0,mod.fill2@react_num)
    }else{
      if( verbose ) cat("\nTry to gapfill", bm.met.name[i],"\n")
      invisible(capture.output( 
        mod.fill2.lst <- gapfill4(mod.orig = mod.fill2, 
                                  mod.full = mod,
                                  rxn.weights = copy(rxn.weights), 
                                  min.gr = min.obj.val,
                                  bcore = bcore,
                                  dummy.weight = dummy.weight,
                                  script.dir = script.dir,
                                  core.only = use.core,
                                  verbose=verbose,
                                  gs.origin = 2,
                                  rXg.tab = rXg.tab,
                                  env = env) ))
      new.reactions <- mod.fill2.lst$rxns.added
      if( length(new.reactions) > 0 ){
        if( verbose ) cat("Added reactions:", new.reactions, "\n")
        mod.fill2 <- mod.fill2.lst$model
        mod.fill2.counter <- mod.fill2.counter + 1
        mod.fill2.names <- c(mod.fill2.names, bm.met.name[i])
      }
      mod.fill2@obj_coef <- rep(0,mod.fill2@react_num)
    }
    if( rm.sink )
      mod.fill2 <- rmReact(mod.fill2, react=paste0("EX_",target.new,"_c0"))
  }
  options(warn=0)
  
  mod.fill2 <- changeObjFunc(mod.fill2, react=paste0("EX_",target.met,"_c0"))
  mod.fill2 <- constrain.model(mod.fill2, media.file = media.file, scaling.fac = 1)
  mod.out <- mod.fill2
  
  cat("\rGapfill summary:\n")
  cat("Filled components:    ",mod.fill2.counter, "(",paste(mod.fill2.names, collapse = ","),")\n")
  cat("Added reactions:      ",length(mod.fill2@react_id)-length(mod.fill1@react_id),"\n")
  cat("Final growth rate:    ",optimizeProb(mod.fill2, retOptSol=F)$obj,"\n")
  
  
  
  media2 <- fread(paste0(script.dir,"/../dat/media/MM_glu.csv")) # load minimal medium and add available carbon sources
  if( nrow(fread(media.file, header=F)[V1=="cpd00007" & V3!=0]) ){
    cat("\n\n2b. Anaerobic biomass gapfilling using core reactions only\n")
    media2 <- media2[name!="O2"] # remove oxygen
  }else{
    cat("\n\n2b. Aerobic biomass gapfilling using core reactions only\n")
  }
  
  mod.orig2 <- mod.out
  
  src.met <- carbon.source[group %in% c("Carbohydrates", "Polymers", "Carboxylic acids", "Amino acids") & seed %in% mod.orig2@react_id, .(seed,name,group)]
  if( nrow(src.met) == 0){
    warnings("No carbon source exchange reactions found in model, considering all available.")
    src.met <- carbon.source[seed %in% mod.orig2@react_id, .(seed,name,group)]
  }
  # if glucose is not usable then add other carbon source(s)
  if( !any(grepl("alpha-D-Glucose", src.met$name)) ){
    src.carbo <- src.met[group=="Carbohydrates"]
    if( nrow(src.carbo)>0 )
      src.add <- src.carbo # if no glucose is there, then add all other available carbohydrates
    else
      src.add <- src.met # if no carbohydrates is avaiable, then take everything else (probably amino acid biosynthesis is not papfilled because amino acids are part of the medium)
    media2 <- rbind(media2, data.table(compounds=gsub("\\[.0\\]","",src.add$seed), name=src.add$name, maxFlux=100))  
  }
  
  
  # constrain model  
  mod.orig2 <- constrain.model(mod.orig2, media = media2)
  mod.orig2@obj_coef <- rep(0,mod.orig2@react_num)
  
  bm.ind      <- which(mod.orig2@react_id == "bio1")
  bm.met.inds <- which(mod.orig2@S[,bm.ind]<0)
  bm.met      <- gsub("\\[.0\\]","",mod.orig2@met_id[bm.met.inds])
  bm.met.name <- mod.orig2@met_name[bm.met.inds]
  mod.fill2    <- mod.orig2
  mod.fill2.counter <- 0
  mod.fill2.names <- c()
  
  if( !verbose ) options(warn=-1)
  for( i in seq_along(bm.met) ){
    cat("\r",i,"/",length(bm.met))
    target.new <- bm.met[i]
    
    # add metabolite objective + sink
    rm.sink = TRUE
    if( paste0("EX_",target.new,"_c0") %in% react_id(mod.fill2) )
      rm.sink = FALSE
    mod.fill2  <- add_met_sink(mod.fill2, target.new, obj = 1)
    
    sol <- optimizeProb(mod.fill2, retOptSol=F)
    
    if(sol$stat == ok & sol$obj >= 1e-6){
      mod.fill2@obj_coef <- rep(0,mod.fill2@react_num)
    }else{
      if( verbose ) cat("\nTry to gapfill", bm.met.name[i],"\n")
      invisible(capture.output( 
        mod.fill2.lst <- gapfill4(mod.orig = mod.fill2, 
                                  mod.full = mod,
                                  rxn.weights = copy(rxn.weights),  
                                  min.gr = min.obj.val,
                                  bcore = bcore,
                                  dummy.weight = dummy.weight,
                                  script.dir = script.dir,
                                  core.only = use.core,
                                  verbose=verbose,
                                  gs.origin = 2,
                                  rXg.tab = rXg.tab,
                                  env = env) ))
      new.reactions <- mod.fill2.lst$rxns.added
      if( length(new.reactions) > 0 ){
        if( verbose ) cat("Added reactions:", new.reactions, "\n")
        mod.fill2 <- mod.fill2.lst$model
        mod.fill2.counter <- mod.fill2.counter + 1
        mod.fill2.names <- c(mod.fill2.names, bm.met.name[i])
      }
      mod.fill2@obj_coef <- rep(0,mod.fill2@react_num)
    }
    if( rm.sink )
      mod.fill2 <- rmReact(mod.fill2, react=paste0("EX_",target.new,"_c0"))
  }
  options(warn=0)
  
  mod.fill2 <- changeObjFunc(mod.fill2, react=paste0("EX_",target.met,"_c0"))
  mod.fill2 <- constrain.model(mod.fill2, media.file = media.file, scaling.fac = 1)
  mod.out <- mod.fill2
  
  cat("\rGapfill summary:\n")
  cat("Filled components:    ",mod.fill2.counter, "(",paste(mod.fill2.names, collapse = ","),")\n")
  cat("Added reactions:      ",length(mod.fill2@react_id)-length(mod.orig2@react_id),"\n")
  cat("Final growth rate:    ",optimizeProb(mod.fill2, retOptSol=F)$obj,"\n")
  
  
  
  # Add list of exchange reactions for step 3 and 4 in order to check for a wide range of carbon sources or fermentation products
  # (Unused exchanges will be deleted afterwards)
  mod.out <- add_missing_exchanges(mod.out)
  carbon.source <- carbon.source[!is.na(name) & !is.na(seed) & seed!=""]
  idx <- which( !carbon.source$seed %in% mod.out@react_id )
  exchanges.new.met  <- str_replace(str_remove(carbon.source$seed[idx], "EX_"),"_e0","\\[e0\\]")
  exchanges.new.name <- mod@met_name[match(exchanges.new.met,mod@met_id)]
  exchanges.new.ids  <- carbon.source$seed[idx]
  exchanges.new.used  <- rep(FALSE, length(exchanges.new.ids))  # delete unused addionally added exchange reactions later
  mod.out       <- add_exchanges(mod.out, exchanges.new.met, metname=exchanges.new.name)

  if ( !quick.gf ){
    cat("\n\n3. Energy source gapfilling with core reactions only\n")
    
    mod.orig3 <- mod.out
    media.org <- fread(paste0(script.dir,"/../dat/media/MM_glu.csv")) # use minimal medium
    #media.org <- fread(paste0(script.dir,"/../dat/media/Mineral_salt.csv")) # use minimal medium
    
    ex          <- findExchReact(mod.orig3)
    ex.ind      <- ex@react_pos
    ex.id       <- ex@react_id
    ex.met      <- ex@met_id
    ex.met.name <- mod.orig3@met_name[ex@met_pos]
    if ( length(met.limit) > 0 ){ # # potentially limit carbon.source
      ex.idx <- match(intersect(ex.id, carbon.source$seed), ex.id)
      ex.ind      <- ex.ind[ex.idx]
      ex.id       <- ex.id[ex.idx]
      ex.met      <- ex.met[ex.idx]
      ex.met.name <- ex.met.name[ex.idx]
    }
    # Exchange reactions to be ignored (metals etc.)
    ignore <- c("EX_cpd17041_e0", "EX_cpd17042_e0", "EX_cpd17043_e0", "EX_cpd11416_e0", "rxn13782_c0", "rxn13783_c0", "rxn13783_c0", "EX_cpd00001_e0","EX_cpd00007_e0", "EX_cpd00009_e0", "EX_cpd00011_e0" ,"EX_cpd00012_e0", "EX_cpd00030_e0", "EX_cpd00034_e0", "EX_cpd00058_e0", "EX_cpd00063_e0", "EX_cpd00067_e0", "EX_cpd00075_e0","EX_cpd00099_e0", "EX_cpd00149_e0", "EX_cpd00205_e0", "EX_cpd00254_e0", "EX_cpd10515_e0", "EX_cpd00971_e0", "EX_cpd01012_e0", "EX_cpd10516_e0", "EX_cpd11574_e0")
    
    # add metabolite objective + sink
    mod.fill3    <- mod.orig3
    mod.fill3@obj_coef <- rep(0,mod.fill3@react_num)
    
    # add biolog like test
    mql  <- "cpd15499[c0]"; mqn   <- "cpd15500[c0]" # menaquinone
    uql  <- "cpd15561[c0]"; uqn   <- "cpd15560[c0]" # ubiquinone
    h    <- "cpd00067[c0]"
    nad  <- "cpd00003[c0]"; nadh  <- "cpd00004[c0]"
    fdox <- "cpd11621[c0]"; fdred <- "cpd11620[c0]" # ferredoxin
    pql  <- "cpd27796[c0]"; pqn   <- "cpd27797[c0]" # plastoquinone
    mod.fill3 <- addReact(mod.fill3, "ESP1", met=c(mql,h,mqn), Scoef=c(-1,2,1), lb=0, ub=1000, metComp = rep(1,3))
    mod.fill3 <- addReact(mod.fill3, "ESP2", met=c(uql,h,uqn), Scoef=c(-1,2,1), lb=0, ub=1000, metComp = rep(1,3))
    mod.fill3 <- addReact(mod.fill3, "ESP3", met=c(nadh,h,nad), Scoef=c(-1,1,1), lb=0, ub=1000, metComp = rep(1,3))
    mod.fill3 <- addReact(mod.fill3, "ESP4", met=c(fdred,fdox), Scoef=c(-1,1), lb=0, ub=1000, metComp = rep(1,2))
    mod.fill3 <- addReact(mod.fill3, "ESP5", met=c(pql,h,pqn), Scoef=c(-1,2,1), lb=0, ub=1000, metComp = rep(1,3))
    mod.fill3 <- changeObjFunc(mod.fill3, react=c("ESP1", "ESP2", "ESP3", "ESP4", "ESP5"), obj_coef=c(1,1,1,1,1))
    mod.fill3.counter <- 0
    mod.fill3.names <- c()
    
    if( !verbose ) options(warn=-1)
    cs.dt <- data.table()
    for( i in seq_along(ex.met) ){
      cat("\r",i,"/",length(ex.met))
      if( ex.id[i] %in% ignore ) 
        next
      if( !checkIfGapfillIsPossible(ex.met[i]) )
        next
      
      src.met      <- ex.met[i]
      src.met.name <- ex.met.name[i]
      media <- media.org[name!="D-Glucose"]
      #media <- media.org[!name %in% c("Benzoate", "co2")]
      media <- rbind(media, data.table(compounds=gsub("\\[.0\\]","",src.met), name=src.met.name, maxFlux=100))
      
      # constrain model
      mod.fill3 <- constrain.model(mod.fill3, media = media)
      
      sol <- optimizeProb(mod.fill3, retOptSol=F)
      
      if(sol$stat == ok & sol$obj >= 1e-7){
        #mod.fill3@obj_coef <- rep(0,mod.fill3@react_num)
        src.status <- TRUE
      }else{
        if( verbose ) cat("\nTry to gapfill", src.met.name, ex.id[i], "\n")
        invisible(capture.output( mod.fill3.lst <- gapfill4(mod.orig = mod.fill3, 
                                                            mod.full = mod, 
                                                            rxn.weights = copy(rxn.weights),  
                                                            min.gr = min.obj.val,
                                                            bcore = bcore,
                                                            dummy.weight = dummy.weight,
                                                            script.dir = script.dir,
                                                            core.only = use.core,
                                                            verbose=verbose,
                                                            gs.origin = 3,
                                                            rXg.tab = rXg.tab,
                                                            env = env) ))
        src.status <- mod.fill3.lst$growth.rate >= 1e-7
        new.reactions <- mod.fill3.lst$rxns.added
        if( length(new.reactions) > 0 ){
          if( verbose ) cat("Added reactions:", new.reactions, "\n")
          mod.fill3 <- mod.fill3.lst$model
          mod.fill3.counter <- mod.fill3.counter + 1
          mod.fill3.names <- c(mod.fill3.names, src.met.name)
          if( ex.id[i] %in% exchanges.new.ids) # delete unused addionally added exchange reactions later
            exchanges.new.used[match(ex.id[i], exchanges.new.ids)] <- TRUE
        }
      }
      cs.dt <- rbind(cs.dt, data.table(id=str_extract(src.met,"cpd[0-9]+"), name=src.met.name, status=src.status))
    }
    options(warn=0)
    
    mod.fill3 <- rmReact(mod.fill3, react=c("ESP1","ESP2", "ESP3", "ESP4", "ESP5"))
    mod.fill3 <- changeObjFunc(mod.fill3, react=paste0("EX_",target.met,"_c0"))
    mod.fill3 <- constrain.model(mod.fill3, media.file = media.file, scaling.fac = 1)
    mod.out <- mod.fill3
    cat("\rGapfill summary:\n")
    cat("Filled components:    ",mod.fill3.counter, "(",paste(mod.fill3.names, collapse = ","),")\n")
    cat("Added reactions:      ",length(mod.fill3@react_id)-length(mod.fill2@react_id),"\n")
    cat("Final growth rate:    ",optimizeProb(mod.fill3, retOptSol=F)$obj,"\n")
  }
  
  
  if ( !quick.gf ){
    cat("\n\n4. Checking for potential metabolic products with core reactions only\n")
    
    mod.orig4 <- mod.out
    
    ex          <- findExchReact(mod.orig4)
    ex.ind      <- ex@react_pos
    ex.id       <- ex@react_id
    ex.met      <- ex@met_id
    ex.met.name <- mod.orig4@met_name[ex@met_pos]
    if ( length(met.limit) > 0 ){ # # potentially limit carbon.source
      ex.idx <- match(intersect(ex.id, carbon.source$seed), ex.id)
      ex.ind      <- ex.ind[ex.idx]
      ex.id       <- ex.id[ex.idx]
      ex.met      <- ex.met[ex.idx]
      ex.met.name <- ex.met.name[ex.idx]
    }
    # Exchange reactions to be ignored (metals etc.)
    ignore <- c("EX_cpd17041_e0", "EX_cpd17042_e0", "EX_cpd17043_e0", "EX_cpd11416_e0", "rxn13782_c0", "rxn13783_c0", "rxn13783_c0", "EX_cpd00001_e0","EX_cpd00007_e0", "EX_cpd00009_e0", "EX_cpd00011_e0" ,"EX_cpd00012_e0", "EX_cpd00030_e0", "EX_cpd00034_e0", "EX_cpd00058_e0", "EX_cpd00063_e0", "EX_cpd00067_e0", "EX_cpd00075_e0","EX_cpd00099_e0", "EX_cpd00149_e0", "EX_cpd00205_e0", "EX_cpd00254_e0", "EX_cpd10515_e0", "EX_cpd00971_e0", "EX_cpd01012_e0", "EX_cpd10516_e0", "EX_cpd11574_e0")
    
    # add metabolite objective + sink
    mod.fill4    <- mod.orig4
    mod.fill4 <- constrain.model(mod.fill4, media.file = media.file)
    
    mod.fill4.counter <- 0
    mod.fill4.names <- c()
    
    ferm.dt <- data.table()
    if( !verbose ) options(warn=-1)
    for( i in seq_along(ex.met) ){
      cat("\r",i,"/",length(ex.met))
      if( ex.id[i] %in% ignore ) 
        next
      if( !checkIfGapfillIsPossible(ex.met[i]) )
        next
      
      src.met      <- ex.met[i]
      src.met.name <- ex.met.name[i]
      src.id       <- ex.id[i]
      mod.fill4@obj_coef <- rep(0,mod.fill4@react_num)
      mod.fill4 <- changeObjFunc(mod.fill4, react=src.id, obj_coef=1)
      
      sol <- optimizeProb(mod.fill4, retOptSol=F)
      
      if(sol$stat == ok & sol$obj >= 1e-7){
        #mod.fill4@obj_coef <- rep(0,mod.fill4@react_num)
        src.status <- TRUE
      }else{
        if( verbose ) cat("\nTry to gapfill", src.met.name, src.id, "\n")
        invisible(capture.output( mod.fill4.lst <- gapfill4(mod.orig = mod.fill4, 
                                                            mod.full = mod, 
                                                            rxn.weights = copy(rxn.weights), 
                                                            min.gr = min.obj.val,
                                                            bcore = bcore,
                                                            dummy.weight = dummy.weight,
                                                            script.dir = script.dir,
                                                            core.only = use.core,
                                                            verbose=verbose,
                                                            gs.origin = 4,
                                                            rXg.tab = rXg.tab,
                                                            env = env) ))
        src.status <- mod.fill4.lst$growth.rate >= 1e-7
        new.reactions <- mod.fill4.lst$rxns.added
        if( length(new.reactions) > 0 ){
          if( verbose ) cat("Added reactions:", new.reactions, "\n")
          mod.fill4 <- mod.fill4.lst$model
          mod.fill4.counter <- mod.fill4.counter + 1
          mod.fill4.names <- c(mod.fill4.names, src.met.name)
          if( ex@react_id[i] %in% exchanges.new.ids) # delete unused addionally added exchange reactions later
            exchanges.new.used[match(ex@react_id[i], exchanges.new.ids)] <- TRUE
        }
      }
      ferm.dt <- rbind(ferm.dt, data.table(id=str_extract(src.met,"cpd[0-9]+"), name=src.met.name, status=src.status))
    }
    options(warn=0)
    
    mod.fill4 <- changeObjFunc(mod.fill4, react=paste0("EX_",target.met,"_c0"))
    mod.fill4 <- constrain.model(mod.fill4, media.file = media.file, scaling.fac = 1)
    mod.out <- mod.fill4
    
    mod.fill4.sol <- optimizeProb(mod.fill4, retOptSol=F, algorithm = "mtf")
    dt.sol        <- data.table(rxn = mod.fill4@react_id, 
                                flux = mod.fill4.sol$fluxes[1:mod.fill4@react_num], 
                                lb = mod.fill4@lowbnd,
                                met.name = gsub("-e0 Exchange","",mod.fill4@react_name))
    dt.sol[, met.name := gsub(" Exchange","", met.name)]
    dt.sol.u      <- copy(dt.sol[flux < 0 & grepl("^EX_", rxn) & flux <= lb*0.999])
    dt.sol.p      <- copy(dt.sol[flux > 0 & grepl("^EX_", rxn)][order(-flux)][1:min(c(10,.N))])
    
    cat("\rGapfill summary:\n")
    cat("Filled components:    ",mod.fill4.counter, "(",paste(mod.fill4.names, collapse = ","),")\n")
    cat("Added reactions:      ",length(mod.fill4@react_id)-length(mod.fill3@react_id),"\n")
    cat("Final growth rate:    ",mod.fill4.sol$fluxes[which(mod.fill4@obj_coef==1)],"\n\n")
    
    cat("Uptake at limit:\n")
    cat(paste0(paste(dt.sol.u$met.name, round(-dt.sol.u$flux, digits = 3), sep = ":"), collapse = ", "),"\n\n")
    
    cat("Top 10 produced metabolites [mmol / (gDW * hr)]:\n")
    cat(paste0(paste(dt.sol.p$met.name, round(dt.sol.p$flux, digits = 3), sep = ":"), collapse = ", "),"\n")
  }
  
  
}
mod.out <- add_missing_exchanges(mod.out)

# delete unused addionally added exchange reactions later
exchanges.rm <- exchanges.new.ids[!exchanges.new.used]
if( length(exchanges.rm) > 0 )
mod.out <- rmReact(mod.out, react=exchanges.rm)
mod.out <- rm_unused_exchanges(mod.out)


# add metabolite-, reaction-, and model attributes
mod.out <- addMetAttr(mod.out, seed_x_mets = seed_x_mets)
if("ec" %in% colnames(mod.out@react_attr)) mod.out <- addReactAttr(mod.out)
mod.out@mod_attr <- bu_mod_attr

out.id <- gsub("\\.xml$|\\.RDS$|\\.rds$|\\.xml\\.gz$","",gsub("-draft","",basename(mod.file)))

if( verbose ){
  mod.out.rxns.added <- setdiff(mod.out@react_id, mod.orig@react_id)
  cat(mod.out.rxns.added, file = paste0(output.dir,"/",out.id,"-gapfilled.rxnlst"))
  
  mod.out.rxns.added.without.seq <- setdiff(gsub("_.0","",mod.out.rxns.added), rxn.weights[bitscore>bcore, seed])
  cat(mod.out.rxns.added.without.seq, file = paste0(output.dir,"/",out.id,"-gapfilled.without.seq.rxnlst"))  
}

# add gapseq version info to model object
gapseq_version <- system(paste0(script.dir,"/.././gapseq -v"), intern = T)
seqdb_version  <- str_match(mod.orig@mod_desc, "Sequence DB md5sum: .*")
if( !is.na(seqdb_version) ){
    mod.out@mod_desc <- paste0(gapseq_version, "; ", seqdb_version)
} else mod.out@mod_desc <- gapseq_version

out.rds <- paste0(output.dir,"/",out.id,".RDS")

if(file.exists(out.rds)) warning("Model file already exists and will be overwritten!")
saveRDS(mod.out, file = out.rds)
# Write SBML
if(!opt$sbml.no.output){
  source(paste0(script.dir,"/sbml_write.R"))
  write_gapseq_sbml(mod.out, paste0(output.dir,"/",out.id))
}

# Save additionally an unconstrained version of the model if desired
if(relaxed.constraints) {
  mod.out@lowbnd[grep("^EX_",mod.out@react_id)] <- -sybil::SYBIL_SETTINGS("MAXIMUM")
  saveRDS(mod.out, file = paste0(output.dir,"/",out.id,"-unconstrained.RDS"))
  if(!opt$sbml.no.output){
    source(paste0(script.dir,"/sbml_write.R"))
    write_gapseq_sbml(mod.out, paste0(output.dir,"/",out.id,"-unconstrained"))
  }
}

# Save found carbon sources and fermentation products
if(write.cs.ferm){
  fwrite(ferm.dt, paste0(output.dir,"/",out.id, "-ferm.tbl"), sep="\t")
  fwrite(cs.dt,   paste0(output.dir,"/",out.id, "-cs.tbl"),   sep="\t")
}

q(status=0)