sidechain_steps.py

from ost import io,mol
from promod3 import sidechain

# load a protein
prot = io.LoadPDB('data/1CRN.pdb')
# load rotamer library
library = sidechain.LoadDunbrackLib()
# we need a rotamer constructor to create any rotamers or 
# frame residues
rot_constructor = sidechain.SCWRLRotamerConstructor()

# 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], 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])

    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")