-
Notifications
You must be signed in to change notification settings - Fork 9
/
gim3e.py
148 lines (122 loc) · 5.38 KB
/
gim3e.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
# -*- coding: utf-8 -*-
# Copyright 2018 Synchon Mandal.
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""Implement the GIM3E method."""
from __future__ import absolute_import
import cobra
import six
from cobra.util import fix_objective_as_constraint
from optlang.symbolics import Zero
from driven.data_sets import ExpressionProfile
def gim3e(model, expression_profile, metabolite_profile, media,
optimum_fraction=0.9, penalty_fraction=1.01, condition=0):
"""
Apply GIM3E [1]_ to a metabolic model.
GIM3E stands for 'Gene Inactivation Moderated by Metabolism, Metabolomics
and Expression'.
Parameters
----------
model: cobra.Model
A constraint-based model to perform GIMME on.
expression_profile: ExpressionProfile
An expression profile to integrate in the model.
metabolite_profile: dict
A dictionary of medium condition as keys and list of
detected metabolites as values.
media: str
The key(media) for the metabolite_profile being used.
optimum_fraction: float, optional (default 0.9)
The fraction of the Required Metabolic Functionalities.
penalty_fraction: float, optional (default 1.01)
The fraction of the penalty bound.
condition: str or int, optional (default 0)
The condition from the expression profile.
If None, the first condition is used.
Returns
-------
context-specfic model: cobra.Model
solution: cobra.Solution
References:
----------
.. [1] Brian J. Schmidt, Ali Ebrahim, Thomas O. Metz, Joshua N. Adkins,
Bernhard Ø. Palsson, Daniel R. Hyduke;
GIM3E: condition-specific models of cellular metabolism developed
from metabolomics and expression data, Bioinformatics, Volume 29,
Issue 22, 15 November 2013, Pages 2900–2908,
https://doi.org/10.1093/bioinformatics/btt493
"""
assert isinstance(model, cobra.Model)
assert isinstance(expression_profile, ExpressionProfile)
assert isinstance(metabolite_profile, dict)
assert isinstance(media, str)
assert isinstance(optimum_fraction, float)
assert isinstance(penalty_fraction, float)
assert isinstance(condition, (str, int))
exp_max, _ = expression_profile.minmax()
tolerance = model.solver.configuration.tolerances.feasibility
objective = model.objective
objective_dir = model.objective_direction
rxn_profile = expression_profile.to_reaction_dict(condition, model)
media_metabolites = metabolite_profile[media]
prob = model.problem
# Fix objective
fix_objective_as_constraint(model, fraction=optimum_fraction)
# Add turnover metabolites
turnover_mets = []
for met in model.metabolites:
turnover_met = cobra.Metabolite("TM_{}".format(met.id))
turnover_mets.append(turnover_met)
model.add_metabolites(turnover_mets)
# Add turnover metabolites to reactions
for rxn in model.reactions:
mets_to_add = {}
for key, value in six.iteritems(rxn.metabolites):
turnover_met = model.metabolites.get_by_id("TM_{}".format(key.id))
mets_to_add[turnover_met] = abs(value)
rxn.add_metabolites(mets_to_add)
# Add sink reactions
for met in model.metabolites.query("^TM_"):
sink_rxn = cobra.Reaction("SK_{}".format(met.id))
# TAKEN FROM THE ORIGINAL IMPLEMENTATION:
# Since both creation and consumption of
# the real metabolite drive the turnover,
# we require 2 units of metabolite per
# 1 unit of flux sink so this matches
# the turnover through the real metabolite.
# METABOLOMICS DATA:
sink_rxn.add_metabolites({met: -2})
if met.id[3:] in media_metabolites:
sink_rxn.lower_bound = penalty_fraction * tolerance
else:
sink_rxn.lower_bound = 0.0
model.add_reaction(sink_rxn)
# TRANSCRIPTOMICS DATA:
# Add penalty coefficients
coefficients = {rxn_id: exp_max - expression for rxn_id, expression in
six.iteritems(rxn_profile)}
# Determine penalty bounds
obj_vars = []
for rxn_id, coefficient in six.iteritems(coefficients):
rxn = model.reactions.get_by_id(rxn_id)
obj_vars.append((rxn.forward_variable, coefficient))
obj_vars.append((rxn.reverse_variable, coefficient))
model.objective = prob.Objective(Zero, sloppy=True, direction="min")
model.objective.set_linear_coefficients({v: c for v, c in obj_vars})
# Penalty bound constraint
penalty_obj = model.slim_optimize()
penalty_bound_const = prob.Constraint(model.objective.expression,
name="penalty_objective_bound",
ub=penalty_fraction * penalty_obj)
model.add_cons_vars(penalty_bound_const)
model.objective = objective
model.objective_direction = objective_dir
sol = model.optimize()
return (model, sol)