9003003376

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

Build # 9003003376

Build Type

Pull #38

travis-ci

Committed by

web-flow

Commit Message

Merge ed0261725 into 2308e2409

Pull Request Pull Request #38: init draft aromatic

Run Details

228 of 245 branches covered (93.06%)

Branch coverage included in aggregate %.

41 of 43 new or added lines in 3 files covered. (95.35%)

6 existing lines in 1 file now uncovered.

423 of 435 relevant lines covered (97.24%)

3.89 hits per line

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

94.44

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

pckroon / pysmiles / 9003003376

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