Select Git revision
init_iplt.py
_ring_punches.py 9.27 KiB
'''Helper functions to deal with ring punchings.'''
import ost
from ost import geom
from promod3 import sidechain, core
from collections import namedtuple
def _AddRing(rings, res, atom_names):
'''Try to add ring for given atoms in residue res to rings.'''
prof = core.StaticRuntimeProfiler.StartScoped("ring_punches::_AddRing")
# get exisiting atom positions
N = len(atom_names)
assert(N >= 3)
pos = list()
center = geom.Vec3(0, 0, 0)
for atom_name in atom_names:
a = res.FindAtom(atom_name)
if a.IsValid():
pos.append(a.GetPos())
center += a.GetPos()
# all good?
allgood = False
if len(pos) == N:
# all the ring there -> good
center = center/N
plane = geom.Plane(pos[0], pos[1], center)
allgood = True
elif res.one_letter_code == 'P' and len(pos) == 3:
# incomplete proline: 3 BB-pos there
# get ring-plane from 3 pos
plane = geom.Plane(pos[0], pos[1], pos[2])
# fit circle in plane through 3 points:
# - build 2 planes normal to ring-plane and vector connecting 2 points
# - intersect 2 planes and resulting line with ring-plane -> DONE
plane1 = geom.Plane((pos[0] + pos[1])/2, pos[0] - pos[1])
plane2 = geom.Plane((pos[1] + pos[2])/2, pos[1] - pos[2])
i_line = geom.IntersectionLine(plane1, plane2)
center = geom.IntersectionPoint(i_line, plane)
allgood = True
# add ring
if allgood:
# get ring radius
radius = 0
for p in pos:
radius = max(radius, geom.Length(p - center))
# append to list
Ring = namedtuple('Ring', 'center plane radius residue')
rings.append(Ring(center, plane, radius, res.handle))
def _IgnoreAtom(ring, atom):
'''Return true if atom is part of ring-residue.'''
return ring.residue == atom.residue
def _CheckAtomVsRing(ring, atom):
'''Check all bonds of atom for punches through ring.
Atom is ignored if it's in ring-residue or direct neighbors.
Return true if punch found.'''
# check this atom (can be view or handle!)
a1 = atom.handle
if _IgnoreAtom(ring, a1):
return False
p1 = a1.GetPos()
d1 = geom.Dot(p1 - ring.center, ring.plane.normal)
# check all bonded partners
for a_other in atom.GetBondPartners():
a2 = a_other.handle
if not a_other.IsValid() or _IgnoreAtom(ring, a2):
continue
# two bonded atoms on diff. sides of the plane?
p2 = a2.GetPos()
d2 = geom.Dot(p2 - ring.center, ring.plane.normal)
if d1*d2 < 0:
# get intersect
line = geom.Line3(p1, p2)
pi = geom.IntersectionPoint(line, ring.plane)
if geom.Length(pi - ring.center) < ring.radius:
return True
return False
def GetRings(ent):
'''Get rings for a protein structure.
A ring is only added if all ring-atoms exist or if it is a proline and
three of the atoms exist (center and radii are estimated then).
:param ent: Structure for which to detect rings.
:type ent: :class:`~ost.mol.EntityHandle` or :class:`~ost.mol.EntityView`
:return: :class:`list` of rings to perform ring checks. Each ring is a named
tuple with:
center (:class:`~ost.geom.Vec3`),
plane (:class:`~ost.geom.Plane`),
radius (:class:`float`),
residue (:class:`~ost.mol.ResidueHandle`).
'''
prof_name = 'ring_punches::GetRings'
prof = core.StaticRuntimeProfiler.StartScoped(prof_name)
rings = list()
for r in ent.residues:
if r.one_letter_code in 'YF':
_AddRing(rings, r, ["CG", "CD1", "CD2", "CE1", "CE2", "CZ"])
elif r.one_letter_code == 'W':
_AddRing(rings, r, ["CG", "CD1", "NE1", "CD2", "CE2"])
_AddRing(rings, r, ["CD2", "CE2", "CE3", "CZ2", "CZ3", "CH2"])
elif r.one_letter_code == 'H':
_AddRing(rings, r, ["CG", "CD2", "ND1", "CE1", "NE2"])
elif r.one_letter_code == 'P':
_AddRing(rings, r, ["N", "CA", "CB", "CD", "CG"])
return rings
def GetRingPunches(rings, ent):
'''Get list of residues with rings that are punched by the given structure.
:param rings: List of rings as provided by :func:`GetRings`.
:param ent: Structure for which to detect punches.
:type ent: :class:`~ost.mol.EntityHandle` or :class:`~ost.mol.EntityView`
:return: :class:`list` of residues (:class:`~ost.mol.ResidueHandle`) which
have a punched ring.
'''
prof_name = 'ring_punches::GetRingPunches'
prof = core.StaticRuntimeProfiler.StartScoped(prof_name)
ring_punches = list()
for ring in rings:
# we don't need to add residues multiple times
if ring.residue in ring_punches:
continue
# check neighborhood (3A should be enough)
for atom in ent.FindWithin(ring.center, 3):
if _CheckAtomVsRing(ring, atom):
ring_punches.append(ring.residue)
break
return ring_punches
def HasRingPunches(rings, ent):
'''Check if any ring is punched by the given structure.
This check is faster than using :func:`GetRingPunches`.
:param rings: List of rings as provided by :func:`GetRings`.
:param ent: Structure for which to detect punches.
:type ent: :class:`~ost.mol.EntityHandle` or :class:`~ost.mol.EntityView`
:return: True, iff any ring is punched
:rtype: :class:`bool`
'''
prof_name = 'ring_punches::HasRingPunches'
prof = core.StaticRuntimeProfiler.StartScoped(prof_name)
for ring in rings:
# check neighborhood (3A should be enough)
for atom in ent.FindWithin(ring.center, 3):
if _CheckAtomVsRing(ring, atom):
return True
return False
def FilterCandidates(candidates, model, gap, orig_indices=[]):
'''Remove loop candidates if they cause ring punches.
:param candidates: Loop candidates meant to fill *gap* within *model*.
Offending candidates are removed from this list.
:type candidates: :class:`~promod3.loop.LoopCandidates`
:param model: Model for which loop is to be filled.
:type model: :class:`~ost.mol.EntityHandle`
:param gap: Gap for which loop is to be filled.
:type gap: :class:`StructuralGap`.
:param orig_indices: Mapping to old indexing of candidates. If given, it
must have as many elements as *candidates*.
:type orig_indices: :class:`list`
'''
prof_name = 'ring_punches::FilterCandidates'
prof = core.StaticRuntimeProfiler.StartScoped(prof_name)
start_resnum = gap.before.GetNumber()
chain_idx = gap.GetChainIndex()
# precompute rings and range of rings to replace
chain_name = gap.GetChainName()
end_resnum = gap.after.GetNumber()
myqueryin = "cname=" + chain_name + " and rnum=" +\
str(start_resnum) + ":" + str(end_resnum)
myqueryout = "cname!=" + chain_name + " or rnum<" +\
str(start_resnum) + " or rnum>" + str(end_resnum)
model_out = model.Select(myqueryout)
rings_out = GetRings(model_out)
# filter loop candidates
lc_idx = 0
while lc_idx < len(candidates):
# get loop-entity for checks
new_loop = candidates[lc_idx].ToEntity()
rings_new = GetRings(new_loop)
check_punches = HasRingPunches(rings_out+rings_new, new_loop) or\
HasRingPunches(rings_new, model_out)
if check_punches:
candidates.Remove(lc_idx)
if orig_indices:
del orig_indices[lc_idx]
else:
lc_idx += 1
def FilterCandidatesWithSC(candidates, model, gap, orig_indices=[]):
'''Remove loop candidates if they (with sidechain) cause ring punches.
See :func:`FilterCandidates`.
'''
prof_name = 'ring_punches::FilterCandidatesWithSC'
prof = core.StaticRuntimeProfiler.StartScoped(prof_name)
start_resnum = gap.before.GetNumber()
chain_idx = gap.GetChainIndex()
cur_model = model.Copy()
# precompute rings and range of rings to replace
chain_name = gap.GetChainName()
end_resnum = gap.after.GetNumber()
myqueryin = "cname=" + chain_name + " and rnum=" +\
str(start_resnum) + ":" + str(end_resnum)
myqueryout = "cname!=" + chain_name + " or rnum<" +\
str(start_resnum) + " or rnum>" + str(end_resnum)
rings_out = GetRings(cur_model.Select(myqueryout))
# filter loop candidates
lc_idx = 0
while lc_idx < len(candidates):
# insert loop into model-copy
bb_list = candidates[lc_idx]
bb_list.InsertInto(cur_model.chains[chain_idx], start_resnum)
# add sidechains and check for clashes
sidechain.Reconstruct(cur_model, keep_sidechains=True)
models_new = cur_model.Select(myqueryin)
rings_new = GetRings(models_new)
check_punches = HasRingPunches(rings_out, models_new) or\
HasRingPunches(rings_new, cur_model)
if check_punches:
candidates.Remove(lc_idx)
if orig_indices:
del orig_indices[lc_idx]
else:
lc_idx += 1
# these methods will be exported into module
__all__ = ('GetRings', 'GetRingPunches', 'HasRingPunches', 'FilterCandidates',
'FilterCandidatesWithSC')