9079427458

Committed 14 May 2024 12:34PM UTC coverage: 95.894% (-0.4%) from 96.296%

Build # 9079427458

Build Type

Pull #41

travis-ci

Committed by

web-flow

Commit Message

Merge 2f1b53645 into a1684afb9

Pull Request Pull Request #41: Aromatics, round 2

Run Details

232 of 248 branches covered (93.55%)

Branch coverage included in aggregate %.

43 of 43 new or added lines in 2 files covered. (100.0%)

2 existing lines in 1 file now uncovered.

422 of 434 relevant lines covered (97.24%)

3.88 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

94.81

/pysmiles/smiles_helper.py

# -*- coding: utf-8 -*-
# Copyright 2018 Peter C Kroon

# 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.

"""
Contains helper functions for parsing and writing SMILES strings, as well as
some convenience functions for adding hydrogens, and detecting aromaticity.
"""

import logging
import re
import operator

import networkx as nx

LOGGER = logging.getLogger(__name__)

ISOTOPE_PATTERN = r'(?P<isotope>[\d]+)?'
ELEMENT_PATTERN = r'(?P<element>b|c|n|o|s|p|as|se|\*|[A-Z][a-z]{0,2})'
STEREO_PATTERN = r'(?P<stereo>@|@@|@TH[1-2]|@AL[1-2]|@SP[1-3]|@OH[\d]{1,2}|'\
                  r'@TB[\d]{1,2})?'
HCOUNT_PATTERN = r'(?P<hcount>H[\d]?)?'
CHARGE_PATTERN = r'(?P<charge>(-|\+)(\++|-+|[\d]{1,2})?)?'
CLASS_PATTERN = r'(?::(?P<class>[\d]+))?'
ATOM_PATTERN = re.compile(r'^\[' + ISOTOPE_PATTERN + ELEMENT_PATTERN +
                          STEREO_PATTERN + HCOUNT_PATTERN + CHARGE_PATTERN +
                          CLASS_PATTERN + r'\]$')

VALENCES = {"B": (3,), "C": (4,), "N": (3, 5), "O": (2,), "P": (3, 5),
            "S": (2, 4, 6), "F": (1,), "Cl": (1,), "Br": (1,), "I": (1,)}

AROMATIC_ATOMS = "B C N O P S Se As *".split()


def parse_atom(atom):
    """
    Parses a SMILES atom token, and returns a dict with the information.

    Note
    ----
    Can not deal with stereochemical information yet. This gets discarded.

    Parameters
    ----------
    atom : str
        The atom string to interpret. Looks something like one of the
        following: "C", "c", "[13CH3-1:2]"

    Returns
    -------
    dict
        A dictionary containing at least 'element', 'aromatic', and 'charge'. If
        present, will also contain 'hcount', 'isotope', and 'class'.
    """
    defaults = {'charge': 0, 'hcount': 0, 'aromatic': False}
    if not atom.startswith('[') and not atom.endswith(']'):
        if atom != '*':
            # Don't specify hcount to signal we don't actually know anything
            # about it
            return {'element': atom.capitalize(), 'charge': 0,
                    'aromatic': atom.islower()}
        else:
            return defaults.copy()
    match = ATOM_PATTERN.match(atom)
    if match is None:
        raise ValueError('The atom {} is malformatted'.format(atom))
    out = defaults.copy()
    out.update({k: v for k, v in match.groupdict().items() if v is not None})

    if out.get('element', 'X').islower():
        out['aromatic'] = True

    parse_helpers = {
        'isotope': int,
        'element': str.capitalize,
        'stereo': lambda x: x,
        'hcount': parse_hcount,
        'charge': parse_charge,
        'class': int,
        'aromatic': lambda x: x,
    }

    for attr, val_str in out.items():
        out[attr] = parse_helpers[attr](val_str)

    if out['element'] == '*':
        del out['element']

    if out.get('element') == 'H' and out.get('hcount', 0):
        raise ValueError("A hydrogen atom can't have hydrogens")

    if 'stereo' in out:
        LOGGER.warning('Atom "%s" contains stereochemical information that will be discarded.', atom)

    return out


def format_atom(molecule, node_key, default_element='*'):
    """
    Formats a node following SMILES conventions. Uses the attributes `element`,
    `charge`, `hcount`, `stereo`, `isotope` and `class`.

    Parameters
    ----------
    molecule : nx.Graph
        The molecule containing the atom.
    node_key : hashable
        The node key of the atom in `molecule`.
    default_element : str
        The element to use if the attribute is not present in the node.

    Returns
    -------
    str
        The atom as SMILES string.
    """
    node = molecule.nodes[node_key]
    name = node.get('element', default_element)
    charge = node.get('charge', 0)
    hcount = node.get('hcount', 0)
    stereo = node.get('stereo', None)
    isotope = node.get('isotope', '')
    class_ = node.get('class', '')
    aromatic = node.get('aromatic', False)
    default_h = has_default_h_count(molecule, node_key)

    if stereo is not None:
        raise NotImplementedError

    if aromatic:
        name = name.lower()

    if (stereo is None and isotope == '' and charge == 0 and default_h and class_ == '' and
            (name.lower() in 'b c n o p s *'.split() or name in 'F Cl Br I'.split())):
        return name

    if hcount:
        hcountstr = 'H'
        if hcount > 1:
            hcountstr += str(hcount)
    else:
        hcountstr = ''

    if charge > 0:
        chargestr = '+'
        if charge > 1:
            chargestr += str(charge)
    elif charge < 0:
        chargestr = '-'
        if charge < -1:
            chargestr += str(-charge)
    else:
        chargestr = ''

    if class_ != '':
        class_ = ':{}'.format(class_)
    fmt = '[{isotope}{name}{stereo}{hcount}{charge}{class_}]'
    return fmt.format(isotope=isotope, name=name, stereo='', hcount=hcountstr,
                      charge=chargestr, class_=class_)


def parse_hcount(hcount_str):
    """
    Parses a SMILES hydrogen count specifications.

    Parameters
    ----------
    hcount_str : str
        The hydrogen count specification to parse.

    Returns
    -------
    int
        The number of hydrogens specified.
    """
    if not hcount_str:
        return 0
    if hcount_str == 'H':
        return 1
    return int(hcount_str[1:])


def parse_charge(charge_str):
    """
    Parses a SMILES charge specification.

    Parameters
    ----------
    charge_str : str
        The charge specification to parse.

    Returns
    -------
    int
        The charge.
    """
    if not charge_str:
        return 0
    signs = {'-': -1, '+': 1}
    sign = signs[charge_str[0]]
    if len(charge_str) > 1 and charge_str[1].isdigit():
        charge = sign * int(charge_str[1:])
    else:
        charge = sign * charge_str.count(charge_str[0])
    return charge


def add_explicit_hydrogens(mol):
    """
    Adds explicit hydrogen nodes to `mol`, the amount is determined by the node
    attribute 'hcount'. Will remove the 'hcount' attribute.

    Parameters
    ----------
    mol : nx.Graph
        The molecule to which explicit hydrogens should be added. Is modified
        in-place.

    Returns
    -------
    None
        `mol` is modified in-place.
    """
    h_atom = parse_atom('[H]')
    if 'hcount' in h_atom:
        del h_atom['hcount']
    for n_idx in list(mol.nodes):
        hcount = mol.nodes[n_idx].get('hcount', 0)
        idxs = range(max(mol) + 1, max(mol) + hcount + 1)
        # Get the defaults from parse_atom.
        mol.add_nodes_from(idxs, **h_atom.copy())
        mol.add_edges_from([(n_idx, jdx) for jdx in idxs], order=1)
        if 'hcount' in mol.nodes[n_idx]:
            del mol.nodes[n_idx]['hcount']


def remove_explicit_hydrogens(mol):
    """
    Removes all explicit, simple hydrogens from `mol`. Simple means it is
    identical to the SMILES string "[H]", and has exactly one bond. Increments
    'hcount' where appropriate.

    Parameters
    ----------
    mol : nx.Graph
        The molecule whose explicit hydrogens should be removed. Is modified
        in-place.

    Returns
    -------
    None
        `mol` is modified in-place.
    """
    to_remove = set()
#    defaults = parse_atom('[H]')
    for n_idx in mol.nodes:
        node = mol.nodes[n_idx]
        neighbors = list(mol[n_idx])
        # TODO: get these defaults from parsing [H]. But do something smart
        #       with the hcount attribute.
        if (node.get('charge', 0) == 0 and node.get('element', '') == 'H' and
                'isotope' not in node and node.get('class', 0) == 0 and
                len(neighbors) == 1):
            neighbor = neighbors[0]
            if (mol.nodes[neighbor].get('element', '') == 'H' or
                    mol.edges[n_idx, neighbor].get('order', 1) != 1):
                # The molecule is H2, or the bond order is not 1.
                continue
            to_remove.add(n_idx)
            mol.nodes[neighbor]['hcount'] = mol.nodes[neighbor].get('hcount', 0) + 1
    mol.remove_nodes_from(to_remove)
    for n_idx in mol.nodes:
        if 'hcount' not in mol.nodes[n_idx]:
            mol.nodes[n_idx]['hcount'] = 0


def fill_valence(mol, respect_hcount=True, respect_bond_order=True,
                 max_bond_order=3):
    """
    Sets the attribute 'hcount' on all nodes in `mol` that don't have it yet.
    The value to which it is set is based on the node's 'element', and the
    number of bonds it has. Default valences are as specified by the global
    variable VALENCES.

    Parameters
    ----------
    mol : nx.Graph
        The molecule whose nodes should get a 'hcount'. Is modified in-place.
    respect_hcount : bool
        If True, don't change the hcount on nodes that already have it set.
    respect_bond_order : bool
        If False, first try to fill the valence by increasing bond orders, and
        add hydrogens after.
    max_bond_order : number
        Only meaningful if respect_bond_order is False. This is the highest
        bond order that will be set.

    Returns
    -------
    None
        `mol` is modified in-place.
    """
    if not respect_bond_order:
        increment_bond_orders(mol, max_bond_order=max_bond_order)
    for n_idx in mol:
        node = mol.nodes[n_idx]
        if 'hcount' in node and respect_hcount:
            continue
        missing = max(bonds_missing(mol, n_idx), 0)
        node['hcount'] = node.get('hcount', 0) + missing


def bonds_missing(mol, node_idx, use_order=True):
    """
    Returns how much the specified node is under valence. If use_order is
    False, treat all bonds as if they are order 1.

    Parameters
    ----------
    mol : nx.Graph
        The molecule.
    node_idx : hashable
        The node to look at. Should be in mol.
    use_order : bool
        If False, treat all bonds as single.

    Returns
    -------
    int
        The number of missing bonds.
    """
    bonds = _bonds(mol, node_idx, use_order)
    bonds += mol.nodes[node_idx].get('hcount', 0)
    valence = _valence(mol, node_idx, bonds)
    return int(valence - bonds)


def _valence(mol, node_idx, minimum=0):
    """
    Returns the valence of the specified node. Since some elements can have
    multiple valences, give the smallest one that is more than `minimum`.

    Parameters
    ----------
    mol : nx.Graph
        The molecule.
    node_idx : hashable
        The node to look at. Should be in mol.
    minimum : int
        The minimum value of valence.

    Returns
    -------
    int
        The smallest valence of node more than `minimum`.
    """
    element = mol.nodes[node_idx].get('element', '').capitalize()
    if element not in VALENCES:
        return 0
    val = VALENCES.get(element)
    try:
        val = min(filter(lambda a: a >= minimum, val))
    except ValueError:  # More bonds than possible
        val = max(val)
    return val


def _bonds(mol, node_idx, use_order=True):
    """
    Returns how many explicit bonds the specified node has. If use_order is
    False, treat all bonds as if they are order 1.

    Parameters
    ----------
    mol : nx.Graph
        The molecule.
    node_idx : hashable
        The node to look at. Should be in mol.
    use_order : bool
        If False, treat all bonds as single.

    Returns
    -------
    int
        The number of bonds.
    """
    if use_order:
        bond_orders = map(operator.itemgetter(2),
                          mol.edges(nbunch=node_idx, data='order', default=1))
        bonds = sum(bond_orders)
    else:
        bonds = len(mol[node_idx])
    return bonds


def has_default_h_count(mol, node_idx, use_order=True):
    """
    Returns whether the hydrogen count for this atom is non-standard.

    Parameters
    ----------
    mol : nx.Graph
        The molecule.
    node_idx : hashable
        The node to look at. Should be in mol.
    use_order : bool
        If False, treat all bonds as single.

    Returns
    -------
    bool
    """
    bonds = _bonds(mol, node_idx, use_order)
    valence = _valence(mol, node_idx, bonds)
    hcount = mol.nodes[node_idx].get('hcount', 0)
    return valence - bonds == hcount


def _hydrogen_neighbours(mol, n_idx):
    neighbours = mol[n_idx]
    h_neighbours = 0
    for n_jdx in neighbours:
        if (mol.nodes[n_jdx].get('element', '*') == 'H' and
                mol.edges[n_idx, n_jdx].get('order', 1) == 1):
            h_neighbours += 1
    return h_neighbours

def _prune_nodes(nodes, mol):
    new_nodes = []
    for node in nodes:
        # all wild card nodes are eligible
        if mol.nodes[node].get('element', '*') == '*':
            new_nodes.append(node)
            continue
        missing = bonds_missing(mol, node, use_order=True) + mol.nodes[node].get('charge', 0)
        if missing > 0:
            new_nodes.append(node)
    return mol.subgraph(new_nodes)


def mark_aromatic_atoms(mol, strict=True):
    """
    Properly kekulizes molecules and sets the aromatic attribute.

    Parameters
    ----------
    mol : nx.Graph
        The molecule.
    strict: bool
        Whether to raise a SyntaxError if the aromatic region of the molecule
        cannot be kekulized.

    Returns
    -------
    None
        `mol` is modified in-place.

    Raises
    ------
    SyntaxError
        If ``strict`` is True and the aromatic region of the molecule cannot
        be kekulized.
    """
    # prune all nodes from molecule that are eligible and have
    # full valency
    arom_atoms = [node for node, aromatic in mol.nodes(data='aromatic') if aromatic]
    ds_graph = _prune_nodes(arom_atoms, mol)

    # set the aromatic attribute to False for all nodes
    # as a precaution
    nx.set_node_attributes(mol, False, 'aromatic')

    sub_ds_graph = mol.subgraph(ds_graph)
    max_match = nx.max_weight_matching(sub_ds_graph)
    # we check if a maximum matching exists and
    # if it is perfect. if it is not perfect,
    # this graph originates from a completely invalid
    # smiles and we raise an error
    if strict and not nx.is_perfect_matching(sub_ds_graph, max_match):
        msg = "Your molecule is invalid and cannot be kekulized."
        raise SyntaxError(msg)

    # From here we want to do 2 things: 1) dekekulize aromatic rings; and 2)
    # kekulize "aromatic" fragments that are not in a cycle. We only consider
    # nodes to be aromatic if they can participate in DIME, which means they
    # must be in a cycle.
    # 1) Dekekulize aromatic rings
    # Add all double bonds to the perfect matching we have, and then for every
    # simple cycle check if the resulting matching is still perfect. You need
    # to take simple cycles (rather than the cycle basis), for cases like the
    # following: Take Naphthalene, as entered by a madman: `c12ccccc1=CC=CC=2`
    # with the following matching produced for the aromatic region:
    # {(0, 1), (2, 3), (4, 5)}, which will result in non-perfect matchings for
    # the cycle basis.
    double_bonds = {e for e in mol.edges if mol.edges[e].get('order', 1) == 2}
    extended_matching = max_match | double_bonds
    cycles = list(nx.simple_cycles(mol))
    for cycle in cycles:
        is_perfect = nx.is_perfect_matching(mol.subgraph(cycle),
                                            {e for e in extended_matching if
                                             all(v in cycle for v in e)})
        if extended_matching and is_perfect:
            for node in cycle:
                mol.nodes[node]['aromatic'] = True

    # 2) Kekulize "aromatic" fragments that do not show DIME
    nodes_in_cycles = {n for cycle in cycles for n in cycle}
    for node in mol:
        if node not in nodes_in_cycles:
            mol.nodes[node]['aromatic'] = False
    for edge in max_match:
        if not all(mol.nodes[v].get('aromatic', False) for v in edge):
            mol.edges[edge]['order'] = 2

def mark_aromatic_edges(mol):
    """
    Set all bonds between aromatic atoms (attribute 'aromatic' is `True`) to
    1.5. Gives all other bonds that don't have an order yet an order of 1.

    Parameters
    ----------
    mol : nx.Graph
        The molecule.

    Returns
    -------
    None
        `mol` is modified in-place.
    """
    for edge in mol.edges:
        if all(mol.nodes[node].get('aromatic', 'False') for node in edge):
            mol.edges[edge]['order'] = 1.5
        elif 'order' not in mol.edges[edge]:
            mol.edges[edge]['order'] = 1


def correct_aromatic_rings(mol):
    """
    Sets hcount for all atoms, marks aromaticity for all atoms, and the order of
    all aromatic bonds to 1.5.

    Parameters
    ----------
    mol : nx.Graph
        The molecule.

    Returns
    -------
    None
        `mol` is modified in-place.
    """
    nx.set_node_attributes(mol, True, 'aromatic')
    fill_valence(mol)
    mark_aromatic_atoms(mol)
    mark_aromatic_edges(mol)


def increment_bond_orders(molecule, max_bond_order=3):
    """
    Increments bond orders up to what the atom's valence allows.

    Parameters
    ----------
    molecule : nx.Graph
        The molecule to process.
    max_bond_order : number
        The highest bond order allowed to make.

    Returns
    -------
    None
        molecule is modified in-place.
    """
    # Gather the number of open spots for all atoms beforehand, since some
    # might have multiple oxidation states (e.g. S). We don't want to change
    # oxidation state halfway through for some funny reason. It shouldn't be
    # nescessary, but it can't hurt.
    missing_bonds = {}
    for idx in molecule:
        missing_bonds[idx] = max(bonds_missing(molecule, idx), 0)

    for idx, jdx in molecule.edges:
        missing_idx = missing_bonds[idx]
        missing_jdx = missing_bonds[jdx]
        edge_missing = min(missing_idx, missing_jdx)
        current_order = molecule.edges[idx, jdx].get("order", 1)
        if current_order == 1.5:
            continue
        new_order = edge_missing + current_order
        new_order = min(new_order, max_bond_order)
        molecule.edges[idx, jdx]['order'] = new_order
        missing_bonds[idx] -= edge_missing
        missing_bonds[jdx] -= edge_missing

1	# -- coding: utf-8 --
2	# Copyright 2018 Peter C Kroon
3
4	# Licensed under the Apache License, Version 2.0 (the "License");
5	# you may not use this file except in compliance with the License.
6	# You may obtain a copy of the License at
7	#
8	# http://www.apache.org/licenses/LICENSE-2.0
9	#
10	# Unless required by applicable law or agreed to in writing, software
11	# distributed under the License is distributed on an "AS IS" BASIS,
12	# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13	# See the License for the specific language governing permissions and
14	# limitations under the License.
15
16	"""	4✔
17	Contains helper functions for parsing and writing SMILES strings, as well as
18	some convenience functions for adding hydrogens, and detecting aromaticity.
19	"""
20
21	import logging	4✔
22	import re	4✔
23	import operator	4✔
24
25	import networkx as nx	4✔
26
27	LOGGER = logging.getLogger(__name__)	4✔
28
29	ISOTOPE_PATTERN = r'(?P<isotope>[\d]+)?'	4✔
30	ELEMENT_PATTERN = r'(?P<element>b\|c\|n\|o\|s\|p\|as\|se\|\*\|[A-Z][a-z]{0,2})'	4✔
31	STEREO_PATTERN = r'(?P<stereo>@\|@@\|@TH[1-2]\|@AL[1-2]\|@SP[1-3]\|@OH[\d]{1,2}\|'\	4✔
32	r'@TB[\d]{1,2})?'
33	HCOUNT_PATTERN = r'(?P<hcount>H[\d]?)?'	4✔
34	CHARGE_PATTERN = r'(?P<charge>(-\|\+)(\++\|-+\|[\d]{1,2})?)?'	4✔
35	CLASS_PATTERN = r'(?::(?P<class>[\d]+))?'	4✔
36	ATOM_PATTERN = re.compile(r'^\[' + ISOTOPE_PATTERN + ELEMENT_PATTERN +	4✔
37	STEREO_PATTERN + HCOUNT_PATTERN + CHARGE_PATTERN +
38	CLASS_PATTERN + r'\]$')
39
40	VALENCES = {"B": (3,), "C": (4,), "N": (3, 5), "O": (2,), "P": (3, 5),	4✔
41	"S": (2, 4, 6), "F": (1,), "Cl": (1,), "Br": (1,), "I": (1,)}
42
43	AROMATIC_ATOMS = "B C N O P S Se As *".split()	4✔
44
45
46	def parse_atom(atom):	4✔
47	"""
48	Parses a SMILES atom token, and returns a dict with the information.
49
50	Note
51	----
52	Can not deal with stereochemical information yet. This gets discarded.
53
54	Parameters
55	----------
56	atom : str
57	The atom string to interpret. Looks something like one of the
58	following: "C", "c", "[13CH3-1:2]"
59
60	Returns
61	-------
62	dict
63	A dictionary containing at least 'element', 'aromatic', and 'charge'. If
64	present, will also contain 'hcount', 'isotope', and 'class'.
65	"""
66	defaults = {'charge': 0, 'hcount': 0, 'aromatic': False}	4✔
67	if not atom.startswith('[') and not atom.endswith(']'):	4✔
68	if atom != '*':	4✔
69	# Don't specify hcount to signal we don't actually know anything
70	# about it
71	return {'element': atom.capitalize(), 'charge': 0,	4✔
72	'aromatic': atom.islower()}
73	else:
74	return defaults.copy()	4✔
75	match = ATOM_PATTERN.match(atom)	4✔
76	if match is None:	4✔
77	raise ValueError('The atom {} is malformatted'.format(atom))	4✔
78	out = defaults.copy()	4✔
79	out.update({k: v for k, v in match.groupdict().items() if v is not None})	4✔
80
81	if out.get('element', 'X').islower():	4✔
82	out['aromatic'] = True	4✔
83
84	parse_helpers = {	4✔
85	'isotope': int,
86	'element': str.capitalize,
87	'stereo': lambda x: x,
88	'hcount': parse_hcount,
89	'charge': parse_charge,
90	'class': int,
91	'aromatic': lambda x: x,
92	}
93
94	for attr, val_str in out.items():	4✔
95	out[attr] = parse_helpers[attr](val_str)	4✔
96
97	if out['element'] == '*':	4✔
98	del out['element']	4✔
99
100	if out.get('element') == 'H' and out.get('hcount', 0):	4✔
101	raise ValueError("A hydrogen atom can't have hydrogens")	4✔
102
103	if 'stereo' in out:	4✔
104	LOGGER.warning('Atom "%s" contains stereochemical information that will be discarded.', atom)	4✔
105
106	return out	4✔
107
108
109	def format_atom(molecule, node_key, default_element='*'):	4✔
110	"""
111	Formats a node following SMILES conventions. Uses the attributes `element`,
112	`charge`, `hcount`, `stereo`, `isotope` and `class`.
113
114	Parameters
115	----------
116	molecule : nx.Graph
117	The molecule containing the atom.
118	node_key : hashable
119	The node key of the atom in `molecule`.
120	default_element : str
121	The element to use if the attribute is not present in the node.
122
123	Returns
124	-------
125	str
126	The atom as SMILES string.
127	"""
128	node = molecule.nodes[node_key]	4✔
129	name = node.get('element', default_element)	4✔
130	charge = node.get('charge', 0)	4✔
131	hcount = node.get('hcount', 0)	4✔
132	stereo = node.get('stereo', None)	4✔
133	isotope = node.get('isotope', '')	4✔
134	class_ = node.get('class', '')	4✔
135	aromatic = node.get('aromatic', False)	4✔
136	default_h = has_default_h_count(molecule, node_key)	4✔
137
138	if stereo is not None:	4!
139	raise NotImplementedError	×
140
141	if aromatic:	4✔
142	name = name.lower()	4✔
143
144	if (stereo is None and isotope == '' and charge == 0 and default_h and class_ == '' and	4✔
145	(name.lower() in 'b c n o p s *'.split() or name in 'F Cl Br I'.split())):
146	return name	4✔
147
148	if hcount:	4✔
149	hcountstr = 'H'	4✔
150	if hcount > 1:	4✔
151	hcountstr += str(hcount)	4✔
152	else:
153	hcountstr = ''	4✔
154
155	if charge > 0:	4✔
156	chargestr = '+'	4✔
157	if charge > 1:	4✔
158	chargestr += str(charge)	4✔
159	elif charge < 0:	4✔
160	chargestr = '-'	4✔
161	if charge < -1:	4✔
162	chargestr += str(-charge)	4✔
163	else:
164	chargestr = ''	4✔
165
166	if class_ != '':	4✔
167	class_ = ':{}'.format(class_)	4✔
168	fmt = '[{isotope}{name}{stereo}{hcount}{charge}{class_}]'	4✔
169	return fmt.format(isotope=isotope, name=name, stereo='', hcount=hcountstr,	4✔
170	charge=chargestr, class_=class_)
171
172
173	def parse_hcount(hcount_str):	4✔
174	"""
175	Parses a SMILES hydrogen count specifications.
176
177	Parameters
178	----------
179	hcount_str : str
180	The hydrogen count specification to parse.
181
182	Returns
183	-------
184	int
185	The number of hydrogens specified.
186	"""
187	if not hcount_str:	4✔
188	return 0	4✔
189	if hcount_str == 'H':	4✔
190	return 1	4✔
191	return int(hcount_str[1:])	4✔
192
193
194	def parse_charge(charge_str):	4✔
195	"""
196	Parses a SMILES charge specification.
197
198	Parameters
199	----------
200	charge_str : str
201	The charge specification to parse.
202
203	Returns
204	-------
205	int
206	The charge.
207	"""
208	if not charge_str:	4✔
209	return 0	4✔
210	signs = {'-': -1, '+': 1}	4✔
211	sign = signs[charge_str[0]]	4✔
212	if len(charge_str) > 1 and charge_str[1].isdigit():	4✔
213	charge = sign * int(charge_str[1:])	4✔
214	else:
215	charge = sign * charge_str.count(charge_str[0])	4✔
216	return charge	4✔
217
218
219	def add_explicit_hydrogens(mol):	4✔
220	"""
221	Adds explicit hydrogen nodes to `mol`, the amount is determined by the node
222	attribute 'hcount'. Will remove the 'hcount' attribute.
223
224	Parameters
225	----------
226	mol : nx.Graph
227	The molecule to which explicit hydrogens should be added. Is modified
228	in-place.
229
230	Returns
231	-------
232	None
233	`mol` is modified in-place.
234	"""
235	h_atom = parse_atom('[H]')	4✔
236	if 'hcount' in h_atom:	4!
237	del h_atom['hcount']	4✔
238	for n_idx in list(mol.nodes):	4✔
239	hcount = mol.nodes[n_idx].get('hcount', 0)	4✔
240	idxs = range(max(mol) + 1, max(mol) + hcount + 1)	4✔
241	# Get the defaults from parse_atom.
242	mol.add_nodes_from(idxs, **h_atom.copy())	4✔
243	mol.add_edges_from([(n_idx, jdx) for jdx in idxs], order=1)	4✔
244	if 'hcount' in mol.nodes[n_idx]:	4✔
245	del mol.nodes[n_idx]['hcount']	4✔
246
247
248	def remove_explicit_hydrogens(mol):	4✔
249	"""
250	Removes all explicit, simple hydrogens from `mol`. Simple means it is
251	identical to the SMILES string "[H]", and has exactly one bond. Increments
252	'hcount' where appropriate.
253
254	Parameters
255	----------
256	mol : nx.Graph
257	The molecule whose explicit hydrogens should be removed. Is modified
258	in-place.
259
260	Returns
261	-------
262	None
263	`mol` is modified in-place.
264	"""
265	to_remove = set()	4✔
266	# defaults = parse_atom('[H]')
267	for n_idx in mol.nodes:	4✔
268	node = mol.nodes[n_idx]	4✔
269	neighbors = list(mol[n_idx])	4✔
270	# TODO: get these defaults from parsing [H]. But do something smart
271	# with the hcount attribute.
272	if (node.get('charge', 0) == 0 and node.get('element', '') == 'H' and	4✔
273	'isotope' not in node and node.get('class', 0) == 0 and
274	len(neighbors) == 1):
275	neighbor = neighbors[0]	4✔
276	if (mol.nodes[neighbor].get('element', '') == 'H' or	4✔
277	mol.edges[n_idx, neighbor].get('order', 1) != 1):
278	# The molecule is H2, or the bond order is not 1.
279	continue	4✔
280	to_remove.add(n_idx)	4✔
281	mol.nodes[neighbor]['hcount'] = mol.nodes[neighbor].get('hcount', 0) + 1	4✔
282	mol.remove_nodes_from(to_remove)	4✔
283	for n_idx in mol.nodes:	4✔
284	if 'hcount' not in mol.nodes[n_idx]:	4✔
285	mol.nodes[n_idx]['hcount'] = 0	4✔
286
287
288	def fill_valence(mol, respect_hcount=True, respect_bond_order=True,	4✔
289	max_bond_order=3):
290	"""
291	Sets the attribute 'hcount' on all nodes in `mol` that don't have it yet.
292	The value to which it is set is based on the node's 'element', and the
293	number of bonds it has. Default valences are as specified by the global
294	variable VALENCES.
295
296	Parameters
297	----------
298	mol : nx.Graph
299	The molecule whose nodes should get a 'hcount'. Is modified in-place.
300	respect_hcount : bool
301	If True, don't change the hcount on nodes that already have it set.
302	respect_bond_order : bool
303	If False, first try to fill the valence by increasing bond orders, and
304	add hydrogens after.
305	max_bond_order : number
306	Only meaningful if respect_bond_order is False. This is the highest
307	bond order that will be set.
308
309	Returns
310	-------
311	None
312	`mol` is modified in-place.
313	"""
314	if not respect_bond_order:	4✔
315	increment_bond_orders(mol, max_bond_order=max_bond_order)	4✔
316	for n_idx in mol:	4✔
317	node = mol.nodes[n_idx]	4✔
318	if 'hcount' in node and respect_hcount:	4✔
319	continue	4✔
320	missing = max(bonds_missing(mol, n_idx), 0)	4✔
321	node['hcount'] = node.get('hcount', 0) + missing	4✔
322
323
324	def bonds_missing(mol, node_idx, use_order=True):	4✔
325	"""
326	Returns how much the specified node is under valence. If use_order is
327	False, treat all bonds as if they are order 1.
328
329	Parameters
330	----------
331	mol : nx.Graph
332	The molecule.
333	node_idx : hashable
334	The node to look at. Should be in mol.
335	use_order : bool
336	If False, treat all bonds as single.
337
338	Returns
339	-------
340	int
341	The number of missing bonds.
342	"""
343	bonds = _bonds(mol, node_idx, use_order)	4✔
344	bonds += mol.nodes[node_idx].get('hcount', 0)	4✔
345	valence = _valence(mol, node_idx, bonds)	4✔
346	return int(valence - bonds)	4✔
347
348
349	def _valence(mol, node_idx, minimum=0):	4✔
350	"""
351	Returns the valence of the specified node. Since some elements can have
352	multiple valences, give the smallest one that is more than `minimum`.
353
354	Parameters
355	----------
356	mol : nx.Graph
357	The molecule.
358	node_idx : hashable
359	The node to look at. Should be in mol.
360	minimum : int
361	The minimum value of valence.
362
363	Returns
364	-------
365	int
366	The smallest valence of node more than `minimum`.
367	"""
368	element = mol.nodes[node_idx].get('element', '').capitalize()	4✔
369	if element not in VALENCES:	4✔
370	return 0	4✔
371	val = VALENCES.get(element)	4✔
372	try:	4✔
373	val = min(filter(lambda a: a >= minimum, val))	4✔
374	except ValueError: # More bonds than possible	4✔
375	val = max(val)	4✔
376	return val	4✔
377
378
379	def _bonds(mol, node_idx, use_order=True):	4✔
380	"""
381	Returns how many explicit bonds the specified node has. If use_order is
382	False, treat all bonds as if they are order 1.
383
384	Parameters
385	----------
386	mol : nx.Graph
387	The molecule.
388	node_idx : hashable
389	The node to look at. Should be in mol.
390	use_order : bool
391	If False, treat all bonds as single.
392
393	Returns
394	-------
395	int
396	The number of bonds.
397	"""
398	if use_order:	4!
399	bond_orders = map(operator.itemgetter(2),	4✔
400	mol.edges(nbunch=node_idx, data='order', default=1))
401	bonds = sum(bond_orders)	4✔
402	else:
403	bonds = len(mol[node_idx])	×
404	return bonds	4✔
405
406
407	def has_default_h_count(mol, node_idx, use_order=True):	4✔
408	"""
409	Returns whether the hydrogen count for this atom is non-standard.
410
411	Parameters
412	----------
413	mol : nx.Graph
414	The molecule.
415	node_idx : hashable
416	The node to look at. Should be in mol.
417	use_order : bool
418	If False, treat all bonds as single.
419
420	Returns
421	-------
422	bool
423	"""
424	bonds = _bonds(mol, node_idx, use_order)	4✔
425	valence = _valence(mol, node_idx, bonds)	4✔
426	hcount = mol.nodes[node_idx].get('hcount', 0)	4✔
427	return valence - bonds == hcount	4✔
428
429
430	def _hydrogen_neighbours(mol, n_idx):	4✔
431	neighbours = mol[n_idx]	×
432	h_neighbours = 0	×
433	for n_jdx in neighbours:	×
434	if (mol.nodes[n_jdx].get('element', '*') == 'H' and	×
435	mol.edges[n_idx, n_jdx].get('order', 1) == 1):
436	h_neighbours += 1	×
437	return h_neighbours	×
438
439	def _prune_nodes(nodes, mol):	4✔
440	new_nodes = []	4✔
441	for node in nodes:	4✔
442	# all wild card nodes are eligible
443	if mol.nodes[node].get('element', '') == '':	4!
UNCOV 444	new_nodes.append(node)	×
UNCOV 445	continue	×
446	missing = bonds_missing(mol, node, use_order=True) + mol.nodes[node].get('charge', 0)	4✔
447	if missing > 0:	4✔
448	new_nodes.append(node)	4✔
449	return mol.subgraph(new_nodes)	4✔
450
451
452	def mark_aromatic_atoms(mol, strict=True):	4✔
453	"""
454	Properly kekulizes molecules and sets the aromatic attribute.
455
456	Parameters
457	----------
458	mol : nx.Graph
459	The molecule.
460	strict: bool
461	Whether to raise a SyntaxError if the aromatic region of the molecule
462	cannot be kekulized.
463
464	Returns
465	-------
466	None
467	`mol` is modified in-place.
468
469	Raises
470	------
471	SyntaxError
472	If ``strict`` is True and the aromatic region of the molecule cannot
473	be kekulized.
474	"""
475	# prune all nodes from molecule that are eligible and have
476	# full valency
477	arom_atoms = [node for node, aromatic in mol.nodes(data='aromatic') if aromatic]	4✔
478	ds_graph = _prune_nodes(arom_atoms, mol)	4✔
479
480	# set the aromatic attribute to False for all nodes
481	# as a precaution
482	nx.set_node_attributes(mol, False, 'aromatic')	4✔
483
484	sub_ds_graph = mol.subgraph(ds_graph)	4✔
485	max_match = nx.max_weight_matching(sub_ds_graph)	4✔
486	# we check if a maximum matching exists and
487	# if it is perfect. if it is not perfect,
488	# this graph originates from a completely invalid
489	# smiles and we raise an error
490	if strict and not nx.is_perfect_matching(sub_ds_graph, max_match):	4✔
491	msg = "Your molecule is invalid and cannot be kekulized."	4✔
492	raise SyntaxError(msg)	4✔
493
494	# From here we want to do 2 things: 1) dekekulize aromatic rings; and 2)
495	# kekulize "aromatic" fragments that are not in a cycle. We only consider
496	# nodes to be aromatic if they can participate in DIME, which means they
497	# must be in a cycle.
498	# 1) Dekekulize aromatic rings
499	# Add all double bonds to the perfect matching we have, and then for every
500	# simple cycle check if the resulting matching is still perfect. You need
501	# to take simple cycles (rather than the cycle basis), for cases like the
502	# following: Take Naphthalene, as entered by a madman: `c12ccccc1=CC=CC=2`
503	# with the following matching produced for the aromatic region:
504	# {(0, 1), (2, 3), (4, 5)}, which will result in non-perfect matchings for
505	# the cycle basis.
506	double_bonds = {e for e in mol.edges if mol.edges[e].get('order', 1) == 2}	4✔
507	extended_matching = max_match \| double_bonds	4✔
508	cycles = list(nx.simple_cycles(mol))	4✔
509	for cycle in cycles:	4✔
510	is_perfect = nx.is_perfect_matching(mol.subgraph(cycle),	4✔
511	{e for e in extended_matching if
512	all(v in cycle for v in e)})
513	if extended_matching and is_perfect:	4✔
514	for node in cycle:	4✔
515	mol.nodes[node]['aromatic'] = True	4✔
516
517	# 2) Kekulize "aromatic" fragments that do not show DIME
518	nodes_in_cycles = {n for cycle in cycles for n in cycle}	4✔
519	for node in mol:	4✔
520	if node not in nodes_in_cycles:	4✔
521	mol.nodes[node]['aromatic'] = False	4✔
522	for edge in max_match:	4✔
523	if not all(mol.nodes[v].get('aromatic', False) for v in edge):	4✔
524	mol.edges[edge]['order'] = 2	4✔
525
526	def mark_aromatic_edges(mol):	4✔
527	"""
528	Set all bonds between aromatic atoms (attribute 'aromatic' is `True`) to
529	1.5. Gives all other bonds that don't have an order yet an order of 1.
530
531	Parameters
532	----------
533	mol : nx.Graph
534	The molecule.
535
536	Returns
537	-------
538	None
539	`mol` is modified in-place.
540	"""
541	for edge in mol.edges:	4✔
542	if all(mol.nodes[node].get('aromatic', 'False') for node in edge):	4✔
543	mol.edges[edge]['order'] = 1.5	4✔
544	elif 'order' not in mol.edges[edge]:	4✔
545	mol.edges[edge]['order'] = 1	4✔
546
547
548	def correct_aromatic_rings(mol):	4✔
549	"""
550	Sets hcount for all atoms, marks aromaticity for all atoms, and the order of
551	all aromatic bonds to 1.5.
552
553	Parameters
554	----------
555	mol : nx.Graph
556	The molecule.
557
558	Returns
559	-------
560	None
561	`mol` is modified in-place.
562	"""
563	nx.set_node_attributes(mol, True, 'aromatic')	4✔
564	fill_valence(mol)	4✔
565	mark_aromatic_atoms(mol)	4✔
566	mark_aromatic_edges(mol)	4✔
567
568
569	def increment_bond_orders(molecule, max_bond_order=3):	4✔
570	"""
571	Increments bond orders up to what the atom's valence allows.
572
573	Parameters
574	----------
575	molecule : nx.Graph
576	The molecule to process.
577	max_bond_order : number
578	The highest bond order allowed to make.
579
580	Returns
581	-------
582	None
583	molecule is modified in-place.
584	"""
585	# Gather the number of open spots for all atoms beforehand, since some
586	# might have multiple oxidation states (e.g. S). We don't want to change
587	# oxidation state halfway through for some funny reason. It shouldn't be
588	# nescessary, but it can't hurt.
589	missing_bonds = {}	4✔
590	for idx in molecule:	4✔
591	missing_bonds[idx] = max(bonds_missing(molecule, idx), 0)	4✔
592
593	for idx, jdx in molecule.edges:	4✔
594	missing_idx = missing_bonds[idx]	4✔
595	missing_jdx = missing_bonds[jdx]	4✔
596	edge_missing = min(missing_idx, missing_jdx)	4✔
597	current_order = molecule.edges[idx, jdx].get("order", 1)	4✔
598	if current_order == 1.5:	4✔
599	continue	4✔
600	new_order = edge_missing + current_order	4✔
601	new_order = min(new_order, max_bond_order)	4✔
602	molecule.edges[idx, jdx]['order'] = new_order	4✔
603	missing_bonds[idx] -= edge_missing	4✔
604	missing_bonds[jdx] -= edge_missing	4✔

pckroon / pysmiles / 9079427458

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous