Skip to content
This repository was archived by the owner on Feb 9, 2024. It is now read-only.

Commit 549aaa3

Browse files
cshenrycshenry
committed
Checking in the multiomics notebook
1 parent df7c9f1 commit 549aaa3

File tree

1 file changed

+177
-28
lines changed

1 file changed

+177
-28
lines changed

‎PLP/MultiomicsDataIntegration.ipynb

+177-28
Original file line numberDiff line numberDiff line change
@@ -29,61 +29,210 @@
2929
" sys.path.append(path)\n",
3030
"import cobrakbase\n",
3131
"from IPython.core.display import HTML\n",
32-
"from modelseedpy import MSPackageManager"
32+
"from modelseedpy import MSMedia, MSPackageManager, MSExpression"
3333
]
3434
},
3535
{
3636
"cell_type": "code",
37-
"execution_count": 3,
37+
"execution_count": 2,
3838
"metadata": {},
3939
"outputs": [
4040
{
4141
"name": "stdout",
4242
"output_type": "stream",
4343
"text": [
44-
"4HTHRK\n",
45-
"FOUND!\n"
46-
]
47-
},
48-
{
49-
"ename": "KeyError",
50-
"evalue": "'bgoalc'",
51-
"output_type": "error",
52-
"traceback": [
53-
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
54-
"\u001b[0;31mKeyError\u001b[0m Traceback (most recent call last)",
55-
"\u001b[0;32m<ipython-input-3-97902581d7e4>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 13\u001b[0m pkgmgr.getpkg(\"ChangeOptPkg\").build_package({\n\u001b[1;32m 14\u001b[0m \u001b[0;34m\"4HTHRK\"\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m{\u001b[0m\u001b[0;34m\"flux\"\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0msolution\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfluxes\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"4HTHRK\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;36m1.2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\"direction\"\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\"high\"\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\"objcoef\"\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 15\u001b[0;31m \u001b[0;34m\"THRS\"\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m{\u001b[0m\u001b[0;34m\"flux\"\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0msolution\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfluxes\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"THRS\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;36m1.2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\"direction\"\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\"high\"\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\"objcoef\"\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 16\u001b[0m })\n\u001b[1;32m 17\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"multiomics.lp\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'w'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mout\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
56-
"\u001b[0;32m~/code/ModelSEEDpy/modelseedpy/fbapkg/changeoptpkg.py\u001b[0m in \u001b[0;36mbuild_package\u001b[0;34m(self, target_values, build_objective)\u001b[0m\n\u001b[1;32m 26\u001b[0m \u001b[0mvar\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuild_variable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"bgoal\"\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mrxn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 27\u001b[0m \u001b[0mobj_coef\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mvar\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtarget_values\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mrxnid\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"objcoef\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 28\u001b[0;31m \u001b[0mconst\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuild_constraint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"bgoalc\"\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mrxn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 29\u001b[0m \u001b[0;31m#Now setting the goal objective if build_objective is true\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 30\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mbuild_objective\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
57-
"\u001b[0;32m~/code/ModelSEEDpy/modelseedpy/fbapkg/changeoptpkg.py\u001b[0m in \u001b[0;36mbuild_constraint\u001b[0;34m(self, type, object)\u001b[0m\n\u001b[1;32m 56\u001b[0m \u001b[0mcoef\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mfluxvar\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 57\u001b[0m \u001b[0mlb\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0mgoal\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"flux\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 58\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mBaseFBAPkg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuild_constraint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mtype\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mlb\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcoef\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mobject\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 59\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
58-
"\u001b[0;32m~/code/ModelSEEDpy/modelseedpy/fbapkg/basefbapkg.py\u001b[0m in \u001b[0;36mbuild_constraint\u001b[0;34m(self, type, lower_bound, upper_bound, coef, object)\u001b[0m\n\u001b[1;32m 77\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mbuild_constraint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mtype\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mlower_bound\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mupper_bound\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcoef\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mobject\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 78\u001b[0m \u001b[0mname\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 79\u001b[0;31m \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconstraint_types\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mtype\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m\"none\"\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 80\u001b[0m \u001b[0mcount\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconstraints\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mtype\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 81\u001b[0m \u001b[0mname\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mstr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcount\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
59-
"\u001b[0;31mKeyError\u001b[0m: 'bgoalc'"
44+
"4.447533484208366\n",
45+
"0.3768000000000051\n"
6046
]
6147
}
6248
],
6349
"source": [
64-
"kbase_api = cobrakbase.KBaseAPI()\n",
6550
"model= cobra.io.read_sbml_model(\"iML1515.xml\")\n",
6651
"model.solver = 'optlang-cplex'\n",
67-
"gmm = kbase_api.get_from_ws(\"Carbon-D-Glucose\",\"KBaseMedia\")\n",
52+
"media_dict = {\n",
53+
" \"glc__D\":22.5,\n",
54+
" \"nh4\":100,\n",
55+
" \"k\":100,\n",
56+
" \"pydxn\":0.00002,\n",
57+
" 'so4':100,\n",
58+
" 'pi':100,\n",
59+
" 'mn2': 100,\n",
60+
" 'fe2': 100,\n",
61+
" 'na1': 100,\n",
62+
" 'mg2': 100,\n",
63+
" 'o2': 20,\n",
64+
" 'h2o': 100,\n",
65+
" 'h': 100,\n",
66+
" 'fe3': 100.0,\n",
67+
" 'zn2': 100.0,\n",
68+
" 'ca2': 100.0,\n",
69+
" 'ni2': 100.0,\n",
70+
" 'cu2': 100.0,\n",
71+
" 'sel': 100.0,\n",
72+
" 'cobalt2': 100.0,\n",
73+
" 'mobd': 100.0,\n",
74+
" 'cl': 100.0,\n",
75+
" 'cbl1':100,\n",
76+
" 'tungs': 100.0,\n",
77+
" 'slnt': 100.0\n",
78+
"}\n",
79+
"media = MSMedia.from_dict(media_dict)\n",
6880
"pkgmgr = MSPackageManager.get_pkg_mgr(model)\n",
69-
"pkgmgr.getpkg(\"KBaseMediaPkg\").build_package(gmm)\n",
81+
"pkgmgr.getpkg(\"KBaseMediaPkg\").build_package(media)\n",
82+
"#Plug in growth rate and minimize nutrient uptake\n",
83+
"model.reactions.BIOMASS_Ec_iML1515_WT_75p37M.lower_bound = 0.3768\n",
84+
"model.reactions.BIOMASS_Ec_iML1515_WT_75p37M.upper_bound = 0.3768\n",
85+
"glucose_uptake_objective = model.problem.Objective(\n",
86+
" 1 * model.reactions.EX_glc__D_e.flux_expression,\n",
87+
" direction='max')\n",
88+
"model.objective = glucose_uptake_objective\n",
89+
"solution = model.optimize()\n",
90+
"#Constraining media to limit glucose to hit measured growth rate\n",
91+
"media_dict[\"glc__D\"] = -1*solution.fluxes[\"EX_glc__D_e\"]\n",
92+
"print(media_dict[\"glc__D\"])\n",
93+
"media = MSMedia.from_dict(media_dict)\n",
94+
"pkgmgr.getpkg(\"KBaseMediaPkg\").build_package(media)\n",
95+
"#Confirming that we now get the correct growth rate\n",
96+
"model.reactions.BIOMASS_Ec_iML1515_WT_75p37M.lower_bound = 0\n",
97+
"model.reactions.BIOMASS_Ec_iML1515_WT_75p37M.upper_bound = 100\n",
7098
"biomass_objective = model.problem.Objective(\n",
7199
" 1 * model.reactions.BIOMASS_Ec_iML1515_WT_75p37M.flux_expression,\n",
72100
" direction='max')\n",
73101
"model.objective = biomass_objective\n",
74102
"solution = cobra.flux_analysis.pfba(model)\n",
103+
"print(solution.fluxes[\"BIOMASS_Ec_iML1515_WT_75p37M\"])"
104+
]
105+
},
106+
{
107+
"cell_type": "code",
108+
"execution_count": 3,
109+
"metadata": {},
110+
"outputs": [
111+
{
112+
"name": "stdout",
113+
"output_type": "stream",
114+
"text": [
115+
"binary_variable_count: 7\n",
116+
"4HTHRK\n",
117+
"FOUND!\n",
118+
"THRS\n",
119+
"FOUND!\n"
120+
]
121+
},
122+
{
123+
"data": {
124+
"text/html": [
125+
"<strong><em>Optimal</em> solution with objective value 0.000</strong><br><div>\n",
126+
"<style scoped>\n",
127+
" .dataframe tbody tr th:only-of-type {\n",
128+
" vertical-align: middle;\n",
129+
" }\n",
130+
"\n",
131+
" .dataframe tbody tr th {\n",
132+
" vertical-align: top;\n",
133+
" }\n",
134+
"\n",
135+
" .dataframe thead th {\n",
136+
" text-align: right;\n",
137+
" }\n",
138+
"</style>\n",
139+
"<table border=\"1\" class=\"dataframe\">\n",
140+
" <thead>\n",
141+
" <tr style=\"text-align: right;\">\n",
142+
" <th></th>\n",
143+
" <th>fluxes</th>\n",
144+
" <th>reduced_costs</th>\n",
145+
" </tr>\n",
146+
" </thead>\n",
147+
" <tbody>\n",
148+
" <tr>\n",
149+
" <th>CYTDK2</th>\n",
150+
" <td>0.000000</td>\n",
151+
" <td>NaN</td>\n",
152+
" </tr>\n",
153+
" <tr>\n",
154+
" <th>XPPT</th>\n",
155+
" <td>0.000000</td>\n",
156+
" <td>NaN</td>\n",
157+
" </tr>\n",
158+
" <tr>\n",
159+
" <th>HXPRT</th>\n",
160+
" <td>0.000000</td>\n",
161+
" <td>NaN</td>\n",
162+
" </tr>\n",
163+
" <tr>\n",
164+
" <th>NDPK5</th>\n",
165+
" <td>-0.009347</td>\n",
166+
" <td>NaN</td>\n",
167+
" </tr>\n",
168+
" <tr>\n",
169+
" <th>SHK3Dr</th>\n",
170+
" <td>0.140141</td>\n",
171+
" <td>NaN</td>\n",
172+
" </tr>\n",
173+
" <tr>\n",
174+
" <th>...</th>\n",
175+
" <td>...</td>\n",
176+
" <td>...</td>\n",
177+
" </tr>\n",
178+
" <tr>\n",
179+
" <th>MPTS</th>\n",
180+
" <td>0.000005</td>\n",
181+
" <td>NaN</td>\n",
182+
" </tr>\n",
183+
" <tr>\n",
184+
" <th>MOCOS</th>\n",
185+
" <td>0.000005</td>\n",
186+
" <td>NaN</td>\n",
187+
" </tr>\n",
188+
" <tr>\n",
189+
" <th>BMOGDS2</th>\n",
190+
" <td>0.000000</td>\n",
191+
" <td>NaN</td>\n",
192+
" </tr>\n",
193+
" <tr>\n",
194+
" <th>FESD2s</th>\n",
195+
" <td>0.000000</td>\n",
196+
" <td>NaN</td>\n",
197+
" </tr>\n",
198+
" <tr>\n",
199+
" <th>OCTNLL</th>\n",
200+
" <td>0.000001</td>\n",
201+
" <td>NaN</td>\n",
202+
" </tr>\n",
203+
" </tbody>\n",
204+
"</table>\n",
205+
"<p>2712 rows × 2 columns</p>\n",
206+
"</div>"
207+
],
208+
"text/plain": [
209+
"<Solution 0.000 at 0x7fe1b903a990>"
210+
]
211+
},
212+
"execution_count": 3,
213+
"metadata": {},
214+
"output_type": "execute_result"
215+
}
216+
],
217+
"source": [
75218
"#Add the gene expression profile to the project notebook folder (Fatime)\n",
76-
"expression_data = MSExpression.from_gene_feature_file(\"DifferentialExpression.tsv\",None,False)\n",
219+
"expression_data = MSExpression.from_gene_feature_file(\"differentialexpression_gene.tsv\",None,False)\n",
77220
"#Covert differential gene expression profile into differential reaction expression profile (Chris)\n",
78-
"#Print this reaction profile and spot check values against Jill's manual numbers (Fatima,Jill)\n",
79-
"#Add table of metabolomic differential expression to project notebook (Fatima)\n",
221+
"reaction_expression = expression_data.build_reaction_expression(model,0)\n",
80222
"#Need to load the metabolite differential profile data in terms of iML1515 IDs\n",
81-
"#metabolite\\tfoldchange\n",
82-
"#Translate data into goals (Chris)\n",
223+
"metabolite_data = MSExpression.from_gene_feature_file(\"Metabolomics.tsv\",None,False)\n",
224+
"#Translate data into goals\n",
225+
"goalflux = {}\n",
226+
"reactions = [\"4HTHRK\",\"THRS\"]\n",
227+
"for rxn in reactions:\n",
228+
" if solution.fluxes[rxn] == 0:\n",
229+
" goalflux[rxn] = 1\n",
230+
" else:\n",
231+
" goalflux[rxn] = 1.5*solution.fluxes[rxn]\n",
83232
"pkgmgr.getpkg(\"BilevelPkg\").build_package(None,7)\n",
84233
"pkgmgr.getpkg(\"ChangeOptPkg\").build_package({\n",
85-
" \"4HTHRK\":{\"flux\":solution.fluxes[\"4HTHRK\"]*1.2,\"direction\":\"high\",\"objcoef\":1},\n",
86-
" \"THRS\":{\"flux\":solution.fluxes[\"THRS\"]*1.2,\"direction\":\"high\",\"objcoef\":1}\n",
234+
" \"4HTHRK\":{\"flux\":goalflux[\"4HTHRK\"]*1.2,\"direction\":\"high\",\"objcoef\":1},\n",
235+
" \"THRS\":{\"flux\":goalflux[\"THRS\"]*1.2,\"direction\":\"high\",\"objcoef\":1}\n",
87236
"})\n",
88237
"with open(\"multiomics.lp\", 'w') as out:\n",
89238
" out.write(str(model.solver))\n",

0 commit comments

Comments
 (0)