How to calculate the pairwise rotamer-to-rotamer energy table between mutable residues allowed to change to all amino acid types

[Anonymous]

How can I calculate the self-energies and the table of pairwise energies between rotamers of mutable (molten or changeable) residues in a redesign problem where the mutable residues can change to all amino acid types? I basically am trying to figure out how to get the tables of pairwise rotamer-to-rotamer energies between mutable residues. I want the pairwise rotamer-to-rotamer energy table between mutable residues where the residue can change to all possible amino acid types and thus all possible rotamers of all possible amino acid types. I want to include all rotamers of all possible amino acid types in the mutable residues because I am trying to extend the side chain packing (only rotamers of a single amino acid type allowed) to redesign (all rotamers of all amino acid types allowed). I am modifying code from Jason Pacheco of Brown University (I changed the rotamer_set_for_moltenresidue function to FixbbRotamerSets function to try to include all rotamers of all amino acid types but it is giving the following error – AttributeError: ‘RotamerSets’ object has no attribute ‘FixbbRotamerSets’):

import numpy as np
import optparse # for option sorting
import rosetta
import rosetta.core.pack.task

from rosetta import *
from scipy.sparse import *
from rosetta.core.pack.rotamer_set import *
from rosetta.core.pack import *
from rosetta.core.pack.interaction_graph import *
from rosetta.core.pack.task import * # for using resfiles

init()

## FUNCTIONS ############################################

def get_rotamer_etable( ig, rotsets, moltenresI, moltenresJ ):
“””Computes rotamer-rotamer energy table between two residues.”””

# check edge exists
if not ig.get_edge_exists(moltenresI, moltenresJ):
return None

# get etable
rotsI = rotsets.FixbbRotamerSets()
rotsJ = rotsets.FixbbRotamerSets()
n_rots_I = rotsI.num_rotamers()
n_rots_J = rotsJ.num_rotamers()
E = np.zeros((n_rots_I,n_rots_J))
for s_i in range(1, n_rots_I+1):
for s_j in range(1, n_rots_J+1):
E[s_i-1,s_j-1] = ig.get_two_body_energy_for_edge(moltenresI, moltenresJ, s_i, s_j)
return E

Please kindly guide me.
Thank you so much.

[Anonymous]

The rotsets object you’re getting is an instance of an object somewhere in the RotamerSets hierarchy. The FixbbRotamerSets() call would be to the constructor of the FixbbRotamerSets class to create a new object instance. You don’t need to call it, as rotsets should already be a FixbbRotamerSets object already. (If it’s not of the correct class for some reason, you’ll need to look at how you’re preparing the rotsets object before you pass it to your get_rotamer_etable() function call. You should be able to call the rotamer_set_for_moltenresidue() function directly.

So how do you get all rotamers of all amino acids? That’s all up to how you’re creating the rotsets object. If you’re creating the object with a packertask which allows for design at a given position, the RotamerSet for that position will contain the rotamers for all of the amino acids allowed at that position. If you’ve set the packertask to be repack only at that position, then the rotamerset will only contain rotamers for a single identity.

So what’s a molten residue? To save on storage space, the packer does a bit of index compression. Instead of indexing across the entire pose, it subsets just those residues for which it is considering rotamer changes. So any residue set to NATRO in the packertask won’t be considered a molten residue. Molten residues contain both design and repacking residues. The important thing to remember is that the parameter passed to rotamer_set_for_moltenresidue() is the moltenresID, which is not equivalent to the residue index in the pose. The translation is packertask dependent, but is stored in the RotamerSets object, and can be accessed with the resid_2_moltenres() and moltenres_2_resid() functions.

So you get your RotamerSet object for your desired position (after converting it from absolute pose numbering to the moltenresID). Each of the rotamers in the rotamerset has a particular index for that particular rotamerset object. (The exact numbering depends on the particular packer task you use. Different levels of design and different packer task settings will result in different indexing.) This RotamerSet-specific index is what’s used internally in the interaction graph and the packer to index and label the rotamers and the rotamer interaction energies. (This is the reason for iterating from 1 to rotamerset.num_rotamers() inclusive in the code you provide.) This index is not amino-acid specific. In fact, all possible amino acids allowed at that position by the packer task are represented. For example, indices 1-20 may be for serine, 21 is for alanine, 22 for glycine, and 23-74 for lysine.

Note that this rotamerset index is almost completely unrelated to the value obtained from core.pack.dunbrack.rotamer_from_chi(). Not only does the rotamerset combine multiple residue types into a single ordinal scheme, if you looked you’ll likely notice that there’s multiple rotamerset indexes which correspond to the same amino acid type and core.pack.dunbrack.rotamer_from_chi() vector. The reason for that is the extra rotamer sampling settings like -ex1 and -ex2. These options add subsamples of each of the rotamer bins. so you’ll get several rotamer samples (each with a different rotamerset index) which all fall at different places within the same dunbrack.rotamer_from_chi() bin.

The essential problem you’re running up against is how to take a multidimensional data source (aa identity by chi identity by chi value) and collapse it into a single dimension which you can put as one side of your 2×2 interaction table. There’s no simple way to do it. The packer and the interaction graph choose to do it on an ad hoc basis, specialized for the particular residue and packer task under current consideration. If you’re looking for a more general indexing scheme, I believe you’ll have to come up with one that fits the criteria you need for your application.

[Anonymous]

How can I get the chi angles or rotamer vector of a particular rotamer for a particular molten residue in an instantiation of the RotamerSets class? Is there a method or function with which I can get the chi angles or rotamer vector of a rotamer for any molten residue in the core::pack::rotamer_set::RotamerSets class? I want to get the rotamer vector or chi angles for a rotamer for a molten residue from a public member function in the core::pack::rotamer_set::RotamerSets class. I want to get the rotamer vector or chi angles for a rotamer for a molten residue in the RotamerSets class so that I can tabulate the energy for each of a molten residue’s rotamers defined by its set of chi angles or defined by the rotamer vector.

core::conformation::ResidueCOP core::pack::rotamer_set::RotamerSets::rotamer (uint rotid)
core::conformation::ResidueCOP core::pack::rotamer_set::RotamerSets::rotamer_for_moltenres (uint moltenres_id, uint rotamerid)

In the core::pack::rotamer_set::RotamerSets class, the public member functions rotamer (uint rotid) and rotamer_for_moltenres (uint moltenres_id, uint rotamerid) give the amino acid type and residue number followed by the pdb coordinates of the atoms, like follows (but they don’t give the rotamer vector or chi angles. Is it possible to calculate the rotamer vector or chi angles from the atom coordinates with a PyRosetta function):

Rotamer 11 of residue 1 is: ASP 3:
N : 27.727 28.474 33.178
CA : 26.925 27.439 32.581
C : 27.786 26.384 31.926
O : 28.925 26.623 31.562
CB : 25.9648 28.0374 31.55
CG : 24.9239 28.9573 32.1737
OD1: 25.031 29.2359 33.3445
OD2: 24.0317 29.3726 31.4734
H : 28.2632 29.087 32.5807
HA : 26.333 26.9649 33.3643
1HB : 26.5316 28.603 30.8101
2HB : 25.4493 27.2338 31.0233

The rotamer sets are created with the command:
rotsets = RotamerSets()
(see in code below).

I modified Jason Pacheco’s code to allow the packer task to do design:

pose = pose_from_pdb(in_file)
nRes = pose.n_residue()
scorefxn = get_fa_scorefxn() # create_score_function()

task_pack = TaskFactory.create_packer_task(pose)
parse_resfile(pose, task_pack, resfile)
print ‘Created packer task from resfile’

print ‘Task pack is ‘,task_pack

pack_scorefxn_pose_handshake( pose, scorefxn )
pose.update_residue_neighbors()
scorefxn.setup_for_packing( pose, task_pack.designing_residues(), task_pack.designing_residues() )
packer_neighbor_graph = create_packer_graph( pose, scorefxn, task_pack )
rotsets = RotamerSets()
rotsets.set_task( task_pack )
rotsets.build_rotamers( pose, scorefxn, packer_neighbor_graph )
rotsets.prepare_sets_for_packing( pose, scorefxn )
ig = InteractionGraphFactory.create_interaction_graph( task_pack, rotsets, pose, scorefxn )
print “built”, rotsets.nrotamers(), “rotamers at”, rotsets.nmoltenres(), “positions.”

rotsets.compute_energies( pose, scorefxn, packer_neighbor_graph, ig )

return pose, task_pack, scorefxn, rotsets, ig

The output spits out the packer task:

Task pack is #Packer_Task

resid pack? design? allowed_aas
1 FALSE FALSE
2 FALSE FALSE
3 TRUE TRUE ALA,CYS,ASP,GLU,PHE,GLY,HIS,HIS_D,ILE,LYS,LEU,MET,ASN,PRO,GLN,ARG,SER,THR,VAL,TRP,TYR
4 FALSE FALSE
5 TRUE TRUE ALA,CYS,ASP,GLU,PHE,GLY,HIS,HIS_D,ILE,LYS,LEU,MET,ASN,PRO,GLN,ARG,SER,THR,VAL,TRP,TYR
6 FALSE FALSE
7 TRUE TRUE ALA,CYS,ASP,GLU,PHE,GLY,HIS,HIS_D,ILE,LYS,LEU,MET,ASN,PRO,GLN,ARG,SER,THR,VAL,TRP,TYR
8 FALSE FALSE
9 FALSE FALSE
10 FALSE FALSE
11 FALSE FALSE
12 FALSE FALSE
13 FALSE FALSE
14 FALSE FALSE
15 FALSE FALSE
16 FALSE FALSE
17 FALSE FALSE
18 FALSE FALSE
19 FALSE FALSE
20 TRUE TRUE ALA,CYS,ASP,GLU,PHE,GLY,HIS,HIS_D,ILE,LYS,LEU,MET,ASN,PRO,GLN,ARG,SER,THR,VAL,TRP,TYR
21 FALSE FALSE
22 FALSE FALSE
23 FALSE FALSE
24 FALSE FALSE
25 FALSE FALSE
26 TRUE TRUE ALA,CYS,ASP,GLU,PHE,GLY,HIS,HIS_D,ILE,LYS,LEU,MET,ASN,PRO,GLN,ARG,SER,THR,VAL,TRP,TYR
27 FALSE FALSE
28 FALSE FALSE
29 FALSE FALSE
30 TRUE TRUE ALA,CYS,ASP,GLU,PHE,GLY,HIS,HIS_D,ILE,LYS,LEU,MET,ASN,PRO,GLN,ARG,SER,THR,VAL,TRP,TYR
31 FALSE FALSE
32 FALSE FALSE
33 FALSE FALSE
34 TRUE TRUE ALA,CYS,ASP,GLU,PHE,GLY,HIS,HIS_D,ILE,LYS,LEU,MET,ASN,PRO,GLN,ARG,SER,THR,VAL,TRP,TYR
35 FALSE FALSE
36 FALSE FALSE
37 FALSE FALSE
38 FALSE FALSE
39 TRUE TRUE ALA,CYS,ASP,GLU,PHE,GLY,HIS,HIS_D,ILE,LYS,LEU,MET,ASN,PRO,GLN,ARG,SER,THR,VAL,TRP,TYR
40 FALSE FALSE
41 FALSE FALSE
42 FALSE FALSE
43 TRUE TRUE ALA,CYS,ASP,GLU,PHE,GLY,HIS,HIS_D,ILE,LYS,LEU,MET,ASN,PRO,GLN,ARG,SER,THR,VAL,TRP,TYR
44 FALSE FALSE
45 FALSE FALSE
46 FALSE FALSE
47 FALSE FALSE
48 FALSE FALSE
49 FALSE FALSE
50 FALSE FALSE
51 FALSE FALSE
52 TRUE TRUE ALA,CYS,ASP,GLU,PHE,GLY,HIS,HIS_D,ILE,LYS,LEU,MET,ASN,PRO,GLN,ARG,SER,THR,VAL,TRP,TYR
53 FALSE FALSE
54 TRUE TRUE ALA,CYS,ASP,GLU,PHE,GLY,HIS,HIS_D,ILE,LYS,LEU,MET,ASN,PRO,GLN,ARG,SER,THR,VAL,TRP,TYR
55 FALSE FALSE
56 FALSE FALSE

[Anonymous]

Thank you so much. Yes I can use the chi() method and it works on the rotamer object returned by the rotamer_for_moltenres (uint moltenres_id, uint rotamerid) RotamerSets method. The chi() method is giving angles back to me, however, sometimes they are empty. If the chi angles is empty, is it [0 0 0 0]?

[Anonymous]

Yes this amino acid was ALA. So are GLY and ALA always only one single rotamer? Is that why they don’t have any chi rotamers?

Rotamer 1 of residue 1 is: ALA 3:
N : 27.727 28.474 33.178
CA : 26.925 27.439 32.581
C : 27.786 26.384 31.926
O : 28.925 26.623 31.562
CB : 25.9599 28.0395 31.5693
H : 28.2632 29.087 32.5807
HA : 26.333 26.9649 33.3643
1HB : 25.3554 27.2469 31.1281
2HB : 25.3088 28.7569 32.0689
3HB : 26.5224 28.5449 30.786

Chi angles are []

[Anonymous]

If you have a RotamerSets instance, I’d first recommend obtaining the appropriate RotamerSet instance from it for whichever particular residue you want. While the RotamerSets (plural) class has some convenience functions, internally they just obtain the RotamerSet (singular) instance for the residue of interest, and then call the appropriate function on that RotamerSet. Or if they don’t take a residue/moltenres id, they’ll simply report numbers across *all* residue positions. You’ll get more flexibility if you obtain the RotamerSet instance yourself, and then you can call any number of functions on that RotamerSet

For any particular rotamer in the RotamerSet you can get the Residue object which represents that rotamer with the rotamer(i) member function of the rotamer set. The returned Rotamer object is the same class as the Rotamer object one gets from the Pose, so you can use the same methods on it, for example the nchi() and chi() methods. (BTW, the rotamer (uint rotid) and rotamer_for_moltenres (uint moltenres_id, uint rotamerid) RotamerSets methods are giving you the same Residue objects.)

[Anonymous]

What’s the amino acid identity for those residues which give you empty chi() values? (Print out the results of the name() method on your Residue object.) I’m guessing they’re going to be something like Gly or Ala, which don’t have any chi rotamers. Those should return an empty list/vector1 for their chis, not [0 0 0 0] (which would imply four chi angles, each set to zero.)

[Anonymous]

Well, glycine doesn’t have a sidechain, so there are no sidechain rotamers. While alanine has a sidechain, it only has a single heavy atom in the sidechain, so there are no heavy atom chis. Also, even though there’s technically rotational freedom in the methyl group, Rosetta typically does not model that rotational freedom, as the methyl group is rather symmetric and for most purposes that rotation doesn’t matter.

So the upshot is that for glycine and alanine there isn’t any conformational flexibility in the sidechain which Rosetta models, so it doesn’t have any chis, and only has a single “rotamer”.

[Author]

Posts