From 2451a2167b500183255717f26c25ab7f373e709c Mon Sep 17 00:00:00 2001
From: Gerardo Tauriello <gerardo.tauriello@unibas.ch>
Date: Wed, 4 May 2016 15:30:39 +0200
Subject: [PATCH] SCHWED-876: Refactored reconstruct_sidechain to add
consider_ligands flag
---
sidechain/doc/index.rst | 38 +-
sidechain/pymod/_reconstruct_sidechains.py | 622 ++++++++++++---------
2 files changed, 361 insertions(+), 299 deletions(-)
diff --git a/sidechain/doc/index.rst b/sidechain/doc/index.rst
index bfca8940..98c570f3 100644
--- a/sidechain/doc/index.rst
+++ b/sidechain/doc/index.rst
@@ -20,43 +20,7 @@ function:
.. literalinclude:: ../../../tests/doc/scripts/sidechain_reconstruct.py
-.. method:: Reconstruct(prot, keep_sidechains=False, build_disulfids, \
- rotamer_model="frm", consider_hbonds=True, \
- rotamer_library=None)
-
- The function takes a structure and reconstructs its sidechains given the input
- parameters.
-
- :param prot: Structure for sidechain reconstruction. Note, that the
- sidechain reconstruction gets directly applied on the
- structure itself.
-
- :param keep_sidechains: Flag, whether complete sidechains in *prot* (i.e.
- containing all required atoms) should be kept rigid
- and directly be added to the frame.
-
- :param build_disulfids: Flag, whether possible disulfid bonds should be
- searched. If a disulfid bond is found, the two
- participating cysteins are fixed and added to
- the frame.
-
- :param rotamer_model: Rotamer model to be used, can either be "frm" or "rrm"
-
- :param consider_hbonds: Flag, whether hbonds should be evaluated in the energy
- function. If set to False, no hydrogens will be built
- when building rotamers and frame and the **add_polar_hydrogens**
- flag won't have any consequences.
-
- :param rotamer_library: A rotamer library to extract the
- rotamers from.
-
-
- :type prot: :class:`ost.mol.EntityHandle`
- :type keep_sidechains: :class:`bool`
- :type build_disulfids: :class:`bool`
- :type rotamer_model: :class:`str`
- :type consider_hbonds: :class:`bool`
- :type rotamer_library: :class:`BBDepRotamerLib` / :class:`RotamerLib`
+.. autofunction:: Reconstruct
Sidechain Module Functionality
diff --git a/sidechain/pymod/_reconstruct_sidechains.py b/sidechain/pymod/_reconstruct_sidechains.py
index c385aac2..bb8f5410 100644
--- a/sidechain/pymod/_reconstruct_sidechains.py
+++ b/sidechain/pymod/_reconstruct_sidechains.py
@@ -2,54 +2,36 @@ from . import _sidechain as sidechain
from ost import geom
from ost import mol
-def Reconstruct(ent, keep_sidechains = False, build_disulfids = True,
- rotamer_model = "frm", consider_hbonds = True,
- rotamer_library = None):
-
- if rotamer_model.lower() == "frm":
- use_frm = True
- elif rotamer_model.lower() == "rrm":
- use_frm = False
- else:
- raise RuntimeError("Only \"rrm\" and \"frm\" allowed for rotamer_model!")
-
-
- rotamer_settings = sidechain.RotamerSettings()
- rotamer_settings.consider_hbonds = consider_hbonds
- if rotamer_library == None:
- rotamer_library = sidechain.LoadDunbrackLib()
- bbdep = False
- if type(rotamer_library) is sidechain.BBDepRotamerLib:
- bbdep = True
- residues_with_rotamer_group = list()
- rotamer_groups = list()
- # ignore ligand chain and any non-peptides
- prot = ent.Select("peptide=true and cname!='_'")
- incomplete_sidechains = list()
-
- #extract rotamer ids
- rotamer_ids = [sidechain.ALA] * len(prot.residues)
- for i,r in enumerate(prot.residues):
+###############################################################################
+# helper functions
+def _GetRotamerIDs(res_list):
+ '''Return list (length = len(res_list)) of rotamer IDs for all residues.'''
+ rotamer_ids = [sidechain.ALA] * len(res_list)
+ for i,r in enumerate(res_list):
rot_id = sidechain.TLCToRotID(r.GetName())
if rot_id == sidechain.XXX:
continue # no idea what it is, so we stick with ALA
# => don't model sidechain
rotamer_ids[i] = rot_id
-
- #extract dihedral angles
- phi_angles = [0.0] * len(prot.residues)
- psi_angles = [0.0] * len(prot.residues)
-
- for i,r in enumerate(prot.residues):
-
+ return rotamer_ids
+
+def _GetDihedrals(res_list):
+ '''Extract dihedral angles for all residues.
+ Returns phi and psi angles as 2 lists with length = len(res_list).
+ '''
+ phi_angles = [0.0] * len(res_list)
+ psi_angles = [0.0] * len(res_list)
+ for i,r in enumerate(res_list):
+ # def. fallback = helix
phi = -1.0472
psi = -0.7854
-
+
+ # try to get phi from torsion angles
tor = r.GetPhiTorsion()
if tor.IsValid():
phi = tor.GetAngle()
elif i > 0:
- r_prev = prot.residues[i-1]
+ r_prev = res_list[i-1]
if mol.InSequence(r_prev.handle,r.handle):
c_prev = r_prev.FindAtom("C")
n = r.FindAtom("N")
@@ -58,12 +40,13 @@ def Reconstruct(ent, keep_sidechains = False, build_disulfids = True,
if c_prev.IsValid() and n.IsValid() and ca.IsValid() and c.IsValid():
phi = geom.DihedralAngle(c_prev.GetPos(),n.GetPos(),
ca.GetPos(),c.GetPos())
-
+
+ # try to get psi from torsion angles
tor = r.GetPsiTorsion()
if tor.IsValid():
psi = tor.GetAngle()
- elif i < (len(prot.residues) - 1):
- r_next = prot.residues[i+1]
+ elif i < (len(res_list) - 1):
+ r_next = res_list[i+1]
if mol.InSequence(r.handle,r_next.handle):
n = r.FindAtom("N")
ca = r.FindAtom("CA")
@@ -72,222 +55,154 @@ def Reconstruct(ent, keep_sidechains = False, build_disulfids = True,
if n.IsValid() and ca.IsValid() and c.IsValid() and n_next.IsValid():
psi = geom.DihedralAngle(n.GetPos(), ca.GetPos(),
c.GetPos(), n_next.GetPos())
-
-
+ # store
phi_angles[i] = phi
psi_angles[i] = psi
+ return phi_angles, psi_angles
-
- #set nter and cter
- for c in prot.chains:
- c.residues[0].SetIntProp("n_ter",1)
- c.residues[-1].SetIntProp("c_ter",1)
-
-
- #build up frame and cysteins
- frame_residues = list()
-
- #build up backbone frame
- for i,r in enumerate(prot.residues):
+def _AddBackboneFrameResidues(frame_residues, res_list, rotamer_ids,
+ rotamer_settings, phi_angles):
+ '''Update frame_residues (list) with BackboneFrameResidues for res_list.'''
+ for i,r in enumerate(res_list):
try:
- frame_residue = sidechain.ConstructBackboneFrameResidue(r.handle,
- rotamer_ids[i],
- i,rotamer_settings,
- phi_angles[i],
- r.HasProp("n_ter"),
- r.HasProp("c_ter"))
+ frame_residue = sidechain.ConstructBackboneFrameResidue(\
+ r.handle, rotamer_ids[i], i, rotamer_settings,
+ phi_angles[i], r.HasProp("n_ter"), r.HasProp("c_ter"))
frame_residues.append(frame_residue)
except:
continue
+def _AddLigandFrameResidues(frame_residues, ent_lig, rotamer_settings, offset):
+ '''Update frame_residues (list) with FrameResidues for res. in ent_lig.
+ Set offset >= number of non-ligand residues (used for residue_index).
+ '''
+ # TODO? Split peptides and non-peptides?
+ # for now: add as static carbons
+ for i,r in enumerate(ent_lig.residues):
+ try:
+ frame_residue = sidechain.ConstructFrameResidue(r.handle, offset+i)
+ frame_residues.append(frame_residue)
+ except:
+ continue
- frame = None
-
- if build_disulfids:
-
- cystein_indices = list()
-
- if keep_sidechains:
- #let's add all complete sidechains except the cysteins to the frame
- for i,r in enumerate(prot.residues):
-
- if rotamer_ids[i] == sidechain.CYS:
- cystein_indices.append(i)
- continue
-
- if rotamer_ids[i] == sidechain.ALA or rotamer_ids[i] == sidechain.GLY:
- continue #no sidechain to model
-
- try:
- frame_residue = sidechain.ConstructSidechainFrameResidue(r.handle,
- rotamer_ids[i],
- i,rotamer_settings)
- frame_residues.append(frame_residue)
- except:
- incomplete_sidechains.append(i)
- else:
- #let's add everything except cysteins to the incomplete sidechains
- for i,r in enumerate(prot.residues):
-
- if rotamer_ids[i] == sidechain.CYS:
- cystein_indices.append(i)
- continue
+def _AddSidechainFrameResidues(frame_residues, incomplete_sidechains,
+ keep_sidechains, res_list, rotamer_ids,
+ rotamer_settings, cystein_indices=None):
+ '''Update frame_residues (list) with SidechainFrameResidues for res_list,
+ incomplete_sidechains (list of indices) with sidechains to be constructed,
+ and (if given) cystein_indices (list of indices) with all CYS (appended).
+ Each residue can only end up in one of the 3 lists.
+ '''
+ if keep_sidechains:
+ # try to generate frame residues for all existing side chains
+ # skip non-existing sidechains and CYS (if cystein_indices) and update
+ # incomplete_sidechains and cystein_indices
+ for i,r in enumerate(res_list):
+
+ if cystein_indices is not None and rotamer_ids[i] == sidechain.CYS:
+ cystein_indices.append(i)
+ continue
- if rotamer_ids[i] == sidechain.ALA or rotamer_ids[i] == sidechain.GLY:
- continue #no sidechain to model
+ if rotamer_ids[i] == sidechain.ALA or rotamer_ids[i] == sidechain.GLY:
+ continue # no sidechain to model
+ try:
+ frame_residue = sidechain.ConstructSidechainFrameResidue(\
+ r.handle, rotamer_ids[i], i, rotamer_settings)
+ frame_residues.append(frame_residue)
+ except:
incomplete_sidechains.append(i)
+ else:
+ # no frame residues to create, just update incomplete_sidechains
+ # and cystein_indices if needed
+ for i,r in enumerate(res_list):
- #let's generate the frame without any cystein sidechains
- #this is required for the disulfid score evaluation
- frame = sidechain.Frame(frame_residues)
-
- #some info we have to keep track of when evaluating disulfid bonds
- cystein_rotamers = list()
- disulfid_indices = list()
- cys_ca_positions = list()
- cys_cb_positions = list()
-
- for i in cystein_indices:
-
- rot_grop = None
- r = prot.residues[i]
- ca = r.FindAtom("CA")
- cb = r.FindAtom("CB")
-
- if not (ca.IsValid() and cb.IsValid()):
+ if cystein_indices is not None and rotamer_ids[i] == sidechain.CYS:
+ cystein_indices.append(i)
continue
- cys_ca_positions.append(ca.GetPos())
- cys_cb_positions.append(cb.GetPos())
-
- if use_frm:
- if bbdep:
- rot_group = sidechain.ConstructFRMRotamerGroup(r.handle,
- sidechain.CYD,i,
- rotamer_library,
- rotamer_settings,
- phi_angles[i],
- psi_angles[i])
- else:
- rot_group = sidechain.ConstructFRMRotamerGroup(r.handle,
- sidechain.CYD,i,
- rotamer_library,
- rotamer_settings)
- else:
- if bbdep:
- rot_group = sidechain.ConstructRRMRotamerGroup(r.handle,
- sidechain.CYD,i,
- rotamer_library,
- rotamer_settings,
- phi_angles[i],
- psi_angles[i])
- else:
- rot_group = sidechain.ConstructRRMRotamerGroup(r.handle,
- sidechain.CYD,i,
- rotamer_library,
- rotamer_settings)
-
- frame.AddFrameEnergy(rot_group)
- cystein_rotamers.append(rot_group)
-
- for i in range(len(cystein_rotamers)):
- for j in range(i+1, len(cystein_rotamers)):
-
- if geom.Distance(cys_ca_positions[i], cys_ca_positions[j]) > 8.0:
- continue #they're too far for a disulfid bond
-
- if cystein_indices[i] in disulfid_indices or cystein_indices[j] in disulfid_indices:
- continue #one of them already participates in a disulfid bond!
- #that one might be bether but right now is first come
- #first served principle
-
- min_score = float("inf")
- min_index_k = -1
- min_index_l = -1
-
- for k in range(len(cystein_rotamers[i])):
- for l in range(len(cystein_rotamers[j])):
- score = sidechain.DisulfidScore(cystein_rotamers[i][k],
- cystein_rotamers[j][l],
- cys_ca_positions[i],
- cys_cb_positions[i],
- cys_ca_positions[j],
- cys_cb_positions[j])
- if score < min_score:
- min_index_k = k
- min_index_l = l
- min_score = score
-
- if min_score < 45.0:
- disulfid_indices.append(cystein_indices[i])
- disulfid_indices.append(cystein_indices[j])
- particle_list = list()
- particle_list.append(cystein_rotamers[i][min_index_k][0])
- new_frame_residue_i = sidechain.FrameResidue(particle_list,cystein_indices[i])
- particle_list = list()
- particle_list.append(cystein_rotamers[j][min_index_l][0])
- new_frame_residue_j = sidechain.FrameResidue(particle_list,cystein_indices[j])
- frame_residues.append(new_frame_residue_i)
- frame_residues.append(new_frame_residue_j)
- #set the position in the proteins residues
- cystein_rotamers[i][min_index_k].ApplyOnResidue(prot.residues[cystein_indices[i]].handle,
- consider_hydrogens=False)
- cystein_rotamers[j][min_index_l].ApplyOnResidue(prot.residues[cystein_indices[j]].handle,
- consider_hydrogens=False)
- sidechain.ConnectSidechain(prot.residues[cystein_indices[i]].handle,sidechain.CYS)
- sidechain.ConnectSidechain(prot.residues[cystein_indices[j]].handle,sidechain.CYS)
-
-
- #All cysteins participating in a disulfid bond have been applied to the
- #structure and added to the frame.
- #All remaining ones have to be handled according the given flags.
- for idx in cystein_indices:
- if idx in disulfid_indices:
- continue #it's all fine
- if keep_sidechains:
- try:
- frame_residue = sidechain.ConstructSidechainFrameResidue(prot.residues[idx].handle,
- rotamer_ids[idx],
- idx,rotamer_settings)
- frame_residues.append(frame_residue)
- except:
- incomplete_sidechains.append(idx)
- else:
- incomplete_sidechains.append(idx)
-
- #we finally have to rebuild the frame if any disulfid bonds have been built
- if len(disulfid_indices) > 0:
- frame = sidechain.Frame(frame_residues)
-
-
- else:
+ if rotamer_ids[i] == sidechain.ALA or rotamer_ids[i] == sidechain.GLY:
+ continue # no sidechain to model
+
+ incomplete_sidechains.append(i)
+
+def _AddCysteinFrameResidues(frame_residues, incomplete_sidechains,
+ keep_sidechains, res_list, rotamer_ids,
+ rotamer_settings, cystein_indices,
+ disulfid_indices, disulfid_rotamers):
+ '''Update frame_residues (list) with cysteins.
+ Parameters as in _AddSidechainFrameResidues.
+ Some cysteins (in disulfid_indices) get special treatment as disulfid
+ bridges (disulfid_indices, disulfid_rotamers from _GetDisulfidBridges).
+ '''
+ # handle cysteins participating in a disulfid bond
+ for cys_idx, cys_rot in zip(disulfid_indices, disulfid_rotamers):
+ # add FrameResidue
+ frame_residue = sidechain.FrameResidue([cys_rot[0]], cys_idx)
+ frame_residues.append(frame_residue)
+ # set the position in the proteins residues
+ cys_rot.ApplyOnResidue(res_list[cys_idx].handle,
+ consider_hydrogens=False)
+ sidechain.ConnectSidechain(res_list[cys_idx].handle, sidechain.CYS)
+
+ # add remaining ones according the given flags
+ for idx in cystein_indices:
+ if idx in disulfid_indices:
+ continue # already handled
if keep_sidechains:
- for i,r in enumerate(prot.residues):
- if rotamer_ids[i] == sidechain.ALA or rotamer_ids[i] == sidechain.GLY:
- continue #no sidechain to model
- try:
- frame_residue = sidechain.ConstructSidechainFrameResidue(r.handle,
- rotamer_ids[i],
- i,rotamer_settings)
- frame_residues.append(frame_residue)
- except:
- incomplete_sidechains.append(i)
+ try:
+ frame_residue = sidechain.ConstructSidechainFrameResidue(\
+ res_list[idx].handle, rotamer_ids[idx],
+ idx, rotamer_settings)
+ frame_residues.append(frame_residue)
+ except:
+ incomplete_sidechains.append(idx)
else:
- for i,r in enumerate(prot.residues):
- if rotamer_ids[i] not in [sidechain.GLY, sidechain.ALA]:
- incomplete_sidechains.append(i)
-
- frame = sidechain.Frame(frame_residues)
+ incomplete_sidechains.append(idx)
-
- #build rotamers for incomplete sidechains
- for i in incomplete_sidechains:
-
- r = prot.residues[i]
- rot_id = rotamer_ids[i]
+def _GetRotamerGroup(res_handle, rot_id, res_idx, rot_lib, rot_settings,
+ phi, psi, use_frm, bbdep):
+ '''Get RotamerGroup for res_handle according to settings.'''
+ if use_frm:
+ if bbdep:
+ return sidechain.ConstructFRMRotamerGroup(res_handle, rot_id,
+ res_idx, rot_lib,
+ rot_settings, phi, psi)
+ else:
+ return sidechain.ConstructFRMRotamerGroup(res_handle, rot_id,
+ res_idx, rot_lib,
+ rot_settings)
+ else:
+ if bbdep:
+ return sidechain.ConstructRRMRotamerGroup(res_handle, rot_id,
+ res_idx, rot_lib,
+ rot_settings, phi, psi)
+ else:
+ return sidechain.ConstructRRMRotamerGroup(res_handle, rot_id,
+ res_idx, rot_lib,
+ rot_settings)
+
+def _GetRotamerGroups(res_list, rot_ids, indices, rot_lib, rot_settings,
+ phi_angles, psi_angles, use_frm, bbdep, frame_residues):
+ '''Get list of rotamer groups from subset of res_list.
+ Residues are chosen as res_list[i] for i in indices and only if a rotamer
+ group can be created (e.g. no ALA, GLY).
+ Rotamer groups are filtered to keep only best ones (given frame).
+ Returns list of rotamer groups and list of res. indices they belong to.
+ '''
+ # res.index (res_list[i]) for each modelled sc
+ residues_with_rotamer_group = list()
+ # linked to residue in residues_with_rotamer_group
+ rotamer_groups = list()
+ # get frame for score evaluation
+ frame = sidechain.Frame(frame_residues)
+ # build rotamers for chosen sidechains
+ for i in indices:
+ # get rotamer ID
+ r = res_list[i]
+ rot_id = rot_ids[i]
- if(rot_id == sidechain.ALA or rot_id == sidechain.GLY):
+ if rot_id == sidechain.ALA or rot_id == sidechain.GLY:
continue
if rot_id == sidechain.CYS:
@@ -300,7 +215,7 @@ def Reconstruct(ent, keep_sidechains = False, build_disulfids = True,
omega = tor.GetAngle()
elif i > 0:
# fallback computation of omega as in OST-code
- prev = prot.residues[i-1]
+ prev = res_list[i-1]
if prev.IsValid() and prev.IsPeptideLinking():
ca_prev = prev.FindAtom("CA")
c_prev = prev.FindAtom("C")
@@ -319,49 +234,232 @@ def Reconstruct(ent, keep_sidechains = False, build_disulfids = True,
else:
rot_id = sidechain.TPR
+ # get RotamerGroup
try:
- if use_frm:
- if bbdep:
- rot_group = sidechain.ConstructFRMRotamerGroup(r.handle,rot_id,i,
- rotamer_library,rotamer_settings,
- phi_angles[i],psi_angles[i])
- else:
- rot_group = sidechain.ConstructFRMRotamerGroup(r.handle,rot_id,i,
- rotamer_library,rotamer_settings)
- else:
- if bbdep:
- rot_group = sidechain.ConstructRRMRotamerGroup(r.handle,rot_id,i,
- rotamer_library,rotamer_settings,
- phi_angles[i],psi_angles[i])
- else:
- rot_group = sidechain.ConstructRRMRotamerGroup(r.handle,rot_id,i,
- rotamer_library,rotamer_settings)
-
+ rot_group = _GetRotamerGroup(r.handle, rot_id, i, rot_lib,
+ rot_settings, phi_angles[i],
+ psi_angles[i], use_frm, bbdep)
except:
continue
+ # keep best ones
rot_group.CalculateInternalEnergies()
frame.SetFrameEnergy(rot_group)
rot_group.ApplySelfEnergyThresh()
rotamer_groups.append(rot_group)
- residues_with_rotamer_group.append(i)
+ residues_with_rotamer_group.append(i)
+
+ return rotamer_groups, residues_with_rotamer_group
+
+def _GetDisulfidBridges(frame_residues, cystein_indices, res_list, rotamer_library,
+ use_frm, bbdep, rotamer_settings, phi_angles, psi_angles):
+ '''Get disulfid bridges for CYS and according rotamers.
+ CYS are identified by by items in cystein_indices (into res_list).
+ Returns: disulfid_indices: list of res. index in bridge,
+ disulfid_rotamers: list of rotamers (best one for bridge).
+ '''
+ # this is required for the disulfid score evaluation
+ frame = sidechain.Frame(frame_residues)
+
+ # some info we have to keep track of when evaluating disulfid bonds
+ cystein_rotamers = list()
+ cys_ca_positions = list()
+ cys_cb_positions = list()
+
+ for i in cystein_indices:
+ # check ca, cb
+ r = res_list[i]
+ ca = r.FindAtom("CA")
+ cb = r.FindAtom("CB")
+ if not (ca.IsValid() and cb.IsValid()):
+ continue
+ cys_ca_positions.append(ca.GetPos())
+ cys_cb_positions.append(cb.GetPos())
+ # get RotamerGroup
+ rot_group = _GetRotamerGroup(r.handle, sidechain.CYD, i, rotamer_library,
+ rotamer_settings, phi_angles[i],
+ psi_angles[i], use_frm, bbdep)
+ frame.AddFrameEnergy(rot_group)
+ cystein_rotamers.append(rot_group)
+
+ # get CYS with disulfid bonds and the chosen rotamers
+ disulfid_indices = list()
+ disulfid_rotamers = list()
+ for i in range(len(cystein_rotamers)):
+ for j in range(i+1, len(cystein_rotamers)):
+
+ # too far for a disulfid bond?
+ if geom.Distance(cys_ca_positions[i], cys_ca_positions[j]) > 8.0:
+ continue
+
+ # already done? NOTE: new one might be better, but here we do
+ # first come, first served
+ if cystein_indices[i] in disulfid_indices \
+ or cystein_indices[j] in disulfid_indices:
+ continue
+ min_score = float("inf")
+ min_index_k = -1
+ min_index_l = -1
+
+ for k in range(len(cystein_rotamers[i])):
+ for l in range(len(cystein_rotamers[j])):
+ score = sidechain.DisulfidScore(cystein_rotamers[i][k],
+ cystein_rotamers[j][l],
+ cys_ca_positions[i],
+ cys_cb_positions[i],
+ cys_ca_positions[j],
+ cys_cb_positions[j])
+ if score < min_score:
+ min_index_k = k
+ min_index_l = l
+ min_score = score
+
+ if min_score < 45.0:
+ # update indices
+ cys_idx_i = cystein_indices[i]
+ cys_idx_j = cystein_indices[j]
+ cys_rot_i = cystein_rotamers[i][min_index_k]
+ cys_rot_j = cystein_rotamers[j][min_index_l]
+ disulfid_indices.append(cys_idx_i)
+ disulfid_indices.append(cys_idx_j)
+ disulfid_rotamers.append(cys_rot_i)
+ disulfid_rotamers.append(cys_rot_j)
+
+ return disulfid_indices, disulfid_rotamers
+
+###############################################################################
+
+def Reconstruct(ent, keep_sidechains=False, build_disulfids=True,
+ rotamer_model="frm", consider_hbonds=True,
+ consider_ligands=True, rotamer_library=None):
+ '''Reconstruct sidechains for the given structure.
+
+ :param ent: Structure for sidechain reconstruction. Note, that the
+ sidechain reconstruction gets directly applied on the
+ structure itself.
+
+ :param keep_sidechains: Flag, whether complete sidechains in *ent* (i.e.
+ containing all required atoms) should be kept rigid
+ and directly be added to the frame.
+
+ :param build_disulfids: Flag, whether possible disulfid bonds should be
+ searched. If a disulfid bond is found, the two
+ participating cysteins are fixed and added to
+ the frame.
+
+ :param rotamer_model: Rotamer model to be used, can either be "frm" or "rrm"
+
+ :param consider_hbonds: Flag, whether hbonds should be evaluated in the
+ energy function. If set to False, no hydrogens will
+ be built when building rotamers and frame and the
+ **add_polar_hydrogens** flag won't have any
+ consequences.
+
+ :param consider_ligands: Flag, whether to add ligands (anything in chain
+ '_') as static objects.
+
+ :param rotamer_library: A rotamer library to extract the rotamers from. The
+ default is the :meth:`Dunbrack <LoadDunbrackLib>`
+ library.
+
+
+ :type ent: :class:`ost.mol.EntityHandle`
+ :type keep_sidechains: :class:`bool`
+ :type build_disulfids: :class:`bool`
+ :type rotamer_model: :class:`str`
+ :type consider_hbonds: :class:`bool`
+ :type consider_ligands: :class:`bool`
+ :type rotamer_library: :class:`BBDepRotamerLib` / :class:`RotamerLib`
+ '''
+
+ # setup settings
+ if rotamer_model.lower() == "frm":
+ use_frm = True
+ elif rotamer_model.lower() == "rrm":
+ use_frm = False
+ else:
+ raise RuntimeError("Only \"rrm\" and \"frm\" allowed for rotamer_model!")
+ rotamer_settings = sidechain.RotamerSettings()
+ rotamer_settings.consider_hbonds = consider_hbonds
+ if rotamer_library == None:
+ rotamer_library = sidechain.LoadDunbrackLib()
+ bbdep = False
+ if type(rotamer_library) is sidechain.BBDepRotamerLib:
+ bbdep = True
+
+ # take out ligand chain and any non-peptides
+ prot = ent.Select("peptide=true and cname!='_'")
+
+ # parse residues (all lists of length len(prot.residues))
+ rotamer_ids = _GetRotamerIDs(prot.residues)
+ phi_angles, psi_angles = _GetDihedrals(prot.residues)
+
+ # set nter and cter (needed in _AddBackboneFrameResidues)
+ for c in prot.chains:
+ c.residues[0].SetIntProp("n_ter",1)
+ c.residues[-1].SetIntProp("c_ter",1)
+
+ # build up frame
+ frame_residues = list() # list of frame residues connected to frame
+ incomplete_sidechains = list() # residue indices
+ _AddBackboneFrameResidues(frame_residues, prot.residues, rotamer_ids,
+ rotamer_settings, phi_angles)
+
+ # add ligands?
+ if consider_ligands:
+ ligs = ent.Select("cname='_'")
+ offset = len(prot.residues)
+ _AddLigandFrameResidues(frame_residues, ligs, rotamer_settings, offset)
+
+ # check special handling of cysteins
+ if build_disulfids:
+ # residue indices of cysteins
+ cystein_indices = list()
+ # update frame_residues, incomplete_sidechains, cystein_indices
+ _AddSidechainFrameResidues(frame_residues, incomplete_sidechains,
+ keep_sidechains, prot.residues, rotamer_ids,
+ rotamer_settings, cystein_indices)
+ # update frame_residues, incomplete_sidechains with cysteins (if needed)
+ if len(cystein_indices) > 0:
+ # get disulfid bridges and according rotamers
+ disulfid_indices, disulfid_rotamers = \
+ _GetDisulfidBridges(frame_residues, cystein_indices, prot.residues,
+ rotamer_library, use_frm, bbdep,
+ rotamer_settings, phi_angles, psi_angles)
+ # update frame_residues, incomplete_sidechains
+ _AddCysteinFrameResidues(frame_residues, incomplete_sidechains,
+ keep_sidechains, prot.residues, rotamer_ids,
+ rotamer_settings, cystein_indices,
+ disulfid_indices, disulfid_rotamers)
+ else:
+ # update frame_residues, incomplete_sidechains
+ _AddSidechainFrameResidues(frame_residues, incomplete_sidechains,
+ keep_sidechains, prot.residues, rotamer_ids,
+ rotamer_settings)
+
+ # get rotamer groups and residues they're linked to
+ rotamer_groups, residues_with_rotamer_group = \
+ _GetRotamerGroups(prot.residues, rotamer_ids, incomplete_sidechains,
+ rotamer_library, rotamer_settings, phi_angles,
+ psi_angles, use_frm, bbdep, frame_residues)
+
+ # set up graph and solve to get best rotamers
if use_frm:
graph = sidechain.Graph.CreateFromFRMList(rotamer_groups)
else:
graph = sidechain.Graph.CreateFromRRMList(rotamer_groups)
-
solution = graph.Solve(100000000,0.02)
-
+ # update structure
for i,rot_group,sol in zip(residues_with_rotamer_group,rotamer_groups,solution):
try:
- rot_group[sol].ApplyOnResidue(prot.residues[i].handle,consider_hydrogens=False)
- sidechain.ConnectSidechain(prot.residues[i].handle,rotamer_ids[i])
+ res_handle = prot.residues[i].handle
+ rot_group[sol].ApplyOnResidue(res_handle, consider_hydrogens=False)
+ sidechain.ConnectSidechain(res_handle, rotamer_ids[i])
except:
- print "there is a backbone atom missing... ",prot.residues[i].GetQualifiedName()
-
+ print "there is a backbone atom missing... ", res_handle.GetQualifiedName()
# these methods will be exported into module
__all__ = ('Reconstruct',)
--
GitLab