Source code for pysb.importers.bngl

from pysb.core import MonomerPattern, ComplexPattern, RuleExpression, \
    ReactionPattern, ANY, WILD, MultiState
from pysb.builder import Builder
from pysb.bng import BngFileInterface, parse_bngl_expr
import xml.etree.ElementTree
import re
import sympy
import warnings
import pysb.logging
import collections
import numbers
import os


def _ns(tag_string):
    """
    Shorthand to apply XML namespace
    """
    return tag_string.format('{http://www.sbml.org/sbml/level3}')


class BnglImportError(Exception):
    pass


class BnglBuilder(Builder):
    """
    Assemble a Model from a .bngl file.

    See :py:func:`model_from_bngl` for further details.
    """
    def __init__(self, filename, force=False, cleanup=True):
        super(BnglBuilder, self).__init__()
        filename = os.path.abspath(filename)
        with BngFileInterface(model=None, cleanup=cleanup) as con:
            con.action('readFile', file=filename, skip_actions=1)
            con.action('writeXML', evaluate_expressions=0)
            con.execute()
            self._x = xml.etree.ElementTree.parse('%s.xml' %
                                                  con.base_filename)\
                                           .getroot().find(_ns('{0}model'))

        self.model.name = os.path.splitext(os.path.basename(filename))[0]
        self._force = force
        self._model_env = {}
        self._renamed_states = collections.defaultdict(dict)
        self._renamed_observables = {}
        self._log = pysb.logging.get_logger(__name__)
        self._parse_bng_xml()

    def _warn_or_except(self, msg):
        """
        Raises a warning or Exception, depending of the value of self._force
        """
        if self._force:
            warnings.warn(msg)
        else:
            raise BnglImportError(msg)

    def _eval_in_model_env(self, expression):
        """
        Evaluates an expression string in the model environment
        """
        components = dict((c.name, c) for c in self.model.all_components())
        self._model_env.update(components)

        # Quick security check on the expression
        if re.match(r'^[\w\s()/+\-._*^]*$', expression):
            return parse_bngl_expr(expression, local_dict=self._model_env)
        else:
            self._warn_or_except('Security check on expression "%s" failed' %
                                 expression)
            return None

    def _parse_species(self, species_xml):
        # Species may be empty for synthesis/degradation reaction patterns
        if species_xml is None:
            return []
        # Give the bonds unique numeric IDs
        bond_ids = {}
        next_bond_id = 1
        for bond in species_xml.iterfind(_ns('{0}ListOfBonds/{0}Bond')):
            for bond_attr in bond.attrib.keys():
                if bond_attr.startswith('site'):
                    bond_ids.setdefault(bond.attrib[bond_attr], [])\
                                        .append(next_bond_id)
            next_bond_id += 1

        # Create a list of monomer patterns
        mon_pats = []
        for mon in species_xml.iterfind(_ns('{0}ListOfMolecules/{0}Molecule')):
            mon_name = mon.get('name')
            mon_obj = self.model.monomers[mon_name]
            mon_states = collections.defaultdict(list)
            for comp in mon.iterfind(_ns('{0}ListOfComponents/{0}Component')):
                # BioNetGen component (state) labels are not supported yet
                if 'label' in comp.attrib.keys():
                    self._warn_or_except('BioNetGen component/state labels '
                                         'are not yet supported in PySB')

                state_nm = comp.get('name')
                bonds = comp.get('numberOfBonds')
                if bonds == "0":
                    last_bond = None
                elif bonds == "?":
                    last_bond = WILD
                elif bonds == "+":
                    last_bond = ANY
                else:
                    bond_list = bond_ids[comp.get('id')]
                    assert int(bonds) == len(bond_list)
                    if len(bond_list) == 1:
                        last_bond = bond_list[0]
                    else:
                        last_bond = bond_list
                state = comp.get('state')
                if state and state != '?':
                    if state in ('PLUS', 'MINUS'):
                        self._warn_or_except(
                            'PLUS/MINUS state values are non-standard BNGL '
                            'used to increment or decrement a numeric state '
                            'value. They are not supported in PySB.'
                        )
                    # If we changed the state string, use the updated version
                    state = self._renamed_states.get(mon_name, {}).get(
                        state, state)
                    if last_bond:
                        # Site has a bond and a state
                        mon_states[state_nm].append((state, last_bond))
                    else:
                        # Site only has a state, no bond
                        mon_states[state_nm].append(state)
                else:
                    mon_states[state_nm].append(last_bond)

            mon_states = {k: MultiState(*v) if len(v) > 1 else v[0]
                          for k, v in mon_states.items()}

            mon_cpt = self.model.compartments.get(mon.get('compartment'))
            mon_pats.append(MonomerPattern(mon_obj, mon_states, mon_cpt))
            if 'label' in mon.attrib.keys():
                try:
                    tag = self.model.tags[mon.get('label')]
                except KeyError:
                    tag = self.tag(mon.get('label'))
                mon_pats[-1]._tag = tag
        return mon_pats

    def _parse_monomers(self):
        for m in self._x.iterfind(_ns('{0}ListOfMoleculeTypes/'
                                      '{0}MoleculeType')):
            mon_name = m.get('id')
            sites = []
            states = collections.defaultdict(dict)
            for ctype in m.iterfind(_ns('{0}ListOfComponentTypes/'
                                        '{0}ComponentType')):
                c_name = ctype.get('id')
                sites.append(c_name)
                states_list = ctype.find(_ns('{}ListOfAllowedStates'))
                if states_list is not None:
                    for s in states_list.iterfind(_ns('{}AllowedState')):
                        state = s.get('id')
                        if re.match('[0-9]*$', state):
                            if state not in states[c_name]:
                                new_state = '_' + state
                                while new_state in states[c_name].values():
                                    new_state = '_' + new_state
                                states[c_name][state] = new_state
                                self._renamed_states[mon_name][state] = \
                                    new_state
                        else:
                            states[c_name][state] = state
                    if self._renamed_states[mon_name]:
                        self._log.info('Monomer "{}" states were renamed as '
                                       'follows: {}'.format(
                            mon_name,
                            self._renamed_states[mon_name])
                        )
            try:
                self.monomer(mon_name, sites,
                             {c_name: list(statedict.values())
                              for c_name, statedict in states.items()})
            except Exception as e:
                if str(e).startswith('Duplicate sites specified'):
                    self._warn_or_except('Molecule %s has multiple '
                                         'sites with the same name. '
                                         'This is not supported in PySB.' %
                                         mon_name)
                else:
                    raise BnglImportError(str(e))

    def _parse_parameters(self):
        for p in self._x.iterfind(_ns('{0}ListOfParameters/{0}Parameter')):
            p_name = p.get('id')
            if p.get('type') == 'Constant':
                p_value = p.get('value').replace('10^', '1e')
                try:
                    self.parameter(name=p_name, value=p_value)
                except ValueError:
                    # Despite the "Constant" label, some constant expressions
                    # appear here e.g. ln(2)/120 in BNG's Repressilator model
                    self.expression(name=p_name,
                                    expr=self._eval_in_model_env(
                                        p.get('value')))
            elif p.get('type') == 'ConstantExpression':
                self.expression(name=p_name,
                                expr=self._eval_in_model_env(p.get('value')))
            else:
                self._warn_or_except('Parameter %s has unknown type: %s' %
                                     (p_name, p.get('type')))

    def _parse_observables(self):
        for o in self._x.iterfind(_ns('{0}ListOfObservables/{0}Observable')):
            o_name = o.get('name')
            match_mode = o.get('type').lower()
            # Some BNG observables have same name as a monomer, but in PySB
            # these must be unique
            if o_name in self.model.monomers.keys():
                o_name_old = o_name
                o_name = 'Obs_{}'.format(o_name)
                self._renamed_observables[o_name_old] = o_name
            cplx_pats = []
            for mp in o.iterfind(_ns('{0}ListOfPatterns/{0}Pattern')):
                cpt = self.model.compartments.get(mp.get('compartment'))
                match_once = mp.get('matchOnce')
                match_once = True if match_once == "1" and \
                    match_mode != 'species' else False
                cplx_pats.append(ComplexPattern(self._parse_species(mp),
                                                compartment=cpt,
                                                match_once=match_once))
            self.observable(o_name,
                            ReactionPattern(cplx_pats),
                            match=match_mode)

    def _parse_initials(self):
        for i in self._x.iterfind(_ns('{0}ListOfSpecies/{0}Species')):
            value_param = i.get('concentration')
            try:
                value = float(value_param)
                # Need to create a new parameter for the initial conc. literal
                name = re.sub(r'[^\w]+', '_', i.get('name').replace(
                    ')', '').replace('(', '')) + '_0'
                try:
                    value_param = self.parameter(name, value)
                except ValueError as ve:
                    raise BnglImportError(str(ve))
            except ValueError:
                # Retrieve existing parameter or (constant) expression
                try:
                    value_param = self.model.parameters[i.get('concentration')]
                except KeyError:
                    value_param = self.model.expressions[i.get(
                        'concentration')]
            mon_pats = self._parse_species(i)
            species_cpt = self.model.compartments.get(i.get('compartment'))
            cp = ComplexPattern(mon_pats, species_cpt)
            fixed = i.get('Fixed') == "1"
            self.initial(cp, value_param, fixed)

    def _parse_compartments(self):
        for c in self._x.iterfind(_ns('{0}ListOfCompartments/{0}compartment')):
            cpt_size = None
            if c.get('size'):
                cpt_size = self.parameter('%s_size' % c.get('id'),
                                          c.get('size'))
            cpt_parent = None
            if c.get('outside'):
                cpt_parent = self.model.compartments[c.get('outside')]
            self.compartment(name=c.get('id'),
                             parent=cpt_parent,
                             dimension=int(c.get('spatialDimensions')),
                             size=cpt_size)

    def _parse_rate_law(self, rl):
        if rl.get('type') == 'Ele':
            rate_law_name = rl.find(_ns('{0}ListOfRateConstants/'
                                        '{0}RateConstant')).get('value')
            try:
                return self.model.parameters[rate_law_name]
            except KeyError:
                return self.model.expressions[rate_law_name]
        elif rl.get('type') == 'Function':
            try:
                return self.model.expressions[rl.get('name')]
            except KeyError:
                return self.model.parameters[rl.get('name')]
        else:
            self._warn_or_except('Rate law %s has unknown type %s' %
                                 (rl.get('id'), rl.get('type')))
            return None

    def _parse_rules(self):
        # Store reversible rates for post-processing (we don't know if we'll
        # encounter fwd or rev rule first)
        rev_rates = {}
        for r in self._x.iterfind(_ns('{0}ListOfReactionRules/'
                                      '{0}ReactionRule')):
            r_name = r.get('name')

            r_rate_xml = r.find(_ns('{}RateLaw'))
            if r_rate_xml is None:
                raise BnglImportError('Rate law missing for rule %s' % r_name)
            r_rate = self._parse_rate_law(r_rate_xml)

            # Store reverse rates for post-processing
            if r_name.startswith('_reverse_'):
                r_name = r_name[9:]
                rev_rates[r_name] = r_rate
                continue

            # Compile reactant and product patterns
            reactant_pats = []
            for rp in r.iterfind(_ns('{0}ListOfReactantPatterns/'
                                     '{0}ReactantPattern')):
                cpt = self.model.compartments.get(rp.get('compartment'))
                reactant_pats.append(ComplexPattern(self._parse_species(rp),
                                                    cpt))
                if 'label' in rp.attrib:
                    reactant_pats[-1]._tag = self.model.tags[rp.get('label')]
            product_pats = []
            for pp in r.iterfind(_ns('{0}ListOfProductPatterns/'
                                     '{0}ProductPattern')):
                cpt = self.model.compartments.get(pp.get('compartment'))
                product_pats.append(ComplexPattern(self._parse_species(pp),
                                                   cpt))

                if 'label' in pp.attrib:
                    product_pats[-1]._tag = self.model.tags[pp.get('label')]
            rule_exp = RuleExpression(ReactionPattern(reactant_pats),
                                      ReactionPattern(product_pats),
                                      is_reversible=False)

            # Process any DeleteMolecules declaration
            delete_molecules = False
            for del_operations in r.iterfind(_ns('{0}ListOfOperations/'
                                                 '{0}Delete')):
                if del_operations.get('DeleteMolecules') == "1":
                    delete_molecules = True
                    break

            # Process any MoveConnected declaration
            move_connected = False
            for change_cpt_ops in r.iterfind(_ns('{0}ListOfOperations/'
                                                 '{0}ChangeCompartment')):
                if change_cpt_ops.get('moveConnected') == "1":
                    move_connected = True
                    break

            # Give warning/error if ListOfExcludeReactants or
            # ListOfExcludeProducts is present
            if r.find(_ns('{}ListOfExcludeReactants')) is not None or \
               r.find(_ns('{}ListOfExcludeProducts')) is not None:
                self._warn_or_except('ListOfExcludeReactants and/or '
                                     'ListOfExcludeProducts declarations '
                                     'are deprecated in BNG, and not supported '
                                     'in PySB.')

            # Give warning/error if ListOfIncludeReactants or
            # ListOfIncludeProducts is present
            if r.find(_ns('{}ListOfIncludeReactants')) is not None or \
               r.find(_ns('{}ListOfIncludeProducts')) is not None:
                self._warn_or_except('ListOfIncludeReactants and/or '
                                     'ListOfIncludeProducts declarations '
                                     'are deprecated in BNG, and not supported '
                                     'in PySB.')

            self.rule(r_name, rule_exp, r_rate,
                      delete_molecules=delete_molecules,
                      move_connected=move_connected)

        # Set the reverse rates
        for r_name, rev_rate in rev_rates.items():
            rule = self.model.rules[r_name]
            rule.rule_expression.is_reversible = True
            rule.is_reversible = True
            rule.rate_reverse = rev_rate

    def _parse_expressions(self):
        expr_namespace = (self.model.parameters | self.model.expressions)
        expr_symbols = {e.name: e for e in expr_namespace}

        for e in self._x.iterfind(_ns('{0}ListOfFunctions/{0}Function')):
            for arg in e.iterfind(_ns('{0}ListOfArguments/{0}Argument')):
                tag_name = arg.get('id')
                try:
                    self.model.tags[tag_name]
                except KeyError:
                    tag = self.tag(tag_name)
                    expr_symbols[tag_name] = tag
            expr_name = e.get('id')
            expr_text = e.find(_ns('{0}Expression')).text
            expr_val = 0
            try:
                expr_val = parse_bngl_expr(expr_text, local_dict=expr_symbols)
            except Exception as ex:
                self._warn_or_except('Could not parse expression %s: '
                                     '%s\n\nError: %s' % (expr_name,
                                                          expr_text,
                                                          str(ex)))
            # Replace observables now, so they get expanded by .expand_expr()
            # Doing this as part of expr_symbols breaks local functions!
            observables = {
                sympy.Symbol(o.name): o for o in self.model.observables
            }
            # Add renamed observables
            observables.update({
                sympy.Symbol(obs_old): self.model.observables[obs_new]
                for obs_old, obs_new in self._renamed_observables.items()
            })
            expr_val = expr_val.xreplace(observables)

            if isinstance(expr_val, numbers.Number):
                self.parameter(expr_name, expr_val)
            else:
                self.expression(expr_name, expr_val)

    def _parse_bng_xml(self):
        self._parse_monomers()
        self._parse_parameters()
        self._parse_compartments()
        self._parse_initials()
        self._parse_observables()
        self._parse_expressions()
        self._parse_rules()


[docs]def model_from_bngl(filename, force=False, cleanup=True): """ Convert a BioNetGen .bngl model file into a PySB Model. Notes ----- The following features are not supported in PySB and will cause an error if present in a .bngl file: * Fixed species (with a ``$`` prefix, like ``$Null``) * BNG excluded or included reaction patterns (deprecated in BNG) * BNG local functions * Molecules with identically named sites, such as ``M(l,l)`` * BNG's custom rate law functions, such as ``MM`` and ``Sat`` (deprecated in BNG) Parameters ---------- filename : string A BioNetGen .bngl file force : bool, optional The default, False, will raise an Exception if there are any errors importing the model to PySB, e.g. due to unsupported features. Setting to True will attempt to ignore any import errors, which may lead to a model that only poorly represents the original. Use at own risk! cleanup : bool Delete temporary directory on completion if True. Set to False for debugging purposes. """ bb = BnglBuilder(filename, force=force, cleanup=cleanup) return bb.model