from promod3 import sidechain
from ost import geom
[docs]def Reconstruct(ent, keep_sidechains = False, build_disulfids = True,
rotamer_model = "frm", consider_hbonds = True,
rotamer_library = None,add_all_polar_hydrogens=False):
name_id_mapper = dict()
name_id_mapper["ARG"] = sidechain.ARG
name_id_mapper["ASN"] = sidechain.ASN
name_id_mapper["ASP"] = sidechain.ASP
name_id_mapper["GLN"] = sidechain.GLN
name_id_mapper["GLU"] = sidechain.GLU
name_id_mapper["LYS"] = sidechain.LYS
name_id_mapper["SER"] = sidechain.SER
name_id_mapper["CYS"] = sidechain.CYS
name_id_mapper["MET"] = sidechain.MET
name_id_mapper["TRP"] = sidechain.TRP
name_id_mapper["TYR"] = sidechain.TYR
name_id_mapper["THR"] = sidechain.THR
name_id_mapper["VAL"] = sidechain.VAL
name_id_mapper["ILE"] = sidechain.ILE
name_id_mapper["LEU"] = sidechain.LEU
name_id_mapper["PRO"] = sidechain.PRO
name_id_mapper["HSD"] = sidechain.HSD
name_id_mapper["HSE"] = sidechain.HSE
name_id_mapper["HIS"] = sidechain.HIS
name_id_mapper["PHE"] = sidechain.PHE
name_id_mapper["ALA"] = sidechain.ALA
name_id_mapper["GLY"] = sidechain.GLY
name_id_mapper["MSE"] = sidechain.MET
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()
residues_with_rotamer_group = list()
rotamer_groups = list()
prot = ent.Select("peptide=true")
incomplete_sidechains = list()
#extract dihedral angles
phi_angles = [0.0] * len(prot.residues)
psi_angles = [0.0] * len(prot.residues)
for i,r in enumerate(prot.residues):
phi = -1.0472
psi = -0.7854
tor = r.GetPhiTorsion()
if tor.IsValid():
phi = tor.GetAngle()
tor = r.GetPsiTorsion()
if tor.IsValid():
psi = tor.GetAngle()
phi_angles[i] = phi
psi_angles[i] = psi
#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):
try:
rot_id = name_id_mapper[r.GetName()]
except:
rot_id = sidechain.ALA
try:
frame_residue = sidechain.ConstructBackboneFrameResidue(r.GetHandle(),rot_id,
i,rotamer_settings,
phi_angles[i],
r.HasProp("n_ter"),
r.HasProp("c_ter"))
frame_residues.append(frame_residue)
except:
continue
if keep_sidechains:
for i,r in enumerate(prot.residues):
try:
rot_id = name_id_mapper[r.GetName()]
except:
continue
if rot_id == sidechain.CYS:
continue #cysteins will be handled seperately
if rot_id == sidechain.ALA or rot_id == sidechain.GLY:
continue
try:
frame_residue = sidechain.ConstructSidechainFrameResidue(r.GetHandle(),rot_id,
i,rotamer_settings)
frame_residues.append(frame_residue)
except:
incomplete_sidechains.append(i)
else:
for i,r in enumerate(prot.residues):
try:
rot_id = name_id_mapper[r.GetName()]
except:
continue
if rot_id not in [sidechain.GLY, sidechain.ALA]:
incomplete_sidechains.append(i)
#look for cysteins and evaluate potential disulfid bonds
final_disulfid_indices = list()
if build_disulfids:
cystein_info = list()
for i,r in enumerate(prot.residues):
try:
rot_id = name_id_mapper[r.GetName()]
except:
continue
if rot_id == sidechain.CYS:
ca = r.FindAtom("CA")
cb = r.FindAtom("CB")
if ca.IsValid() and cb.IsValid():
cystein_info.append((ca.GetPos(),cb.GetPos(),i))
frame = sidechain.Frame(frame_residues)
for i in range(len(cystein_info)):
if cystein_info[i][2] in final_disulfid_indices:
continue
complete_i = keep_sidechains and prot.residues[cystein_info[i][2]].FindAtom("SG").IsValid()
i_ca_pos = cystein_info[i][0]
i_cb_pos = cystein_info[i][1]
i_index = cystein_info[i][2]
for j in range(i+1,len(cystein_info)):
if cystein_info[j][2] in final_disulfid_indices:
continue
complete_j = keep_sidechains and prot.residues[cystein_info[j][2]].FindAtom("SG").IsValid()
if complete_i and complete_j:
continue #they're already there
if geom.Distance(i_ca_pos,cystein_info[j][0]) < 8:
j_ca_pos = cystein_info[j][0]
j_cb_pos = cystein_info[j][1]
j_index = cystein_info[j][2]
if use_frm:
rot_group_one = sidechain.ConstructFRMRotamerGroup(prot.residues[cystein_info[i][2]].GetHandle(),
sidechain.CYD,cystein_info[i][2],
rotamer_library,rotamer_settings,
phi_angles[i_index],psi_angles[i_index])
rot_group_two = sidechain.ConstructFRMRotamerGroup(prot.residues[cystein_info[j][2]].GetHandle(),
sidechain.CYD,cystein_info[j][2],
rotamer_library,rotamer_settings,
phi_angles[j_index],psi_angles[j_index])
else:
rot_group_one = sidechain.ConstructRRMRotamerGroup(prot.residues[cystein_info[i][2]].GetHandle(),
sidechain.CYD,cystein_info[i][2],
rotamer_library,rotamer_settings,
phi_angles[i_index],psi_angles[i_index])
rot_group_two = sidechain.ConstructRRMRotamerGroup(prot.residues[cystein_info[j][2]].GetHandle(),
sidechain.CYD,cystein_info[j][2],
rotamer_library,rotamer_settings,
phi_angles[j_index],psi_angles[j_index])
frame.AddFrameEnergy(rot_group_one)
frame.AddFrameEnergy(rot_group_two)
min_index_k = -1
min_index_l = -1
min_score = 10000000
for k in range(len(rot_group_one)):
for l in range(len(rot_group_two)):
score = sidechain.DisulfidScore(rot_group_one[k],rot_group_two[l],
i_ca_pos,i_cb_pos,j_ca_pos,j_cb_pos)
if score < min_score:
min_index_k = k
min_index_l = l
min_score = score
if min_score < 45:
final_disulfid_indices.append(i_index)
final_disulfid_indices.append(j_index)
#create new frame residue...
particle_list = list()
particle_list.append(rot_group_one[min_index_k][0])
particle_list.append(rot_group_two[min_index_l][0])
new_frame_residue = sidechain.FrameResidue(particle_list,100000)
frame_residues.append(new_frame_residue)
#set the position in the proteins residues
rot_group_one[min_index_k].ApplyOnResidue(prot.residues[i_index].GetHandle())
rot_group_two[min_index_l].ApplyOnResidue(prot.residues[j_index].GetHandle())
#we still have to add the cysteins, that are complete, even though
#they're not participating in a disulfid bond
if keep_sidechains:
for i,r in enumerate(prot.residues):
try:
rot_id = name_id_mapper[r.GetName()]
except:
continue
if rot_id != sidechain.CYS:
continue #all other sidechains are already handled
if i in final_disulfid_indices:
continue #already added to the frame residues
try:
frame_residue = sidechain.ConstructSidechainFrameResidue(r.GetHandle(),rot_id,
i,rotamer_settings)
frame_residues.append(frame_residue)
except:
incomplete_sidechains.append(i)
else:
#The disulfid cysteins must be removed from incomplete sidechains
for i in final_disulfid_indices:
try:
incomplete_sidechains.remove(i)
except:
pass #it's not there anyway, so we don't have to remove it
elif keep_sidechains:
for i,r in enumerate(prot.residues):
try:
rot_id = name_id_mapper[r.GetName()]
except:
continue
if rot_id != CYS:
continue #all other sidechains are already handled
try:
frame_residue = sidechain.ConstructSidechainFrameResidue(r.GetHandle(),rot_id,
i,rotamer_settings)
frame_residues.append(frame_residue)
except:
incomplete_sidechains.append(i)
#finally build complete frame
frame = sidechain.Frame(frame_residues)
#build rotamers
for i in incomplete_sidechains:
r = prot.residues[i]
try:
rot_id = name_id_mapper[r.GetName()]
except:
continue
if(rot_id == sidechain.ALA or rot_id == sidechain.GLY):
continue
if rot_id == sidechain.CYS:
if i in final_disulfid_indices:
continue
rot_id = sidechain.CYH
if rot_id == sidechain.PRO:
tor = r.GetOmegaTorsion()
if tor.IsValid():
if abs(tor.GetAngle()) < 1.57:
rot_id = sidechain.CPR
else:
rot_id = sidechain.TPR
try:
if use_frm:
rot_group = sidechain.ConstructFRMRotamerGroup(r.GetHandle(),rot_id,i,
rotamer_library,rotamer_settings,
phi_angles[i],psi_angles[i])
else:
rot_group = sidechain.ConstructRRMRotamerGroup(r.GetHandle(),rot_id,i,
rotamer_library,rotamer_settings,
phi_angles[i],psi_angles[i])
except:
continue
rot_group.CalculateInternalEnergies()
frame.SetFrameEnergy(rot_group)
rot_group.ApplySelfEnergyThresh()
rotamer_groups.append(rot_group)
residues_with_rotamer_group.append(i)
if use_frm:
graph = sidechain.Graph.CreateFromFRMList(rotamer_groups)
else:
graph = sidechain.Graph.CreateFromRRMList(rotamer_groups)
solution = graph.Solve(100000000,0.02)
for i,rot_group,sol in zip(residues_with_rotamer_group,rotamer_groups,solution):
try:
rot_group[sol].ApplyOnResidue(prot.residues[i].GetHandle(),consider_hydrogens=True)
except:
print "there is a backbone atom missing... ",prot.residues[i].GetQualifiedName()
if add_all_polar_hydrogens:
for r in frame_residues:
i=r.GetResidueIndex()
if i==100000:continue#These are cysteins that were already handled
rh=prot.residues[i].GetHandle()
r.ApplyOnResidue(rh,True)