8940101947

Committed 03 May 2024 01:55PM UTC coverage: 95.256% (-2.5%) from 97.708%

Build # 8940101947

Build Type

Pull #39

travis-ci

Committed by

web-flow

Commit Message

Merge bdbee7b10 into 2308e2409

Pull Request Pull Request #39: EZ isomerism and chirality

Run Details

264 of 287 branches covered (91.99%)

Branch coverage included in aggregate %.

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

2 existing lines in 1 file now uncovered.

479 of 493 relevant lines covered (97.16%)

3.88 hits per line

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

93.77

/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 mark_aromatic_atoms(mol, atoms=None):
    """
    Sets the 'aromatic' attribute for all nodes in `mol`. Requires that
    the 'hcount' on atoms is correct.

    Parameters
    ----------
    mol : nx.Graph
        The molecule.
    atoms: collections.abc.Iterable
        The atoms to act on. Will still analyse the full molecule.

    Returns
    -------
    None
        `mol` is modified in-place.
    """
    if atoms is None:
        atoms = set(mol.nodes)
    aromatic = set()
    # Only cycles can be aromatic
    for cycle in nx.cycle_basis(mol):
        # All atoms should be sp2, so each contributes an electron. We make
        # sure they are later.
        electrons = len(cycle)
        maybe_aromatic = True

        for node_idx in cycle:
            node = mol.nodes[node_idx]
            element = node.get('element', '*').capitalize()
            hcount = node.get('hcount', 0)
            degree = mol.degree(node_idx) + hcount
            hcount += _hydrogen_neighbours(mol, node_idx)
            # Make sure they are possibly aromatic, and are sp2 hybridized
            if element not in AROMATIC_ATOMS or degree not in (2, 3):
                maybe_aromatic = False
                break
            # Some of the special cases per group. N and O type atoms can
            # donate an additional electron from a lone pair.
            # missing cases:
            #   extracyclic sp2 heteroatom (e.g. =O)
            #   some charged cases
            if element in 'N P As'.split() and hcount == 1:
                electrons += 1
            elif element in 'O S Se'.split():
                electrons += 1
            if node.get('charge', 0) == +1 and not (element == 'C' and hcount == 0):
                electrons -= 1
        if maybe_aromatic and int(electrons) % 2 == 0:
            # definitely (anti) aromatic
            aromatic.update(cycle)
    for node_idx in atoms:
        node = mol.nodes[node_idx]
        if node_idx not in aromatic:
            node['aromatic'] = False
        else:
            node['aromatic'] = True


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 cycle in nx.cycle_basis(mol):
        for idx, jdx in mol.edges(nbunch=cycle):
            if idx not in cycle or jdx not in cycle:
                continue
            if (mol.nodes[idx].get('aromatic', False)
                    and mol.nodes[jdx].get('aromatic', False)):
                mol.edges[idx, jdx]['order'] = 1.5
    for idx, jdx in mol.edges:
        if 'order' not in mol.edges[idx, jdx]:
            mol.edges[idx, jdx]['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)
    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

def mark_chiral_atoms(molecule):
    """
    For all nodes tagged as chiral, figure out the three
    substiuents and annotate the node with a tuple that
    has the order in which to rotate. This essentially
    corresponds to the definition of an improper dihedral
    angle centered on the chiral atom.
    """
    chiral_nodes = nx.get_node_attributes(molecule, 'stereo')
    for node, direction in chiral_nodes.items():
        # discard the node with the lowest index because we
        # are looking down that node; it is not part of the
        # stereo projection
        neighbors = sorted(molecule.neighbors(node))[1:]
        if len(neighbors) != 3:
            n = len(neighbors) + 1
            msg = (f"Chiral node {node} has {n} neighbors, which "
                    "is different than the four expected for "
                    "tetrahedral chirality. If the chiral center "
                    "contains hydrogen you need to set explicit_hydrogen "
                    "to True.")
            raise ValueError(msg)
        # the default is anti-clockwise sorting indicated by '@'
        # in this case the nodes are sorted with increasing
        # node index; however @@ means clockwise and the
        # order of nodes is reversed (i.e. with decreasing index)
        if direction == '@@':
            neighbors = [neighbors[0], neighbors[2], neighbors[1]]
        molecule.nodes[node]['stereo'] = (node, *neighbors)

def annotate_ez_isomers(molecule, ez_pairs):
    for first, second in ez_pairs:
        ligand_first, anchor_first, ez_first = first
        ligand_second, anchor_second, ez_second = second
        # here come all the ridcilous cases of EZ in smiles
        if ligand_first < anchor_first:
            # case 1 F/C=C/F is trans
            if ez_first == '/' and ez_second == '/':
                ez_isomer = 'trans'
            # case 2 F\C=C/F
            elif ez_first == '\\' and ez_second == '/':
                ez_isomer = 'cis'
            # case 3 F/C=C\F
            elif ez_first == '/' and ez_second == '\\':
                ez_isomer = 'cis'
            # case 4 F\C=C\F
            elif ez_fist == '\\' and ez_second == '\\':
                ez_isomer = 'trans'
        # as in C(/F)
        elif ligand_first > anchor_first:
            # case 5 C(\F)=C/F is trans
            if ez_first == '\\' and ez_second == '/':
                ez_isomer = 'trans'
            # case 6 C(/F)=C/F
            elif ez_first == '/' and ez_second == '/':
                ez_isomer = 'cis'
            # case 7 C(/F)=C\F
            elif ez_first == '/' and ez_second == '\\':
                ez_isomer = 'trans'
            # case 4 C(\F)=C\F
            elif ez_fist == '\\' and ez_second == '\\':
                ez_isomer = 'cis'
        # annotate ligands
        print(ligand_first, ligand_second)
        molecule.nodes[ligand_first]['ez_isomer'] = (ligand_first,
                                                     anchor_first,
                                                     anchor_second,
                                                     ligand_second,
                                                     ez_isomer)
        molecule.nodes[ligand_second]['ez_isomer'] = (ligand_second,
                                                      anchor_second,
                                                      anchor_first,
                                                      ligand_first,
                                                      ez_isomer)

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	# if 'stereo' in out:
103	# LOGGER.warning('Atom "%s" contains stereochemical information that will be discarded.', atom)
104
105	return out	4✔
106
107
108	def format_atom(molecule, node_key, default_element='*'):	4✔
109	"""
110	Formats a node following SMILES conventions. Uses the attributes `element`,
111	`charge`, `hcount`, `stereo`, `isotope` and `class`.
112
113	Parameters
114	----------
115	molecule : nx.Graph
116	The molecule containing the atom.
117	node_key : hashable
118	The node key of the atom in `molecule`.
119	default_element : str
120	The element to use if the attribute is not present in the node.
121
122	Returns
123	-------
124	str
125	The atom as SMILES string.
126	"""
127	node = molecule.nodes[node_key]	4✔
128	name = node.get('element', default_element)	4✔
129	charge = node.get('charge', 0)	4✔
130	hcount = node.get('hcount', 0)	4✔
131	stereo = node.get('stereo', None)	4✔
132	isotope = node.get('isotope', '')	4✔
133	class_ = node.get('class', '')	4✔
134	aromatic = node.get('aromatic', False)	4✔
135	default_h = has_default_h_count(molecule, node_key)	4✔
136
137	if stereo is not None:	4!
UNCOV 138	raise NotImplementedError	×
139
140	if aromatic:	4✔
141	name = name.lower()	4✔
142
143	if (stereo is None and isotope == '' and charge == 0 and default_h and class_ == '' and	4✔
144	(name.lower() in 'b c n o p s *'.split() or name in 'F Cl Br I'.split())):
145	return name	4✔
146
147	if hcount:	4✔
148	hcountstr = 'H'	4✔
149	if hcount > 1:	4✔
150	hcountstr += str(hcount)	4✔
151	else:
152	hcountstr = ''	4✔
153
154	if charge > 0:	4✔
155	chargestr = '+'	4✔
156	if charge > 1:	4✔
157	chargestr += str(charge)	4✔
158	elif charge < 0:	4✔
159	chargestr = '-'	4✔
160	if charge < -1:	4✔
161	chargestr += str(-charge)	4✔
162	else:
163	chargestr = ''	4✔
164
165	if class_ != '':	4✔
166	class_ = ':{}'.format(class_)	4✔
167	fmt = '[{isotope}{name}{stereo}{hcount}{charge}{class_}]'	4✔
168	return fmt.format(isotope=isotope, name=name, stereo='', hcount=hcountstr,	4✔
169	charge=chargestr, class_=class_)
170
171
172	def parse_hcount(hcount_str):	4✔
173	"""
174	Parses a SMILES hydrogen count specifications.
175
176	Parameters
177	----------
178	hcount_str : str
179	The hydrogen count specification to parse.
180
181	Returns
182	-------
183	int
184	The number of hydrogens specified.
185	"""
186	if not hcount_str:	4✔
187	return 0	4✔
188	if hcount_str == 'H':	4✔
189	return 1	4✔
190	return int(hcount_str[1:])	4✔
191
192
193	def parse_charge(charge_str):	4✔
194	"""
195	Parses a SMILES charge specification.
196
197	Parameters
198	----------
199	charge_str : str
200	The charge specification to parse.
201
202	Returns
203	-------
204	int
205	The charge.
206	"""
207	if not charge_str:	4✔
208	return 0	4✔
209	signs = {'-': -1, '+': 1}	4✔
210	sign = signs[charge_str[0]]	4✔
211	if len(charge_str) > 1 and charge_str[1].isdigit():	4✔
212	charge = sign * int(charge_str[1:])	4✔
213	else:
214	charge = sign * charge_str.count(charge_str[0])	4✔
215	return charge	4✔
216
217
218	def add_explicit_hydrogens(mol):	4✔
219	"""
220	Adds explicit hydrogen nodes to `mol`, the amount is determined by the node
221	attribute 'hcount'. Will remove the 'hcount' attribute.
222
223	Parameters
224	----------
225	mol : nx.Graph
226	The molecule to which explicit hydrogens should be added. Is modified
227	in-place.
228
229	Returns
230	-------
231	None
232	`mol` is modified in-place.
233	"""
234	h_atom = parse_atom('[H]')	4✔
235	if 'hcount' in h_atom:	4!
236	del h_atom['hcount']	4✔
237	for n_idx in list(mol.nodes):	4✔
238	hcount = mol.nodes[n_idx].get('hcount', 0)	4✔
239	idxs = range(max(mol) + 1, max(mol) + hcount + 1)	4✔
240	# Get the defaults from parse_atom.
241	mol.add_nodes_from(idxs, **h_atom.copy())	4✔
242	mol.add_edges_from([(n_idx, jdx) for jdx in idxs], order=1)	4✔
243	if 'hcount' in mol.nodes[n_idx]:	4✔
244	del mol.nodes[n_idx]['hcount']	4✔
245
246
247	def remove_explicit_hydrogens(mol):	4✔
248	"""
249	Removes all explicit, simple hydrogens from `mol`. Simple means it is
250	identical to the SMILES string "[H]", and has exactly one bond. Increments
251	'hcount' where appropriate.
252
253	Parameters
254	----------
255	mol : nx.Graph
256	The molecule whose explicit hydrogens should be removed. Is modified
257	in-place.
258
259	Returns
260	-------
261	None
262	`mol` is modified in-place.
263	"""
264	to_remove = set()	4✔
265	# defaults = parse_atom('[H]')
266	for n_idx in mol.nodes:	4✔
267	node = mol.nodes[n_idx]	4✔
268	neighbors = list(mol[n_idx])	4✔
269	# TODO: get these defaults from parsing [H]. But do something smart
270	# with the hcount attribute.
271	if (node.get('charge', 0) == 0 and node.get('element', '') == 'H' and	4✔
272	'isotope' not in node and node.get('class', 0) == 0 and
273	len(neighbors) == 1):
274	neighbor = neighbors[0]	4✔
275	if (mol.nodes[neighbor].get('element', '') == 'H' or	4✔
276	mol.edges[n_idx, neighbor].get('order', 1) != 1):
277	# The molecule is H2, or the bond order is not 1.
278	continue	4✔
279	to_remove.add(n_idx)	4✔
280	mol.nodes[neighbor]['hcount'] = mol.nodes[neighbor].get('hcount', 0) + 1	4✔
281	mol.remove_nodes_from(to_remove)	4✔
282	for n_idx in mol.nodes:	4✔
283	if 'hcount' not in mol.nodes[n_idx]:	4✔
284	mol.nodes[n_idx]['hcount'] = 0	4✔
285
286
287	def fill_valence(mol, respect_hcount=True, respect_bond_order=True,	4✔
288	max_bond_order=3):
289	"""
290	Sets the attribute 'hcount' on all nodes in `mol` that don't have it yet.
291	The value to which it is set is based on the node's 'element', and the
292	number of bonds it has. Default valences are as specified by the global
293	variable VALENCES.
294
295	Parameters
296	----------
297	mol : nx.Graph
298	The molecule whose nodes should get a 'hcount'. Is modified in-place.
299	respect_hcount : bool
300	If True, don't change the hcount on nodes that already have it set.
301	respect_bond_order : bool
302	If False, first try to fill the valence by increasing bond orders, and
303	add hydrogens after.
304	max_bond_order : number
305	Only meaningful if respect_bond_order is False. This is the highest
306	bond order that will be set.
307
308	Returns
309	-------
310	None
311	`mol` is modified in-place.
312	"""
313	if not respect_bond_order:	4✔
314	increment_bond_orders(mol, max_bond_order=max_bond_order)	4✔
315	for n_idx in mol:	4✔
316	node = mol.nodes[n_idx]	4✔
317	if 'hcount' in node and respect_hcount:	4✔
318	continue	4✔
319	missing = max(bonds_missing(mol, n_idx), 0)	4✔
320	node['hcount'] = node.get('hcount', 0) + missing	4✔
321
322
323	def bonds_missing(mol, node_idx, use_order=True):	4✔
324	"""
325	Returns how much the specified node is under valence. If use_order is
326	False, treat all bonds as if they are order 1.
327
328	Parameters
329	----------
330	mol : nx.Graph
331	The molecule.
332	node_idx : hashable
333	The node to look at. Should be in mol.
334	use_order : bool
335	If False, treat all bonds as single.
336
337	Returns
338	-------
339	int
340	The number of missing bonds.
341	"""
342	bonds = _bonds(mol, node_idx, use_order)	4✔
343	bonds += mol.nodes[node_idx].get('hcount', 0)	4✔
344	valence = _valence(mol, node_idx, bonds)	4✔
345	return int(valence - bonds)	4✔
346
347
348	def _valence(mol, node_idx, minimum=0):	4✔
349	"""
350	Returns the valence of the specified node. Since some elements can have
351	multiple valences, give the smallest one that is more than `minimum`.
352
353	Parameters
354	----------
355	mol : nx.Graph
356	The molecule.
357	node_idx : hashable
358	The node to look at. Should be in mol.
359	minimum : int
360	The minimum value of valence.
361
362	Returns
363	-------
364	int
365	The smallest valence of node more than `minimum`.
366	"""
367	element = mol.nodes[node_idx].get('element', '').capitalize()	4✔
368	if element not in VALENCES:	4✔
369	return 0	4✔
370	val = VALENCES.get(element)	4✔
371	try:	4✔
372	val = min(filter(lambda a: a >= minimum, val))	4✔
373	except ValueError: # More bonds than possible	4✔
374	val = max(val)	4✔
375	return val	4✔
376
377
378	def _bonds(mol, node_idx, use_order=True):	4✔
379	"""
380	Returns how many explicit bonds the specified node has. If use_order is
381	False, treat all bonds as if they are order 1.
382
383	Parameters
384	----------
385	mol : nx.Graph
386	The molecule.
387	node_idx : hashable
388	The node to look at. Should be in mol.
389	use_order : bool
390	If False, treat all bonds as single.
391
392	Returns
393	-------
394	int
395	The number of bonds.
396	"""
397	if use_order:	4!
398	bond_orders = map(operator.itemgetter(2),	4✔
399	mol.edges(nbunch=node_idx, data='order', default=1))
400	bonds = sum(bond_orders)	4✔
401	else:
UNCOV 402	bonds = len(mol[node_idx])	×
403	return bonds	4✔
404
405
406	def has_default_h_count(mol, node_idx, use_order=True):	4✔
407	"""
408	Returns whether the hydrogen count for this atom is non-standard.
409
410	Parameters
411	----------
412	mol : nx.Graph
413	The molecule.
414	node_idx : hashable
415	The node to look at. Should be in mol.
416	use_order : bool
417	If False, treat all bonds as single.
418
419	Returns
420	-------
421	bool
422	"""
423	bonds = _bonds(mol, node_idx, use_order)	4✔
424	valence = _valence(mol, node_idx, bonds)	4✔
425	hcount = mol.nodes[node_idx].get('hcount', 0)	4✔
426	return valence - bonds == hcount	4✔
427
428
429	def _hydrogen_neighbours(mol, n_idx):	4✔
430	neighbours = mol[n_idx]	4✔
431	h_neighbours = 0	4✔
432	for n_jdx in neighbours:	4✔
433	if (mol.nodes[n_jdx].get('element', '*') == 'H' and	4✔
434	mol.edges[n_idx, n_jdx].get('order', 1) == 1):
435	h_neighbours += 1	4✔
436	return h_neighbours	4✔
437
438
439	def mark_aromatic_atoms(mol, atoms=None):	4✔
440	"""
441	Sets the 'aromatic' attribute for all nodes in `mol`. Requires that
442	the 'hcount' on atoms is correct.
443
444	Parameters
445	----------
446	mol : nx.Graph
447	The molecule.
448	atoms: collections.abc.Iterable
449	The atoms to act on. Will still analyse the full molecule.
450
451	Returns
452	-------
453	None
454	`mol` is modified in-place.
455	"""
456	if atoms is None:	4!
457	atoms = set(mol.nodes)	4✔
458	aromatic = set()	4✔
459	# Only cycles can be aromatic
460	for cycle in nx.cycle_basis(mol):	4✔
461	# All atoms should be sp2, so each contributes an electron. We make
462	# sure they are later.
463	electrons = len(cycle)	4✔
464	maybe_aromatic = True	4✔
465
466	for node_idx in cycle:	4✔
467	node = mol.nodes[node_idx]	4✔
468	element = node.get('element', '*').capitalize()	4✔
469	hcount = node.get('hcount', 0)	4✔
470	degree = mol.degree(node_idx) + hcount	4✔
471	hcount += _hydrogen_neighbours(mol, node_idx)	4✔
472	# Make sure they are possibly aromatic, and are sp2 hybridized
473	if element not in AROMATIC_ATOMS or degree not in (2, 3):	4✔
474	maybe_aromatic = False	4✔
475	break	4✔
476	# Some of the special cases per group. N and O type atoms can
477	# donate an additional electron from a lone pair.
478	# missing cases:
479	# extracyclic sp2 heteroatom (e.g. =O)
480	# some charged cases
481	if element in 'N P As'.split() and hcount == 1:	4✔
482	electrons += 1	4✔
483	elif element in 'O S Se'.split():	4✔
484	electrons += 1	4✔
485	if node.get('charge', 0) == +1 and not (element == 'C' and hcount == 0):	4✔
486	electrons -= 1	4✔
487	if maybe_aromatic and int(electrons) % 2 == 0:	4✔
488	# definitely (anti) aromatic
489	aromatic.update(cycle)	4✔
490	for node_idx in atoms:	4✔
491	node = mol.nodes[node_idx]	4✔
492	if node_idx not in aromatic:	4✔
493	node['aromatic'] = False	4✔
494	else:
495	node['aromatic'] = True	4✔
496
497
498	def mark_aromatic_edges(mol):	4✔
499	"""
500	Set all bonds between aromatic atoms (attribute 'aromatic' is `True`) to
501	1.5. Gives all other bonds that don't have an order yet an order of 1.
502
503	Parameters
504	----------
505	mol : nx.Graph
506	The molecule.
507
508	Returns
509	-------
510	None
511	`mol` is modified in-place.
512	"""
513	for cycle in nx.cycle_basis(mol):	4✔
514	for idx, jdx in mol.edges(nbunch=cycle):	4✔
515	if idx not in cycle or jdx not in cycle:	4✔
516	continue	4✔
517	if (mol.nodes[idx].get('aromatic', False)	4✔
518	and mol.nodes[jdx].get('aromatic', False)):
519	mol.edges[idx, jdx]['order'] = 1.5	4✔
520	for idx, jdx in mol.edges:	4✔
521	if 'order' not in mol.edges[idx, jdx]:	4✔
522	mol.edges[idx, jdx]['order'] = 1	4✔
523
524
525	def correct_aromatic_rings(mol):	4✔
526	"""
527	Sets hcount for all atoms, marks aromaticity for all atoms, and the order of
528	all aromatic bonds to 1.5.
529
530	Parameters
531	----------
532	mol : nx.Graph
533	The molecule.
534
535	Returns
536	-------
537	None
538	`mol` is modified in-place.
539	"""
540	fill_valence(mol)	4✔
541	mark_aromatic_atoms(mol)	4✔
542	mark_aromatic_edges(mol)	4✔
543
544
545	def increment_bond_orders(molecule, max_bond_order=3):	4✔
546	"""
547	Increments bond orders up to what the atom's valence allows.
548
549	Parameters
550	----------
551	molecule : nx.Graph
552	The molecule to process.
553	max_bond_order : number
554	The highest bond order allowed to make.
555
556	Returns
557	-------
558	None
559	molecule is modified in-place.
560	"""
561	# Gather the number of open spots for all atoms beforehand, since some
562	# might have multiple oxidation states (e.g. S). We don't want to change
563	# oxidation state halfway through for some funny reason. It shouldn't be
564	# nescessary, but it can't hurt.
565	missing_bonds = {}	4✔
566	for idx in molecule:	4✔
567	missing_bonds[idx] = max(bonds_missing(molecule, idx), 0)	4✔
568
569	for idx, jdx in molecule.edges:	4✔
570	missing_idx = missing_bonds[idx]	4✔
571	missing_jdx = missing_bonds[jdx]	4✔
572	edge_missing = min(missing_idx, missing_jdx)	4✔
573	current_order = molecule.edges[idx, jdx].get("order", 1)	4✔
574	if current_order == 1.5:	4✔
575	continue	4✔
576	new_order = edge_missing + current_order	4✔
577	new_order = min(new_order, max_bond_order)	4✔
578	molecule.edges[idx, jdx]['order'] = new_order	4✔
579	missing_bonds[idx] -= edge_missing	4✔
580	missing_bonds[jdx] -= edge_missing	4✔
581
582	def mark_chiral_atoms(molecule):	4✔
583	"""
584	For all nodes tagged as chiral, figure out the three
585	substiuents and annotate the node with a tuple that
586	has the order in which to rotate. This essentially
587	corresponds to the definition of an improper dihedral
588	angle centered on the chiral atom.
589	"""
590	chiral_nodes = nx.get_node_attributes(molecule, 'stereo')	4✔
591	for node, direction in chiral_nodes.items():	4✔
592	# discard the node with the lowest index because we
593	# are looking down that node; it is not part of the
594	# stereo projection
595	neighbors = sorted(molecule.neighbors(node))[1:]	4✔
596	if len(neighbors) != 3:	4!
NEW 597	n = len(neighbors) + 1	×
NEW 598	msg = (f"Chiral node {node} has {n} neighbors, which "	×
599	"is different than the four expected for "
600	"tetrahedral chirality. If the chiral center "
601	"contains hydrogen you need to set explicit_hydrogen "
602	"to True.")
NEW 603	raise ValueError(msg)	×
604	# the default is anti-clockwise sorting indicated by '@'
605	# in this case the nodes are sorted with increasing
606	# node index; however @@ means clockwise and the
607	# order of nodes is reversed (i.e. with decreasing index)
608	if direction == '@@':	4✔
609	neighbors = [neighbors[0], neighbors[2], neighbors[1]]	4✔
610	molecule.nodes[node]['stereo'] = (node, *neighbors)	4✔
611
612	def annotate_ez_isomers(molecule, ez_pairs):	4✔
613	for first, second in ez_pairs:	4✔
614	ligand_first, anchor_first, ez_first = first	4✔
615	ligand_second, anchor_second, ez_second = second	4✔
616	# here come all the ridcilous cases of EZ in smiles
617	if ligand_first < anchor_first:	4✔
618	# case 1 F/C=C/F is trans
619	if ez_first == '/' and ez_second == '/':	4✔
620	ez_isomer = 'trans'	4✔
621	# case 2 F\C=C/F
622	elif ez_first == '\\' and ez_second == '/':	4✔
623	ez_isomer = 'cis'	4✔
624	# case 3 F/C=C\F
625	elif ez_first == '/' and ez_second == '\\':	4!
626	ez_isomer = 'cis'	4✔
627	# case 4 F\C=C\F
NEW 628	elif ez_fist == '\\' and ez_second == '\\':	×
NEW 629	ez_isomer = 'trans'	×
630	# as in C(/F)
631	elif ligand_first > anchor_first:	4!
632	# case 5 C(\F)=C/F is trans
633	if ez_first == '\\' and ez_second == '/':	4✔
634	ez_isomer = 'trans'	4✔
635	# case 6 C(/F)=C/F
636	elif ez_first == '/' and ez_second == '/':	4!
637	ez_isomer = 'cis'	4✔
638	# case 7 C(/F)=C\F
NEW 639	elif ez_first == '/' and ez_second == '\\':	×
NEW 640	ez_isomer = 'trans'	×
641	# case 4 C(\F)=C\F
NEW 642	elif ez_fist == '\\' and ez_second == '\\':	×
NEW 643	ez_isomer = 'cis'	×
644	# annotate ligands
645	print(ligand_first, ligand_second)	4✔
646	molecule.nodes[ligand_first]['ez_isomer'] = (ligand_first,	4✔
647	anchor_first,
648	anchor_second,
649	ligand_second,
650	ez_isomer)
651	molecule.nodes[ligand_second]['ez_isomer'] = (ligand_second,	4✔
652	anchor_second,
653	anchor_first,
654	ligand_first,
655	ez_isomer)

pckroon / pysmiles / 8940101947

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