diff --git a/networkComplexityBigg/qc/compGrowth.R b/networkComplexityBigg/qc/compGrowth.R new file mode 100644 index 0000000000000000000000000000000000000000..f50fd4c2d2d6325bf111f5df6931dad81b2f1f00 --- /dev/null +++ b/networkComplexityBigg/qc/compGrowth.R @@ -0,0 +1,149 @@ +#!/usr/bin/Rscript +library(sybil) +library(sybilSBML) +library(sybilSWITCH) +library(parallel) + +eco <- list() +print(load("~/workspace/models/iAF1260_the-seed-adaption.Rdata")) +eco$seed <- upgradeModelorg(iAFseed) + +load("../models/universalBiGG.ver1.1.Rdata") +eco$bigg <- changeObjFunc(rmReact(uni1.1, setdiff(react_id(uni1.1), modelReactMap$iAF1260)), modelBiomassMapSelection["iAF1260"]) + +#load("../modelsCheck/checkedModels.Rdata") +#eco$biggChecked <- allModels$iAF1260 + + +print(load("../sourceData/seed.mobd.media.Rdata")) +print(load("../sourceData/balazs.media.Rdata")) + + +print(load("../sourceData/seed.mobd.media.Rdata")) +print(load("../sourceData/balazs.media.Rdata")) +media <- list() +media$seed <- c(seed.mobd.media, balazs.media) +media$bigg <- get(load("../convertMedia/mediaBiGG.Rdata")) + +stopifnot(setequal(names(media$seed), names(media$bigg))) +mediaNames <- names(media$seed) +media <- lapply(media, function(x) x[mediaNames]) + +namespace <- c(bigg="bigg", seed="seed")#, biggChecked="bigg") +mnames <- names(eco) +names(mnames) <- names(eco) + + + +#mediaArm <- function(id="", medium){ +# m <- eco[[id]] +# me <- media[[namespace[id]]] +# +# ex <- findExchReact(m) +# lowbnd(m)[react_pos(ex)] <- -10 +# exMap <- gsub("\\[.+$", "", met_id(ex)) +# +# +# sp <- suggestedArmSolverSettings() +# sp$CPX_PARAM_EPINT <- 1e-9 +# sp$CPX_PARAM_EPRHS <- 1e-9 +# +# ar <- setdiff(react_pos(ex), react_pos(ex)[exMap %in% me[[medium]]]) +# +# arm <- optimizeProb(m, +# additionalReact=ar, +# biomassThreshold=0.1, +# algorithm="arm", +# solverParm=sp) +# +# fluxes <- getArmReactionFluxes(m, arm)[react_pos(ex)] +# names(fluxes) <- met_id(ex) +# browser() +# fluxes[fluxes < -1e-9] +#} + +#add <- mediaArm("bigg", "C_cpd00751") + + + +df <- as.data.frame(lapply(mnames, function(e){ + m <- eco[[e]] + + ex <- findExchReact(m) + + react <- list(react_pos(ex)) + mid <- namespace[e] + r <- rep(list(react_pos(ex)), length(media[[mid]])) + + exMap <- gsub("\\[.+$", "", met_id(ex)) + lb <- lapply(media[[mid]], function(x){ + v <- rep(0, length(exMap)) + v[exMap %in% x] <- -10 + print(x[!x %in% exMap]) + v + }) + + ub <- rep(list(uppbnd(m)[react_pos(ex)]), length(media[[mid]])) + opt <- optimizer(m, react=r, lb=lb, ub=ub, verboseMode=1) + obj <- opt$obj + obj[checkSolStat(opt$stat)] <- 0 + return(round(obj, digits=6)) +})) + +rownames(df) <- mediaNames + +growth <- cbind(name=rownames(df), as.data.frame(df > 1e-6)) + + +dfg1 <- data.frame(expand.grid(seed=c(T, F), bigg=c(T, F))) +dfg1$sums <- apply(dfg1, 1, function(x) sum(growth[,"seed"] == x[1] & growth[,"bigg"] == x[2] )) +print(dfg1) + +stopifnot(sum(dfg1[2:3, 3])==0) + +#dfg1 <- data.frame(expand.grid(seed=c(T, F), bigg=c(T, F))) +#dfg1$sums <- apply(dfg1, 1, function(x) sum(growth[,"seed"] == x[1] & growth[,"biggChecked"] == x[2] )) +#print(dfg1) + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +