1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
|
from ost import io,mol
from promod3 import sidechain
# load a protein
prot = io.LoadPDB('data/1CRN.pdb')
# load rotamer library
library = sidechain.LoadBBDepLib()
# we need a rotamer constructor to create any rotamers or
# frame residues
rot_constructor = sidechain.SCWRL4RotamerConstructor(False)
# create new entity from protein only containing the amino acids
prot = mol.CreateEntityFromView(prot.Select("peptide=true"), True)
# gather some data, the rotamer ids and backbone torsion angles
torsion_angles = list()
rotamer_ids = list()
for r in prot.residues:
rotamer_ids.append(sidechain.TLCToRotID(r.GetName()))
phi_handle = r.GetPhiTorsion()
psi_handle = r.GetPsiTorsion()
# set typical default values for an alpha helix...
phi = -1.0472
psi = -0.7854
if phi_handle.IsValid():
phi = phi_handle.GetAngle()
if psi_handle.IsValid():
psi = psi_handle.GetAngle()
torsion_angles.append((phi,psi))
# first build a frame representing the rigid parts including
# cystein sidechains
frame_residues = list()
for i,r in enumerate(prot.residues):
frame_residue = rot_constructor.ConstructBackboneFrameResidue(
r, rotamer_ids[i], i,
torsion_angles[i][0],
torsion_angles[i][1],
i == 0,
i == (len(rotamer_ids)-1))
frame_residues.append(frame_residue)
frame = sidechain.Frame(frame_residues)
# let's build up rotamer groups
rotamer_groups = list()
aa_with_rotamers = list()
for i,r in enumerate(prot.residues):
if r.GetName() == "ALA" or r.GetName() == "GLY":
continue
rot_group = rot_constructor.ConstructFRMRotamerGroup(
r, rotamer_ids[i], i, library,
torsion_angles[i][0], torsion_angles[i][1],
i == 0, i == (len(rotamer_ids)-1), 0.98)
# calculate pairwise energies towards the rigid frame
# the values will be stored and used by the solving algorithm
# later on when considering self energies
rot_group.SetFrameEnergy(frame)
# remove super unlikely rotamer in rotamer group
# e.g. those who clash with the frame
rot_group.ApplySelfEnergyThresh()
# finally add it and keep track of the indices
rotamer_groups.append(rot_group)
aa_with_rotamers.append(i)
# buildup a graph given the rotamer groups
graph = sidechain.RotamerGraph.CreateFromFRMList(rotamer_groups)
# and get a solution out of it using a minimal width tree approach
solution = graph.TreeSolve()[0]
# let's finally apply the solution to the residues
for (i, rot_group, sol) in zip(aa_with_rotamers,
rotamer_groups,
solution):
rot_group[sol].ApplyOnResidue(prot.residues[i])
io.SavePDB(prot, "example_reconstruction.pdb")
|
