findMDCFlux: find Manhaten-Distance-Closest Flux

In sybilEFBA: Using Gene Expression Data to Improve Flux Balance Analysis Predictions

Description

Given a wildtype flux (some fluxes can be absent i.e NA), find FBA solution of network given by model, such that Sum(|v_i-v_wt|) is minimal (like lMOMA but with option to exclude some reactions)

Usage

findMDCFlux(model, wtflux, objVal = NA, pct_objective=100, 
        lpdir = SYBIL_SETTINGS("OPT_DIRECTION"), 
       solver = SYBIL_SETTINGS("SOLVER"), method = SYBIL_SETTINGS("METHOD"), 
       solverParm = data.frame(CPX_PARAM_EPRHS = 1e-06), verboseMode = 2)

Arguments

model	An object of class modelorg.
wtflux	desired flux distribution, when wtflux is NA its constraint won't be included in LP.
pct_objective	Biomass will be garnteed to be at least this value multiplied by max biomass calculated using standard FBA. Values are 0 to 100.
objVal	lower bound of objective function, if not set it will be set to maximum biomass.
lpdir	Character value, direction of optimisation. Can be set to "min" or "max". Default: SYBIL_SETTINGS("OPT_DIRECTION").
solver	Single character string giving the solver package to use. See SYBIL_SETTINGS for possible values. Default: SYBIL_SETTINGS("SOLVER").
method	Single character string giving the method the desired solver has to use. SYBIL_SETTINGS for possible values. Default: SYBIL_SETTINGS("METHOD").
solverParm	A named data frame or list containing parameters for the specified solver. Parameters can be set as data frame or list: solverParm = list(parm1 = val1, parm2 = val2) with parm1 and parm2 being the names of two different parameters and val1 and val2 the corresponding values. For possible parameters and values see the documentation of the used solver package (e.g. glpkAPI). Default: SYBIL_SETTINGS("SOLVER_CTRL_PARM").
verboseMode	An integer value indicating the amount of output to stdout: 0: nothing, 1: status messages, 2: like 1 plus with more details, 3: generates files of the LP problem. Default: 2.

Value

return list with slot mdcflx containing the new calculated fluxes and status returned from solver.

Author(s)

Abdelmoneim Amer Desouki

See Also

modelorg, lmoma

Examples

## Not run: 
## The function is currently defined as
function (model, wtflux, objVal = NA, lpdir = SYBIL_SETTINGS("OPT_DIRECTION"), 
    solver = SYBIL_SETTINGS("SOLVER"), method = SYBIL_SETTINGS("METHOD"), 
    solverParm = data.frame(CPX_PARAM_EPRHS = 1e-06), verboseMode = 2) 
{
    if (!is(model, "modelorg")) {
        stop("needs an object of class modelorg!")
    }
    if (is.na(objVal)) {
        sol = optimizeProb(model, solver = solver, method = method, 
            solverParm = solverParm)
        objVal = lp_obj(sol)
    }
    out <- FALSE
    nc <- react_num(model)
    nr <- met_num(model)
    nd <- sum(!is.na(wtflux))
    absMAX <- SYBIL_SETTINGS("MAXIMUM")
    nRows = nr + 2 * nd + 1
    nCols = nc + 2 * nd
    LHS <- Matrix::Matrix(0, nrow = nRows, ncol = nCols, sparse = TRUE)
    LHS[1:nr, 1:(nc)] <- S(model)
    ii = matrix(c((nr + 1):(nr + nd), which(!is.na(wtflux))), 
        ncol = 2)
    LHS[ii] <- 1
    if (nd == 1) {
        LHS[(nr + 1):(nr + nd), (nc + 1):(nc + nd)] <- 1
    }
    else {
        diag(LHS[(nr + 1):(nr + nd), (nc + 1):(nc + nd)]) <- 1
    }
    ii = matrix(c((nr + nd + 1):(nr + 2 * nd), which(!is.na(wtflux))), 
        ncol = 2)
    LHS[ii] <- -1
    if (nd == 1) {
        LHS[(nr + nd + 1):(nr + 2 * nd), (nc + nd + 1):(nc + 
            2 * nd)] <- 1
    }
    else {
        diag(LHS[(nr + nd + 1):(nr + 2 * nd), (nc + nd + 1):(nc + 
            2 * nd)]) <- 1
    }
    LHS[(nr + 2 * nd + 1), 1:nc] <- obj_coef(model)
    if (verboseMode > 2) 
        print(sprintf("nrows:%d, ncols=%d, nd=%d,nc=%d,nr=%d", 
            nRows, nCols, nd, nc, nr))
    lower <- c(lowbnd(model), rep(0, 2 * nd))
    upper <- c(uppbnd(model), rep(absMAX, 2 * nd))
    rlower <- c(rep(0, nr), wtflux[!is.na(wtflux)], -wtflux[!is.na(wtflux)], 
        objVal)
    rupper <- c(rep(0, nr), rep(absMAX, 2 * nd + 1))
    cobj <- c(rep(0, nc), rep(1, 2 * nd))
    cNames = paste(c(rep("x", nc), rep("dp", nd), rep("dn", nd)), 
        c(1:nc, which(!is.na(wtflux)), which(!is.na(wtflux))), 
        sep = "_")
    switch(solver, glpkAPI = {
        out <- vector(mode = "list", length = 5)
        prob <- glpkAPI::initProbGLPK()
        rtype <- c(rep(glpkAPI::GLP_FX, nr), rep(glpkAPI::GLP_LO, 
            2 * nd))
        if (lpdir == "max") {
            rtype <- c(rtype, glpkAPI::GLP_LO)
        } else {
            rtype <- c(rtype, glpkAPI::GLP_UP)
        }
        TMPmat <- as(LHS, "TsparseMatrix")
        out[[1]] <- glpkAPI::addRowsGLPK(prob, nrows = nRows)
        outj <- glpkAPI::addColsGLPK(prob, ncols = nCols)
        mapply(setColNameGLPK, j = c(1:nCols), cname = cNames, 
            MoreArgs = list(lp = prob))
        glpkAPI::setObjDirGLPK(prob, glpkAPI::GLP_MIN)
        out[[2]] <- glpkAPI::loadMatrixGLPK(prob, length(TMPmat@x), 
            TMPmat@i + 1, TMPmat@j + 1, TMPmat@x)
        out[[3]] <- glpkAPI::setColsBndsObjCoefsGLPK(prob, c(1:nCols), 
            lower, upper, cobj)
        out[[4]] <- glpkAPI::setRowsBndsGLPK(prob, c(1:nRows), 
            rlower, rupper, rtype)
        parm <- sapply(dimnames(solverParm)[[2]], function(x) eval(parse(text = x)))
        val <- solverParm[1, ]
        if (method == "interior") {
            glpkAPI::setInteriorParmGLPK(parm, val)
            out[[5]] <- TRUE
        } else {
            glpkAPI::setSimplexParmGLPK(parm, val)
            out[[5]] <- TRUE
        }
        if (verboseMode > 2) {
            fname = format(Sys.time(), "glpk_ManhatDist_%Y%m%d_%H%M.lp")
            print(sprintf("write problem: %s/%s", getwd(), fname))
            writeLPGLPK(prob, fname)
            print("Solving...")
        }
        lp_ok <- glpkAPI::solveSimplexGLPK(prob)
        lp_obj <- glpkAPI::getObjValGLPK(prob)
        lp_stat <- glpkAPI::getSolStatGLPK(prob)
        if (is.na(lp_stat)) {
            lp_stat <- lp_ok
        }
        lp_fluxes <- glpkAPI::getColsPrimGLPK(prob)
    }, cplexAPI = {
        out <- vector(mode = "list", length = 4)
        prob <- openProbCPLEX()
        out <- setIntParmCPLEX(prob$env, CPX_PARAM_SCRIND, CPX_OFF)
        chgProbNameCPLEX(prob$env, prob$lp, "ManhatenDist cplex")
        setObjDirCPLEX(prob$env, prob$lp, CPX_MIN)
        rtype <- c(rep("E", nr), rep("G", 2 * nd), "G")
        TMPmat <- as(LHS, "TsparseMatrix")
        out[[1]] <- newRowsCPLEX(prob$env, prob$lp, nRows, rlower, 
            rtype)
        out[[2]] <- newColsCPLEX(prob$env, prob$lp, nCols, cobj, 
            lower, upper)
        out[[3]] <- chgCoefListCPLEX(prob$env, prob$lp, length(TMPmat@x), 
            TMPmat@i, TMPmat@j, TMPmat@x)
        parm <- sapply(dimnames(solverParm)[[2]], function(x) eval(parse(text = x)))
        out[[4]] <- setDblParmCPLEX(prob$env, parm, solverParm)
        if (verboseMode > 2) {
            fname = format(Sys.time(), "Cplex_ManhatDist_%Y%m%d_%H%M.lp")
            print(sprintf("write problem: %s/%s", getwd(), fname))
            writeProbCPLEX(prob$env, prob$lp, fname)
            print("Solving...")
        }
        lp_ok <- lpoptCPLEX(prob$env, prob$lp)
        lp_obj <- getObjValCPLEX(prob$env, prob$lp)
        lp_stat <- getStatCPLEX(prob$env, prob$lp)
        if (is.na(lp_stat)) {
            lp_stat <- lp_ok
        }
        lp_fluxes <- getProbVarCPLEX(prob$env, prob$lp, 0, nCols - 
            1)
    }, {
        wrong_type_msg(solver)
    })
    optsol <- list(ok = lp_ok, obj = lp_obj, stat = lp_stat, 
        fluxes = lp_fluxes, mdcflx = lp_fluxes[1:nc], wtflx = wtflux)
    return(optsol)
  }

## End(Not run)

sybilEFBA documentation built on May 29, 2017, 9:35 a.m.