You are here

Rotamers one- and two-body energies

4 posts / 0 new
Last post
Rotamers one- and two-body energies

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!

Post Situation: 
Thu, 2024-04-04 00:48

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:

  • a pose
  • a score function
  • a packer task -- create as with any PyRosetta packing approach - note that if you have a TaskFactory, you'll need to call `packer_factory.create_task_and_apply_taskoperations(pose)` to get the packer task.
  • a RotamerSets object


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

Thu, 2024-04-04 09:44


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:




Any idea what may be wrong?
Again, thank you so much!!

Mon, 2024-04-08 13:04


Using rmoretti's inputs (thank you!) and taking lots of inspiration from this code, I was able to get the 1- and 2-body rotamer energy matrices for a given pose (built from an input pdb file). It's been a while, but I thought it was a good idea to post the code here, maybe it's useful for someone. 

from import tqdm
import sys
import pyrosetta
from pyrosetta.rosetta.core.pack import create_packer_graph, pack_rotamers_setup
from pyrosetta.rosetta.core.pack.task import TaskFactory
from pyrosetta.rosetta.core.pack.rotamer_set import RotamerSets
from pyrosetta.rosetta.core.pack.interaction_graph import InteractionGraphFactory

def compute_1_and_2_body_rotamer_energies(
        pdb_file_path, init_options, 
        only_save_E2_above_cutoff=True, cutoff_E2=0.001
    computes the 1- and 2-body rotamer energies for a pose specified by the pdb file.
    returns and saves the data in the form of dictionaries.

    - pdb_file_path (str): path of pdb file whose energies will be calculated;
    - init_options (str or None): commandline rosetta initialization options;
    - only_save_E2_above_cutoff (bool). if False, save everything.
    - cutoff_E2 (float): cutoff 2-body energy value to save.
    # ======================================================
    # setup and ig construction
    # ======================================================

    pyrosetta.init(f"-mute all {init_options}")

    pose = pyrosetta.pose_from_pdb(pdb_file_path)

    print("\nBuilding InteractionGraph...")

    sfxn = pyrosetta.get_fa_scorefxn()

    task_design = TaskFactory.create_packer_task(pose)
    if init_options:

    # this is important to get a PDInteractionGraph object -- it's the only one i could make work
    # without this, the ig returned is a LinearMemoryInteractionGraph, which doesn't have 
    # the method `get_two_body_energy` implemented! (tested for PyRosetta-4 2024)


    png = create_packer_graph(pose, sfxn, task_design)

    rotsets = RotamerSets()

    ig = pack_rotamers_setup(pose, sfxn, task_design, rotsets)
    ig = InteractionGraphFactory.create_and_initialize_two_body_interaction_graph(
        task_design, rotsets, pose, sfxn, png


    # ======================================================
    # dicts for saving data
    # ======================================================

    # data format: `{pos: n_rotamers}`
    dict_n_rots_at_position = {}

    # data format: `{pos : {rot_index: (name1, [chi_angles])}}`
    # with this, it's possible to retrieve the actual side chain identity and chi angles
    # for a given sequence of rotamers -- for getting final answer, for instance
    dict_pos_rotamer_angles_id = {}

    # data format: `{pos: {rot_index: E1}}`
    # this stores the one-body energies
    dict_E1 = {}

    # data format: `{(pos_i, pos_j): {(rot_index_i, rot_index_j): E2}}`
    # this stores the two-body energies
    dict_E2 = {}

    # ======================================================
    # one-body energies
    # ======================================================

    print("\nComputing 1-body energies...")

    for pos_i in tqdm(range(1, ig.get_num_nodes()+1), desc="position"):

        n_rots_i = rotsets.rotamer_set_for_moltenresidue(pos_i).num_rotamers()
        dict_n_rots_at_position[pos_i] = n_rots_i

        dict_E1[pos_i] = {}
        dict_pos_rotamer_angles_id[pos_i] = {}

        for rot_index_i in range(1, n_rots_i+1):

            # one-body energy
            E1 = ig.get_one_body_energy_for_node_state(pos_i, rot_index_i)
            dict_E1[pos_i][rot_index_i] = E1

            # rotamer res id and chi angles
            rotamer = rotsets.rotamer_set_for_moltenresidue(pos_i).rotamer(rot_index_i)
            name1 = rotamer.name1()
            chi_angles = list(rotamer.chi())
            dict_pos_rotamer_angles_id[pos_i][rot_index_i] = (name1, chi_angles)


    # ======================================================
    # two-body energies
    # ======================================================

    print("\nComputing 2-body energies...")

    for pos_i in tqdm(range(1,ig.get_num_nodes()+1), desc="all pairs for position"):
        for pos_j in range(pos_i+1,ig.get_num_nodes()+1):

            # calculate energy and save energies only if edge exists
            if (ig.get_edge_exists(pos_i, pos_j)):

                n_rots_i = rotsets.rotamer_set_for_moltenresidue(pos_i).num_rotamers()
                n_rots_j = rotsets.rotamer_set_for_moltenresidue(pos_j).num_rotamers()
                dict_E2[(pos_i, pos_j)] = {}

                for r_i in range(1, n_rots_i+1):
                    for r_j in range(1, n_rots_j+1):

                        # two-body energy
                        E2 = ig.get_two_body_energy_for_edge(pos_i, pos_j, r_i, r_j)

                        # won't store small energies
                        if only_save_E2_above_cutoff and np.abs(E2) <= cutoff_E2:

                        dict_E2[(pos_i, pos_j)][(r_i, r_j)] = E2


    # saving pickle with dicts
    if output_file_path:

        init_options_str = init_options.replace(" ", "_") if init_options else "None"

        pdb_code = pdb_file_path.split("/")[-1].strip(".pdb")
        with open(os.path.join(output_file_path, f'dicts_energy_matrices_{pdb_code}_cutoff_E2={cutoff_E2}_init_options={init_options_str}.pkl'), 'wb') as f:
                (dict_E1, dict_E2, dict_n_rots_at_position, dict_pos_rotamer_angles_id), 
        print("\nPickle with dictionaries saved successfully!")

    return dict_E1, dict_E2, dict_n_rots_at_position, dict_pos_rotamer_angles_id

Thanks for the help!

Wed, 2024-05-29 05:39