Use taxonomic domain information to exclude specific reactions

Mainly meant to prevent gapfilling with domain-specific reactions in
reconstructions for different domains.

dat/biomass/biomass_Gram_neg.json

@@ -2,6 +2,7 @@
   "name" : "Bacterial Gram-negative biomass reaction",
   "ref"  : "derived from ModelSEED and https://doi.org/10.1038/msb.2011.65",
   "energy_GAM" :  40,
+  "domain" : "Bacteria",
   "met_groups" : [
     {
       "group_name" : "DNA",

dat/biomass/biomass_Gram_pos.json

@@ -2,6 +2,7 @@
   "name" : "Bacterial Gram-positive biomass reaction",
   "ref" : "derived from ModelSEED and https://doi.org/10.1074/jbc.m703759200", 
   "energy_GAM" :  41.257,
+  "domain" : "Bacteria",
   "met_groups" : [
     {
       "group_name" : "DNA",

dat/biomass/biomass_archaea.json

@@ -2,6 +2,7 @@
   "name" : "Archaeal biomass reaction",
   "ref" : "derived from https://doi.org/10.1038/msb4100046", 
   "energy_GAM" :  30,
+  "domain" : "Archaea",
   "met_groups" : [
     {
       "group_name" : "DNA",

dat/biomass/excluded_reactions.tsv

@@ -0,0 +1,5 @@
+domain	exclude.reaction	exclusion.reason
+Bacteria	rxn15962	involves Archaea-specific co-factor
+Bacteria	rxn15964	involves Archaea-specific co-factor
+Bacteria	rxn90052	involves Archaea-specific co-factor
+Bacteria	rxn03127	involves Archaea-specific co-factor

dat/media/MM_anaerobic_Acetate.csv

@@ -5,7 +5,7 @@ cpd11640,H2,0
 cpd00009,Phosphate,100
 cpd00030,Mn2+,100
 cpd00034,Zn2+,100
-cpd00048,Sulfate,100
+cpd00048,Sulfate,2.5
 cpd00058,Cu2+,100
 cpd00063,Ca2+,100
 cpd00067,H+,100

src/generate_GSdraft.R

@@ -261,25 +261,22 @@ build_draft_model_from_blast_results <- function(blast.res, transporter.res, bio
     mod <- add_reaction_from_db(mod, react = "rxn05683", gs.origin = 1)
   }
 
-  #mod@genes_table <- copy(dt[bitscore>0])
   cat("\n")
   warnings()
   # Adding Biomass reaction
   if(biomass == "neg" | biomass == "Gram_neg"){
     ls.bm <- parse_BMjson(paste0(script.dir, "/../dat/biomass/biomass_Gram_neg.json"), seed_x_mets)
     dt.bm <- ls.bm$bmS
-    #dt.bm <- fread(paste0(script.dir, "/../dat/biomass/seed_biomass.DT_gramNeg.tsv"))
   }
   if(biomass == "pos" | biomass == "Gram_pos"){
     ls.bm <- parse_BMjson(paste0(script.dir, "/../dat/biomass/biomass_Gram_pos.json"), seed_x_mets)
     dt.bm <- ls.bm$bmS
-    #dt.bm <- fread(paste0(script.dir, "/../dat/biomass/seed_biomass.DT_gramPos.tsv"))
   }
   if(biomass == "archaea" | biomass == "Archaea"){
     ls.bm <- parse_BMjson(paste0(script.dir, "/../dat/biomass/biomass_archaea.json"), seed_x_mets)
     dt.bm <- ls.bm$bmS
-    #dt.bm <- fread(paste0(script.dir, "/../dat/biomass/seed_biomass.DT_archaea.tsv"))
   }
+  mod@mod_attr <- data.frame(annotation = paste0("tax_domain:",ls.bm$domain))
 
 
   if(biomass %in% c("neg","pos","Gram_neg","Gram_pos")) {

src/gf.suite.R

@@ -146,10 +146,39 @@ if ( toupper(file_ext(mod.file)) == "RDS" ){
 }else{ 
   mod.orig <- readSBMLmod(mod.file)}
 
+bu_mod_attr <- mod.orig@mod_attr
+
 # 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),
@@ -624,9 +653,10 @@ if(nrow(mseed.t)>0) { # Skip steps 2,2b,3, and 4 if core-reaction list does not
 
 mod.out <- add_missing_exchanges(mod.out)
 
-# add metabolite & reactions attributes
+# add metabolite-, reaction-, and model attributes
 mod.out <- addMetAttr(mod.out, seed_x_mets = seed_x_mets)
 mod.out <- addReactAttr(mod.out)
+mod.out@mod_attr <- bu_mod_attr
 
 if(!dir.exists(output.dir))
   system(paste0("mkdir ",output.dir))

src/parse_BMjson.R

@@ -98,9 +98,10 @@ parse_BMjson <- function(json.file, gs_mets) {
   grp_tab[, id := paste0(id,"[",comp,"0]")]
 
   res <- list()
-  res$id   <- bmjs$id
-  res$name <- bmjs$name
-  res$bmS  <- grp_tab
+  res$id     <- bmjs$id
+  res$name   <- bmjs$name
+  res$domain <- bmjs$domain
+  res$bmS    <- grp_tab
 
   return(res)
 }