Source code for promod3.modelling._ring_punches

'''Helper functions to deal with ring punchings.'''
import ost
from ost import geom
from promod3 import core
from _reconstruct_sidechains import ReconstructSidechains
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

[docs]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
[docs]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
[docs]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
[docs]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
[docs]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 ReconstructSidechains(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')