import collections
from collections.abc import Iterable
from .core import ComplexPattern, MonomerPattern, Monomer, \
ReactionPattern, ANY, as_complex_pattern, DanglingBondError, \
ReusedBondError, Rule
import networkx as nx
from networkx.algorithms.isomorphism.vf2userfunc import GraphMatcher
from networkx.algorithms.isomorphism import categorical_node_match
import numpy as np
from abc import ABCMeta, abstractmethod
import re
[docs]class FilterPredicate(object):
"""
Base class for building predicates
For use with :func:`pysb.core.ComponentSet.filter`.
"""
__metaclass__ = ABCMeta
@abstractmethod
def __call__(self, component):
pass
def __and__(self, other):
return _AndPredicate(self, other)
def __or__(self, other):
return _OrPredicate(self, other)
def __invert__(self):
return _NotPredicate(self)
class _AndPredicate(FilterPredicate):
"""
Logical conjunction of two predicates.
This class is not generally used directly -- it exists to support operator overloading
in the base FilterPredicate class.
"""
def __init__(self, pred1, pred2):
self.pred1 = pred1
self.pred2 = pred2
def __call__(self, component):
return self.pred1(component) and self.pred2(component)
def __repr__(self):
return '{} & {}'.format(self.pred1.__repr__(), self.pred2.__repr__())
class _OrPredicate(FilterPredicate):
"""
Logical disjunction of two predicates.
This class is not generally used directly -- it exists to support operator overloading
in the base FilterPredicate class.
"""
def __init__(self, pred1, pred2):
self.pred1 = pred1
self.pred2 = pred2
def __call__(self, component):
return self.pred1(component) or self.pred2(component)
def __repr__(self):
return '{} | {}'.format(self.pred1.__repr__(), self.pred2.__repr__())
class _NotPredicate(FilterPredicate):
"""
Logical negation of a predicate.
This class is not generally used directly -- it exists to support operator overloading
in the base FilterPredicate class.
"""
def __init__(self, pred):
self.pred = pred
def __call__(self, component):
return not self.pred(component)
def __repr__(self):
return '~{}'.format(self.pred.__repr__())
[docs]class Name(FilterPredicate):
"""
Predicate to filter a ComponentSet by regular expression name search
Note that this uses re.search which matches anywhere in the component name.
Use ^ to explicitly anchor the match to the beginning.
"""
def __init__(self, regex):
self.regex = regex
def __call__(self, component):
return re.search(self.regex, component.name) is not None
def __repr__(self):
return '{}(\'{}\')'.format(self.__class__.__name__, self.regex)
[docs]class Pattern(FilterPredicate):
"""
Predicate to filter a ComponentSet by matching a ComplexPattern
See :func:`pysb.core.ComponentSet.filter` for examples.
"""
def __init__(self, pattern):
self.pattern = as_complex_pattern(pattern)
def __call__(self, component):
if isinstance(component, (Monomer, MonomerPattern, ComplexPattern)):
return match_complex_pattern(self.pattern,
as_complex_pattern(component))
elif isinstance(component, Rule):
if not isinstance(self.pattern, ComplexPattern):
raise ValueError(
'Cannot currently compare {} to rules'.format(
self.pattern.__class__.__name__))
return any(match_complex_pattern(self.pattern, cp) for cp in
(component.reactant_pattern.complex_patterns +
component.product_pattern.complex_patterns))
else:
raise ValueError('Cannot apply pattern to {}'.format(
component.__class__.__name__))
def __repr__(self):
return '{}({})'.format(self.__class__.__name__, self.pattern)
[docs]class Module(FilterPredicate):
"""
Predicate to filter a ComponentSet by module where components are defined
"""
def __init__(self, regex):
self.regex = regex
def __call__(self, component):
return any(re.search(self.regex, m) for m in component._modules)
[docs]class Function(FilterPredicate):
"""
Predicate to filter a ComponentSet by function where components are defined
"""
def __init__(self, regex):
self.regex = regex
def __call__(self, component):
return re.search(self.regex, component._function)
[docs]def get_half_bonds_in_pattern(pat):
"""
Return the list of integer bond numbers used in a pattern
To return as a set, use :func:`get_bonds_in_pattern`.
Parameters
----------
pat : ComplexPattern, MonomerPattern, or None
A pattern from which bond numberings are extracted
Returns
-------
list
Bond numbers used in the supplied pattern
Examples
--------
>>> A = Monomer('A', ['b1', 'b2'], _export=False)
>>> get_half_bonds_in_pattern(A(b1=None, b2=None))
[]
>>> get_half_bonds_in_pattern(A(b1=1) % A(b2=1))
[1, 1]
"""
bonds_used = list()
def _get_bonds_in_monomer_pattern(mp):
for sc in mp.site_conditions.values():
if isinstance(sc, int):
bonds_used.append(sc)
elif not isinstance(sc, str) and \
isinstance(sc, Iterable):
[bonds_used.append(b) for b in sc if isinstance(b, int)]
if pat is None:
return bonds_used
if isinstance(pat, MonomerPattern):
_get_bonds_in_monomer_pattern(pat)
elif isinstance(pat, ComplexPattern):
for mp in pat.monomer_patterns:
_get_bonds_in_monomer_pattern(mp)
else:
raise ValueError('Unknown pattern type: %s' % type(pat))
return bonds_used
[docs]def get_bonds_in_pattern(pat):
"""
Return the set of integer bond numbers used in a pattern
To return as a list (with duplicates), use
:func:`get_half_bonds_in_pattern`
Parameters
----------
pat : ComplexPattern, MonomerPattern, or None
A pattern from which bond numberings are extracted
Returns
-------
set
Bond numbers used in the supplied pattern
Examples
--------
>>> A = Monomer('A', ['b1', 'b2'], _export=False)
>>> get_bonds_in_pattern(A(b1=None, b2=None)) == set()
True
>>> get_bonds_in_pattern(A(b1=1) % A(b2=1)) == {1}
True
>>> get_bonds_in_pattern(A(b1=1) % A(b1=2, b2=1) % A(b1=2)) == {1, 2}
True
"""
return set(get_half_bonds_in_pattern(pat))
[docs]def check_dangling_bonds(pattern):
"""
Check for dangling bonds in a PySB ComplexPattern/ReactionPattern
Raises a DanglingBondError if a dangling bond is found
"""
if isinstance(pattern, ReactionPattern):
for cp in pattern.complex_patterns:
check_dangling_bonds(cp)
return
bond_counts = collections.Counter(get_half_bonds_in_pattern(pattern))
dangling_bonds = [bond for bond, count in bond_counts.items()
if count == 1]
if dangling_bonds:
raise DanglingBondError('Dangling bond(s) {} in {}'
.format(dangling_bonds, pattern))
reused_bonds = [bond for bond, count in bond_counts.items()
if count > 2]
if reused_bonds:
raise ReusedBondError('Bond(s) {} used more than twice in pattern '
'{}'.format(reused_bonds, pattern))
def _match_graphs(pattern, candidate, exact, count):
""" Compare two pattern graphs for isomorphism """
node_matcher = categorical_node_match('id', default=None)
if exact:
match = nx.is_isomorphic(pattern._as_graph(),
candidate._as_graph(),
node_match=node_matcher)
return 1 if count else match
else:
gm = GraphMatcher(
candidate._as_graph(), pattern._as_graph(),
node_match=node_matcher
)
if count:
if pattern.match_once:
return 1 if gm.subgraph_is_isomorphic() else 0
return sum(1 for _ in gm.subgraph_isomorphisms_iter())
else:
return gm.subgraph_is_isomorphic()
[docs]def match_complex_pattern(pattern, candidate, exact=False, count=False):
"""
Compare two ComplexPatterns against each other
Parameters
----------
pattern: pysb.ComplexPattern
candidate: pysb.Complex.Pattern
exact: bool
Set to True for exact matches (i.e. species equivalence,
or exact graph isomorphism). Set to False to compare as a
pattern (i.e. subgraph isomorphism).
count: bool
Provide match counts for pattern in candidate
Returns
-------
True if pattern matches candidate, False otherwise
"""
if exact:
if not pattern.is_concrete():
raise ValueError('Pattern must be concrete for '
'exact matching: {}'.format(pattern))
if not candidate.is_concrete():
raise ValueError('Candidate must be concrete for '
'exact matching: {}'.format(candidate))
if exact and len(pattern.monomer_patterns) != len(
candidate.monomer_patterns):
return False
# Compare the monomer counts in the patterns so we can fail fast
# without having to compare bonds using graph isomorphism checks, which
# are more computationally expensive
mons_pat = collections.Counter([mp.monomer for mp in \
pattern.monomer_patterns])
mons_cand = collections.Counter([mp.monomer for mp in \
candidate.monomer_patterns])
for mon, mon_count_cand in mons_cand.items():
mon_count_pat = mons_pat.get(mon, 0)
if exact and mon_count_cand != mon_count_pat:
return False
if mon_count_pat > mon_count_cand:
return False
# If we've got this far, we'll need to do a full pattern match
# by searching for a graph isomorphism
return _match_graphs(pattern, candidate, exact=exact, count=count)
[docs]def match_reaction_pattern(pattern, candidate):
"""
Compare two ReactionPatterns against each other
This function tests that every ComplexPattern in pattern has a
matching ComplexPattern in candidate. If there's a one-to-one
mapping of ComplexPattern matches, this is straightforward.
Otherwise, we need to check for a maximum matching - a graph theory
term referring to the maximum number of edges possible in a
bipartite graph (representing ComplexPattern compatibility between
pattern and candidate) without overlapping nodes. If every
ComplexPattern in pattern has a match, then return True, otherwise
return False. This algorithm is polynomial time (although the
ComplexPattern isomorphism comparisons using match_complex_pattern are
not).
Parameters
----------
pattern: pysb.ReactionPattern
candidate: pysb.ReactionPattern
Returns
-------
True if pattern matches candidate, False otherwise.
"""
if len(pattern.complex_patterns) > len(candidate.complex_patterns):
return False
matches = []
for cplx_pat in pattern.complex_patterns:
matches_this = [cand_cplx_pat.matches(
cplx_pat) for cand_cplx_pat in
candidate.complex_patterns]
matches_this = set(np.where(matches_this)[0])
if len(matches_this) == 0:
return False
matches.append(matches_this)
# If a unique 1:1 mapping exists, a match is assured
if len(set.intersection(*matches)) == 0:
return True
# Find the maximum matching in a bipartite graph representing the
# two sets of ComplexPatterns
g = nx.Graph()
g.add_nodes_from(['p%d' % n for n in
range(len(pattern.complex_patterns))], bipartite=0)
g.add_nodes_from(['c%d' % n for n in
range(len(candidate.complex_patterns))], bipartite=1)
for src_pat_id, src_pat_matches in enumerate(matches):
g.add_edges_from([('p%d' % src_pat_id, 'c%d' % cand_pat_id) for
cand_pat_id in src_pat_matches])
return (len(nx.bipartite.maximum_matching(g)) // 2) == len(
pattern.complex_patterns)
[docs]def monomers_from_pattern(pattern):
""" Return the set of monomers used in a pattern """
if isinstance(pattern, ReactionPattern):
return set.union(*[monomers_from_pattern(cp)
for cp in pattern.complex_patterns])
if isinstance(pattern, ComplexPattern):
return set([mp.monomer for mp in pattern.monomer_patterns])
elif isinstance(pattern, MonomerPattern):
return {pattern.monomer}
elif isinstance(pattern, Monomer):
return {pattern}
else:
raise Exception('Unsupported pattern type: %s' % type(pattern))
[docs]class SpeciesPatternMatcher(object):
"""
Match a pattern against a model's species list
Examples
--------
Create a PatternMatcher for the EARM 1.0 model
>>> from pysb.examples.earm_1_0 import model
>>> from pysb.bng import generate_equations
>>> from pysb.pattern import SpeciesPatternMatcher
>>> from pysb import ANY, WILD, Model, Monomer, as_complex_pattern
>>> generate_equations(model)
>>> spm = SpeciesPatternMatcher(model)
Assign two monomers to variables (only needed when importing a model
instead of defining one interactively)
>>> Bax4 = model.monomers['Bax4']
>>> Bcl2 = model.monomers['Bcl2']
Search using a Monomer
>>> spm.match(Bax4)
[Bax4(b=None), Bax4(b=1) % Bcl2(b=1), Bax4(b=1) % Mito(b=1)]
>>> spm.match(Bcl2) # doctest:+NORMALIZE_WHITESPACE
[Bax2(b=1) % Bcl2(b=1),
Bax4(b=1) % Bcl2(b=1),
Bcl2(b=None),
Bcl2(b=1) % MBax(b=1)]
Search using a MonomerPattern (ANY and WILD keywords can be used)
>>> spm.match(Bax4(b=WILD))
[Bax4(b=None), Bax4(b=1) % Bcl2(b=1), Bax4(b=1) % Mito(b=1)]
>>> spm.match(Bcl2(b=ANY))
[Bax2(b=1) % Bcl2(b=1), Bax4(b=1) % Bcl2(b=1), Bcl2(b=1) % MBax(b=1)]
Search using a ComplexPattern
>>> spm.match(Bax4(b=1) % Bcl2(b=1))
[Bax4(b=1) % Bcl2(b=1)]
>>> spm.match(Bax4() % Bcl2())
[Bax4(b=1) % Bcl2(b=1)]
Contrived example to test a site with both a bond and state defined
>>> model = Model(_export=False)
>>> A = Monomer('A', ['a'], {'a': ['u', 'p']}, _export=False)
>>> model.add_component(A)
>>> species = [ \
A(a='u'), \
A(a=1) % A(a=1), \
A(a=('u', 1)) % A(a=('u', 1)), \
A(a=('p', 1)) % A(a=('p', 1)) \
]
>>> model.species = [as_complex_pattern(sp) for sp in species]
>>> spm2 = SpeciesPatternMatcher(model)
>>> spm2.match(A()) # doctest:+NORMALIZE_WHITESPACE
[A(a='u'), A(a=1) % A(a=1), A(a=('u', 1)) % A(a=('u', 1)),
A(a=('p', 1)) % A(a=('p', 1))]
>>> spm2.match(A(a='u'))
[A(a='u')]
>>> spm2.match(A(a=('u', ANY)))
[A(a=('u', 1)) % A(a=('u', 1))]
>>> spm2.match(A(a=('u', WILD)))
[A(a='u'), A(a=('u', 1)) % A(a=('u', 1))]
"""
def __init__(self, model, species=None):
self.model = model
if not species and not model.species:
raise Exception('Model needs species list - run '
'generate_equations() first')
if not species:
species = model.species
self.species = species
self._species_cache = collections.defaultdict(set)
for idx, sp in enumerate(species):
self._add_species(idx, sp)
def _add_species(self, idx, sp):
if sp.compartment:
raise NotImplementedError
for mp in sp.monomer_patterns:
if mp.compartment:
raise NotImplementedError
self._species_cache[mp.monomer].add(idx)
[docs] def add_species(self, species, check_duplicate=True):
"""
Add a species to the search list without adding to the model
Parameters
----------
species : ComplexPattern
A concrete ComplexPattern (molecular species) to add to the
search list
check_duplicate : bool, optional
If True, check the species list to make sure the new species
is not already in the list
"""
if check_duplicate and self.match(species, exact=True):
return
self.species.append(species)
self._add_species(len(self.species) - 1, species)
[docs] def match(self, pattern, index=False, exact=False, counts=False):
"""
Match a pattern against the list of species
Parameters
----------
pattern: pysb.Monomer or pysb.MonomerPattern or pysb.ComplexPattern
index: bool
If True, return species numerical index, rather than species itself
exact: bool
Treat Match as exact equivalence, not a pattern match (i.e. must be
concrete if a MonomerPattern or ComplexPattern)
counts: bool
If True, return match counts for the search pattern within each
species.
Returns
-------
list of pysb.ComplexPattern or list of int
A list of species matching the pattern is returned, unless
index=True, in which case a list of the numerical indices of
matching species is returned instead
Examples
--------
>>> from pysb.examples import earm_1_0
>>> from pysb.bng import generate_equations
>>> model = earm_1_0.model
>>> generate_equations(model)
>>> spm = SpeciesPatternMatcher(model)
>>> L = model.monomers['L']
>>> spm.match(L())
[L(b=None), L(b=1) % pR(b=1)]
"""
if isinstance(pattern, ReactionPattern):
if len(pattern.complex_patterns) == 1:
pattern = pattern.complex_patterns[0]
else:
raise NotImplementedError()
if not isinstance(pattern, (Monomer, MonomerPattern, ComplexPattern)):
raise ValueError('A Monomer, MonomerPattern or ComplexPattern is '
'required to match species')
monomers = monomers_from_pattern(pattern)
if exact:
if isinstance(pattern, (Monomer, MonomerPattern)):
num_mon_pats = 1
else:
num_mon_pats = len(pattern.monomer_patterns)
else:
# Don't check the number of monomer patterns in search
# candidates if we're not doing an exact match of the species
num_mon_pats = None
shortlist, shortlist_indexes = self._species_containing_monomers(
monomers, num_mon_pats)
# If pattern is a Monomer and we don't need match counts, we're done
if isinstance(pattern, Monomer) and not counts:
return shortlist_indexes if index else shortlist
matches = collections.OrderedDict() if counts else []
for idx, sp in enumerate(shortlist):
val = match_complex_pattern(
as_complex_pattern(pattern), sp, exact=exact, count=counts
)
if val:
key = shortlist_indexes[idx] if index else sp
if counts:
matches[key] = val
else:
matches.append(key)
return matches
def _species_containing_monomers(self, monomer_list, num_mon_pats=None):
"""
Identifies species containing a list of monomers
Parameters
----------
monomer_list: list of Monomers
A list of monomers with which to search the model's species
num_mon_pats: int or None
Restrict matches to species with exactly the specified number of
MonomerPatterns
Returns
-------
Model species containing all of the monomers in the list
"""
sp_indexes = set.intersection(*[self._species_cache[mon] for mon in
monomer_list])
if num_mon_pats:
retval = zip(*[(self.species[sp], sp) for sp in sp_indexes
if len(self.species[sp].monomer_patterns)
== num_mon_pats])
return retval if retval else ((), ())
else:
return [self.species[sp] for sp in sp_indexes], list(sp_indexes)
[docs] def rule_firing_species(self, rules_to_consider=None,
include_reverse=True):
"""
Return the species which match the reactants of a set of rules
Parameters
----------
rules_to_consider: list of pysb.Rule or None
A list of rules to use. If None, use all rules in the model.
include_reverse: bool, optional
For reversible rules, include species triggering the rule in
reverse
Returns
-------
collections.OrderedDict
Dictionary of PySB rules whose reactants contain at least one of
the species in the model. Keys are PySB rules, values are a list
of lists. Each outer list corresponding to each
ComplexPattern in the ReactantPattern (or ReactantPattern and
ProductPattern, if rule is reversible). Each inner list contains
the list of species matching the corresponding ComplexPattern.
Examples
--------
>>> from pysb.examples import robertson
>>> from pysb.bng import generate_equations
>>> model = robertson.model
>>> generate_equations(model)
>>> spm = SpeciesPatternMatcher(model)
Get a list of species which fire each rule:
>>> spm.rule_firing_species() \
#doctest: +NORMALIZE_WHITESPACE
OrderedDict([(Rule('A_to_B', A() >> B(), k1), [[A()]]),
(Rule('BB_to_BC', B() + B() >> B() + C(), k2), [[B()], [B()]]),
(Rule('BC_to_AC', B() + C() >> A() + C(), k3), [[B()], [C()]])])
"""
if rules_to_consider is None:
rules_to_consider = self.model.rules
rules_fired = collections.OrderedDict()
for r in rules_to_consider:
rp = r.reactant_pattern
if len(rp.complex_patterns) == 0:
# Synthesis rules are always fired
rules_fired[r] = []
else:
species_fired = self.species_fired_by_reactant_pattern(rp)
if include_reverse and r.is_reversible:
species_fired += self.species_fired_by_reactant_pattern(
r.product_pattern)
if species_fired:
rules_fired[r] = species_fired
return rules_fired
[docs] def species_fired_by_reactant_pattern(self, reaction_pattern):
"""
Get list of species matching a reactant pattern
Parameters
----------
reaction_pattern: pysb.ReactionPattern
Returns
-------
list of lists of pysb.ComplexPattern
List of lists of species matching each ComplexPattern in the
ReactantPattern.
Examples
--------
>>> from pysb.examples import bax_pore
>>> from pysb.bng import generate_equations
>>> model = bax_pore.model
>>> generate_equations(model)
>>> spm = SpeciesPatternMatcher(model)
Get a list of species which fire each rule:
>>> rxn_pat = model.rules['bax_dim'].reactant_pattern
>>> print(rxn_pat)
BAX(t1=None, t2=None) + BAX(t1=None, t2=None)
>>> spm.species_fired_by_reactant_pattern(rxn_pat) \
#doctest: +NORMALIZE_WHITESPACE
[[BAX(t1=None, t2=None, inh=None),
BAX(t1=None, t2=None, inh=1) % MCL1(b=1)],
[BAX(t1=None, t2=None, inh=None),
BAX(t1=None, t2=None, inh=1) % MCL1(b=1)]]
"""
species_fired = []
for i, cp in enumerate(reaction_pattern.complex_patterns):
species_fired_this_cp = self.match(cp)
if not species_fired_this_cp:
return []
else:
species_fired.append(species_fired_this_cp)
return species_fired
[docs]class RulePatternMatcher(object):
"""
Match a pattern against a model's species list
Methods are provided to match against rule reactants, products or both.
Searches can be Monomers, MonomerPatterns, ComplexPatterns or
ReactionPatterns.
Examples
--------
Create a PatternMatcher for the EARM 1.0 model
>>> from pysb.examples.earm_1_0 import model
>>> from pysb.pattern import RulePatternMatcher
>>> rpm = RulePatternMatcher(model)
Assign some monomers to variables (only needed when importing a model
instead of defining one interactively)
>>> AMito, mCytoC, mSmac, cSmac = [model.monomers[m] for m in \
('AMito', 'mCytoC', 'mSmac', 'cSmac')]
Search using a Monomer
>>> rpm.match_reactants(AMito) # doctest:+NORMALIZE_WHITESPACE
[Rule('bind_mCytoC_AMito', AMito(b=None) + mCytoC(b=None) |
AMito(b=1) % mCytoC(b=1), kf20, kr20),
Rule('produce_ACytoC_via_AMito', AMito(b=1) % mCytoC(b=1) >>
AMito(b=None) + ACytoC(b=None), kc20),
Rule('bind_mSmac_AMito', AMito(b=None) + mSmac(b=None) |
AMito(b=1) % mSmac(b=1), kf21, kr21),
Rule('produce_ASmac_via_AMito', AMito(b=1) % mSmac(b=1) >>
AMito(b=None) + ASmac(b=None), kc21)]
>>> rpm.match_products(mSmac) # doctest:+NORMALIZE_WHITESPACE
[Rule('bind_mSmac_AMito', AMito(b=None) + mSmac(b=None) |
AMito(b=1) % mSmac(b=1), kf21, kr21)]
Search using a MonomerPattern
>>> rpm.match_reactants(AMito(b=1)) # doctest:+NORMALIZE_WHITESPACE
[Rule('produce_ACytoC_via_AMito', AMito(b=1) % mCytoC(b=1) >>
AMito(b=None) + ACytoC(b=None), kc20),
Rule('produce_ASmac_via_AMito', AMito(b=1) % mSmac(b=1) >>
AMito(b=None) + ASmac(b=None), kc21)]
>>> rpm.match_rules(cSmac(b=1)) # doctest:+NORMALIZE_WHITESPACE
[Rule('inhibit_cSmac_by_XIAP', cSmac(b=None) + XIAP(b=None) |
cSmac(b=1) % XIAP(b=1), kf28, kr28)]
Search using a ComplexPattern
>>> rpm.match_reactants(AMito() % mSmac()) # doctest:+NORMALIZE_WHITESPACE
[Rule('produce_ASmac_via_AMito', AMito(b=1) % mSmac(b=1) >>
AMito(b=None) + ASmac(b=None), kc21)]
>>> rpm.match_rules(AMito(b=1) % mCytoC(b=1)) \
# doctest:+NORMALIZE_WHITESPACE
[Rule('bind_mCytoC_AMito', AMito(b=None) + mCytoC(b=None) |
AMito(b=1) % mCytoC(b=1), kf20, kr20),
Rule('produce_ACytoC_via_AMito', AMito(b=1) % mCytoC(b=1) >>
AMito(b=None) + ACytoC(b=None), kc20)]
Search using a ReactionPattern
>>> rpm.match_reactants(mCytoC() + mSmac())
[]
>>> rpm.match_reactants(AMito() + mCytoC()) # doctest:+NORMALIZE_WHITESPACE
[Rule('bind_mCytoC_AMito', AMito(b=None) + mCytoC(b=None) |
AMito(b=1) % mCytoC(b=1), kf20, kr20)]
"""
def __init__(self, model):
self.model = model
self._reactant_cache = collections.defaultdict(set)
self._product_cache = collections.defaultdict(set)
for rule in model.rules:
for cache, rp in ((self._reactant_cache, rule.reactant_pattern),
(self._product_cache, rule.product_pattern)):
for cp in rp.complex_patterns:
if cp.compartment:
raise NotImplementedError
for mp in cp.monomer_patterns:
if mp.compartment:
raise NotImplementedError
cache[mp.monomer].add(rule.name)
def match_reactants(self, pattern):
return self._match_reaction_patterns(pattern, 'reactant')
def match_products(self, pattern):
return self._match_reaction_patterns(pattern, 'product')
def match_rules(self, pattern):
return [r for r in self.model.rules if
r in self.match_reactants(pattern) or
r in self.match_products(pattern)]
def _match_reaction_patterns(self, pattern, reaction_side):
if not isinstance(pattern, (Monomer, MonomerPattern, ComplexPattern,
ReactionPattern)):
raise ValueError('A Monomer, MonomerPattern, ComplexPattern or '
'ReactionPattern required to match rules')
monomers = monomers_from_pattern(pattern)
if reaction_side == 'reactant':
cache = self._reactant_cache
def pat_fn(r):
return r.reactant_pattern
elif reaction_side == 'product':
cache = self._product_cache
def pat_fn(r):
return r.product_pattern
else:
raise Exception('reaction_side must be "reactant" or "product"')
shortlist = self._cache_containing_monomers(cache, monomers)
# If pattern is a Monomer, we're done
if isinstance(pattern, Monomer):
return shortlist
if isinstance(pattern, (MonomerPattern, ComplexPattern)):
new_shortlist = []
for rule in shortlist:
reaction_pattern = pat_fn(rule)
if self._match_complex_pattern_to_reaction_pattern(
as_complex_pattern(pattern), reaction_pattern):
new_shortlist.append(rule)
return new_shortlist
else:
return [rule for rule in shortlist if
pat_fn(rule).matches(pattern)]
@classmethod
def _match_complex_pattern_to_reaction_pattern(cls, pattern, test_pattern):
for cp in test_pattern.complex_patterns:
if match_complex_pattern(pattern, cp):
return True
return False
def _cache_containing_monomers(self, cache, monomer_list):
"""
Identifies rules containing a list of monomers
Parameters
----------
monomer_list: list of Monomers
A list of monomers with which to search the model's rules
Returns
-------
Model rules containing all of the monomers in the list
"""
rule_names = set.intersection(*[cache[mon] for mon in
monomer_list])
return [r for r in self.model.rules if r.name in rule_names]
[docs]class ReactionPatternMatcher(object):
"""
Match a pattern against a model's reactions list
Methods are provided to match against reaction reactants, products or
both. Searches can be Monomers, MonomerPatterns, ComplexPatterns or
ReactionPatterns.
Examples
--------
Create a PatternMatcher for the EARM 1.0 model
>>> from pysb.examples.earm_1_0 import model
>>> from pysb.bng import generate_equations
>>> from pysb.pattern import ReactionPatternMatcher
>>> generate_equations(model)
>>> rpm = ReactionPatternMatcher(model)
Assign some monomers to variables (only needed when importing a model
instead of defining one interactively)
>>> AMito, mCytoC, mSmac, cSmac = [model.monomers[m] for m in \
('AMito', 'mCytoC', 'mSmac', 'cSmac')]
Search using a Monomer
>>> rpm.match_products(mSmac) # doctest:+NORMALIZE_WHITESPACE
[Rxn (reversible):
Reactants: {'__s15': mSmac(b=None), '__s45': AMito(b=None)}
Products: {'__s47': AMito(b=1) % mSmac(b=1)}
Rate: kf21*__s15*__s45 - kr21*__s47
Rules: [Rule('bind_mSmac_AMito', AMito(b=None) + mSmac(b=None) |
AMito(b=1) % mSmac(b=1), kf21, kr21)]]
Search using a MonomerPattern
>>> rpm.match_reactants(AMito(b=ANY)) # doctest:+NORMALIZE_WHITESPACE
[Rxn (one-way):
Reactants: {'__s46': AMito(b=1) % mCytoC(b=1)}
Products: {'__s45': AMito(b=None), '__s48': ACytoC(b=None)}
Rate: kc20*__s46
Rules: [Rule('produce_ACytoC_via_AMito', AMito(b=1) % mCytoC(b=1) >>
AMito(b=None) + ACytoC(b=None), kc20)],
Rxn (one-way):
Reactants: {'__s47': AMito(b=1) % mSmac(b=1)}
Products: {'__s45': AMito(b=None), '__s49': ASmac(b=None)}
Rate: kc21*__s47
Rules: [Rule('produce_ASmac_via_AMito', AMito(b=1) % mSmac(b=1) >>
AMito(b=None) + ASmac(b=None), kc21)]]
>>> rpm.match_products(cSmac(b=ANY)) # doctest:+NORMALIZE_WHITESPACE
[Rxn (reversible):
Reactants: {'__s7': XIAP(b=None), '__s51': cSmac(b=None)}
Products: {'__s53': XIAP(b=1) % cSmac(b=1)}
Rate: kf28*__s51*__s7 - kr28*__s53
Rules: [Rule('inhibit_cSmac_by_XIAP', cSmac(b=None) + XIAP(b=None) |
cSmac(b=1) % XIAP(b=1), kf28, kr28)]]
Search using a ComplexPattern
>>> rpm.match_reactants(AMito() % mSmac()) # doctest:+NORMALIZE_WHITESPACE
[Rxn (one-way):
Reactants: {'__s47': AMito(b=1) % mSmac(b=1)}
Products: {'__s45': AMito(b=None), '__s49': ASmac(b=None)}
Rate: kc21*__s47
Rules: [Rule('produce_ASmac_via_AMito', AMito(b=1) % mSmac(b=1) >>
AMito(b=None) + ASmac(b=None), kc21)]]
>>> rpm.match_reactions(AMito(b=3) % mCytoC(b=3)) \
# doctest:+NORMALIZE_WHITESPACE
[Rxn (reversible):
Reactants: {'__s14': mCytoC(b=None), '__s45': AMito(b=None)}
Products: {'__s46': AMito(b=1) % mCytoC(b=1)}
Rate: kf20*__s14*__s45 - kr20*__s46
Rules: [Rule('bind_mCytoC_AMito', AMito(b=None) + mCytoC(b=None) |
AMito(b=1) % mCytoC(b=1), kf20, kr20)],
Rxn (one-way):
Reactants: {'__s46': AMito(b=1) % mCytoC(b=1)}
Products: {'__s45': AMito(b=None), '__s48': ACytoC(b=None)}
Rate: kc20*__s46
Rules: [Rule('produce_ACytoC_via_AMito', AMito(b=1) % mCytoC(b=1) >>
AMito(b=None) + ACytoC(b=None), kc20)]]
"""
def __init__(self, model, species_pattern_matcher=None):
self.model = model
# In this cache, our caches map species to reactions
self._reactant_cache = collections.defaultdict(set)
self._product_cache = collections.defaultdict(set)
if not species_pattern_matcher:
self.spm = SpeciesPatternMatcher(model)
for r_id, rxn in enumerate(model.reactions_bidirectional):
for cache, species_ids in (
(self._reactant_cache, rxn['reactants']),
(self._product_cache, rxn['products'])):
for sp_id in species_ids:
sp = model.species[sp_id]
if sp.compartment:
raise NotImplementedError
cache[sp].add(r_id)
def match_reactants(self, pattern):
return self._match_reactions_against_cache(pattern, 'reactant')
def match_products(self, pattern):
return self._match_reactions_against_cache(pattern, 'product')
def match_reactions(self, pattern):
return self._match_reactions_against_cache(pattern, 'both')
def _match_reactions_against_cache(self, pattern, reaction_side):
species = self.spm.match(pattern)
rxn_ids = set()
if reaction_side in ['reactant', 'both']:
rxn_ids.update(*[self._reactant_cache[sp] for sp in species])
if reaction_side in ['product', 'both']:
rxn_ids.update(*[self._product_cache[sp] for sp in species])
rxn_ids = list(rxn_ids)
rxn_ids.sort()
return [_Reaction(rxn_dict=self.model.reactions_bidirectional[rxn_id],
model=self.model) for rxn_id in rxn_ids]
class _Reaction(object):
__slots__ = ['_rxn_dict', 'reactants', 'model', 'products', 'species']
"""
Store reactions in object form for pretty-printing
"""
def __init__(self, rxn_dict=None, model=None, species=None):
self._rxn_dict = rxn_dict
if model is None:
raise ValueError('Must specify model or species list')
self.model = model
if species:
self.species = species
else:
self.species = model.species
self.reactants = collections.defaultdict(int)
self.products = collections.defaultdict(int)
for r_id in rxn_dict['reactants']:
self.reactants[r_id] += 1
for p_id in rxn_dict['products']:
self.products[p_id] += 1
@property
def reversible(self):
return self._rxn_dict.get('reversible', None)
@property
def reverse(self):
return self._rxn_dict.get('reverse', None)
@property
def rules(self):
return [self.model.rules[r] for r in self._rxn_dict['rule']]
def add_rule(self, rule):
if rule.name not in self._rxn_dict['rule']:
self._rxn_dict['rule'].append(rule.name)
@property
def rate(self):
return self._rxn_dict['rate']
def __repr__(self):
return 'Rxn (%s): \n Reactants: %s\n Products: %s\n ' \
'Rate: %s\n Rules: %s' % (
'reversible' if self.reversible else
('one-way [reverse]' if self.reverse else
'one-way'),
self._repr_species_dict(self.reactants),
self._repr_species_dict(self.products),
self.rate,
self.rules
)
def __cmp__(self, other):
try:
return self._rxn_dict == other._rxn_dict
except AttributeError:
return False
def _repr_species_dict(self, species_dict):
return '{%s}' % ', '.join(["'__s%d': %s" % (k, self.species[k])
for k in sorted(species_dict.keys())])
# doctest fixtures
def setup_module(module):
# Patch OrderedDict repr to match Python 3.11 and earlier.
global _old_od
_old_od = collections.OrderedDict
class OrderedDict(_old_od):
def __repr__(self):
if not self:
return '%s()' % (self.__class__.__name__,)
return '%s(%r)' % (self.__class__.__name__, list(self.items()))
collections.OrderedDict = OrderedDict
setup_module.__test__ = False
def teardown_module(module):
# Restore OrderedDict.
global _old_od
collections.OrderedDict = _old_od
del _old_od
teardown_module.__test__ = False