Upload New File

networkComplexityBigg/qc/missingForGrowth.R

+#!/usr/bin/Rscript
+library(methods)
+library(sybil)
+library(cplexAPI)
+library(sybilSWITCH)
+library(parallel)
+source("../envirDist/runFBA.R")
+source("../envirDist/initData.R")
+
+
+uptakeLowBnd <- -10
+optRounds <- 3
+bioThreshold <- 1e-3
+
+print(load("universalBiGG.ver2.Rdata"))
+uni <- uni2
+
+print(load("mediaBiGGwithRandom.Rdata"))
+
+#mediaBiGG <- lapply(mediaBiGG, function(x) union(x, "ni2"))
+
+ex <- findExchReact(uni)
+exOnly <- ex[grep("^EX_", react_id(ex))]
+exMap <- react_pos(exOnly)
+names(exMap) <- gsub("\\[\\w\\]$", "", met_id(exOnly))
+
+# setting media constraints
+mediaLowBounds <- lapply(mediaBiGG, function(m){
+	lb <- lowbnd(uni)
+	lb[exMap] <- 0
+	lb[exMap[names(exMap) %in% m]] <- uptakeLowBnd
+	return(lb)
+})
+
+modelReactMapById <- lapply(modelReactMap, function(x){
+	union(match(x, react_id(uni)), exMap)
+})
+
+modelBiomassMapSelectionById <- match(modelBiomassMapSelection, react_id(uni))
+names(modelBiomassMapSelectionById) <- names(modelBiomassMapSelection)
+
+modelRmEnergyCyclesById <- lapply(modelRmEnergyCycles, function(x){
+	lapply(x, function(y) match(y, react_id(uni)))
+})
+
+
+
+mn <- "iAPECO1_1312"
+#mn <- "iB21_1397"
+
+mediaName <- "ArgonneNMSMedia"
+
+
+mod <- changeObjFunc(uni, modelBiomassMapSelection[mn])
+
+uppbnd(mod)[modelRmEnergyCyclesById[[c(mn, "fwd")]]] <- 0
+lowbnd(mod)[modelRmEnergyCyclesById[[c(mn, "bwd")]]] <- 0
+
+lowbnd(mod)[react_pos(ex)[grep("^EX_", react_id(ex))]] <- -10
+
+sp <- suggestedArmSolverSettings(threads=7, timelimit=120, workMemLimit=10, treeMemLimit=10)
+sp$CPX_PARAM_EPINT <- 1e-9
+sp$CPX_PARAM_EPRHS <- 1e-9
+sp$CPX_PARAM_MIPEMPHASIS <- CPX_MIPEMPHASIS_FEASIBILITY
+
+
+cat(">>> missing for full model growth\n")
+opt <- optimizeProb(mod, algorithm="arm",
+						biomassThreshold=0.001,
+						additionalReact=exMap[setdiff(names(exMap), mediaBiGG[[mediaName]])])
+print(getMeanStatus(lp_stat(opt)))
+fluxes <- getArmReactionFluxes(mod, res=opt)
+fl <- rep(F, length(fluxes))
+fl[fluxes <= -1e-6] <- T
+fl[-exMap] <- F
+fl[lowbnd(mod)==0] <- F
+fl[exMap[mediaBiGG[[mediaName]]]] <- F
+print(rbind(react_id(uni)[fl], fluxes[fl]))
+cat(">>> missing for sub model growth\n")
+modfull <- mod
+lowbnd(mod)[!react_id(mod) %in% modelReactMap[[mn]]] <- 0
+uppbnd(mod)[!react_id(mod) %in% modelReactMap[[mn]]] <- 0
+opt <- optimizeProb(mod, algorithm="arm",
+						biomassThreshold=0.001,
+						additionalReact=exMap[setdiff(names(exMap), mediaBiGG[[mediaName]])])
+print(getMeanStatus(lp_stat(opt)))
+fluxes <- getArmReactionFluxes(mod, res=opt)
+fl <- rep(F, length(fluxes))
+fl[fluxes <= -1e-6] <- T
+fl[-exMap] <- F
+fl[lowbnd(mod)==0] <- F
+fl[exMap[mediaBiGG[[mediaName]]]] <- F
+print(rbind(react_id(uni)[fl], fluxes[fl]))
+
+
+
+#lb <- mediaLowBounds[[mediaName]]
+#lb[modelRmEnergyCyclesById[[c(mn, "bwd")]]] <- 0
+#ub <- uppbnd(mod)
+#ub[modelRmEnergyCyclesById[[c(mn, "fwd")]]] <- 0
+
+#print(optimizeProb(mod, react=1:react_num(mod), lb=lb, ub=ub))
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+