R/addReact.R

In sybil: Efficient Constrained Based Modelling

#  addReact.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: addReact
#
#
# The function addReact() is inspired by the function
# addReaction() contained in the COBRA Toolbox.
# The algorithm is (more or less) the same.


setMethod("addReact", signature(model = "modelorg"),
	function(model,
			id,
			met,
			Scoef,
			reversible = FALSE,
			lb = 0,
			ub = SYBIL_SETTINGS("MAXIMUM"),
			obj = 0,
			subSystem = NA,
			gprAssoc = NA,
			reactName = NA,
			metName = NA,
			metComp = NA) {

  
    # ------------------------------------------------------------------------ #
    # check arguments
    # ------------------------------------------------------------------------ #

    if (!is(model, "modelorg")) {
        stop("needs an object of class modelorg!")
    }
	
	stopifnot(checkVersion(model))

    if (length(met) != length(Scoef)) {
        stop("arguments 'met' and 'Scoef' must have the same length")
    }

    if (length(id) > 1) {
        stop("add/change one reaction")
    }

    if ( ( (ub > 0) && (lb < 0) ) && (!isTRUE(reversible)) ) {
        Crev <- TRUE
        warning(paste("'lb' and 'ub' are signed different,",
                      "set reversible to 'TRUE'"))
    }
    else {
       Crev <- reversible
    }


    # ------------------------------------------------------------------------ #
    # check, if we need to add columns and/or rows
    # ------------------------------------------------------------------------ #

    # reaction
    colInd <- match(id, react_id(model))
    addCol <- FALSE
    nCols  <- react_num(model)
    if (is.na(colInd)) {
        # new reaction
        colInd <- react_num(model) + 1
        addCol <- TRUE
        nCols  <- nCols + 1
    }


    # metabolites
    rowInd <- match(met, met_id(model))
    
    newM     <- which(is.na(rowInd))
    nRows    <- met_num(model)          # number of rows in the model
    nNewRows <- length(newM)            # number of new rows
    addRow   <- FALSE
    
    for (i in seq(along = newM)) {
        addRow   <- TRUE
        nRows    <- nRows + 1
        rowInd[newM[i]] <- nRows
    }
    

    if ( (isTRUE(addCol)) || (isTRUE(addRow)) ) {    

        # -------------------------------------------------------------------- #
        # make a new model
        # -------------------------------------------------------------------- #

        # -------------------------------------------------------------------- #
        # data structures

        newmet_num      <- met_num(model)
        newmet_id       <- met_id(model)
        newmet_name     <- met_name(model)
        newmet_comp     <- met_comp(model)
        newmet_single   <- met_single(model)
        newmet_de       <- met_de(model)

        newreact_num    <- react_num(model)
        newreact_rev    <- react_rev(model)
        newreact_id     <- react_id(model)
        newreact_name   <- react_name(model)
        newreact_single <- react_single(model)
        newreact_de     <- react_de(model)
        newlowbnd       <- lowbnd(model)
        newuppbnd       <- uppbnd(model)
        newobj_coef     <- obj_coef(model)

        newgprRules     <- gprRules(model)
        newgenes        <- genes(model)
        newgpr          <- gpr(model)
        newallGenes     <- allGenes(model)
        newrxnGeneMat   <- rxnGeneMat(model)
        newsubSys       <- subSys(model)

        newS            <- S(model)
        
        newMetAttr <- met_attr(model)
        newReactAttr <- react_attr(model)
        newCompAttr <- comp_attr(model)
        newModAttr <- mod_attr(model)

    
        if (isTRUE(addRow)) {
            
            # new number of metabolites
            newmet_num  <- nRows
            
            # new metabolite id's
            newmet_id   <- append(met_id(model), met[newM])

            # new metabolite names
            if (any(is.na(metName))) {
                newmet_name <- append(met_name(model), met[newM])
            }
            else {
                newmet_name <- append(met_name(model), metName[newM])
            }

            # new metabolite compartments
            if (any(is.na(metComp))) {
                newmet_comp <- append(met_comp(model), rep(NA, nNewRows))
            }
            else {
                if (is(metComp, "numeric")) {
                    newmet_comp <- append(met_comp(model), metComp[newM])
                }
                else {
                    newmet_comp <- append(met_comp(model),
                                          match(metComp[newM],
                                                mod_compart(model)))
                }
            }

            # singleton and dead end metabolites (not checked!)
            newmet_single <- append(met_single(model), rep(NA, nNewRows))
            newmet_de     <- append(met_de(model),     rep(NA, nNewRows))

            # new rows in stoichiometric matrix
            newRows <- Matrix::Matrix(0,
                                      nrow = nNewRows,
                                      ncol = react_num(model))
            newS <- rbind(newS, newRows)
            
            # new met attrs
            if(ncol(newMetAttr) > 0){
            	newMetAttr[nrow(newMetAttr)+1:nNewRows, ] <- NA
            }
        }
    
        if (isTRUE(addCol)) {                        # we add at most one column
            # new number of reactions
            newreact_num  <- nCols
            
            # new reaction id
            newreact_id   <- append(react_id(model), id)

            # new reaction name
            if (is.na(reactName)) {
                newreact_name <- append(react_name(model), id)
            }
            else {
                newreact_name <- append(react_name(model), reactName)
            }

            # reaction contains singleton or dead end metabolites (not checked!)
            newreact_single <- append(react_single(model), NA)
            newreact_de     <- append(react_de(model),     NA)

            # reversibility, lower and upper bounds, objective coefficient
            newreact_rev <- append(react_rev(model), Crev)
            newlowbnd    <- append(lowbnd(model),    lb)
            newuppbnd    <- append(uppbnd(model),    ub)
            newobj_coef  <- append(obj_coef(model),  obj)

            # new column in stoichiometric matrix
            newS <- cbind(newS, rep(0, nrow(newS)))
            
            # new react Attr
            # only one new row, /bc we can only add one reaction a time.
            if(ncol(newReactAttr) > 0){
            	newReactAttr[nrow(newReactAttr)+1, ] <- NA
            }
            
            # subsystems
            if (any(is.na(subSystem))) {
            	ss <- subSys(model)
            	if(ncol(ss)==0){ # if no subSys defined, rbind (see else) failed
            		dim(ss) <- c(nrow(ss)+1, ncol(ss))
            		newsubSys <- ss
            	}
            	else {
            		newsubSys <- rbind(ss, rep(FALSE, ncol(subSys(model))))
            	}
            }
            else {
                if (is(subSystem, "logical")) {
                    newsubSys <- rbind(subSys(model), subSystem)
                }
                else {
                    nSubsRow  <- colnames(subSys(model)) %in% subSystem
                    newsubSys <- rbind(subSys(model), nSubsRow)
                }
            }


            # gpr association
            if (ncol(rxnGeneMat(model)) > 0) {
                newrxnGeneMat   <- rbind(rxnGeneMat(model),
                                         rep(FALSE, ncol(rxnGeneMat(model))))
            }
            else { #if (nrow(rxnGeneMat(model)) > 0) {
                newrxnGeneMat <- rxnGeneMat(model)
                dim(newrxnGeneMat) <- c(nrow(newrxnGeneMat)+1,
                                        ncol(newrxnGeneMat))
            }
            # do above else always.

            if ( (is.na(gprAssoc)) || (gprAssoc == "") ) {
                if ((length(gprRules(model)) > 0)) {
                    newgprRules     <- append(gprRules(model), "")
                    newgenes        <- append(genes(model), "")
                    newgpr          <- append(gpr(model), "")
                }
            }
            else {
                gene_rule <- .parseBoolean(gprAssoc)
            
                geneInd <- match(gene_rule$gene, allGenes(model))
            
                # indices of new genes
                new_gene <- which(is.na(geneInd))

                # if we have new gene(s), add a column in rxnGeneMat and
                # gene name(s) to allGenes
                if (length(new_gene) > 0) {
                    newallGenes <- append(allGenes(model),
                                          gene_rule[["gene"]][new_gene])

                    # update geneInd
                    geneInd <- match(gene_rule[["gene"]], newallGenes)

                    # if we have an empty modelorg object, we need to
                    # initialize rxnGeneMat
                    if (ncol(newrxnGeneMat) == 0) {
                        newrxnGeneMat <- Matrix::Matrix(FALSE,
                                                        nCols, max(geneInd))
                    }
                    else {
                        for (i in seq(along = gene_rule[["gene"]][new_gene])) {
	    					newrxnGeneMat <- cbind(newrxnGeneMat,
		    								   rep(FALSE, nrow(newrxnGeneMat)))
			    		}
					}
                }

                # rxnGeneMat
                newrxnGeneMat[nCols, geneInd] <- TRUE
 
                # new rule
                newgpr <- append(gpr(model), gprAssoc)

                # genes per reaction
                newgenes <- append(genes(model), list(gene_rule$gene))
                newrule  <- gene_rule$rule
				
				# not needed for modelorg version 2.0
#                for (j in 1 : length(geneInd)) {
#                    pat  <- paste("x(", j, ")", sep = "")
#                    repl <- paste("x[", geneInd[j], "]", sep = "")
#    
#                    newrule <- gsub(pat, repl, newrule, fixed = TRUE)
#                }

                newgprRules <- append(gprRules(model), newrule)
            }
        }
        
        # values for stoichiometric matrix
        newS[ , colInd]      <- 0    
        newS[rowInd, colInd] <- Scoef    
        
#        for (i in seq(along = rowInd)) {
#            newS[rowInd[i], colInd] <- Scoef[i]
#        }
    
        # -------------------------------------------------------------------- #
        # new model
        # -------------------------------------------------------------------- #
        
        if (is(model, "modelorg_irrev")) {
            mod_out <- modelorg_irrev(mod_id(model), mod_name(model))
            irrev(mod_out)     <- TRUE
            matchrev(mod_out)  <- append(matchrev(model), 0L)
            
            revReactId <- as.integer(max(irrev2rev(model))+1)
            irrev2rev(mod_out) <- append(irrev2rev(model), revReactId)
            rev2irrev(mod_out) <- rbind(rev2irrev(model), c(nCols, nCols))
        }
        else {
            mod_out <- modelorg(mod_id(model), mod_name(model))
        }
        
        mod_desc(mod_out)    <- mod_desc(model)
        mod_compart(mod_out) <- mod_compart(model)


        met_num(mod_out)      <- as.integer(newmet_num)
        met_id(mod_out)       <- newmet_id
        met_name(mod_out)     <- newmet_name
        met_comp(mod_out)     <- as.integer(newmet_comp)
        met_single(mod_out)   <- newmet_single
        met_de(mod_out)       <- newmet_de

        react_num(mod_out)    <- as.integer(newreact_num)
        react_rev(mod_out)    <- newreact_rev
        react_id(mod_out)     <- newreact_id
        react_name(mod_out)   <- newreact_name
        react_single(mod_out) <- newreact_single
        react_de(mod_out)     <- newreact_de
        lowbnd(mod_out)       <- newlowbnd
        uppbnd(mod_out)       <- newuppbnd
        obj_coef(mod_out)     <- newobj_coef

        gprRules(mod_out)     <- newgprRules
        genes(mod_out)        <- newgenes
        gpr(mod_out)          <- newgpr
        allGenes(mod_out)     <- newallGenes
        rxnGeneMat(mod_out)   <- newrxnGeneMat
        subSys(mod_out)       <- newsubSys

        S(mod_out)            <- newS
        
        react_attr(mod_out) <- newReactAttr
        met_attr(mod_out) <- newMetAttr
        comp_attr(mod_out) <- newCompAttr
        mod_attr(mod_out) <- newModAttr
        

    }
    else {
    
        # -------------------------------------------------------------------- #
        # modify old model
        # -------------------------------------------------------------------- #
        
        mod_out <- model
        
        react_rev(mod_out)[colInd] <- Crev
        lowbnd(mod_out)[colInd]    <- lb
        uppbnd(mod_out)[colInd]    <- ub
        obj_coef(mod_out)[colInd]  <- obj
        S(mod_out)[ , colInd]      <- 0
        S(mod_out)[rowInd, colInd] <- Scoef

    }
    
    
    check <- validObject(mod_out, test = TRUE)

    if (check != TRUE) {
        msg <- paste("Validity check failed:", check, sep = "\n    ")
        warning(msg)
    }

    return(mod_out)

})