diff --git a/R/generateWT.R b/R/generateWT.R
index 16cf50a677d7273981691af362fc95e06d78d656..595a005faa4eea3f23b854f31431fc0ddd04d4bf 100644
--- a/R/generateWT.R
+++ b/R/generateWT.R
@@ -28,7 +28,8 @@
#
#
-.generateWT <- function(model, react = NULL, lb = NULL, ub = NULL, ...) {
+.generateWT <- function(model, react = NULL, lb = NULL, ub = NULL,
+ algorithm="fba", ...) {
ca <- match.call()
@@ -59,11 +60,11 @@
if (is(react, "list")) {
message("calculating fba solutions ... ", appendLF = FALSE)
suppressMessages(
- tmp <- optimizer(model, algorithm = "fba",
+ tmp <- optimizer(model, algorithm = algorithm,
lpdir = rep("max", length(react)),
react = react, lb = lb, ub = ub, verboseMode = 0,
solver = me[["sol"]], method = me[["met"]],
- solverParm = as.data.frame(NA))
+ solverParm = as.data.frame(NA), ...)
)
message("OK")
}
@@ -71,11 +72,11 @@
tmp <- optimizeProb(model,
react = react, lb = lb, ub = ub,
retOptSol = FALSE,
- algorithm = "fba",
+ algorithm = algorithm,
lpdir = "max",
solver = me[["sol"]],
method = me[["met"]],
- solverParm = as.data.frame(NA))
+ solverParm = as.data.frame(NA), ...)
}
return(tmp)
diff --git a/R/sysBiolAlg_mtfEasyConstraintClass.R b/R/sysBiolAlg_mtfEasyConstraintClass.R
index f76e49a296085d7d4f63d72958bb28904be74c44..eac646e7be9d96b7a38971440c0d1bd873140fed 100644
--- a/R/sysBiolAlg_mtfEasyConstraintClass.R
+++ b/R/sysBiolAlg_mtfEasyConstraintClass.R
@@ -29,7 +29,8 @@
setClass(Class = "sysBiolAlg_mtfEasyConstraint",
representation(
- maxobj = "numeric"
+ maxobj = "numeric",
+ easyConstraint = "list"
),
contains = "sysBiolAlg"
)
@@ -54,12 +55,16 @@ setMethod(f = "initialize",
rnames = NULL,
pname = NULL,
scaling = NULL,
+ easyConstraint = NULL,
writeProbToFileName = NULL, ...) {
if ( ! missing(model) ) {
if (is.null(wtobj)) {
- tmp <- .generateWT(model, react, lb, ub, ...)
+ tmp <- .generateWT(model, react, lb, ub,
+ algorithm="fbaEasyConstraint",
+ easyConstraint=easyConstraint,
+ ...)
wtobj <- tmp[["obj"]]
}
@@ -233,7 +238,57 @@ setMethod(f = "initialize",
rowNames <- NULL
probName <- NULL
}
-
+
+ #add easyConstraints:
+ if(!is.null(easyConstraint)){
+ if( length(easyConstraint$react) != length(easyConstraint$x)
+ | length(easyConstraint$react) != length(easyConstraint$rtype)
+ ){
+ stop("easyConstraint elements have to have equal lengths")
+ }
+ stopifnot(is.list(easyConstraint$react))
+ stopifnot(is.list(easyConstraint$x))
+ stopifnot(all(easyConstraint$rtype %in% c("F", "L", "U", "D", "E")))
+
+ # setting and checking rlb
+ if(is.null(easyConstraint$lb)){
+ rlower <- c(rlower, rep(0, length(easyConstraint$react)))
+ }else{
+ if(length(easyConstraint$react) != length(easyConstraint$lb)){
+ stop("easyConstraint$lb length has to match length of react argument")
+ }else{
+ stopifnot(is.numeric(easyConstraint$lb))
+ rlower <- c(rlower, easyConstraint$lb)
+ }
+ }
+
+ # setting and checking rub
+ if(is.null(easyConstraint$ub)){
+ rupper <- c(rupper, rep(0, length(easyConstraint$react)))
+ }else{
+ if(length(easyConstraint$react) != length(easyConstraint$ub)){
+ stop("easyConstraint$ub length has to match length of react argument")
+ }else{
+ stopifnot(is.numeric(easyConstraint$ub))
+ rupper <- c(rupper, easyConstraint$ub)
+ }
+ }
+
+ m <- Matrix(0, ncol=nCols, nrow=length(easyConstraint$react))
+
+ for(i in 1:length(easyConstraint$react)){
+ m[i, easyConstraint$react[[i]]] <- easyConstraint$x[[i]]
+ }
+
+
+ LHS <- rbind2(LHS, m)
+ rtype <- c(rtype, easyConstraint$rtype)
+ nRows <- nRows + length(easyConstraint$react)
+ if(!is.null(rowNames)){
+ rowNames <- c(rowNames, paste0("easyConstraint", 1:length(easyConstraint$react)))
+ }
+
+ }
# ---------------------------------------------
# build problem object
@@ -264,6 +319,7 @@ setMethod(f = "initialize",
...)
.Object@maxobj <- as.numeric(maxobj)
+ .Object@easyConstraint <- easyConstraint
if (!is.null(writeProbToFileName)) {
writeProb(problem(.Object),
diff --git a/inst/NEWS.Rd b/inst/NEWS.Rd
index 9eb40dc87d2b1237e4718606d8acc07c0d976755..82b09df06cae065a8d2b1d1c2a024eee7625ca5f 100644
--- a/inst/NEWS.Rd
+++ b/inst/NEWS.Rd
@@ -6,7 +6,16 @@
\newcommand{\CRANpkg}{\href{http://CRAN.R-project.org/package=#1}{\pkg{#1}}}
-\section{Changes in version 1.3.1 2015-08-06}{
+\section{Changes in version 1.3.2 2015-10-21}{
+ \itemize{
+ \item New Algorithm \code{fbaEasyConstraint} and \code{mtfEasyConstraint}
+ implemented. With these new \code{sysBiolAlg}-Classes it is easier to add
+ linear constraints to a model. See package vignette or manual for more details.
+ \item \code{switch} needs the first parameter explicitly named.
+ }
+}
+
+\section{Changes in version 1.3.1 2015-10-02}{
\itemize{
\item \code{rmReact} error, if resulting model had only one gene left, corrected.
\item \code{deadEndMetabolite()} could miss deadEndMetabolites if reactions differ in stoichiometry e.g. 1 and 2.
diff --git a/man/sysBiolAlg_fbaEasyConstraint-class.Rd b/man/sysBiolAlg_fbaEasyConstraint-class.Rd
index 93ed0b56c27efeea52c050467d59a234f712c2e9..2f7ba78ffb17ec34dee9ed3f4319768dfcdfebbc 100644
--- a/man/sysBiolAlg_fbaEasyConstraint-class.Rd
+++ b/man/sysBiolAlg_fbaEasyConstraint-class.Rd
@@ -182,13 +182,6 @@
-
-
-
-
-
-
-
}
\keyword{classes}
diff --git a/vignettes/sybil.Rnw b/vignettes/sybil.Rnw
index a05497414dd54965ec313d6005e634ecea4183e2..7078f2de1456589398fa19df1df1095c3c6f0bb8 100644
--- a/vignettes/sybil.Rnw
+++ b/vignettes/sybil.Rnw
@@ -645,6 +645,53 @@ mod_obj(mtf)
which is of course the same value as for the FBA algorithm.
+% ---------------------------------------------------------------------------- %
+\subsection{Example: Limit carbon intake (easyConstraint)} \label{carbonIntake}
+
+This example is meant to demonstrate the usage of the classes
+\Comp{sysBiolAlg\_fbaEasyConstraint} and \Comp{sysBiolAlg\_mtfEasyConstraint}.
+These two algorithms work like the normal FBA and MTF, but also have the
+ability to easily add additional linear constraints to the linear programm. So
+the $i^{th}$ constraint woult look like follows:
+
+$$\gamma_i \leq v_{r_i} * x_i^\mathrm{T} \leq \delta_i$$
+
+With $\gamma$ and $\delta$ beeing the vectors of upper and lower bounds,
+$v$ is flux vector and $r$ and $x$ are sets of vectors indicating the
+affected reactions and the corresponding coefficients, respectively.
+
+In our example we want to limit the carbon intake on basis of the C-atoms in the
+molekules: we only allow uptake for glucose (6 carbon atoms) and fumarate
+(4 carbon atoms):
+<<>>=
+data(Ec_core)
+optimizeProb(Ec_core)
+lowbnd(Ec_core)[react_id(Ec_core) == "EX_co2(e)"] <- 0
+lowbnd(Ec_core)[react_id(Ec_core) == "EX_fum(e)"] <- -1000
+lowbnd(Ec_core)[react_id(Ec_core) == "EX_glc(e)"] <- -1000
+findExchReact(Ec_core)
+@
+Now we define the additional linear constraint to limit the overall carbon
+uptake to be as high as before ($-10 * 6$ C-atoms of glucose):
+<<>>=
+help("fbaEasyConstraint")
+ec <- list(
+ react=list(match(c("EX_glc(e)", "EX_fum(e)"), react_id(Ec_core))),
+ x=list(c(6, 4)),
+ rtype="L",
+ lb=-60)
+@
+Optimize with the constraint:
+<<print=true>>=
+o <- optimizeProb(Ec_core, algorithm="fbaEasyConstraint", easyConstraint=ec)
+@
+The result shows, that for this model it is more efficient to take up glucose
+than fumarate in order to obtain maximal growth:
+<<print=false>>=
+fluxes(o)[match(c("EX_glc(e)", "EX_fum(e)"), react_id(Ec_core))]
+@
+
+
% ---------------------------------------------------------------------------- %
\subsection{Genetic perturbations} \label{komut}