Skip to content
Snippets Groups Projects
Commit 2451a216 authored by Gerardo Tauriello's avatar Gerardo Tauriello
Browse files

SCHWED-876: Refactored reconstruct_sidechain to add consider_ligands flag

parent a76f54cb
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
sidechain/doc/index.rst
+ 1
− 37
View file @ 2451a216
......@@ -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
......
sidechain/pymod/_reconstruct_sidechains.py
+ 360
− 262
View file @ 2451a216
......@@ -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")
......@@ -59,11 +41,12 @@ def Reconstruct(ent, keep_sidechains = False, build_disulfids = True,
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,46 +55,51 @@ 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()
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:
#let's add all complete sidechains except the cysteins to the frame
for i,r in enumerate(prot.residues):
# 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 rotamer_ids[i] == sidechain.CYS:
if cystein_indices is not None and rotamer_ids[i] == sidechain.CYS:
cystein_indices.append(i)
continue
......@@ -119,17 +107,17 @@ def Reconstruct(ent, keep_sidechains = False, build_disulfids = True,
continue # no sidechain to model
try:
frame_residue = sidechain.ConstructSidechainFrameResidue(r.handle,
rotamer_ids[i],
i,rotamer_settings)
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):
# no frame residues to create, just update incomplete_sidechains
# and cystein_indices if needed
for i,r in enumerate(res_list):
if rotamer_ids[i] == sidechain.CYS:
if cystein_indices is not None and rotamer_ids[i] == sidechain.CYS:
cystein_indices.append(i)
continue
......@@ -138,69 +126,176 @@ def Reconstruct(ent, keep_sidechains = False, build_disulfids = True,
incomplete_sidechains.append(i)
#let's generate the frame without any cystein sidechains
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:
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:
incomplete_sidechains.append(idx)
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:
continue
if rot_id == sidechain.CYS:
rot_id = sidechain.CYH
if rot_id == sidechain.PRO:
tor = r.GetOmegaTorsion()
omega = None
if tor.IsValid():
omega = tor.GetAngle()
elif i > 0:
# fallback computation of omega as in OST-code
prev = res_list[i-1]
if prev.IsValid() and prev.IsPeptideLinking():
ca_prev = prev.FindAtom("CA")
c_prev = prev.FindAtom("C")
n = r.FindAtom("N")
ca = r.FindAtom("CA")
valid = ca_prev.IsValid() and c_prev.IsValid() \
and n.IsValid() and ca.IsValid()
if valid and mol.BondExists(c_prev.handle, n.handle):
omega = geom.DihedralAngle(ca_prev.GetPos(),
c_prev.GetPos(),
n.GetPos(), ca.GetPos())
# omega not set if prev. res. missing
if omega is not None:
if abs(omega) < 1.57:
rot_id = sidechain.CPR
else:
rot_id = sidechain.TPR
# get RotamerGroup
try:
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)
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()
disulfid_indices = list()
cys_ca_positions = list()
cys_cb_positions = list()
for i in cystein_indices:
rot_grop = None
r = prot.residues[i]
# 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())
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)
# 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 #they're too far for a disulfid bond
continue
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
# 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
......@@ -220,148 +315,151 @@ def Reconstruct(ent, keep_sidechains = False, build_disulfids = True,
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)
# 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
#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)
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.
else:
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)
else:
for i,r in enumerate(prot.residues):
if rotamer_ids[i] not in [sidechain.GLY, sidechain.ALA]:
incomplete_sidechains.append(i)
: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.
frame = sidechain.Frame(frame_residues)
: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"
#build rotamers for incomplete sidechains
for i in incomplete_sidechains:
: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.
r = prot.residues[i]
rot_id = rotamer_ids[i]
:param consider_ligands: Flag, whether to add ligands (anything in chain
'_') as static objects.
if(rot_id == sidechain.ALA or rot_id == sidechain.GLY):
continue
:param rotamer_library: A rotamer library to extract the rotamers from. The
default is the :meth:`Dunbrack <LoadDunbrackLib>`
library.
if rot_id == sidechain.CYS:
rot_id = sidechain.CYH
if rot_id == sidechain.PRO:
tor = r.GetOmegaTorsion()
omega = None
if tor.IsValid():
omega = tor.GetAngle()
elif i > 0:
# fallback computation of omega as in OST-code
prev = prot.residues[i-1]
if prev.IsValid() and prev.IsPeptideLinking():
ca_prev = prev.FindAtom("CA")
c_prev = prev.FindAtom("C")
n = r.FindAtom("N")
ca = r.FindAtom("CA")
valid = ca_prev.IsValid() and c_prev.IsValid() \
and n.IsValid() and ca.IsValid()
if valid and mol.BondExists(c_prev.handle, n.handle):
omega = geom.DihedralAngle(ca_prev.GetPos(),
c_prev.GetPos(),
n.GetPos(), ca.GetPos())
# omega not set if prev. res. missing
if omega is not None:
if abs(omega) < 1.57:
rot_id = sidechain.CPR
else:
rot_id = sidechain.TPR
: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`
'''
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])
# setup settings
if rotamer_model.lower() == "frm":
use_frm = True
elif rotamer_model.lower() == "rrm":
use_frm = False
else:
rot_group = sidechain.ConstructRRMRotamerGroup(r.handle,rot_id,i,
rotamer_library,rotamer_settings)
raise RuntimeError("Only \"rrm\" and \"frm\" allowed for rotamer_model!")
except:
continue
rot_group.CalculateInternalEnergies()
frame.SetFrameEnergy(rot_group)
rot_group.ApplySelfEnergyThresh()
rotamer_groups.append(rot_group)
residues_with_rotamer_group.append(i)
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',)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment