Reaction flux coupling to biomass formation

- Adds constraints to community models, were the absolute reaction flux
bounds of organism x is propotional to the biomass formation of organism
x.

DESCRIPTION

@@ -44,6 +44,7 @@ Collate:
     'deadEndMetabolites.R'
     'fba.R'
     'fixBMRatios.R'
+    'fluxBMCoupling.R'
     'fluxdist_analysis.R'
     'fva.R'
     'geneDel.R'

NAMESPACE

@@ -19,6 +19,7 @@ export(constraint_num)
 export(deadEndMetabolites)
 export(fba)
 export(fixBMRatios)
+export(fluxBMCoupling)
 export(fva)
 export(geneDel)
 export(gene_num)

R/LPproblem_glpkClass.R

@@ -247,12 +247,14 @@ setMethod("solveLp", signature(lp = "LPproblem_glpk"),
                      out <- solveInterior(lp@ptr)
                    },
                    "exact" = {
+                     invisible(scaleSimplexProb(lp@ptr))
                      out <- solveSimplexExact(lp@ptr)
                    },
                    "mip" = {
                      out <- solveMIP(lp@ptr)
                    },
                    {
+                     scaleSimplexProb(lp@ptr)
                      out <- solveSimplex(lp@ptr)
                    }
             )

R/ModelOrgClass.R

@@ -47,7 +47,6 @@
 #'
 #' @aliases ModelOrg
 #'
-#' @family Object classes
 #' @exportClass ModelOrg
 setClass("ModelOrg",
 
@@ -154,7 +153,6 @@ setClass("ModelOrg",
 #'
 #' @docType methods
 #' @rdname react_num-methods
-#' @family Model characteristics
 #' @export
 setGeneric("react_num", valueClass = "numeric", function(model) {
   standardGeneric("react_num")
@@ -181,7 +179,6 @@ setMethod("react_num", signature(model = "ModelOrg"),
 #'
 #' @docType methods
 #' @rdname react_pos-methods
-#' @family Model characteristics
 #' @export
 setGeneric("react_pos", valueClass = "numeric", function(model, react) {
   standardGeneric("react_pos")
@@ -217,7 +214,6 @@ setMethod("react_pos", signature(model = "ModelOrg", react = "missing"),
 #'
 #' @docType methods
 #' @rdname met_num-methods
-#' @family Model characteristics
 #' @export
 setGeneric("met_num", valueClass = "numeric", function(model) {
   standardGeneric("met_num")
@@ -244,7 +240,6 @@ setMethod("met_num", signature(model = "ModelOrg"),
 #'
 #' @docType methods
 #' @rdname met_pos-methods
-#' @family Model characteristics
 #' @export
 setGeneric("met_pos", valueClass = "numeric", function(model, met) {
   standardGeneric("met_pos")
@@ -279,7 +274,6 @@ setMethod("met_pos", signature(model = "ModelOrg", met = "missing"),
 #'
 #' @docType methods
 #' @rdname gene_num-methods
-#' @family Model characteristics
 #' @export
 setGeneric("gene_num", valueClass = "numeric", function(model) {
   standardGeneric("gene_num")
@@ -307,7 +301,6 @@ setMethod("gene_num", signature(model = "ModelOrg"),
 #'
 #' @docType methods
 #' @rdname gene_pos-methods
-#' @family Model characteristics
 #' @export
 setGeneric("gene_pos", valueClass = "numeric", function(model, gene) {
   standardGeneric("gene_pos")
@@ -342,7 +335,6 @@ setMethod("gene_pos", signature(model = "ModelOrg", gene = "missing"),
 #'
 #' @docType methods
 #' @rdname comp_num-methods
-#' @family Model characteristics
 #' @export
 setGeneric("comp_num", valueClass = "numeric", function(model) {
   standardGeneric("comp_num")
@@ -369,7 +361,6 @@ setMethod("comp_num", signature(model = "ModelOrg"),
 #'
 #' @docType methods
 #' @rdname comp_pos-methods
-#' @family Model characteristics
 #' @export
 setGeneric("comp_pos", valueClass = "numeric", function(model, comp) {
   standardGeneric("comp_pos")
@@ -411,7 +402,6 @@ setMethod("comp_pos", signature(model = "ModelOrg", comp = "logical"),
 #'
 #' @docType methods
 #' @rdname constraint_num-methods
-#' @family Model characteristics
 #' @export
 setGeneric("constraint_num", valueClass = "numeric", function(model) {
   standardGeneric("constraint_num")
@@ -432,7 +422,6 @@ setMethod("constraint_num", signature(model = "ModelOrg"),
 #'
 #' @docType methods
 #' @rdname subsys_num-methods
-#' @family Model characteristics
 #' @export
 setGeneric("subsys_num", valueClass = "numeric", function(model) {
   standardGeneric("subsys_num")
@@ -459,7 +448,6 @@ setMethod("subsys_num", signature(model = "ModelOrg"),
 #'
 #' @docType methods
 #' @rdname subsys_pos-methods
-#' @family Model characteristics
 #' @export
 setGeneric("subsys_pos", valueClass = "numeric", function(model, subsys) {
   standardGeneric("subsys_pos")
@@ -526,16 +514,17 @@ setMethod("printObjFunc", signature(object = "ModelOrg"),
 #'
 #' @param object S4-object of class \link{ModelOrg}.
 #'
-#' @family Model characteristics
 #' @export
 setMethod("show", signature(object = "ModelOrg"),
           function(object) {
             cat("model ID:                  ", object@mod_id, "\n")
             cat("model name:                ", object@mod_name, "\n")
             cat("number of compartments:    ", length(object@mod_compart), "\n")
-            for(i in 1:length(object@mod_compart)) {
-              cat("                           ", object@mod_compart[i], " (",
-                  object@mod_compart_name[i], ")\n")
+            if(comp_num(object) > 0 && comp_num(object) <= 10) {
+              for(i in 1:length(object@mod_compart)) {
+                cat("                           ", object@mod_compart[i], " (",
+                    object@mod_compart_name[i], ")\n")
+              }
             }
             cat("number of reactions:       ", react_num(object), "\n")
             cat("number of metabolites:     ", met_num(object), "\n")

R/RcppExports.R

@@ -53,6 +53,10 @@ solveSimplexExact <- function(xp) {
     .Call(`_cobrar_solveSimplexExact`, xp)
 }
 
+scaleSimplexProb <- function(xp) {
+    .Call(`_cobrar_scaleSimplexProb`, xp)
+}
+
 getObjVal <- function(xp) {
     .Call(`_cobrar_getObjVal`, xp)
 }

R/fluxBMCoupling.R

@@ -0,0 +1,57 @@
+#' Couple reaction flux bounds to biomass production
+#'
+#' Resets the absolute flux bounds of reactions proportional to the flux through
+#' the biomass reaction.
+#'
+#' @param model Community model of class \link{ModelComm}
+#' @param BMreact IDs of individual biomass reactions
+#' @param cpl_c Coupling constraint \code{c}
+#' @param cpl_u Coupling constraint \code{u}
+#'
+#' @details
+#' See (and cite) Heinken et al. (2013) Gut Microbes (doi: 10.4161/gmic.22370).
+#'
+#' @export
+fluxBMCoupling <- function(model, BMreact = guessBMReaction(model),
+                           cpl_c = 400, cpl_u = 0.01) {
+  if(length(BMreact) != nrow(model@community))
+    stop("Length of 'BMreact' not equal to number of organisms in the community.")
+
+  bmIdx <- unlist(regmatches(BMreact,gregexpr("^M[0-9]+_", BMreact)))
+  bmIdx <- gsub("_$","",bmIdx)
+
+  if(any(duplicated(bmIdx)))
+    stop("There are two or more biomass reactions for at least on organism, which is not allowed. Please check 'BMreact'.")
+
+  if(!all(bmIdx %in% model@community$index))
+    stop("At least one provided biomass reaction ID cannot be matched to a community organism.")
+
+  # make sure biomass reactions are in the same order as organisms within the
+  # community
+  names(BMreact) <- bmIdx
+  BMreact <- BMreact[match(bmIdx,model@community$index)]
+  n <- length(BMreact)
+
+  # identify reactions to be coupled (exclude organisms' exchange reactions)
+  mrxns <- model@react_id[grepl("^M[0-9]+_", model@react_id)]
+  mrxns <- mrxns[!grepl("^M[0-9]+_EX_", mrxns)]
+  mrxns <- mrxns[!(mrxns %in% BMreact)]
+  mrxnsIdx <- unlist(regmatches(mrxns,gregexpr("^M[0-9]+_", mrxns)))
+  mrxnsIdx <- gsub("_$","",mrxnsIdx)
+  reaBMpair <- mapply(function(x,y) {c(x,y)}, BMreact[mrxnsIdx], mrxns,
+                      SIMPLIFY = FALSE)
+
+  nr <- length(mrxns)
+  reaCoefL <- mapply(function(x,y) {c(x,y)}, rep(cpl_c,nr), rep(-1,nr),
+                     SIMPLIFY = FALSE)
+  reaCoefU <- mapply(function(x,y) {c(x,y)}, rep(-cpl_c,nr), rep(-1,nr),
+                     SIMPLIFY = FALSE)
+
+  model <- addConstraint(model,
+                         react = c(reaBMpair, reaBMpair),
+                         coeff = c(reaCoefL, reaCoefU),
+                         rtype = c(rep("L",nr),rep("U",nr)),
+                         lb = rep(-cpl_u,nr*2),
+                         ub = rep(cpl_u,nr*2))
+  return(model)
+}

R/io.R

@@ -10,7 +10,6 @@
 #' mod
 #'
 #' @import Matrix
-#' @family Model import and export helpers
 #' @export
 readSBMLmod <- function(file_path) {
 
@@ -142,7 +141,6 @@ sboterm2int <- function(sbo) {
 #'
 #' @returns TRUE if file export was successful.
 #'
-#' @family Model import and export helpers
 #' @export
 writeSBMLmod <- function(model, file_path = NULL) {
 
@@ -287,7 +285,6 @@ writeSBMLmod <- function(model, file_path = NULL) {
 #' objects.
 #'
 #' @import Matrix
-#' @family Model import and export helpers
 #' @export
 readSybilmod <- function(file_path) {
   sybildoc.lst  <- readRDS(normalizePath(file_path))

R/joinModels.R

@@ -118,7 +118,6 @@ joinModels <- function(models, mergeLB = "none", mergeUB = "none",
   })
 
   for(i in 1:length(inds)) {
-    print(i)
     indstmp <- inds[[i]]
     indstmp[,"row"] <- indstmp[,"row"]+metoffset[i]
     indstmp[,"col"] <- indstmp[,"col"]+rxnoffset[i]
@@ -247,7 +246,6 @@ joinModels <- function(models, mergeLB = "none", mergeUB = "none",
   })
 
   for(i in 1:length(inds)) {
-    print(i)
     indstmp <- inds[[i]]
     indstmp[,"row"] <- indstmp[,"row"]+constroffset[i]
     indstmp[,"col"] <- indstmp[,"col"]+rxnoffset[i]

R/model_modifications.R

@@ -228,7 +228,6 @@ setMethod("addConstraint", signature(model = "ModelOrg",
                                      coeff = "list",
                                      rtype = "character"),
           function(model, react, coeff, rtype, lb = NULL, ub = NULL) {
-
             nc <- length(react)
 
             if(is.null(lb))
@@ -259,31 +258,30 @@ setMethod("addConstraint", signature(model = "ModelOrg",
             lb[indtmp] <- ub[indtmp]
             indtmp <- which(rtype == "E" & is.na(ub))
             ub[indtmp] <- lb[indtmp]
-
             if(any(unlist(lapply(react, length)) != unlist(lapply(coeff, length))))
               stop("List elementes of 'react' must have the same length as the corresponding list elements in 'coeff'.")
-
             if(any(unlist(lapply(react, duplicated))))
               stop("'react' IDs cannot be duplicated within a constraint definition.")
-
-            if(any(unlist(lapply(react, function(x) !(x %in% model@react_id)))))
+            if(any(!(unique(unlist(react)) %in% model@react_id)))
               stop("Not all reaction IDs in 'react' are part of the model.")
 
-            I <- matrix(c(rep(1:nc, unlist(lapply(react, length))),
-                          unlist(lapply(react, function(x) match(x, model@react_id)))),
+            lookup_table <- setNames(1:react_num(model), model@react_id)
+            creacts <- unlist(lapply(react, function(x) x))
+            creactsIdx <- lookup_table[creacts]
+            I <- matrix(c(rep(1:nc, lengths(react)), creactsIdx),
                         ncol = 2)
 
-
             out <- Matrix(0, nrow = nc, ncol = react_num(model), sparse = T)
             out[I] <- unlist(coeff)
 
             model@constraints@coeff <- rbind(model@constraints@coeff,
-                                              out)
+                                             out)
+
             model@constraints@lb <- c(model@constraints@lb, lb)
             model@constraints@ub <- c(model@constraints@ub, ub)
             model@constraints@rtype <- c(model@constraints@rtype, rtype)
 
-            return(rmDuplicateConstraints(model))
+            return(model)
           }
 )