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

Member Site Forums PyRosetta PyRosetta – General How to calculate the pairwise rotamer-to-rotamer energy table between mutable residues allowed to change to all amino acid types

Viewing 4 reply threads
  • Author
    • #1808

        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


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

      • #9708

          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.

        • #9724

            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 )
            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
            4 FALSE FALSE
            6 FALSE FALSE
            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
            21 FALSE FALSE
            22 FALSE FALSE
            23 FALSE FALSE
            24 FALSE FALSE
            25 FALSE FALSE
            27 FALSE FALSE
            28 FALSE FALSE
            29 FALSE FALSE
            31 FALSE FALSE
            32 FALSE FALSE
            33 FALSE FALSE
            35 FALSE FALSE
            36 FALSE FALSE
            37 FALSE FALSE
            38 FALSE FALSE
            40 FALSE FALSE
            41 FALSE FALSE
            42 FALSE FALSE
            44 FALSE FALSE
            45 FALSE FALSE
            46 FALSE FALSE
            47 FALSE FALSE
            48 FALSE FALSE
            49 FALSE FALSE
            50 FALSE FALSE
            51 FALSE FALSE
            53 FALSE FALSE
            55 FALSE FALSE
            56 FALSE FALSE

          • #9729

              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]?

            • #9735

                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 []

              • #9728

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

                • #9732

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

                  • #9736

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

                  Viewing 4 reply threads
                  • You must be logged in to reply to this topic.