Hello! I'm quite new to (Py)Rosetta, so my question may have a very obvious answer, but I tried to find it in many places but couldn't find it, so I figured this was a good place to ask.
My goal is to implement "from scratch" some methods to solve the side-chain optimization problem, using a pairwise-decomposable energy function. In order to do so, I need to have the precomputed one- and two-body rotamer energies . Now, although I managed to calculate these values, I believe I didn't do it in the most efficient way (if it's relevant, what I did was: getting 3 different (equally-spaced) angle values for each chi, for each residue type; for each position in an input pose, I mutated the respective residue to each one of the rotamers; I then calculated the one-body residue energy, and did similarly for pairwise mutated residues energies).
I know this approach is not the best one because: a) it doesn't use Dunbrack's library; b) It's very slow.
In my understanding, Rosetta already has these "rotamer energy matrices" precomputed somewhere (or it precomputes them easily from Dunbrack's library, given a pose), so my question is simply: where is it, and how can I use it?
I think that the information I need is in the pyrosetta/database/rotamer folder, but I really don't know how to use the files there. I tried to look in (Py)Rosetta's source code to see how the files are used, but didn't have much luck there.
If anyone can help me with any information/direction, I really appreciate it.
Thank you so much!
You're right that the Rosetta packer already sets up these sorts of one-body and two-body energy tables.
The datastructure you're looking for is the InteractionGraph. In particular, to create one you're looking at the rosetta.core.pack.pack_rotamers_setup() function, which should return the InteractionGraph object.
This takes the following:
There's no initialization needed for the RotamerSets object. (The pack_rotamers_setup will initialize it based on the task) You just have to instantiate either rosetta.core.pack.rotamer_set.RotamerSets or rosetta.core.pack.rotamer_set.rosetta.core.pack.symmetry.SymmetricRotamerSets depending on whether your pose is symmetric or not.
Note you may want to specify a pre-computed interaction graph, with task.or_precompute_ig(True) -- This will likely give you a PDInteractionGraph.
Using an interaction graph is a bit of a trick. It doesn't really work with residues & conformers, but rather node and states. And nodes only go into the graph if they're packable/designable (so node numbers don't . To some extent you have to work back and forth between the RotamerSets object and the InteractionGraph object.
You can get the (weighted summed) one-body energies with ig.get_one_body_energy_for_node_state(nnode, nstate). (Use ig.get_num_nodes() and ig.get_num_states_for_node(nnode) for the number of nodes and states) -- Note that we don't really store the one-body energies on a per-term basis.
For two-body energies, you need to first get the Edge corresponding to the interaction. Use ig.get_edge_exists(nnode1, nnode2) first to check if the edge exists (no edge, no pairwise energy), then ig.find_edge(nnode1, nnode2) to get the interaction edge and edge.get_two_body_energy(nstate1,nstate2) to get the weighted energies of the pairwise interactions. -- IMPORTANT!!! There really isn't any bounds checking here, so you have to be careful that you don't exceed the state size of each of the nodes. Also important is that the order of the node numbers is important: always pass the smaller-indexed node first (and the state for the smaller indexed node first).
For working back and forth with the RotamerSet, I believe you can use the rotsets.moltenres_2_resid(nnode) and the rotsets.moltenres_rotid_2_rotid( nnode, nstate ) to convert from node numbers to residue numbers and state numbers. Or alternatively, you can just get the corresponding rotamer (as a Residue object) with the rotsets.rotamer_for_moltenres( nnode, nstate )
CAVEAT -- I've only played around with this briefly. I haven't done extensive testing here, so there may be some inaccuracies, particularly with indicies, offsets and conversons. I'd highly recommend testing things out with a small test case first. In particular, I'd recommend double checking the results given with sfxn.eval_intrares_energy() and then the sum of the sfxn.eval_cd_2b() and sfxn.eval_ci_2b() energies. -- I'd also recommend double checking the score of a re-constitued pose of your selected rotamers ... or at least, the particular one-body and two-body terms from the Energies object.
Additional caveat -- not all energies are considered in the packing problem. Most are, but there are certain long-range (most notably disulfide potentials) and whole structure energies which are not necessarily considered during packing, and wouldn't show up in the analysis above. (Also, if you aren't changing a position, it won't necessarily show up in energy calculation, depending on what else is going on.)
Hello!
Thank you so much for the really detailed answer, I really appreciate it! I learned a lot about the InteractionGraph object by playing around, so thank you!
Unfortunately, though, I wasn't able to get the actual desired interaction graph, because whenever I try to run the function rosetta.core.pack.pack_rotamers_setup(), my Python crashes... I tried to do this in different ways, but the simpler one is as follows:
```
pyrosetta.rosetta.core.pack.pack_rotamers_setup(
pose,
sfxn,
pyrosetta.rosetta.core.pack.task.PackerTask_(),
pyrosetta.rosetta.core.pack.rotamer_set.RotamerSets(),
pyrosetta.rosetta.core.pack.interaction_graph.AnnealableGraphBase()
)
```
Any idea what may be wrong?
Again, thank you so much!!