I am working with calculating the rotamer-rotamer pairwise energies over the rotamer library for a fixed backbone. To do this I've essentially reimplemented the body of pack_rotamers_setup(), which ultimately calls RoatmersSet.calculate_energy(). I have a few questions on the topic, some of these are code-based others are to clear up my lack of biology background:
* What is the field RotamerSets.nmoltenres? This appears to be the number of rotamer positions, but is strictly less than the number of residues on the backbone. Shouldn't there be rotamers at every residue?
* RotamerSets.rotamer(id): Is this argument expecting the ID of the rotamer in the full sequence (i.e. over *all* rotamers in the set)?
* After computing the InteractionGraph I get an edge with ig.find_edge(i,j). Looking at the RotamerSet for rotamers i & j I see that they have Mi and Mj rotamers, respectively. However, when I calculate the two-body energy with ig.calculate_two_body_energy(i,j,si,sj) it seems I have nonzero energy for states si>Mi and sj>Mj outside the range of rotamer states.
* Is the two-body energy calculated in InteractionGraph.calculate_two_body_energy() a weighted energy? What energy terms are generally involved? Is it only the Dunbrack energy?
Thanks,
Jason
The number of molten residues for packing isn't necessarily the number of total residues. If you're keeping the rotamer at a position fixed, you really don't need to consider it as a "rotamer", and instead can think of it similar to the fixed backbone.
Yes, the RotamerSets.rotamer() is expecting the rotamer ID for the RotamerSets object, which in general will be different from the ID in a particular rotamer set. You can get the RotamerSets ID -> rotamer set ID translation with the RotamerSets.rotid_on_moltenresidue(rotid) method, and the reverse with the RotamerSets.moltenres_rotid_2_rotid( moltenres, moltenresrotid ) method.
I can't seem to find an InteractionGraph.calculate_two_body_energy() method, so I can't really say much about it. (Other than to question where you found the function.) I should point out that there are a number of different interaction graphs with different behaviors (e.g. how much is precomputed, how information is stored, etc.), so how the function behaves would be likely be highly dependent on which particular subclass of interaction graph you were using.
The Dunbrack energy is actually a one-body energy. Two body energies include fa_atr, fa_rep, fa_sol, fa_pair, the hydrogen bond terms, and the disulfide terms. At any rate, any time you get a single total number as the sum of multiple score terms, that should be the weighted score. (Otherwise, how would you apply the weights afterwards?) It's only when you get a per-scoreterm breakdown of energies that you have to question if the scores are weighted or unweighted.
Sorry, the method is called get_two_body_energies_for_edge(). It is implemented in FixedBBInteractionGraph, but I am using a DensePDInteractionGraph, which must subclass the former. It your point regarding an energy being necessarily weighted if it is a single number makes sense.
Regarding nmoltenres being less than the number of residues: I did not supply a resfile and used standard_packer_task() to create the packer task object. I was assuming this would mean all residues are able to be packed, but maybe it infers certain ones which are fixed? Out of curiosity, why are they called "molten residues"? My background is in compsci, not microbiology, so probably it is terminology I am unfamiliar with.
Jason
Yup, DensePDInteractionGraph is a subclass of FixedBBInteractionGraph (via the intermediate class PrecomputedPairEnergiesInteractionGraph). The get_two_body_energies_for_edge() explicitly has the behavior that it will return zero if it can't find the edge corresponding to the two residues i&j passed to it (even if the values are outside the normal range) - but for speed purposes it doesn't have any bounds safety checking that the rotamer indexes are reasonable. You're deep in the internals of Rosetta, and the code sort of expects that you know what you're doing. If you exceed the bounds of the allowed interactions, you'll get nonsense values - whatever happened to be in a random location of memory. The upshot is that it's up to the calling function to make sure that the indices are valid before calling get_two_body_energies_for_edge(). (I believe the get_num_states_for_node(i) method will tell you what the maximum value of Mi can be for a particular residue/node i.)
By default, the packer task object will allow all the residues to be repackable (and designable). But if you modify the packer task to restrict repacking at some locations (either with a resfile, or with the other packer task manipulation code that is available), then the non-default packer task may specify that certain residues might not be packed. (i.e. in your particular case the number of molten residues is the same as the total number of residues, but that's not necessarily going to be always the case, and you can greatly speed things up by reducing the number of molten residues.)
"Molten" is a bit of a metaphor. I'm not exactly sure where this comes from, but I'm guessing is that it's because they're the residues allowed to move and "flow" like a liquid, as opposed to the fixed residues, which are frozen solid. It's also probably a reference to the optimization technique used in the packer, "Simulated annealing", which is a metaphor based on the metallurgical process of annealing, where controlled heating and cooling is used to partially semi-melt metal objects to strengthen them and fix microdefects.