9002551850

Committed 08 May 2024 01:18PM UTC coverage: 95.71% (-2.0%) from 97.708%

Build # 9002551850

Build Type

Pull #38

travis-ci

Committed by

web-flow

Commit Message

Merge 8fed50e3a into 2308e2409

Pull Request Pull Request #38: init draft aromatic

Run Details

226 of 243 branches covered (93.0%)

Branch coverage included in aggregate %.

39 of 41 new or added lines in 3 files covered. (95.12%)

6 existing lines in 1 file now uncovered.

421 of 433 relevant lines covered (97.23%)

3.88 hits per line

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

94.38

/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)
        # if the subgraph is three nodes it might be
        # a triangle, which is the only special case
        # where there is no maximum match but
        is_triangle = (len(sub_ds_graph.nodes) == 3 and nx.cycle_basis(sub_ds_graph))
        if not is_triangle:
            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

        if set(nodes_in_cycles) == set(sub_ds_graph.nodes):
                nx.set_node_attributes(mol, {node: True for node in sub_ds_graph.nodes}, 'aromatic')
        else:
            nx.set_node_attributes(mol, {node: False for node in sub_ds_graph.nodes}, 'aromatic')
            for edge in max_match:
                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	# if the subgraph is three nodes it might be
481	# a triangle, which is the only special case
482	# where there is no maximum match but
483	is_triangle = (len(sub_ds_graph.nodes) == 3 and nx.cycle_basis(sub_ds_graph))	4✔
484	if not is_triangle:	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 not nx.is_perfect_matching(sub_ds_graph, max_match):	4!
NEW 491	msg = "Your molecule is invalid and cannot be kekulized."	×
NEW 492	raise SyntaxError(msg)	×
493
494	# we consider a node aromatic if it can take part in DIME
495	# to do so all nodes in a delocalized subgraph have to be
496	# part of a cycle system
497	cycles = nx.cycle_basis(sub_ds_graph)	4✔
498	nodes_in_cycles = []	4✔
499	for cycle in cycles:	4✔
500	nodes_in_cycles += cycle	4✔
501
502	if set(nodes_in_cycles) == set(sub_ds_graph.nodes):	4✔
503	nx.set_node_attributes(mol, {node: True for node in sub_ds_graph.nodes}, 'aromatic')	4✔
504	else:
505	nx.set_node_attributes(mol, {node: False for node in sub_ds_graph.nodes}, 'aromatic')	4✔
506	for edge in max_match:	4✔
507	mol.edges[edge]['order'] = 2	4✔
508
509	def mark_aromatic_edges(mol):	4✔
510	"""
511	Set all bonds between aromatic atoms (attribute 'aromatic' is `True`) to
512	1.5. Gives all other bonds that don't have an order yet an order of 1.
513
514	Parameters
515	----------
516	mol : nx.Graph
517	The molecule.
518
519	Returns
520	-------
521	None
522	`mol` is modified in-place.
523	"""
524	for edge in mol.edges:	4✔
525	if all(mol.nodes[node].get('aromatic', 'False') for node in edge):	4✔
526	mol.edges[edge]['order'] = 1.5	4✔
527	elif 'order' not in mol.edges[edge]:	4✔
528	mol.edges[edge]['order'] = 1	4✔
529
530	def correct_aromatic_rings(mol):	4✔
531	"""
532	Sets hcount for all atoms, marks aromaticity for all atoms, and the order of
533	all aromatic bonds to 1.5.
534
535	Parameters
536	----------
537	mol : nx.Graph
538	The molecule.
539
540	Returns
541	-------
542	None
543	`mol` is modified in-place.
544	"""
545	fill_valence(mol)	4✔
546	mark_aromatic_atoms(mol, atoms=mol.nodes)	4✔
547	mark_aromatic_edges(mol)	4✔
548
549
550	def increment_bond_orders(molecule, max_bond_order=3):	4✔
551	"""
552	Increments bond orders up to what the atom's valence allows.
553
554	Parameters
555	----------
556	molecule : nx.Graph
557	The molecule to process.
558	max_bond_order : number
559	The highest bond order allowed to make.
560
561	Returns
562	-------
563	None
564	molecule is modified in-place.
565	"""
566	# Gather the number of open spots for all atoms beforehand, since some
567	# might have multiple oxidation states (e.g. S). We don't want to change
568	# oxidation state halfway through for some funny reason. It shouldn't be
569	# nescessary, but it can't hurt.
570	missing_bonds = {}	4✔
571	for idx in molecule:	4✔
572	missing_bonds[idx] = max(bonds_missing(molecule, idx), 0)	4✔
573
574	for idx, jdx in molecule.edges:	4✔
575	missing_idx = missing_bonds[idx]	4✔
576	missing_jdx = missing_bonds[jdx]	4✔
577	edge_missing = min(missing_idx, missing_jdx)	4✔
578	current_order = molecule.edges[idx, jdx].get("order", 1)	4✔
579	if current_order == 1.5:	4✔
580	continue	4✔
581	new_order = edge_missing + current_order	4✔
582	new_order = min(new_order, max_bond_order)	4✔
583	molecule.edges[idx, jdx]['order'] = new_order	4✔
584	missing_bonds[idx] -= edge_missing	4✔
585	missing_bonds[jdx] -= edge_missing	4✔

pckroon / pysmiles / 9002551850

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