"""
Module for R-group enumeration.
"""
import collections
import csv
import io
import itertools
import scipy
from schrodinger.infra import fast3d
from schrodinger.infra import mm
from schrodinger.infra.mmbitset import Bitset
from schrodinger.structutils import analyze
from schrodinger.structutils import rmsd
from schrodinger.ui.qt.filter_dialog_dir import filter_core
from schrodinger.utils import csv_unicode
from schrodinger.utils import ligfilter
from schrodinger.utils import log
from schrodinger.utils.cgxutils import _get_stereo
from schrodinger.utils.cgxutils import _order
from schrodinger.utils.cgxutils import recompute_stereo
logger = log.get_output_logger('rgroup_enumerate')
"""
RGroup properties:
- atom_index: index of the atom bound to the core (aka the "leaving atom")
- source_index: index of the R-group source that will replace this group
- leaving_atoms: list of all the atoms in the leaving group (by index)
- staying_atom: index of the core atom bound to the leaving group
"""
RGroup = collections.namedtuple('RGroup', [
'atom_index', 'source_index', 'leaving_atoms', 'staying_atom', 'bond_order'
])
# TODO: the RGroup namedtuple should probably be renamed LeavingGroup
# to reflect its current usage better.
"""
RGroupSource properties:
- source_index: index of the R-group source that will replace this group
- length: number of atom indices divided by concentration
- rgroups_indices: list of R group indices for each source
"""
RGroupSource = collections.namedtuple(
'RGroupSource', ['source_index', 'length', 'rgroups_indices'])
# Descriptors that this module cares to compute. Used to implement the -filter
# option of r_group_enumerate.py.
DESCRIPTOR_DICT = {
# name: (func, limiter))
'r_rge_Molecular_weight': (ligfilter.Molecular_weight, (0.0, 1000.0)),
'i_rge_Num_rings': (ligfilter.Num_rings, (0, 10)),
'i_rge_Num_aromatic_rings': (ligfilter.Num_aromatic_rings, (0, 10)),
'i_rge_Num_aliphatic_rings': (ligfilter.Num_aliphatic_rings, (0, 10)),
'i_rge_Num_heteroaromatic_rings': (ligfilter.Num_heteroaromatic_rings,
(0, 10)),
'i_rge_Num_rotatable_bonds': (ligfilter.Num_rotatable_bonds, (0, 50)),
'i_rge_Num_atoms': (ligfilter.Num_atoms, (0, 500)),
'i_rge_Num_chiral_centers': (ligfilter.Num_chiral_centers, (0, 10)),
'i_rge_Total_charge': (ligfilter.Total_charge, (-5, 5)),
'i_rge_Num_positive_atoms': (ligfilter.Num_positive_atoms, (0, 10)),
'i_rge_Num_negative_atoms': (ligfilter.Num_negative_atoms, (0, 10)),
'i_rge_Num_heavy_atoms': (ligfilter.Num_heavy_atoms, (0, 300)),
}
DUMMY_ATOMIC_NUMBER = -2
[docs]class RgroupError(Exception):
"""
Exception class for errors specific to this module, which the caller
may want to present to the user as a simple error message, as opposed
to a traceback. This is meant for "user errors", as opposed to bugs;
for example, when an input structure doesn't fulfill the requirements.
"""
[docs]class BondOrderMismatch(RgroupError):
"""
A specific kind of error that we'll ignore to allow libraries that
include R-groups with different bond orders.
"""
[docs]class RgroupEnumerator(object):
"""
Enumerate a structure using R-group sources.
A source is a sequence with an iterable of Structure as its first element,
followed by one or more core atom atom indices where the side chains from
the source should be inserted.
RgroupEnumerator objects are iterable. Example::
sources = [
(StructureReader('r1.maegz'), 4, 12),
(StructureReader('r2.maegz'), 8),
]
for prod_st in RgroupEnumerator(core_st, sources):
...
will use the first reader to replace atoms 4 and 12 in an homo fashion
(meaning that for a given product, the groups attached to atoms 4 and 12 are
always the same), in combination with the structures for the second reader
for atom 8.
The generated structures have the title of the core structure and the title
of each of the R-groups, encoded in CSV format. For ease of parsing, this
information is also stored as separate properties: i_rge_num_r_groups has
the number of R groups, and the title of each is goes in properties
r_rge_R1, r_rge_R2, etc.
As an option, all CT properties from the R groups can be copied to each
product molecule. These properties have the original name prefixed with
<type-char>_rge_R<index>_; for example, r_i_glide_gscore for the first
R group becomes r_rge_R1_r_i_glide_gscore.
The structures in each R-group source should each have one dummy atom
(symbol '', atomic number zero).
The user of the class can request only a slice of the full set of
combinations to be yielded, by providing the optional 'start' and 'stop'
constructor arguments. These follow the standard Python slicing convention.
If the core structure came from an SDF file with R-group labels
("M RGP" lines), the attachment atoms don't need to be specified; the
labels from the file can be used implicitly.
"""
[docs] def __init__(self,
core_st,
sources,
optimize_sidechains=True,
deduplicate=True,
start=0,
stop=None,
copy_properties=False,
enumerate_cistrans=True,
yield_renum_maps=False,
concentrations=None):
"""
Initialize an enumerator for a given core structure and specification of
rgroup sources.
:param core_st: core structure
:type core_st: `schrodinger.structure.Structure`
:param sources: side chain sources. See class description for details.
:type sources: list of list
:param optimize_sidechains: if true, generate 3D coordinates for the
side chain atoms using Fast3D. The input coordinates will only be
used for determining stereochemistry. If false, position the side
chains using rigid rotation and translation (and an arbitrary
torsional angle around the new bond).
:type optimize_sidechains: bool
:param deduplicate: use unique SMILES to identify
and reject duplicate products
:type deduplicate: bool
:param start: beginning of results slice (used by subjobs)
:type start: int
:param stop: end of results slice (used by subjobs)
:type stop: int or None
:param copy_properties: if true, copy all CT properties from each
R-group to the constructed molecule.
:type copy_properties: bool
:param enumerate_cistrans: if True (default), emit both cis and
trans isomers for double-bonded R-groups.
:type enumerate_cistrans: bool
:param yield_renum_maps: if True then on each iteration yield not only the
product structure but also the relevant old-to-new atom index map
:type yield_renum_maps: bool
:param concentrations: List of concentrations for each source. Must have
the same length as sources
:type concentrations: list of float or None
"""
self.core_st = core_st
self.optimize_sidechains = optimize_sidechains
self.deduplicate = deduplicate
self.start = start
self.stop = stop
self.iters = [source[0] for source in sources]
if all(len(source) == 1 for source in sources):
# Attachment atoms were not specified for any of the R-group
# sources, so we'll try to figure them out from SDF atom properties.
try:
sources = get_sources_from_r_labels(core_st, self.iters)
except ValueError as e:
raise RgroupError(str(e))
if concentrations is None:
# Set default concentrations to 1.
concentrations = [1.] * len(self.iters)
if len(concentrations) != len(self.iters):
raise (RgroupError('Concentrations length must be equal to sources '
'length.'))
self.rgroups = []
self.copy_properties = copy_properties
self.enumerate_cistrans = enumerate_cistrans
self.yield_renum_maps = yield_renum_maps
self.rgroup_sources = []
seen_atoms = set()
natoms = core_st.atom_total
for source_index, (source, concentration) \
in enumerate(zip(sources, concentrations)):
rgroups_indices = []
for atom_index in source[1:]:
if atom_index in seen_atoms:
raise RgroupError(
'Core attachment atom %d specified more than once.' %
(atom_index))
if not 1 <= atom_index <= natoms:
raise RgroupError('Invalid core attachment atom index %d. '
'(Core has %d atoms.)' % (atom_index,
natoms))
leaving_atom = core_st.atom[atom_index]
try:
staying_atom = find_staying_atom(core_st, leaving_atom)
except ValueError as e:
raise RgroupError(str(e))
leaving_atoms = list(
core_st.getMovingAtoms(staying_atom, leaving_atom))
bond_order = core_st.getBond(leaving_atom, staying_atom).order
r = RGroup(atom_index, source_index, leaving_atoms,
staying_atom.index, bond_order)
self.rgroups.append(r)
seen_atoms.add(atom_index)
rgroups_indices.append(len(self.rgroups) - 1)
# Length of the sequence, needed for combinations
# Ensure that it is minimum 1
length = int(concentration * len(rgroups_indices)) or 1
self.rgroup_sources.append(
RGroupSource(source_index, length, rgroups_indices))
if self.optimize_sidechains:
self.f3d_engine = fast3d.SingleConformerEngine()
def __iter__(self):
yield from self._enumerate_and_duduplicate() \
if self.deduplicate else self._enumerate()
def _enumerate_and_duduplicate(self):
seen = set()
for prod in self._enumerate():
code = analyze.generate_smiles(prod[0]
if self.yield_renum_maps else prod)
if code in seen:
continue
else:
seen.add(code)
yield prod
def _enumerate(self):
product = itertools.product(*self.iters)
source_combinations = (itertools.combinations(rs.rgroups_indices,
rs.length)
for rs in self.rgroup_sources)
product_of_combinations = itertools.product(product,
*source_combinations)
combinations = itertools.islice(product_of_combinations, self.start,
self.stop)
for data in combinations:
structs, all_indices = data[0], data[1:]
sidechains = [None] * len(self.rgroups)
for indices, source in zip(all_indices, structs):
for idx in indices:
sidechains[idx] = source
try:
for prod, renum_map in self.attachSidechains(sidechains):
self._addMetadata(prod, tuple(structs))
if self.yield_renum_maps:
yield prod, renum_map
else:
yield prod
except BondOrderMismatch as e:
logger.warning(str(e))
[docs] def attachSidechains(self, sidechains):
"""
Attach the sidechains to the core structure and return the resulting
structure and index map.
:param sidechains: list of sidechains. Should have the same length as
the number of attachment atoms in the core.
:type sidechains: list of `schrodinger.structure.Structure`
:yield: product structures and index maps
:ytype: `schrodinger.structure.Structure`, dict
"""
prod = self.core_st.copy()
atoms_to_delete = []
r_atoms = []
ezs = []
for r, sidechain in zip(self.rgroups, sidechains):
if sidechain is None:
continue
try:
to_del, r_atom, s_atom = self._attachSidechain(
prod, sidechain, r)
atoms_to_delete += to_del
r_atoms.append(r_atom)
if r.bond_order == 2 and self.enumerate_cistrans:
ezs.append((r_atom, s_atom))
except (ValueError, RgroupError) as e:
cls = type(e) if isinstance(e, RgroupError) else RgroupError
raise cls('Error attaching sidechain "%s" from source %d '
'to core atom %d: %s' %
(sidechain.title, r.source_index + 1, r.atom_index,
e.args[0]))
renum_map = prod.deleteAtoms(atoms_to_delete, renumber_map=True)
product_core_atoms = [
renum_map[i]
for i in self.core_st.getAtomIndices()
if renum_map[i] is not None
]
input_core_atoms = [
i for i in self.core_st.getAtomIndices() if renum_map[i] is not None
]
product_ezs = set()
for (i, j) in ezs:
_i, _j = renum_map[i], renum_map[j]
if _i is not None and _j is not None:
product_ezs.add(_order(_i, _j))
if len(input_core_atoms) < 3:
# Not enough core atoms are left for perfect alignment, so we'll
# also try to align the first atom from each R group on top of the
# core attachment atom it replaced.
input_core_atoms += [r.atom_index for r in self.rgroups]
product_core_atoms += [renum_map[i] for i in r_atoms]
for prod_iso in _enumerate_ezs(prod, product_core_atoms, product_ezs):
if self.optimize_sidechains:
self._generateCoordinates(prod_iso, product_core_atoms)
self._alignCore(prod_iso, input_core_atoms, product_core_atoms)
yield prod_iso, renum_map
def _attachSidechain(self, prod, sidechain, r):
"""
Attach a sidechain to the core at the specified atom index. Returns the
indices of the two atoms that should be deleted, as well as the index of
the sidechain atom that is now bonded to the core. (This low-level
method doesn't delete them because the caller may be attaching multiple
sidechains, and deleting each time could trigger confusing renumbering
and be less efficient.)
:param prod: core/product, to be modified in place
:type prod: `schrodinger.structure.Structure`
:param sidechain: sidechain structure (not modified)
:type sidechain: `schrodinger.structure.Structure`
:param r: RGroup object, which specifies the atom to bind to and the
atoms to delete.
:type r: RGroup
:return: list of atoms to delete, R-group atom bound to core,
staying core atom
:rtype: 3-tuple (list of int, int, int)
"""
sidechain = sidechain.copy()
dummy_atom = _find_dummy(sidechain)
r_atom, r_order = _find_neighbor(sidechain.atom[dummy_atom])
try:
core_atom, core_order = _find_neighbor(prod.atom[r.atom_index])
except ValueError:
# The reason we don't always use the staying atom is that it might
# have been deleted in the edge case of a 2-atom core.
core_atom = r.staying_atom
core_order = r.bond_order
if core_order != r_order:
raise BondOrderMismatch('Attachment bond order mismatch.')
rmsd.superimpose_bond(prod, (core_atom, r.atom_index), sidechain,
(dummy_atom, r_atom))
new_dummy_atom = prod.atom_total + dummy_atom
old_atom_total = prod.atom_total
new_r_atom = prod.atom_total + r_atom
prod.extend(sidechain)
# Change sidechain atoms residue name and number to match the ones
# for the 'staying' atom
res_num = prod.atom[r.staying_atom].resnum
res_name = prod.atom[r.staying_atom].pdbres
ins_code = prod.atom[r.staying_atom].inscode
chain = prod.atom[r.staying_atom].chain
for atom_index in range(old_atom_total + 1, prod.atom_total + 1):
prod.atom[atom_index].resnum = res_num
prod.atom[atom_index].pdbres = res_name
prod.atom[atom_index].inscode = ins_code
prod.atom[atom_index].chain = chain
prod.addBond(core_atom, new_r_atom, r_order)
prod.deleteBond(core_atom, r.atom_index)
atoms_to_delete = r.leaving_atoms + [new_dummy_atom]
return atoms_to_delete, new_r_atom, core_atom
def _generateCoordinates(self, st, core_atoms):
"""
Generate 3D coordinates for the non-core atoms in the structure using
fast3d.
:param st: structure, to be modified in place
:type st: `schrodinger.structure.Structure`
:param core_atoms: list of core atoms which will be kept frozen (may be
empty)
:type core_atoms: list of int
"""
core_atom_indices = set(core_atoms)
amide_ez_props = add_amide_constraints(st, core_atom_indices)
core_atom_indices |= get_metals_and_neighbors(st)
try:
self.f3d_engine.run(st, list(core_atom_indices))
except RuntimeError:
pass
for prop in amide_ez_props:
del st.property[prop]
def _alignCore(self, st, input_core_atoms, product_core_atoms):
"""
Superimpose the product structure on the input core.
:param st: structure to rotate in place
:type st: `schrodinger.structure.Structure`
:param input_core_atoms: list of core atoms in self.core_st
:type input_core_atoms: list of int
:param product_core_atoms: list of core atoms in 'st
:type product_core_atoms: list of int
"""
if len(input_core_atoms) > 2:
rmsd.superimpose(
self.core_st,
input_core_atoms,
st,
product_core_atoms,
use_symmetry=False,
move_which=rmsd.CT)
elif len(input_core_atoms) == 2:
rmsd.superimpose_bond(self.core_st, input_core_atoms, st,
product_core_atoms)
else:
# This should never happen; cores with < 2 atoms would have
# failed before this point is reached, but leaving the else
# for completeness.
raise RgroupError("Core must have at least two atoms.")
def _addMetadata(self, prod, sidechains):
prod.property['i_rge_num_r_groups'] = len(sidechains)
# Set concentration
for idx, source in enumerate(self.rgroup_sources, start=1):
concentration = source.length / len(source.rgroups_indices)
prod.property['r_rge_occupation_R%d' % idx] = concentration
for i, st in enumerate(sidechains, 1):
rlabel = 's_rge_R%d' % i
prod.property[rlabel] = st.title
if self.copy_properties:
for prop in st.property:
new_prop = '%s_rge_R%d_%s' % (prop[0], i, prop)
prod.property[new_prop] = st.property[prop]
prod.title = list_to_csv([st.title for st in (prod,) + sidechains])
[docs] def combinations(self):
"""
Return the number of combintations that will be generated for each
tuple of R-groups. That is, combinations due to occupancy of the various
attachment points when all the concentrations are not 1.0.
"""
total = 1
for source in self.rgroup_sources:
total *= scipy.special.comb(
len(source.rgroups_indices), source.length, exact=True)
return total
def _flip_bond(st, core, i, j, k, l):
'''
Flips i-j-k-l dihedral 180 degrees. Assumes that i-j, j-k,
and k-l are bonded, and j-k bond is exocyclic.
:param st: Structure
:type st: `schrodinger.structure.Structure`
:param core: Indices of the "core" atoms (that not to be moved).
:type core: set(int)
:param i: Atom index.
:type i: int
:param j: Atom index.
:type j: int
:param k: Atom index.
:type k: int
:param l: Atom index.
:type l: int
:raises: ValueError, mm.MmException
'''
bs = Bitset(size=st.atom_total)
if k not in core:
fixed, moving = j, k
elif j not in core:
fixed, moving = k, j
else:
raise ValueError('no movable atoms')
mm.mmct_atom_get_moving(st, fixed, st, moving, bs)
angle = st.measure(i, j, k, l)
if angle > 0.0:
angle -= 180.0
else:
angle += 180.0
mm.mmct_atom_set_dihedral_angle(
angle, mm.MM_ANGLE_DEG,
st, i, st, j, st, k, st, l, bs) # yapf: disable
def _enumerate_ezs(st, core, ezs):
'''
Generator that enumerates E/Z isomers for the given bonds.
:param st: Original structure.
:type st: `schrodinger.structure.Structure`
:param core: Indices of the "core" atoms (that are not to be moved).
:type core: set(int)
:param ezs: Set of ordered pairs of atom indices for the
E/Z bonds to enumerate.
:type ezs: set(tuple(int, int))
:yield: Structures.
:ytype: `schrodinger.structure.Structure`
'''
recompute_stereo(st)
_, cistrans_bonds = _get_stereo(st, logger)
todo = [ez for ez in ezs if ez in cistrans_bonds]
for comb in itertools.product((False, True), repeat=len(todo)):
st_iso = st.copy()
for (ez, flip) in zip(todo, comb):
if flip:
atoms = cistrans_bonds[ez][1:]
_flip_bond(st_iso, core, *atoms)
recompute_stereo(st_iso)
yield st_iso
def _find_neighbor(atom):
"""
Return the index of the atom bonded to a given atom
and corresponding bond order.
:param atom: atom
:type atom: `schrodinger.structure._StructureAtom`
:return: atom index and bond order
:rtype: (int, int)
:raises: ValueError if the atom does not have exactly one bond.
"""
if atom.bond_total != 1:
raise ValueError("R-group attachment atom must have exactly one bond.")
bond = atom.bond[1] # 1-based indices
return (bond.atom2.index, bond.order)
def _find_dummy(st):
"""
Return the index of the dummy atom in a structure. Atoms with atomic number
-2 (mmat convention) or 0 (RDKit convention) are considered dummies.
:param st: structure to search
:type st: `schrodinger.structure.Structure`
:return: dummy atom index
:rtype: int
:raises: ValueError if the structure does not contain exactly one dummy atom.
"""
dummies = [at.index for at in st.atom if at.atomic_number in (0, -2)]
if len(dummies) != 1:
raise ValueError(
"R-group structure must have exactly one dummy atom.\n"
"R-group libraries should be prepared using the R-Group Creator "
"panel or the r_group_creator.py command-line script.")
return dummies[0]
def _get_attachment_bond_order(st):
'''
Returns order of the bond that joins "dummy" atom to to the rest.
:param st: R-group structure.
:type st: `schrodinger.structure.Structure`
:return: Attachment bond order.
:rtype: int
:raises: ValueError if the structure does not contain exactly one
"dummy" atom or if the "dummy" atom does not have exactly one bond.
'''
_, order = _find_neighbor(st.atom[_find_dummy(st)])
return order
[docs]def find_rgroup_from_smarts(st,
smarts,
leaving_atom_pos,
staying_atom_pos,
bond_order=None):
"""
Find the various ways in which a structure can be split into "R-group" and
"functional group" using a SMARTS pattern.
The SMARTS pattern must consist of at least two atoms. Two of the atoms,
identified by their position in the SMARTS string, are used to define
the bond to be broken between the R group and the "leaving group".
If the two atoms are not directly connected, the bond leading from
the leaving atom to the staying atom is broken.
For example, consider the structure c1ccccc1cC(=O)O and the SMARTS pattern
*C(=O)O. With leaving_atom_pos=2, staying_atom_pos=1, the entire carboxylate
is removed, producing the R-group c1ccccc1*. With leaving_atom_pos=4,
staying_atom_pos=2, only the terminal O is removed, leading to the R-group
c1ccccc1C(=O)*. (The asterisks are shown here only to highlight the bond
that was broken.)
The return value is a list of tuples, where the first element is the
attachment atom index and the second is a list of the indexes of the atoms
comprising the R-group. In the first example above, if we pretend there are
no hydrogens, the return value might be [(7, [1,2,3,4,5,6])].
Notes: 1) ring bonds can't be broken because they don't split the structure
in two; 2) if `bond_order` is not `None`, skip matches having the attachment
bond of different order.
:param st: structure to analyze
:type st: `schrodinger.structure.Structure`
:param smarts: SMARTS pattern describing the functional group
:type smarts: str
:param leaving_atom_pos: position of the leaving atom in the
SMARTS pattern (1-based)
:type leaving_atom_pos: index
:param staying_atom_pos: position of the attachment atom in the
SMARTS pattern (1-based)
:type staying_atom_pos: index
:param bond_order: If `None` (default), has no effect. Otherwise
skip matches having R-group attachment bond of different order.
:type bond_order: int or NoneType
:return: list of tuples (attachment atom, list of R-group atom indexes). If
no matches satisfied all the requirements, the list may be empty.
May include duplicate R-groups (R-group in this context is the
substructure made of the newly found R-group atoms).
:rtype: list
"""
seen = set()
matches = []
for match in analyze.evaluate_smarts_canvas(st, smarts, uniqueFilter=False):
leaving_atom_idx = match[leaving_atom_pos - 1]
staying_atom_idx = match[staying_atom_pos - 1]
if (leaving_atom_idx, staying_atom_idx) in seen:
continue
else:
seen.add((leaving_atom_idx, staying_atom_idx))
try:
leaving_atoms_idx = st.getMovingAtoms(staying_atom_idx,
leaving_atom_idx)
rgroup_atoms_idx = set(st.getAtomIndices()) - leaving_atoms_idx
except ValueError:
continue # bond is in a ring, so skip
leaving_atom = st.atom[leaving_atom_idx]
bond = next(
b for b in leaving_atom.bond if b.atom2.index in rgroup_atoms_idx)
if bond_order is not None and bond.order != bond_order:
continue # undesired bond order
matches.append((leaving_atom_idx, sorted(rgroup_atoms_idx)))
return matches
[docs]def find_staying_atom(st, leaving_atom):
"""
Given a picked "leaving" atom, determine which of the atoms it is bonded
to is part of the larger molecule - the "staying" atom. All other atoms
bound to the leaving atom are considered to be part of the leaving group.
:param leaving_atom: atom which defines the start of the leaving group
:type leaving_atom: `schrodinger.structure._StructureAtom`
:return: "staying atom": the core atom bound to the leaving atom
:rtype leaving_atom: `schrodinger.structure._StructureAtom`
"""
num_bonds = len(leaving_atom.bond)
if num_bonds == 0:
raise ValueError("Selected atom %i is not bound to any other atom" %
leaving_atom.index)
# For all other errors, there is a standard error text:
err_msg = "Not possible to detect the leaving group from selected atom %i." % leaving_atom.index
moving_atoms_dict = {}
for bond in leaving_atom.bond:
at2 = bond.atom2
try:
moving_atoms_dict[at2] = len(st.getMovingAtoms(at2, leaving_atom))
except ValueError:
# these atoms are in a ring; so whole molecule would move
pass
if len(moving_atoms_dict) == 0:
# All atoms are in rings; freezing either will cause the
# whole molecule to move.
raise ValueError(err_msg)
# Size of the smallest group of bonded atoms:
smallest_moving = min(moving_atoms_dict.values())
if smallest_moving > len(leaving_atom.getMolecule().atom) / 2:
# The wrong "side" of the only non-ring bond was selected. The
# leaving group should never be bigger than 1/2 of the molecule.
raise ValueError(err_msg)
# Ideal atom(s) to freeze:
best_bonded_atoms = [
b for b in moving_atoms_dict if moving_atoms_dict[b] == smallest_moving
]
if len(best_bonded_atoms) > 1:
# 2 (or more) neighbors, if fixed, would yield the same sized
# leaving group.
# e.g. carbon in a methane, or middle carbon in a propane.
raise ValueError(err_msg)
return best_bonded_atoms[0]
[docs]def get_dummy_filter():
"""
Return a filter which has as criteria all the descriptors that can be
computed by this module, along with their suggested default limiters (ranges).
:rtype: `schrodinger.ui.qt.filter_dialog_dir.filter_core.Filter`
"""
criteria = []
for name in sorted(DESCRIPTOR_DICT):
_, limiter = DESCRIPTOR_DICT[name]
criteria.append(filter_core.Criterion(name, limiter))
return filter_core.Filter('dummy filters', criteria)
[docs]def add_descriptors(st, filter_obj):
"""
Add the descriptors required by a filter to a given Structure.
:type st: `schrodinger.structure.Structure`
:type filter_obj: `schrodinger.ui.qt.filter_dialog_dir.filter_core.Filter`
"""
for c in filter_obj.criteria:
try:
func = DESCRIPTOR_DICT[c.prop][0]
except KeyError:
# We'll ignore properties that we don't know how to compute,
# which might be the case for properties that are supposed to be
# already present in the structure.
continue
st.property[c.prop] = func(st)
[docs]def list_to_csv(fields):
"""
Convert a list into a CSV string representation.
:param fields: list to convert
:type fields: list
:rtype: str
"""
dialect = csv.excel()
dialect.lineterminator = ''
buf = io.StringIO()
writer = csv.writer(buf, dialect=dialect)
writer.writerow(fields)
return buf.getvalue()
[docs]def add_amide_constraints(st, frozen_set):
"""
For amide bonds which have one atom frozen and the other not, add the
necessary ct properties to tell fast3d to constrain the amides to the trans
conformation.
:param st: structure, to be modified in place
:type st: `schrodinger.structure.Structure`
:param frozen_set: set of frozen atoms, by atom index
:type frozen_set: set of int
:return: names of properties that were added
:rtype: list of str
"""
amide_bonds = analyze.evaluate_smarts_canvas(
st, r'O=C[NH1X3][#6]', uniqueFilter=True)
index = get_last_EZ_property_index(st)
props = []
for i, j, k, l in amide_bonds:
if len({j, k} & frozen_set) == 1:
index += 1
prop = 's_st_EZ_%d' % index
props.append(prop)
# Making the O=CNC dihedral Z makes it a "trans amide".
st.property[prop] = '%d_%d_%d_%d_Z' % (i, j, k, l)
return props
[docs]def get_last_EZ_property_index(st):
"""
Return the maximum index of the s_st_EZ_<index> properties of st.
If there are no such properties, return 0.
:type st: `schrodinger.structure.Structure`
:rtype: int
"""
prefix = 's_st_EZ_'
indexes = [
int(p.split(prefix)[1]) for p in st.property if p.startswith(prefix)
]
return max(indexes) if indexes else 0
[docs]def get_sources_from_r_labels(st, iters, prop='i_rdkit__MolFileRLabel'):
"""
Given a Structure and a list of iterables, return the "sources" data
structure needed by RgroupEnumerator. The structure must have (some) atoms
with the specified property; the values of this property must be in the
range [1, len(iters)] and all the values in that range must be represented
at least once. If this condition is not met, raise a ValueError.
:param st: core Structure
:type st: schrodinger.structure.Structure
:param iters: list of iterables of structures
:type iters: list
:param prop: name of the atom property holding the R-group labels. The
default is what comes from reading an SD file with "M RGP" fields using
RDKit.
:type prop: str
:return: sources data structure. See RgroupEnumerator for details.
:rtype: list of list
"""
n = len(iters)
source_dict = collections.defaultdict(list)
for atom in st.atom:
r_index = atom.property.get(prop)
if r_index is not None:
source_dict[r_index].append(atom.index)
expected_keys = set(range(1, n + 1))
got_keys = set(source_dict.keys())
if got_keys != expected_keys:
raise ValueError('Structure does not have the necessary R-group '
'atom properties; expected {}, got {}'.format(
expected_keys, got_keys))
return [[iters[i]] + source_dict[i + 1] for i in range(n)]
[docs]def convert_attachment_point(struct, bond_order=None):
"""
Converts attachment point from methyl to dummy atom. If r-group fragment
is 'Null' returns False, otherwise returns True. Null r-group has atom
with atomic number -2 and growname 'rpc1'.
:param struct: structure object
:type struct: `structure.Structure`
:param bond_order: If `None`, has no effect. Otherwise
return False if attachment bond is of different order.
:type bond_order: int or NoneType
:return: True if conversion succeeded and False otherwise.
:rtype: bool
"""
# Check whether r-group fragment is a 'Null' group or has more than
# one attachment point. The later case is identified by the
# presence of s_cg_attachment property.
for atom in struct.atom:
if atom.growname == 'rpc1' and \
atom.atomic_number == DUMMY_ATOMIC_NUMBER:
return False
if 's_cg_attachment' in atom.property:
return False
for atom in struct.atom:
if atom.growname == "rpc2":
delete_atoms = [
a for a in atom.bonded_atoms if not a.growname == "rpc1"
]
struct.deleteAtoms(delete_atoms)
atom.atomic_number = DUMMY_ATOMIC_NUMBER
try:
_, order = _find_neighbor(atom)
except ValueError:
return False
if bond_order is not None and order != bond_order:
return False
return True
# No conversion was needed. Here we need to check that fragment contains
# attachment point (a dummy atom).
return check_attachment_point(struct, bond_order)
[docs]def get_attachment_point(st):
"""
Identifies attachment point (dummy atom with a single neighbor)
in the provided structure.
:param st: Structure
:type st: `schrodinger.structure.Structure`
:return: Index of the dummy atom and order of the bond that
joins the dummy to the rest.
:rtype: (int, int)
"""
try:
dummy = _find_dummy(st)
neighbor, order = _find_neighbor(st.atom[dummy])
return dummy, order
except ValueError:
return (None, None)
[docs]def check_attachment_point(struct, bond_order=None):
"""
Checks that provided structure contains attachment point.
:param bond_order: If `None`, has no effect. Otherwise
return False if attachment bond is of different order.
:type bond_order: int or NoneType
:return: True if attachment point was found and False otherwise.
:rtype: bool
"""
dummy, order = get_attachment_point(struct)
if dummy is not None:
if bond_order is not None and order != bond_order:
return False
else:
return True
return False
[docs]def create_fragment_structure(st, rgroup_data):
"""
Creates r-group fragment structures from a given structure and a list
of atoms that should be included in the r-group.
:param st: structure object
:type st: `structure.Structure`
:param rgroups: r-group data that contains index of attachment atom and
indices of R-group fragment atoms.
:type rgroups: tuple
:return: fragment structure
:rtype: `structure.Structure`
"""
attach_atom, rgroup_indices = rgroup_data
frag_st = st.copy()
frag_st.atom[attach_atom].atomic_number = DUMMY_ATOMIC_NUMBER
# r-group structure should include an attachment point, so we
# are adding it here.
rgroup_indices.append(attach_atom)
frag_st = frag_st.extract(rgroup_indices, copy_props=True)
return frag_st