# blockedReact.R # FBA and friends with R. # # Copyright (C) 2010-2014 Gabriel Gelius-Dietrich, Dpt. for Bioinformatics, # Institute for Informatics, Heinrich-Heine-University, Duesseldorf, Germany. # All right reserved. # Email: geliudie@uni-duesseldorf.de # # This file is part of sybil. # # Sybil is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # Sybil is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with sybil. If not, see <http://www.gnu.org/licenses/>. ################################################ # Function: blockedReact # # This function finds blocked reactions in a specified IO-environment. # # blockedReact <- function(model, tol = SYBIL_SETTINGS("TOLERANCE"), exex = TRUE, fld = FALSE, retOptSol = FALSE, verboseMode = 2, ... ) { .Deprecated(new = "fluxVar") if (!is(model, "modelorg")) { stop("needs an object of class modelorg!") } modIrr <- ifelse(is(model, "modelorg_irrev"), TRUE, FALSE) intReact <- NA # remove exchange reactions from analysis if (isTRUE(exex)) { exchReact <- findExchReact(model) ex <- react_pos(exchReact) intReact <- 1:react_num(model) intReact <- intReact[-ex] if (length(intReact) < 1) { stop("model contains no internal reactions!") } } else { intReact <- 1:react_num(model) } # memory for results nObj <- ifelse(isTRUE(modIrr), length(intReact), 2 * length(intReact)) blocked_react <- logical(react_num(model)) if (isTRUE(retOptSol)) { obj <- numeric(nObj) ok <- integer(nObj) stat <- integer(nObj) if (isTRUE(fld)) { flux <- Matrix::Matrix(0, nrow = react_num(model), ncol = nObj) } else { flux <- NA } } obj_mod <- obj_coef(model) obj_coef(model) <- as.integer(rep(0, react_num(model))) lpmod <- sysBiolAlg(model, algorithm = "fv", tol = tol, ...) #------------------------------------------------------------------------------# # finding blocked reactions # #------------------------------------------------------------------------------# if (verboseMode > 0) { message("calculating blocked reactions ...") } obj_max <- 0 obj_min <- 0 if (verboseMode > 1) { progr <- .progressBar() } for (i in seq(along = intReact)) { solpl <- ifelse(isTRUE(modIrr), i+1, 2*i) # if (uppbnd(model)[intReact[i]] > 0) { #print(intReact[i]) sol <- optimizeProb(lpmod, lpdir = "max", react = intReact[i], obj_coef = 1) obj_max <- sol$obj if (isTRUE(retOptSol)) { obj[(solpl-1)] <- sol$obj ok[(solpl-1)] <- sol$ok stat[(solpl-1)] <- sol$stat if (isTRUE(fld)) { flux[,(solpl-1)] <- sol$fluxes } } # } # else { # obj_max <- 0 # if (isTRUE(retOptSol)) { # obj[(i*2-1)] <- 0 # ok[(i*2-1)] <- NA # stat[(i*2-1)] <- NA # if (isTRUE(fld)) { # flux[,(i*2-1)] <- NA # } # } # } # if (lowbnd(model)[intReact[i]] < 0) { if (!is(model, "modelorg_irrev")) { sol <- optimizeProb(lpmod, lpdir = "min", react = intReact[i], obj_coef = 1) obj_min <- sol$obj if (isTRUE(retOptSol)) { obj[solpl] <- sol$obj ok[solpl] <- sol$ok stat[solpl] <- sol$stat if (isTRUE(fld)) { flux[,solpl] <- sol$fluxes } } } # } # else { # obj_min <- 0 # if (isTRUE(retOptSol)) { # obj[(i*2)] <- 0 # ok[(i*2)] <- NA # stat[(i*2)] <- NA # if (isTRUE(fld)) { # flux[,(i*2)] <- NA # } # } # } #print(paste("i", i)) #print(paste("max", obj_max)) #print(paste("min", obj_min)) #print(" ") if (is(model, "modelorg_irrev")) { blocked_react[intReact[i]] <- ifelse(abs(obj_max) < tol, TRUE, FALSE) } else { if ( (abs(obj_max) < tol) && (abs(obj_min) < tol) ) { blocked_react[intReact[i]] <- TRUE } else { blocked_react[intReact[i]] <- FALSE } } if (verboseMode > 1) { progr <- .progressBar(i, length(intReact), progr) } } #print(nObj) #print(obj) #print(ok) #print(blocked_react) if (isTRUE(retOptSol)) { optsol <- new("optsol_blockedReact", mod_id = mod_id(model), mod_key = mod_key(model), solver = solver(problem(lpmod)), method = method(problem(lpmod)), algorithm = algorithm(lpmod), num_of_prob = as.integer(nObj), lp_num_cols = nc(lpmod), lp_num_rows = nr(lpmod), lp_obj = as.numeric(obj), lp_ok = as.integer(ok), lp_stat = as.integer(stat), lp_dir = factor(rep(c("max", "min"), length(intReact))), obj_coef = obj_mod, fldind = fldind(lpmod), fluxdist = fluxDistribution(flux), blocked = blocked_react, react = new("reactId", mod_id = mod_id(model), mod_key = mod_key(model), pnt = intReact, id = react_id(model)[intReact]) ) } else { optsol <- blocked_react } return(optsol) }