Pack vs. Minimize vs. Relax post-mutation

[sean-workman]

Having never had the opportunity (yet) to attend a (py)Rosetta workshop, I've been working my way through some little projects myself but I'm not certain my current workflow makes sense for what I am trying to do.

I am generating a clean version of my input PDB with cleanATOM(), then running a FastRelax() to get my input ready to go for mutations. I am separately making point mutations to all non-native residues at a number of locations with MutateResidue() and these are at an interface between two proteins. The plan is to use Surfaces to look at the perturbations to the interfaces after, but I'm unsure if I should just try packing rotamers, minimizing, or doing a FastRelax() on the mutated structure files before analysis. I had thought FastRelax() but that seems maybe a bit overkill for dealing with the effects of the point mutations.

Any guidance or insight would be greatly appreciated.

P.S. I am also wondering how to directly reproduce the behavior of the default command line minimize in pyRosetta.

[roccomoretti]

How much pre-optimization you do depends a bit on your system, where it's coming from and what sort of score function it was (or was not) optimized into before.

The basic philosophy (especially with design) is that the Rosetta energy function is rugged, and if you have an input structure which has slight (sub-angstrom) movements off from where Rosetta thinks it should be, it could potentially score poorly. This means that Rosetta could put in a mutation (or make a massive structural change) not because the native sequence is bad, but just because the atoms are slightly off and the mutation just happens to fit that "bad" atom position better than the native one does. Pre-optimizing the structure allows you to relieve those slight deficiencies, which should hopefully allow Rosetta to be less aggressive about changing things.

You have to be slightly careful, though, as you can over-optimize things. If you do too much pre-optimization, you can dig the native structure deep into an energy well, which means that any change you make is going to be worse in comparison. The trick is to optimize things just enough to avoid issues with false positives from under-optimizing, while not drilling down too much on the input structure.

Different approaches thus vary based on where you want to strike that balance. Just packing rotamers will keep the backbone fixed. This might work, or small deviations in the backbone might cause different residues to be favored (e.g. the rama terms are residue-identity dependent.). Minimization generally moves everything, but doesn't necessarily sample different rotamers and has the potential to (it also has the potential to move the backbone significantly). Relax combines the two, and adds a more thorough sampling scheme. This has the benefits and drawbacks of both approaches. Which one you choose depends a bit on what you're hoping to optimized, and what sort of sampling you're doing downstream.

One approach is to use Relax, but constraint the ability of relax to move the structure too much. The relax application can take the -constrain_relax_to_start_coords and the -relax:ramp_constraints false flags, which will apply somewhat strong restraints to the backbone coordinates, keeping them from moving too much during relax (generally sub-angstrom movements). (There's also the -coord_constrain_sidechains option if you also want to add the restraints to the sidechain positions. ) This has been shown to greatly improve the energetics of the protein without changing the structure much (see https://doi.org/10.1371/journal.pone.0059004 and https://onlinelibrary.wiley.com/doi/10.1002/pro.4807#pro4807-sec-0801), and would probably be what I'd recommend if you don't have evidence/inclinations that another preparation method might be better.

To do the equivalent of the minimize command line program in PyRosetta, take a look at the rosetta.protocols.minimization_packing.MinMover mover.

[sean-workman]

Hi Rocco,

Thank you so much for your detailed and explanatory answer. I will definitely give your suggestions a shot!

My concern in this case largely stem from the fact that I'm making point mutations at an interface where I know the that certain side chains are interacting. If I mutate a Thr to an Arg, for example, beyond the change in the charge, I'm imagining that I am also introducing some pretty huge steric clashes.

What I tried while waiting for input (and your answer was very fast - thank you again!), was using the MinMover() to do a local minimization around the location of the mutation. I appreciate this wouldn't be as thorough as running a full relaxation protocol, but it seemed an okay compromise in terms of computational load.

Does something like this look reasonable?

pyrosetta.init('-ex1 -ex2 -use_input_sc -ignore_unrecognized_res -ignore_zero_occupancy false -multithreading:total_threads 25' )

# Dictionary mapping one-letter to three-letter amino acid codes
AA_MAP = {
    'A': 'ALA', 'C': 'CYS', 'D': 'ASP', 'E': 'GLU', 'F': 'PHE', 'G': 'GLY',
    'H': 'HIS', 'I': 'ILE', 'K': 'LYS', 'L': 'LEU', 'M': 'MET', 'N': 'ASN',
    'P': 'PRO', 'Q': 'GLN', 'R': 'ARG', 'S': 'SER', 'T': 'THR', 'V': 'VAL',
    'W': 'TRP', 'Y': 'TYR'
}

# Get a score function to use for updating neighbors
scorefxn = pyrosetta.get_fa_scorefxn()
input_pose = pyrosetta.pose_from_pdb('####.pdb')

# Initialize the neighbor graph for the input pose
input_pose.update_residue_neighbors()

# Set up FastRelax
relax = FastRelax()
relax.set_scorefxn(scorefxn)
input_relaxed = pyrosetta.Pose()
input_relaxed.assign(input_pose)
relax.apply(input_relaxed)

input_relaxed.dump_pdb("####_relaxed.pdb")

with open("to_mutate.txt") as f:
    targets = list(set([line.strip() for line in f.readlines()]))
    residues = [res[:-1] for res in targets]
    indices = [int(re.findall(r"\d+", res)[0]) for res in targets]
    paired = list(zip(residues, indices))

# Define minimization radius (in Angstroms)
radius = 10.0

# Make mutations, minimize, save PDBs
for res, idx in paired:
    ridx = input_relaxed.pdb_info().pdb2pose("A", idx)
    native_aa = res[:3].upper()

    for aa in AA_MAP.values():
        if aa != native_aa:

            pose_copy = pyrosetta.Pose()
            pose_copy.assign(input_relaxed)
        
            # Perform mutation
            mutate = MutateResidue()
            mutate.set_target(ridx)
            mutate.set_res_name(aa)
            mutate.apply(pose_copy)
        
            # Save the non-minimized structure if needed
            pose_copy.dump_pdb(f"pdbs/{res}_{aa}.pdb")
        
            # Update the neighbor graph after mutation
            pose_copy.update_residue_neighbors()
        
            # Set up movemap for local minimization
            movemap = MoveMap()
            movemap.set_bb(False)  # Don't move backbone by default
            movemap.set_chi(False)  # Don't move side chains by default
        
            # Find residues within radius of mutation
            nearby_residues = pyrosetta.rosetta.core.select.residue_selector.NeighborhoodResidueSelector()
            nearby_residues.set_focus(str(ridx))
            nearby_residues.set_distance(radius)
            nearby_residues.apply(pose_copy)
        
            # Enable movement for residues in neighborhood
            for i in range(1, pose_copy.total_residue() + 1):
                if nearby_residues.apply(pose_copy)[i]:
                    movemap.set_bb(i, True)
                    movemap.set_chi(i, True)
        
            # Set up and apply minimization
            min_mover = MinMover(movemap, scorefxn, 'lbfgs_armijo_nonmonotone', 0.001, True)
            min_mover.apply(pose_copy)
        
            # Save the minimized structure
            pose_copy.dump_pdb(f"pdbs/{res}_{aa}_minimized.pdb")

[roccomoretti]

Sorry, I realize much of my response to you was based on *pre-*optimization considerations (versus what you were actually asking, which is the optimization after mutation.) -- You said so in your question, but I was reading it through a different lens.

That said, optimization before and after the mutation should be considered together. You'll want to coordinate the two somewhat. That is, you'll want to make sure that your pre-optimization approach matches your post-optimization, such that you (A) don't spend an inordinate amount of time post-optimizing for multiple mutations when you can do it just once with a pre-optimizaiton and (B) such that you don't unfairly penalize the native structure due to a bad starting geometry. (Though the latter can be compensated somewhat by also doing a "null" mutation -- run the same protocol, but include a wt -> wt "mutation" with all the same optimizations as a reference.)

Regarding what to do for post-mutations, it depends a bit on how aggressive you want to get. What we have observed is that the more you deviate from the experimental structure, typically the more chance of error creeping in. As such, you may want to reduce how much movement you have in the structure post-mutation to minimize the "disruptive-ness" of the mutation.

That is, if you pack just the residue, you'll avoid mutations which disrupt the structure of other residues, but you'll miss residues which would be accommodated if the other sidechains moved slightly. If you only repack sidechains, you'll screen against mutations which could mess up backbone structure, but you'll miss mutations which could be accommodated with slight backbone movement. And so on. How much you attempt to minimize things is thus mostly a call on the sort of mutations you want to accommodate experimentally, and also how much trust you have in Rosetta not to mess up the structure and give you false positives.

Your protocol looks reasonable. The one thing to be aware of is lever-arm effects in backbone minimization. That is, if you minimize a backbone torsion, everything downstream of it will move in XYZ space, even if it's not enabled in the MoveMap. Often that's not a problem, as the Rosetta energy function will work to keep things together and folded. But there's a chance a mutation+minimization combo will result in the protein opening up. -- It may be worth adding checks to see how much the structure has moved from the original structure, and flagging/discarding those runs where things change too much.

P.S. As mentioned, I'd skip the if aa != native_aa check -- run the null mutation to get a good baseline energy for the native at each position. (Particularly important as you're changing what's minimized at each position.)

I'd also recommend pose_copy = input_relaxed.clone(). And you may just want to cache the results of the nearby_residues.apply(pose_copy), as the way you're doing it re-computes the neighbors for each position.

[sean-workman]

Ah, I see. I definitely like the idea of the null mutation to account for the potential effects of starting with a relaxed structure that's too far from what I actually started with.

This is all very helpful. In terms of rejecting or flagging structures, would you recommend something as simple as a backbone RMSD calculation? And I know this all gets very subjective and there are no hard answers, but if you were to use backbone RMSD, what type of deviation from the starting structure do you think you would personally see as something to be manually inspected or discarded?

P.S. Small world, it turns out we were actually on a publication together with Vikram looking at NDM-1 a few years back! :)

[roccomoretti]

Backbone RMSD could work. I'd say it all depends on your use case and system, and how much change you expect/are willing to allow.

One caveate with RMSD is that the absolute numbers can get perturbed based on how much you're including in the calculation. The same structural deviation can be a tiny RMSD in a large protein, but a comparatively bigger one in a smaller protein, if you're calculating the RMSD across the entire protein. So the actual threshold value for RMSD being "too big" a change will depend a bit on how you're calculating it.

[sean-workman]

Ah I see. At this risk of asking for too much of your time in this thread (which I honestly feel I've already done), is there an alternative strategy that you find is more robust to changes in protein size and might be a bit more broadly applicable? Or are these things you just tend to start to get a feel for as you gain more hands on experience?

Thanks again for all of your help thus far!