Addiing files to run model building on sequoia

Author: cshenry

CarbonSources/GapfillGenomeModels.py

@@ -0,0 +1,81 @@
+import platform
+import sys
+import json
+import cobra
+import cplex
+import re
+import os
+from os.path import exists
+import logging
+from configparser import ConfigParser
+config = ConfigParser()
+config.read("config.cfg")
+paths = config.get("script", "syspaths").split(";")
+for path in paths:
+    sys.path.append(path)
+import cobrakbase
+from escher import Builder
+from optlang.symbolics import Zero, add
+from modelseedpy import MSPackageManager, MSGapfill, MSGrowthPhenotypes, MSModelUtil, MSATPCorrection
+from cobrakbase.core.kbasefba.newmodeltemplate_builder import NewModelTemplateBuilder
+from annotation_ontology_api.annotation_ontology_apiServiceClient import annotation_ontology_api
+from modelseedpy.helpers import get_template
+import pandas as pd
+
+kbase_api = cobrakbase.KBaseAPI()
+anno_api = annotation_ontology_api()
+glc_o2_atp_media = kbase_api.get_from_ws("Glc.O2.atp",94026)
+gmm = kbase_api.get_from_ws("Carbon-D-Glucose","KBaseMedia")
+template = kbase_api.get_from_ws("GramNegModelTemplateV4","NewKBaseModelTemplates")
+core = NewModelTemplateBuilder.from_dict(get_template('template_core'), None).build()
+types = ["RAST","DRAM","DRAM.RAST"]
+genomeid = sys.argv[1]
+genomews = int(sys.argv[2])
+modelws = int(sys.argv[3])
+#Computing reaction scores
+reaction_genes = {}
+output = anno_api.get_annotation_ontology_events({
+    "input_ref" : str(genomews)+"/"+genomeid,
+})
+events = output["events"]
+for event in events:
+    for gene in event["ontology_terms"]:
+        for term in event["ontology_terms"][gene]:
+            if "modelseed_ids" in term:
+                for rxn in term["modelseed_ids"]:
+                    newrxn = re.sub("^MSRXN:","",rxn)
+                    if newrxn not in reaction_genes:
+                        reaction_genes[newrxn] = {}
+                    if gene not in reaction_genes[newrxn]:
+                        reaction_genes[newrxn][gene] = 0            
+                    reaction_genes[newrxn][gene] += 1
+#Gapfilling RAST, DRAM, and RAST.DRAM models
+for mdltype in types:
+	mdlid = genomeid+"."+mdltype+".mdl"
+	if not exists("GFModels/"+mdltype+"/"+genomeid+".json"):
+		#Getting model, media, templates
+		model = kbase_api.get_from_ws(mdlid,modelws)
+		model.solver = 'optlang-cplex'
+		#Computing ATP to build tests to ensure gapfilling doesn't break ATP
+		atpmethod = MSATPCorrection(model,core,[glc_o2_atp_media])
+		evaluation = atpmethod.evaluate_growth_media()
+		atpmethod.restore_noncore_reactions()
+		threshold = 4*evaluation[glc_o2_atp_media.id]
+		if threshold > 50:
+			threshold = 50
+		tests = [{"media":glc_o2_atp_media,"objective":atpmethod.atp_hydrolysis.id,"is_max_threshold":True,"threshold":threshold}]
+		#Gapfilling
+		msgapfill = MSGapfill(model,[template],[],tests,reaction_genes,[])
+		gfresults = msgapfill.run_gapfilling(gmm,"bio1")
+		if gfresults == None:
+			print("Gapfilling failed on "+mdlid)
+		else:
+			model = msgapfill.integrate_gapfill_solution(gfresults)
+			pkgmgr = MSPackageManager.get_pkg_mgr(model)
+			pkgmgr.getpkg("KBaseMediaPkg").build_package(gmm)
+			solution = model.optimize()
+			print(mdltype+"\t"+genomeid+"\t"+str(solution.objective_value))
+			if solution.objective_value > 0.01:
+				cobra.io.save_json_model(model,"GFModels/"+mdltype+"/"+genomeid+".json")
+			else:
+				print("Model failed validation and will not be saved!")