Upload New File

networkComplexityBigg/qc/testGrowth.R

+#!/usr/bin/Rscript
+
+
+library(sybil)
+SYBIL_SETTINGS("SOLVER", "cplexAPI")
+
+
+versions <- c(v1.1="1.1", v1.2="1.2", v2="2")
+df <- lapply(versions, function(ver){
+
+	print(load(paste0("../models/universalBiGG.ver", ver, ".Rdata")))
+	uni <- get(paste0("uni", ver))
+
+
+	react <- rep(list(1:react_num(uni)), length(modelBiomassMap))
+	ub <- lapply(modelReactMap, function(x){
+		v <- uppbnd(uni)
+		v[-match(x, react_id(uni))] <- 0
+		v
+	})
+	lb <- lapply(modelReactMap, function(x){
+		v <- lowbnd(uni)
+		v[-match(x, react_id(uni))] <- 0
+		v
+	})
+
+	obj <- lapply(modelBiomassMapSelection, function(x){
+		v <- rep(0, react_num(uni))
+		v[match(x, react_id(uni))] <- 1
+		v
+	})
+	
+	
+	if(ver=="2"){
+		for(n in names(modelReactMap)){
+			offF <- match(modelRmEnergyCycles[[n]]$removedFwd, react_id(uni))
+			offB <- match(modelRmEnergyCycles[[n]]$removedBwd, react_id(uni))
+			ub[[n]][offF] <- 0
+			lb[[n]][offB] <- 0
+		}
+	}
+	
+	opt <- optimizer(uni, react=react, lb=lb, ub=ub, obj_coef=obj)
+	
+	# change to generic BOF
+	obj <- lapply(modelBiomassMapSelection, function(x){
+		v <- rep(0, react_num(uni))
+		v[match("GENERAL_BOF", react_id(uni))] <- 1
+		v
+	})
+	optG <- optimizer(uni, react=react, lb=lb, ub=ub, obj_coef=obj)
+	
+	return(data.frame(model=names(modelReactMap), growth=opt$obj, growthGeneral=optG$obj, hasBOF=!is.na(modelBiomassMapSelection)))
+})
+
+
+versionGrowth <- do.call(cbind, df)
+rownames(versionGrowth) <- versionGrowth$v1.1.model
+
+cat(">>> following models growth is blocked after mass balance:\n")
+print(as.character(with(versionGrowth, v1.1.model[v1.1.growth>1e-6 & v1.2.growth<1e-6])))
+
+
+cat(">>> following models have lost their BOF because of mass balancing:\n")
+print(as.character(with(versionGrowth, v1.1.model[v1.1.hasBOF & !v1.2.hasBOF])))
+
+pdf("biomassProduction.pdf", width=9, paper="a4r")
+for(ver in versions){
+	for(bof in c(growth="growth", growthGeneral="growthGeneral")){
+		barplot(versionGrowth[,paste0("v", ver, ".", bof)], las=2, 
+			main=paste0("Biomass growth: ", paste0("v", ver, ".", bof)),
+			names.arg=versionGrowth$v1.1.model,
+			border=NA, cex.names=0.7)
+	}
+}
+dev.off()
+
+
+
+
+
+
+
+
+
+
+
+save(versionGrowth, file="versionGrowth.Rdata")
+
+
+
+
+
+
+
+
+