from collections import defaultdict
import networkx as nx
import schrodinger.structure as structure
import schrodinger.structutils.analyze as analyze
import schrodinger.utils.log as log
from schrodinger.application.desmond import constants
N_HATOM = 1
MAX_NUM_REST_ATOM = 25
MIN_FUSED_RING = 7
LARGE_FRAG_ON_RING = 3
logger = log.get_output_logger(name="r_group_asl")
[docs]def countHeavyAtoms(st):
"""
This function counts the number of heavy atoms in the structure.
:type st: `structure.Structure`
:param st: Structure
:rtype: `int`
:return: number of heavy atoms
"""
n = 0
for a in st.atom:
if a.element != "H":
n += 1
return n
[docs]def getImportantRotatableBonds(st, atom_cutoff=6):
"""
This function returns important rotatable bonds. Important rotatable bonds
are defined when both sides of the bond have either rings or heavy atoms
greater than atom_cutoff.
:type st: `structure.Structure`
:param st: Structure
:type atom_cutoff: `int`
:param atom_cutoff: # of heavy atoms
:return: generator for atom pairs in important rotatable bonds
"""
for (i, j) in analyze.rotatable_bonds_iterator(st):
disconnected_st = st.copy()
disconnected_st.deleteBond(i, j)
# probably macro cycle
if disconnected_st.mol_total == 1:
yield (i, j)
mol_i = disconnected_st.atom[i].getMolecule().extractStructure()
mol_j = disconnected_st.atom[j].getMolecule().extractStructure()
nheavy_i = countHeavyAtoms(mol_i)
nheavy_j = countHeavyAtoms(mol_j)
nring_i = len(mol_i.ring)
nring_j = len(mol_j.ring)
if ((nheavy_i > atom_cutoff or nring_i > 0) and
(nheavy_j > atom_cutoff or nring_j > 0)):
yield (i, j)
[docs]def get_fep_cts(cts):
"""
Find reactant and product ct from the input, and return touple of reference ligand and mutated ligand structures
:param cts: input structures
:return: ref_ct, mut_ct
"""
ref_ct = []
mut_ct = []
for ct in cts:
if "s_fep_fragname" in ct.property:
if ct.property["s_fep_fragname"] == "none":
ref_ct.append(ct)
else:
mut_ct.append(ct)
return ref_ct, mut_ct
[docs]def count_heavy_atoms(ct, atom_index_list):
n = 0
for index in atom_index_list:
if ct.atom[index].element != "H":
n += 1
return n
[docs]def get_heavy_atoms(ct, atom_index_list):
ret = []
for index in atom_index_list:
if ct.atom[index].element != "H":
ret.append(index)
return ret
[docs]def get_fragment(b0, b1, ct):
"""
:param b0: first atom index for bond
:param b1: second atom index for bond, this connects to the fragment returned
:return: list of atom indice for all atom connected to b1 when bond b0-b1 is removed
"""
fragments = ct.copy()
for a in fragments.atom:
a.property["i_org_index"] = a.index
if b0 and b1:
fragments.deleteBond(b0, b1)
if fragments.mol_total == 1:
return []
mol_b1 = fragments.atom[b1].getMolecule().extractStructure()
return [a.property["i_org_index"] for a in mol_b1.atom]
def _get_heavy_atom_fragments(ct, atom1, atom2):
"""
:param ct: input structure
:type ct: structure.Structure
:param atom1: first atom index in a bond
:type atom1: int
:param atom2: second atom index in a bond
:type atom2: int
:return: lists of atom indices of connected groups of heavy atoms after bond,
(atom1, atom2), is removed from ct. The element lists are sorted
in ascending number of atoms.
:rtype: List[List[Int]]
"""
atoms = {a for a in ct.atom if not a.element == "H"}
mol_graph = analyze.create_nx_graph(ct, atoms=atoms)
mol_graph.remove_edge(atom1, atom2)
return sorted(nx.connected_components(mol_graph), key=len)
[docs]def get_attachment_bonds(source_ct, dest_ct):
"""
:param source_ct: reference structure
:type source_ct: `structure.Structure`
:param dest_ct: mutant structure
:type dest_ct: `structure.Structure`
:return: [(int, int), (int, int)] list of attachment bond pairs:
((source_anchor, source_bridge), (dest_anchor, dest_bridge))
source_anchor and dest_anchor are matched pair of core atom numbers;
source_bridge and dest_bridge are atoms connect to the anchor
atom numbers (these could be None).
The algorithm works by looking at each matched pair of atoms (core atoms)
and track chemical changes. If there is a dummy atom attached to the mutant core atom,
an attachment point is generated. If the chemical change is connected to more than
one core atoms, ring or linker atom mutation is detected for attachment point
generation after the dummy atom mutations. If any of the fragment to which
linker atom connects has less than 6 heavy atoms, treat this fragment
as functional group.
In an attachment bond, the second atom connects to a mutating fragment.
If there is a dummy atom, it is always the second atom.
For bonds attached to ring/linker atom mutations, both atoms in one structure
may map to atoms in the other. In case one atom is a core atom and the other
is a mutating linker atom, the ring/linker atom must be the second atom.
This way the calling function can use the second atom to decide which fragment
is mutating across the attachment bond.
"""
d2s_map = {
a: source_ct.atom[a.property[constants.FEP_MAPPING]]
for a in dest_ct.atom
if a.property[constants.FEP_MAPPING]
}
source_core = set(d2s_map.values())
dest_core = set(d2s_map)
source_dummy = set(source_ct.atom) - source_core
dest_dummy = set(dest_ct.atom) - dest_core
ring_list = [r.getAtomIndices() for r in dest_ct.ring]
linker = set()
d_anchor = []
d_anchor_dummy = defaultdict(list)
d_anchor_mapped = defaultdict(set)
linker_as_dummy = set()
def atom_index(atom):
return atom.index
for d in sorted(d2s_map, key=atom_index):
s = d2s_map[d]
dest_bonded_atoms = frozenset(d.bonded_atoms)
source_bonded_atoms = frozenset(s.bonded_atoms)
s_linked_core = source_bonded_atoms.intersection(source_core)
d_linked_core = dest_bonded_atoms.intersection(dest_core)
# linker/ring atoms first, treat small linker atom as functional group
# mutation
if len(s_linked_core) > 1 and len(d_linked_core) > 1:
if (len(dest_bonded_atoms) !=
len(source_bonded_atoms)) or (s.element != d.element):
for a in sorted(d_linked_core, key=atom_index):
if not any(
{d.index, a.index}.issubset(r)
for r in ring_list) and count_heavy_atoms(
dest_ct, get_fragment(d.index, a.index, dest_ct)
) + count_heavy_atoms(
source_ct,
get_fragment(s.index, d2s_map[a].index,
source_ct)) < 12:
d_anchor_dummy[d].append(a)
linker_as_dummy.add((d, a))
else:
linker.add((s, d))
for d_bridge in sorted(dest_bonded_atoms.intersection(dest_dummy),
key=atom_index):
d_anchor_dummy[d].append(d_bridge)
if d not in d_anchor:
d_anchor.append(d)
if source_bonded_atoms.intersection(source_dummy):
if d not in d_anchor:
d_anchor.append(d)
for l_s, l_d in sorted(linker, key=lambda atom_pair: atom_pair[0].index):
for d_bridge in sorted(set(l_d.bonded_atoms).intersection(dest_core),
key=atom_index):
if (l_d, d_bridge) not in linker_as_dummy:
d_anchor_mapped[l_d].add(d_bridge)
# add linker/ring atoms back so that dummy atoms get picked first
if l_d not in d_anchor:
d_anchor.append(l_d)
ret_val = []
s_linker = set()
if linker:
s_linker, _ = zip(*linker)
for b in d_anchor:
s_bridge = 0
s_anchor = d2s_map[b]
for s_b in sorted(set(s_anchor.bonded_atoms).intersection(source_dummy),
key=atom_index):
# pair up core-dummy, if there are more core-dummy in any of the molecules
# use None for the dummy in the other molecule
if s_bridge < len(d_anchor_dummy[b]):
ret_val.append(((s_anchor.index, s_b.index),
(b.index, d_anchor_dummy[b][s_bridge].index)))
s_bridge += 1
else:
ret_val.append(((s_anchor.index, s_b.index), (b.index, None)))
dest_bridge_atoms = d_anchor_dummy[b]
for d_bridge_atom in sorted(dest_bridge_atoms[s_bridge:],
key=atom_index):
# unpaired dummy dest bridge atoms don't have any corresponding
# source bridge dummy but there are matched linker dest bridge atoms
# which have matched source bridge atoms. Need to look up the
# source bridge atoms in atom matches.
s_bridge_index = None
if (b, d_bridge_atom) in linker_as_dummy:
s_bridge_index = d2s_map[d_bridge_atom].index
ret_val.append(((s_anchor.index, s_bridge_index),
(b.index, d_bridge_atom.index)))
for d_mapped in sorted(d_anchor_mapped[b], key=atom_index):
# need to make attachment bond a (core atom, mutating atom)
# pair
s_mapped = d2s_map[d_mapped]
if s_mapped not in s_linker:
# s_mapped is not changing, i.e. core atom
# Checking d_mapped and b yields the same result
att_pair = ((s_mapped.index, s_anchor.index), (d_mapped.index,
b.index))
else:
# This is the case where both atoms in the attachment bond are ring/linker atoms.
# Use the number of heavy atoms in the connected fragments connected to determine the direction,
# i.e. make the second atom always in the smaller fragment (in terms of the number of heavy atoms).
# If the attachment bond connects the same sized fragments in both of the molecules,
# print a warning.
source_frags = _get_heavy_atom_fragments(
source_ct, s_anchor.index, s_mapped.index)
if len(source_frags) == 1:
att_pair = ((s_mapped.index, s_anchor.index),
(d_mapped.index, b.index))
if len(
_get_heavy_atom_fragments(dest_ct, b.index,
d_mapped.index)) > 1:
logger.warning(
"linker atom mapped to ring atoms, core-hopping workflow should be used"
)
elif len(source_frags[0]) == len(source_frags[1]):
dest_frags = _get_heavy_atom_fragments(
dest_ct, b.index, d_mapped.index)
if len(dest_frags) == 1:
logger.warning(
"linker atom mapped to ring atoms, core-hopping workflow should be used"
)
att_pair = ((s_mapped.index, s_anchor.index),
(d_mapped.index, b.index))
elif len(dest_frags[0]) == len(dest_frags[1]):
logger.warning(
f"linker ({s_anchor.index}, {s_mapped.index}) connects two fragments of the same number of heavy atoms"
)
att_pair = ((s_mapped.index, s_anchor.index),
(d_mapped.index, b.index))
else:
if b.index in dest_frags[0]:
att_pair = ((s_mapped.index, s_anchor.index),
(d_mapped.index, b.index))
else:
att_pair = ((s_anchor.index,
d2s_map[d_mapped].index),
(b.index, d_mapped.index))
else:
if s_anchor.index in source_frags[0]:
att_pair = ((s_mapped.index, s_anchor.index),
(d_mapped.index, b.index))
else:
att_pair = ((s_anchor.index, s_mapped.index),
(b.index, d_mapped.index))
ret_val.append(att_pair)
return ret_val
[docs]def get_big_fused_rings(ct):
""" return heavy atoms of fused rings from SSSR set
"""
visited = set([])
all_rings = []
ring_list = [r.getAtomIndices() for r in ct.ring]
for (i, r) in enumerate(ring_list):
if i not in visited:
visited.add(i)
fusable = False
for r0 in all_rings:
if r0.intersection(r):
r0.update(r)
fusable = True
if not fusable:
all_rings.append(set(r))
remove_list = []
for r in all_rings:
if count_heavy_atoms(ct, set(r)) < MIN_FUSED_RING:
remove_list.append(r)
for r in remove_list:
all_rings.remove(r)
return all_rings
[docs]def combine_rings(rings):
""" combine fused rings from SSSR set
:return: list of lists of atoms
"""
visited = set([])
all_rings = []
for (i, r) in enumerate(rings):
if i not in visited:
visited.add(i)
fusable = False
for r0 in all_rings:
if r0.intersection(r):
r0.update(r)
fusable = True
if not fusable:
all_rings.append(set(r))
return all_rings
[docs]def is_ring_atom(ct, atom):
for r in ct.ring:
if atom in r.getAtomIndices():
return True
return False
[docs]def whole_ring_test(ct, atom_list):
"""
test if atom list contains entire rings
"""
alist = set(atom_list)
ret_val = False
for r in ct.ring:
ring_atoms = r.getAtomIndices()
if alist.issuperset(set(ring_atoms)):
ret_val = True
return ret_val
[docs]def find_closest_rotbond(graph, anchor_index, rotable_bonds):
"""
given anchor atom and connection graph of a molecule,
find the rotable bonds that is cloest (in terms of number of bonds to the anchor atom)
:param graph: molecular graph
:type graph: networkx graph
:param anchor_index: anchor atom index
:param rotable_bonds: list of rotable bonds
:type rotable_bonds: list of touples
"""
distance = 100
for bond in rotable_bonds:
#find closest rotable bond
path0 = nx.shortest_path(graph, anchor_index, bond[0])
path1 = nx.shortest_path(graph, anchor_index, bond[1])
if len(path1) < len(path0):
if distance > len(path1):
distance = len(path1)
rot_atom1 = bond[1]
rot_atom0 = bond[0]
#save the bond in original atom order
save_rot_bond = bond
else:
if distance > len(path0):
distance = len(path0)
rot_atom1 = bond[0]
rot_atom0 = bond[1]
#save the bond in original atom order
save_rot_bond = bond
return (distance, (rot_atom0, rot_atom1), save_rot_bond)
[docs]def find_frag_from_closest_rotbond(ct, anchor_index, rotable_bonds):
"""
find rotable bond that is closest to the anchor atom in terms of number of
bonds return the fragment, original rotable bond, bond used to find the
fragment, and list of rotable bonds from which large fragment is cut off
(extra_rot). The last two items are used to find the matching fragment in
the other molecule
:return: (hot_r, save_rot_bond, (rot_atom0, rot_atom1), extra_rot)
where:
- hot_r: hot region found
- save_rot_bond: original rotable bond used to find the fragment
- (rot_atom0, rot_atom1): rotable bond in the same atom order
passed to get_fragment
- extra_rot: list of extra_rotable bond directly connected to the
fragment; the fragments attached by these bonds are removed
already
"""
fuse_rings = get_big_fused_rings(ct)
graph = analyze.create_nx_graph(ct)
distance = 100
rot_atom1 = 0
rot_atom0 = 0
save_rot_bond = (0, 0)
if anchor_index:
(distance, (rot_atom0, rot_atom1),
save_rot_bond) = find_closest_rotbond(graph, anchor_index,
rotable_bonds)
if distance == 100:
return ([], (0, 0), (0, 0), [])
hot_r = get_fragment(rot_atom0, rot_atom1, ct)
extra_rot = remove_extra_rotable_atoms(ct, graph, anchor_index, hot_r,
rotable_bonds)
for r in fuse_rings:
if set(hot_r).intersection(r):
hot_r = set()
break
return (hot_r, save_rot_bond, (rot_atom0, rot_atom1), extra_rot)
[docs]def get_fragments_from_att_points(source_ct, dest_ct, source_att, dest_att,
source_rotable_bonds, dest_rotable_bonds):
""" given attachment pairs, find the rotable fragments containing the attachment bond
return fragments and the rotable bond associated
The bigger one of the fragment pair is used to determine which rotabond to use.
:param source_ct: source structure
:param dest_ct: dest structure
:param source_att: source attachment point
:type source_att: touple of integers
:param dest_att: dest attachment point
:type dest_att: touple of integers
:param source_rotable_bonds: source rotable bond list
:param dest_rotable_bonds: dest rotable bond list
:return: (source_hot_frag, dest_hot_frag, (s_rot0, s_rot1), (d_rot0, d_rot1)). The order of atoms in the rotable bonds returned is in accord to the way the fragments are generated.
"""
dest_hot_frag = []
source_hot_frag = []
(s_rot0, s_rot1) = (0, 0)
(d_rot0, d_rot1) = (0, 0)
s_frag = get_fragment(source_att[0], source_att[1], source_ct)
d_frag = get_fragment(dest_att[0], dest_att[1], dest_ct)
if len(s_frag) > len(d_frag) and source_rotable_bonds:
(source_hot_frag, (s_rot0, s_rot1), bond,
extra_rot) = find_frag_from_closest_rotbond(source_ct, source_att[0],
source_rotable_bonds)
if len(source_hot_frag):
dest_hot_frag = get_fragment(
source_ct.atom[bond[0]].property[constants.FEP_MAPPING],
source_ct.atom[bond[1]].property[constants.FEP_MAPPING],
dest_ct)
for (atom0, atom1) in extra_rot:
extra_frag = get_fragment(
source_ct.atom[atom0].property[constants.FEP_MAPPING],
source_ct.atom[atom1].property[constants.FEP_MAPPING],
dest_ct)
dest_hot_frag = list(set(dest_hot_frag) - set(extra_frag))
elif dest_rotable_bonds:
(dest_hot_frag, (d_rot0, d_rot1), bond,
extra_rot) = find_frag_from_closest_rotbond(dest_ct, dest_att[0],
dest_rotable_bonds)
if len(dest_hot_frag):
source_hot_frag = get_fragment(
dest_ct.atom[bond[0]].property[constants.FEP_MAPPING],
dest_ct.atom[bond[1]].property[constants.FEP_MAPPING],
source_ct)
for (atom0, atom1) in extra_rot:
extra_frag = get_fragment(
dest_ct.atom[atom0].property[constants.FEP_MAPPING],
dest_ct.atom[atom1].property[constants.FEP_MAPPING],
source_ct)
source_hot_frag = list(set(source_hot_frag) - set(extra_frag))
if (s_rot0, s_rot1) == (0, 0) and (d_rot0, d_rot1) != (0, 0):
if (d_rot0 == dest_att[0] and
d_rot1 == dest_att[1]) or (d_rot0 == dest_att[1] and
d_rot1 == dest_att[0]):
(s_rot0, s_rot1) = source_att
else:
s_rot0 = dest_ct.atom[d_rot0].property[constants.FEP_MAPPING]
s_rot1 = dest_ct.atom[d_rot1].property[constants.FEP_MAPPING]
elif (d_rot0, d_rot1) == (0, 0) and (s_rot0, s_rot1) != (0, 0):
if (s_rot0 == source_att[0] and
s_rot1 == source_att[1]) or (s_rot0 == source_att[1] and
s_rot1 == source_att[0]):
(d_rot0, d_rot1) = dest_att
else:
d_rot0 = source_ct.atom[s_rot0].property[constants.FEP_MAPPING]
d_rot1 = source_ct.atom[s_rot1].property[constants.FEP_MAPPING]
return (source_hot_frag, dest_hot_frag, (s_rot0, s_rot1), (d_rot0, d_rot1))
[docs]def is_rotabond(bond, rota_bond_list):
for b in rota_bond_list:
if set(b) == set(bond):
return True
return False
[docs]def apply_exclusion(source_hot_frag, source_excl, dest_hot_frag, dest_excl):
s_excl = set(source_hot_frag) & set(source_excl)
for a in s_excl:
source_hot_frag.remove(a)
d_excl = set(dest_hot_frag) & set(dest_excl)
for a in d_excl:
dest_hot_frag.remove(a)
[docs]def extend_hotregion(total_heavy, source_ct, dest_ct, att_bonds, \
hotregion_source, hotregion_dest, \
source_rot_bonds, dest_rot_bonds, \
source_frag_rot_bonds_visited, dest_frag_rot_bonds_visited,
max_rest_atoms,
source_excl, dest_excl):
"""
extend the hot region by moving to the next rotable bond
:param total_heavy: total heavy atoms in the hot region so far
:param source_ct: source structure
:param dest_ct: dest structure
:param att_bonds: all attachment bonds for mutations
:param source_rot_bonds: rotable bond list of reactant
:param dest_rot_bonds: rotable bond list of product
:param source_frag_rot_bonds_visite: dictionary of rotable bonds already visited for attchment bonds in reactant
:param dest_frag_rot_bonds_visite: dictionary of rotable bonds already visited for attchment bonds in product
"""
while total_heavy < max_rest_atoms:
total_rotbond_left = 0
new_frag = 0
# old_total = total_heavy
for (source_att, dest_att) in att_bonds:
#try to increase the rest region by moving to the next rotable bond for each fragment
source_rot_bond_for_frag = [
bond for bond in source_rot_bonds
if bond not in source_frag_rot_bonds_visited[source_att]
]
dest_rot_bond_for_frag = [
bond for bond in dest_rot_bonds
if bond not in dest_frag_rot_bonds_visited[dest_att]
]
total_rotbond_left += len(source_rot_bond_for_frag) + len(
dest_rot_bond_for_frag)
if not source_rot_bond_for_frag and not dest_rot_bond_for_frag:
continue
(source_hot_frag, dest_hot_frag, source_rot_bond, dest_rot_bond) = \
get_fragments_from_att_points(source_ct, dest_ct,
source_att, dest_att,
source_rot_bond_for_frag,
dest_rot_bond_for_frag)
apply_exclusion(source_hot_frag, source_excl, dest_hot_frag,
dest_excl)
new_frag += len(source_hot_frag) + len(dest_hot_frag)
if source_rot_bond != (0, 0):
source_frag_rot_bonds_visited[source_att].append(
source_rot_bond)
if dest_rot_bond != (0, 0):
dest_frag_rot_bonds_visited[dest_att].append(dest_rot_bond)
tmp_source_hot = hotregion_source.copy()
tmp_source_hot.update(source_hot_frag)
tmp_dest_hot = hotregion_dest.copy()
tmp_dest_hot.update(dest_hot_frag)
if (count_heavy_atoms(source_ct, tmp_source_hot) +
count_heavy_atoms(dest_ct, tmp_dest_hot)) < max_rest_atoms:
hotregion_source.update(source_hot_frag)
hotregion_dest.update(dest_hot_frag)
total_heavy = count_heavy_atoms(
source_ct, hotregion_source) + count_heavy_atoms(
dest_ct, hotregion_dest)
elif whole_ring_test(source_ct, source_hot_frag) and \
whole_ring_test(dest_ct, dest_hot_frag) and \
(count_heavy_atoms(source_ct, tmp_source_hot) + \
count_heavy_atoms(dest_ct, tmp_dest_hot)) < max_rest_atoms + 6:
hotregion_source.update(source_hot_frag)
hotregion_dest.update(dest_hot_frag)
total_heavy = count_heavy_atoms(
source_ct, hotregion_source) + count_heavy_atoms(
dest_ct, hotregion_dest)
else:
continue
if total_rotbond_left == 0 or new_frag == 0:
break
[docs]def get_rotable_bonds(source_ct, dest_ct, min_h_atom):
"""
Find rotable bonds paired for source and dest cts
return paired rotable bonds in core, and all possible rotable bonds for each ct
"""
source_rot_bonds = [
x for x in getImportantRotatableBonds(source_ct, atom_cutoff=min_h_atom)
]
dest_rot_bonds = [
x for x in getImportantRotatableBonds(dest_ct, atom_cutoff=min_h_atom)
]
all_s_rot_bonds = [
x for x in getImportantRotatableBonds(source_ct, atom_cutoff=min_h_atom)
]
all_d_rot_bonds = [
x for x in getImportantRotatableBonds(dest_ct, atom_cutoff=min_h_atom)
]
#remove the rotable bonds consisting of dummpy atoms
source_rot_remove = []
for bond in source_rot_bonds:
if source_ct.atom[bond[0]].property[
constants.FEP_MAPPING] == 0 or source_ct.atom[bond[1]].property[
constants.FEP_MAPPING] == 0:
source_rot_remove.append(bond)
for bond in source_rot_remove:
source_rot_bonds.remove(bond)
dest_rot_remove = []
for bond in dest_rot_bonds:
if dest_ct.atom[bond[0]].property[
constants.FEP_MAPPING] == 0 or dest_ct.atom[bond[1]].property[
constants.FEP_MAPPING] == 0:
dest_rot_remove.append(bond)
for bond in dest_rot_remove:
dest_rot_bonds.remove(bond)
if len(source_rot_bonds) < len(dest_rot_bonds):
source_rot_bonds = [(dest_ct.atom[x[0]].property[constants.FEP_MAPPING],
dest_ct.atom[x[1]].property[constants.FEP_MAPPING])
for x in dest_rot_bonds]
elif len(source_rot_bonds) > len(dest_rot_bonds):
dest_rot_bonds = [(source_ct.atom[x[0]].property[constants.FEP_MAPPING],
source_ct.atom[x[1]].property[constants.FEP_MAPPING])
for x in source_rot_bonds]
return (source_rot_bonds, dest_rot_bonds, all_s_rot_bonds, all_d_rot_bonds)
[docs]def set_hotregion(cts, \
max_rest_atoms=MAX_NUM_REST_ATOM, \
min_h_atom=N_HATOM, \
exclusion_lists=(set(), set()), \
get_fragment=get_fragment, \
get_attachment_bonds=get_attachment_bonds):
"""
:param cts: input structures, s_fep_fragname should be defined for sourec
and dest cts
:type cts: Lists of structure
:param max_rest_atoms: maximum allowed number of heavy atoms in hot region,
soft limit
:param min_h_atom: minimum heavy atom moved by a rotable bond
:param exclusion_lists: touple of sets of core atoms (hydrogen included)
:param get_fragment: callable to compute fragment given (atom_index0,
atom_index1, ct)
:param get_attachment_bonds: callable to compute attathment points given
two cts with i_fep_mapping defined (ct0, ct1)
:return: number of heavy atoms assigned
Set the hotregions in the cts according to mutated fragments and maximum
allowed number of atoms. Core atoms can be exclude via exclusion_lists.
Steps to set REST region:
1. Include all the changed atoms (?99, ?97, ?98). If there are too many
atoms, stop. Linker atoms are ignored (linker connects two fragments
with more than 6 heavy atoms in each of the fragments)
2. If a rotable bond is also an attachment point, include all atoms in
the fragment.
3. For each changed fragment pair (attachment point), try to include all
the atoms up to the closest rotatble bond. How far away the rotable bond
is is measure by the minimum number of bonds separating the fragment and
the rotable bonds. Fused rings are excluded. If the fragment is a ring,
the maximum allowed atoms may be six more than the allowed number.
4. Repeat 3, if maximum number of atoms are not reached and it is still
possible to add more atoms.
"""
ref_ct, mut_ct = get_fep_cts(cts)
source_ct = ref_ct[0]
dest_ct = mut_ct[0]
(source_excl, dest_excl) = exclusion_lists
for s_atom in source_ct.atom:
s_atom.property[constants.FEP_MAPPING] = 0
for d_atom in dest_ct.atom:
s_mapped = d_atom.property[constants.FEP_MAPPING]
if s_mapped:
source_ct.atom[s_mapped].property[
constants.FEP_MAPPING] = d_atom.index
att_bonds = get_attachment_bonds(source_ct, dest_ct)
hotregion_source = set()
hotregion_dest = set()
total_heavy = 0
(source_rot_bonds, dest_rot_bonds, all_s_rot_bonds, all_d_rot_bonds) = \
get_rotable_bonds(source_ct, dest_ct, min_h_atom)
source_frag_rot_bonds_visited = {}
dest_frag_rot_bonds_visited = {}
for (source_att, dest_att) in att_bonds:
source_frag_rot_bonds_visited[source_att] = []
dest_frag_rot_bonds_visited[dest_att] = []
old_total_robond_left = 0
#try to include all the mutated fragment pairs once
for (source_att, dest_att) in att_bonds:
# if the attachment bonds is also a rotable bond, add the fragment
tmp_source_fragment = set()
tmp_dest_fragment = set()
if is_rotabond(source_att, all_s_rot_bonds):
tmp_source_fragment = get_fragment(source_att[0], source_att[1],
source_ct)
tmp_dest_fragment = get_fragment(dest_att[0], dest_att[1], dest_ct)
elif is_rotabond(dest_att, all_d_rot_bonds):
tmp_source_fragment = get_fragment(source_att[0], source_att[1],
source_ct)
tmp_dest_fragment = get_fragment(dest_att[0], dest_att[1], dest_ct)
elif is_ring_atom(source_ct, source_att[0]) or is_ring_atom(source_ct, source_att[1]) or\
is_ring_atom(dest_ct, dest_att[0]) or is_ring_atom(dest_ct, dest_att[0]):
#try to include the whole ring when attachment is part of the ring or
#directly attached to a ring
(source_hot_frag, dest_hot_frag, source_rot_bond, dest_rot_bond) = \
get_fragments_from_att_points(source_ct, dest_ct,
source_att, dest_att,
source_rot_bonds,
dest_rot_bonds)
apply_exclusion(source_hot_frag, source_excl, dest_hot_frag,
dest_excl)
(tmp_source_fragment, tmp_dest_fragment) = (source_hot_frag,
dest_hot_frag)
if source_rot_bond != (0, 0):
source_frag_rot_bonds_visited[source_att].append(
source_rot_bond)
if dest_rot_bond != (0, 0):
dest_frag_rot_bonds_visited[dest_att].append(dest_rot_bond)
if total_heavy + (
count_heavy_atoms(source_ct, tmp_source_fragment) +
count_heavy_atoms(dest_ct, tmp_dest_fragment)) < max_rest_atoms:
hotregion_source.update(tmp_source_fragment)
hotregion_dest.update(tmp_dest_fragment)
total_heavy = count_heavy_atoms(
source_ct, hotregion_source) + count_heavy_atoms(
dest_ct, hotregion_dest)
elif whole_ring_test(source_ct, tmp_source_fragment) and \
whole_ring_test(dest_ct, tmp_dest_fragment) and \
total_heavy + (count_heavy_atoms(source_ct, tmp_source_fragment) + \
count_heavy_atoms(dest_ct, tmp_dest_fragment)) < max_rest_atoms + 6:
hotregion_source.update(tmp_source_fragment)
hotregion_dest.update(tmp_dest_fragment)
total_heavy = count_heavy_atoms(
source_ct, hotregion_source) + count_heavy_atoms(
dest_ct, hotregion_dest)
else:
logger.debug("trying to assign more than %d atoms for purterbation %s->%s " % \
(max_rest_atoms, source_ct.property['s_m_title'], dest_ct.property['s_m_title']))
extend_hotregion(total_heavy, source_ct, dest_ct, att_bonds, \
hotregion_source, hotregion_dest, \
source_rot_bonds, dest_rot_bonds, \
source_frag_rot_bonds_visited, dest_frag_rot_bonds_visited, \
max_rest_atoms, source_excl, dest_excl)
for a in source_ct.atom:
a.property[constants.REST_HOTREGION] = int(int(a) in hotregion_source)
for a in dest_ct.atom:
a.property[constants.REST_HOTREGION] = int(int(a) in hotregion_dest)
return total_heavy
if ("__main__" == __name__):
import sys
if (len(sys.argv) == 1):
print(
"Usage: $SCHRODINGER/run %s <input-mae-file> [<output-mae-file>]" %
sys.argv[0])
sys.exit(0)
mae_fname = sys.argv[1]
cts = []
for ct in structure.StructureReader(mae_fname):
if len(ct.atom) > 500:
continue
ct.deletePropertyFromAllAtoms(constants.REST_HOTREGION)
ct.deletePropertyFromAllAtoms('i_user_hotregion')
cts.append(ct)
n_set = set_hotregion(cts,
max_rest_atoms=MAX_NUM_REST_ATOM,
min_h_atom=N_HATOM)
print("heavy atoms assigned", n_set)
for i, ct in enumerate(cts):
for atom in ct.atom:
if (atom.property.get(constants.REST_HOTREGION) or
atom.property.get("i_user_hotregion")):
print(i, atom.property.get(constants.REST_HOTREGION), atom.property.get("i_user_hotregion"), \
atom.property.get(constants.FEP_MAPPING), atom.element)
if (len(sys.argv) > 2):
out_fname = sys.argv[2]
with structure.StructureWriter(out_fname) as writer:
writer.extend(cts)