Source code for promod3.sidechain.reconstruct_sidechains

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)