From 5e07329bc2b2b05faef8765aeb5299741a91bc0b Mon Sep 17 00:00:00 2001
From: Gabriel Studer <gabriel.studer@unibas.ch>
Date: Tue, 2 Jul 2024 10:36:53 +0200
Subject: [PATCH] avoid usage of n_ter/c_ter generic props when reconstructing
sidechains
There was a collision of properties with same name that are set in
OpenStructure energy minimization. If you did a sidechain reconstruction
where the properties are set, then you attach some stuff to terminal
residues, non-terminal residues still have the n_ter/c_ter property set.
This crashed the parameterization of the energy minimization.
---
modelling/pymod/_reconstruct_sidechains.py | 79 ++++++++++++++--------
1 file changed, 51 insertions(+), 28 deletions(-)
diff --git a/modelling/pymod/_reconstruct_sidechains.py b/modelling/pymod/_reconstruct_sidechains.py
index e10f1948..f50f2860 100644
--- a/modelling/pymod/_reconstruct_sidechains.py
+++ b/modelling/pymod/_reconstruct_sidechains.py
@@ -84,14 +84,17 @@ def _GetDihedrals(res_list):
return phi_angles, psi_angles
def _AddBackboneFrameResidues(frame_residues, res_list, rotamer_ids,
- rotamer_constructor, phi_angles, psi_angles):
+ rotamer_constructor, phi_angles, psi_angles,
+ n_ter_residues, c_ter_residues):
'''Update frame_residues (list) with BackboneFrameResidues for res_list.'''
for i,r in enumerate(res_list):
try:
+ is_n_ter = r.handle.GetHashCode() in n_ter_residues
+ is_c_ter = r.handle.GetHashCode() in c_ter_residues
frame_residue = rotamer_constructor.ConstructBackboneFrameResidue(\
r.handle, rotamer_ids[i], i,
phi_angles[i], psi_angles[i],
- r.HasProp("n_ter"), r.HasProp("c_ter"))
+ is_n_ter, is_c_ter)
frame_residues.append(frame_residue)
except:
continue
@@ -149,7 +152,8 @@ def _AddLigandFrameResidues(frame_residues, ent_lig, rotamer_constructor, offset
def _AddSidechainFrameResidues(frame_residues, incomplete_sidechains,
keep_sidechains, res_list, rotamer_ids,
phi_angles, psi_angles,
- rotamer_constructor, cystein_indices=None):
+ rotamer_constructor, cystein_indices,
+ n_ter_residues, c_ter_residues):
'''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).
@@ -166,11 +170,13 @@ def _AddSidechainFrameResidues(frame_residues, incomplete_sidechains,
continue
try:
+ is_n_ter = r.handle.GetHashCode() in n_ter_residues
+ is_c_ter = r.handle.GetHashCode() in c_ter_residues
frame_residue = rotamer_constructor.ConstructSidechainFrameResidue(\
r.handle, rotamer_ids[i], i, phi_angles[i],
- psi_angles[i], r.HasProp("n_ter"), r.HasProp("c_ter"))
+ psi_angles[i], is_n_ter, is_c_ter)
frame_residues.append(frame_residue)
- except:
+ except Exception as e:
incomplete_sidechains.append(i)
else:
# no frame residues to create, just update incomplete_sidechains
@@ -187,7 +193,8 @@ def _AddCysteinFrameResidues(frame_residues, incomplete_sidechains,
keep_sidechains, res_list, rotamer_ids,
phi_angles, psi_angles,
rot_constructor, cystein_indices,
- disulfid_indices, disulfid_rotamers):
+ disulfid_indices, disulfid_rotamers,
+ n_ter_residues, c_ter_residues):
'''Update frame_residues (list) with cysteins.
Parameters as in _AddSidechainFrameResidues.
Some cysteins (in disulfid_indices) get special treatment as disulfid
@@ -208,11 +215,12 @@ def _AddCysteinFrameResidues(frame_residues, incomplete_sidechains,
continue # already handled
if keep_sidechains:
try:
+ is_n_ter = res_list[idx].handle.GetHashCode() in n_ter_residues
+ is_c_ter = res_list[idx].handle.GetHashCode() in c_ter_residues
frame_residue = rot_constructor.ConstructSidechainFrameResidue(\
res_list[idx].handle, rotamer_ids[idx],
idx, phi_angles[idx], psi_angles[idx],
- res_list[idx].handle.HasProp("n_ter"),
- res_list[idx].handle.HasProp("c_ter"))
+ is_n_ter, is_c_ter)
frame_residues.append(frame_residue)
except:
incomplete_sidechains.append(idx)
@@ -220,27 +228,29 @@ def _AddCysteinFrameResidues(frame_residues, incomplete_sidechains,
incomplete_sidechains.append(idx)
def _GetRotamerGroup(res_handle, rot_id, res_idx, rot_lib, rot_constructor,
- phi, psi, use_frm, bbdep,
- probability_cutoff = 0.98):
+ phi, psi, use_frm, bbdep, probability_cutoff,
+ n_ter_residues, c_ter_residues):
+
+ is_n_ter = res_handle.handle.GetHashCode() in n_ter_residues
+ is_c_ter = res_handle.handle.GetHashCode() in c_ter_residues
if use_frm:
return rot_constructor.ConstructFRMRotamerGroup(res_handle, rot_id,
res_idx, rot_lib,
phi, psi,
- res_handle.HasProp("n_ter"),
- res_handle.HasProp("c_ter"),
+ is_n_ter, is_c_ter,
probability_cutoff)
else:
return rot_constructor.ConstructRRMRotamerGroup(res_handle, rot_id,
res_idx, rot_lib,
phi, psi,
- res_handle.HasProp("n_ter"),
- res_handle.HasProp("c_ter"),
+ is_n_ter, is_c_ter,
probability_cutoff)
def _GetRotamerGroups(res_list, rot_ids, indices, rot_lib, rot_constructor,
- phi_angles, psi_angles, use_frm, bbdep, frame_residues):
+ phi_angles, psi_angles, use_frm, bbdep, frame_residues,
+ n_ter_residues, c_ter_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.
@@ -295,7 +305,8 @@ def _GetRotamerGroups(res_list, rot_ids, indices, rot_lib, rot_constructor,
try:
rot_group = _GetRotamerGroup(r.handle, rot_id, i, rot_lib,
rot_constructor, phi_angles[i],
- psi_angles[i], use_frm, bbdep)
+ psi_angles[i], use_frm, bbdep, 0.98,
+ n_ter_residues, c_ter_residues)
except:
continue
# keep best ones
@@ -308,7 +319,8 @@ def _GetRotamerGroups(res_list, rot_ids, indices, rot_lib, rot_constructor,
def _GetDisulfidBridges(frame_residues, keep_sidechains, cystein_indices,
res_list, rotamer_library, use_frm, bbdep,
- rotamer_constructor, phi_angles, psi_angles):
+ rotamer_constructor, phi_angles, psi_angles,
+ n_ter_residues, c_ter_residues):
'''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,
@@ -346,12 +358,13 @@ def _GetDisulfidBridges(frame_residues, keep_sidechains, cystein_indices,
# residue in the rotamer constructor
phi = _GetPhiAngle(r)
psi = _GetPsiAngle(r)
+ is_n_ter = r.handle.GetHashCode() in n_ter_residues
+ is_c_ter = r.handle.GetHashCode() in c_ter_residues
cys_frame_res = \
rotamer_constructor.ConstructSidechainFrameResidue(r.handle,
sidechain.CYD,
0, phi, psi,
- r.HasProp("n_ter"),
- r.HasProp("c_ter"))
+ is_n_ter, is_c_ter)
for j in range(len(cys_frame_res)):
if cys_frame_res[j].GetName() == "SG":
particle_list = [cys_frame_res[j]]
@@ -372,7 +385,9 @@ def _GetDisulfidBridges(frame_residues, keep_sidechains, cystein_indices,
rotamer_library,
rotamer_constructor,
phi_angles[i],
- psi_angles[i], True, bbdep)
+ psi_angles[i], True, bbdep, 0.98,
+ n_ter_residues, c_ter_residues)
+
rot_group.AddFrameEnergy(frame)
cystein_rot_groups.append(rot_group)
indices_with_rot_groups.append(i)
@@ -513,15 +528,18 @@ def ReconstructSidechains(ent, keep_sidechains=False, build_disulfids=True,
phi_angles, psi_angles = _GetDihedrals(prot.residues)
# set nter and cter (needed in _AddBackboneFrameResidues)
+ n_ter_residues = set()
+ c_ter_residues = set()
for c in prot.chains:
- c.residues[0].SetIntProp("n_ter",1)
- c.residues[-1].SetIntProp("c_ter",1)
+ n_ter_residues.add(c.residues[0].handle.GetHashCode())
+ c_ter_residues.add(c.residues[-1].handle.GetHashCode())
# 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_constructor, phi_angles, psi_angles)
+ rotamer_constructor, phi_angles, psi_angles,
+ n_ter_residues, c_ter_residues)
# add ligands?
if consider_ligands:
@@ -537,7 +555,8 @@ def ReconstructSidechains(ent, keep_sidechains=False, build_disulfids=True,
_AddSidechainFrameResidues(frame_residues, incomplete_sidechains,
keep_sidechains, prot.residues, rotamer_ids,
phi_angles, psi_angles,
- rotamer_constructor, cystein_indices)
+ rotamer_constructor, cystein_indices,
+ n_ter_residues, c_ter_residues)
# update frame_residues, incomplete_sidechains with cysteins (if needed)
if len(cystein_indices) > 0:
# get disulfid bridges and according rotamers
@@ -545,25 +564,29 @@ def ReconstructSidechains(ent, keep_sidechains=False, build_disulfids=True,
_GetDisulfidBridges(frame_residues, keep_sidechains,
cystein_indices, prot.residues,
rotamer_library, use_frm, bbdep,
- rotamer_constructor, phi_angles, psi_angles)
+ rotamer_constructor, phi_angles, psi_angles,
+ n_ter_residues, c_ter_residues)
# update frame_residues, incomplete_sidechains
_AddCysteinFrameResidues(frame_residues, incomplete_sidechains,
keep_sidechains, prot.residues, rotamer_ids,
phi_angles, psi_angles,
rotamer_constructor, cystein_indices,
- disulfid_indices, disulfid_rotamers)
+ disulfid_indices, disulfid_rotamers,
+ n_ter_residues, c_ter_residues)
else:
# update frame_residues, incomplete_sidechains
_AddSidechainFrameResidues(frame_residues, incomplete_sidechains,
keep_sidechains, prot.residues, rotamer_ids,
phi_angles, psi_angles,
- rotamer_constructor)
+ rotamer_constructor, None,
+ n_ter_residues, c_ter_residues)
# 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_constructor, phi_angles,
- psi_angles, use_frm, bbdep, frame_residues)
+ psi_angles, use_frm, bbdep, frame_residues,
+ n_ter_residues, c_ter_residues)
# set up graph and solve to get best rotamers
if use_frm:
--
GitLab