9003330176

Committed 08 May 2024 02:12PM UTC coverage: 96.296% (-1.4%) from 97.708%

Build # 9003330176

Build Type

push

travis-ci

Committed by

web-flow

Commit Message

Merge pull request #38 from gruenewald-lab/aromatic_improved

init draft aromatic

Run Details

228 of 243 branches covered (93.83%)

Branch coverage included in aggregate %.

40 of 40 new or added lines in 3 files covered. (100.0%)

6 existing lines in 1 file now uncovered.

422 of 432 relevant lines covered (97.69%)

3.9 hits per line

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

95.55

/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, atoms=None):
    """
    Properly kekeulizes molecules and sets the aromatic attribute.

    Parameters
    ----------
    mol : nx.Graph
        The molecule.
    atoms: collections.abc.Iterable
        The atoms to act on; all other nodes are pruned

    Returns
    -------
    None
        `mol` is modified in-place.
    """
    if atoms is None:
        atoms = set(mol.nodes)
    # prune all nodes from molecule that are eligible and have
    # full valency
    ds_graph = _prune_nodes(atoms, mol)

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

    for sub_ds in nx.connected_components(ds_graph):
        sub_ds_graph = mol.subgraph(sub_ds)
        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 not nx.is_perfect_matching(sub_ds_graph, max_match):
            msg = "Your molecule is invalid and cannot be kekulized."
            raise SyntaxError(msg)

        # we consider a node aromatic if it can take part in DIME
        # to do so all nodes in a delocalized subgraph have to be
        # part of a cycle system
        cycles = nx.cycle_basis(sub_ds_graph)
        nodes_in_cycles = []
        for cycle in cycles:
            nodes_in_cycles += cycle
            nx.set_node_attributes(mol, {node: True for node in cycle}, 'aromatic')

        for node in sub_ds_graph.nodes:
            if node not in nodes_in_cycles:
                mol.nodes[node]['aromatic'] = False

        for edge in max_match:
            if not edge[0] in nodes_in_cycles and not edge[1] in nodes_in_cycles:
                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.
    """
    fill_valence(mol)
    mark_aromatic_atoms(mol, atoms=mol.nodes)
    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✔
UNCOV 431	neighbours = mol[n_idx]	×
UNCOV 432	h_neighbours = 0	×
UNCOV 433	for n_jdx in neighbours:	×
UNCOV 434	if (mol.nodes[n_jdx].get('element', '*') == 'H' and	×
435	mol.edges[n_idx, n_jdx].get('order', 1) == 1):
UNCOV 436	h_neighbours += 1	×
UNCOV 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✔
444	new_nodes.append(node)	4✔
445	continue	4✔
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	def mark_aromatic_atoms(mol, atoms=None):	4✔
452	"""
453	Properly kekeulizes molecules and sets the aromatic attribute.
454
455	Parameters
456	----------
457	mol : nx.Graph
458	The molecule.
459	atoms: collections.abc.Iterable
460	The atoms to act on; all other nodes are pruned
461
462	Returns
463	-------
464	None
465	`mol` is modified in-place.
466	"""
467	if atoms is None:	4✔
468	atoms = set(mol.nodes)	4✔
469	# prune all nodes from molecule that are eligible and have
470	# full valency
471	ds_graph = _prune_nodes(atoms, mol)	4✔
472
473	# set the aromatic attribute to False for all nodes
474	# as a precaution
475	nx.set_node_attributes(mol, False, 'aromatic')	4✔
476
477	for sub_ds in nx.connected_components(ds_graph):	4✔
478	sub_ds_graph = mol.subgraph(sub_ds)	4✔
479	max_match = nx.max_weight_matching(sub_ds_graph)	4✔
480	# we check if a maximum matching exists and
481	# if it is perfect. if it is not perfect,
482	# this graph originates from a completely invalid
483	# smiles and we raise an error
484	if not nx.is_perfect_matching(sub_ds_graph, max_match):	4✔
485	msg = "Your molecule is invalid and cannot be kekulized."	4✔
486	raise SyntaxError(msg)	4✔
487
488	# we consider a node aromatic if it can take part in DIME
489	# to do so all nodes in a delocalized subgraph have to be
490	# part of a cycle system
491	cycles = nx.cycle_basis(sub_ds_graph)	4✔
492	nodes_in_cycles = []	4✔
493	for cycle in cycles:	4✔
494	nodes_in_cycles += cycle	4✔
495	nx.set_node_attributes(mol, {node: True for node in cycle}, 'aromatic')	4✔
496
497	for node in sub_ds_graph.nodes:	4✔
498	if node not in nodes_in_cycles:	4✔
499	mol.nodes[node]['aromatic'] = False	4✔
500
501	for edge in max_match:	4✔
502	if not edge[0] in nodes_in_cycles and not edge[1] in nodes_in_cycles:	4✔
503	mol.edges[edge]['order'] = 2	4✔
504
505	def mark_aromatic_edges(mol):	4✔
506	"""
507	Set all bonds between aromatic atoms (attribute 'aromatic' is `True`) to
508	1.5. Gives all other bonds that don't have an order yet an order of 1.
509
510	Parameters
511	----------
512	mol : nx.Graph
513	The molecule.
514
515	Returns
516	-------
517	None
518	`mol` is modified in-place.
519	"""
520	for edge in mol.edges:	4✔
521	if all(mol.nodes[node].get('aromatic', 'False') for node in edge):	4✔
522	mol.edges[edge]['order'] = 1.5	4✔
523	elif 'order' not in mol.edges[edge]:	4✔
524	mol.edges[edge]['order'] = 1	4✔
525
526	def correct_aromatic_rings(mol):	4✔
527	"""
528	Sets hcount for all atoms, marks aromaticity for all atoms, and the order of
529	all aromatic bonds to 1.5.
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	fill_valence(mol)	4✔
542	mark_aromatic_atoms(mol, atoms=mol.nodes)	4✔
543	mark_aromatic_edges(mol)	4✔
544
545
546	def increment_bond_orders(molecule, max_bond_order=3):	4✔
547	"""
548	Increments bond orders up to what the atom's valence allows.
549
550	Parameters
551	----------
552	molecule : nx.Graph
553	The molecule to process.
554	max_bond_order : number
555	The highest bond order allowed to make.
556
557	Returns
558	-------
559	None
560	molecule is modified in-place.
561	"""
562	# Gather the number of open spots for all atoms beforehand, since some
563	# might have multiple oxidation states (e.g. S). We don't want to change
564	# oxidation state halfway through for some funny reason. It shouldn't be
565	# nescessary, but it can't hurt.
566	missing_bonds = {}	4✔
567	for idx in molecule:	4✔
568	missing_bonds[idx] = max(bonds_missing(molecule, idx), 0)	4✔
569
570	for idx, jdx in molecule.edges:	4✔
571	missing_idx = missing_bonds[idx]	4✔
572	missing_jdx = missing_bonds[jdx]	4✔
573	edge_missing = min(missing_idx, missing_jdx)	4✔
574	current_order = molecule.edges[idx, jdx].get("order", 1)	4✔
575	if current_order == 1.5:	4✔
576	continue	4✔
577	new_order = edge_missing + current_order	4✔
578	new_order = min(new_order, max_bond_order)	4✔
579	molecule.edges[idx, jdx]['order'] = new_order	4✔
580	missing_bonds[idx] -= edge_missing	4✔
581	missing_bonds[jdx] -= edge_missing	4✔

pckroon / pysmiles / 9003330176

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