You are here

Energy Minimization

25 posts / 0 new
Last post
Energy Minimization
#1

Hi,

I searched for a basic protocol that can be used to minimize my four residue turn structure and I didn't find any. Therefore, I used the following method :-

scorefxn = create_score_function_ws_patch("standard","score12")
mm = MoveMap()
mm.set_bb(True)
pose_move_map.set_chi(True) ## For side chain minimization
mm.set_bb_true_range(1,4) ## Here range means 1 to 4 or the minimzation is applied only to residues 1 and4 ??

minmover = MinMover(mm, scorefxn, 'dfpmin', 10, True) ## I don't know the meaning of dfpmin,10,True. I saw it somewhere and used it

minmover.apply(pose)

After using this protocol the energy of the initial structure which was 20 came down to 6 ... But the energy is still positive. Here can we compare the positive energies or not ??

An another query is regarding the link for the documentation of various functions in pyrosetta, for eg. ( MinMover(mm, scorefxn, 'dfpmin', 10, True))

Post Situation: 
Wed, 2013-07-31 20:45
bharat_46010

Minimization isn't a miracle - it will sometimes get stuck in local minima. Repacking then minimizing is often recommended, and if you really want to reduce the energy, the ramping repulsive FastRelax protocol is suggested (though more for full proteins rather than peptides). Also, if you have an isolated four residue peptide, you might not get the same per-residue energies that a full protein will (because you're missing a number of favorable protein-protein interactions). That's probably not your problem, though (read on).

Regarding what "dfpmin" means, see https://www.rosettacommons.org/manuals/archive/rosetta3.5_user_guide/d7/... It's simply the minimization algorithm that's being used. It's a decent choice, though "dfpmin_armijo_nonmonotone" is likely a better overall choice (though it will depend on your particular case).

The True is a designation regarding neighbor list updating during minimization. This is an implementation detail. There are some score terms which are environmental (neighbor) dependent. Minimization can change the number of neighbors, so constantly updating the neighborlist during minimization gives more accurate energy evaluations, but at the cost of a slower runtime. Leaving it as True will probably be fine. (Looking at the code, there's some complications, and I don't actually think that MinMover setting actually controls much.)

10 is the tolerance, and is probably where your issues are. The default is for relative tolerance, and this is typically a fractional value. Right now you're telling the minimizer to stop when it reaches a value which it estimates to be +- 10x of the final minimum value. E.g., if you expect a final minimum at about -6 REU, you're telling the minimizer you're fine with a final value anywhere from -6 REU up to -6 + -6*-10 = +54 REU. Try using something smaller, like 0.01 or 0.0001. It'll take longer, but you should get better energies.

You can also specify absolute tolerances, by using method names like "dfpmin_atol" and "dfpmin_armijo_nonmonotone_atol", but then you'd with a setting of 10 you'd still be looking at a large range.

Thu, 2013-08-01 11:15
rmoretti

Dear Sir,

Thank you for such a detailed reply. I just want to one thing that how to decide the minimum value for minimizer to stop or a value of 0.01 or 0.001 would be be fine for any system, such as small peptides.

I need to ask one more thing regarding the H-bond detection using Pyrosetta. Is it possible to loosen the criteria of h-bond detection by increasing the bond length from default value (???) to 4-4.5 angstrom

Thu, 2013-08-01 16:33
bharat_46010

The threshold isn't a firm value - it's simply the convergence criteria that is used for the DFP algorithm. Regarding what to set it at, that depends on how stringently you wish to minimize. The tolerances are set up as relative values. So a tolerance of 0.01 means that the final result should be within 1% of the energy value of the local minimum. For a large system giving you total energies in the range of -300 REU or so, that means the end result will likely be within 3 REU of the true local minimum. For a smaller system with total energies around -20 REU, that means the end result will likely be within 0.2 REU of the true local minimum. If you reduce that to 0.0001, the minimizer will execute extra cycles such that the value will be 0.03 REU for the aforementioned large system, or 0.002 REU for the smaller system. Smaller numbers give more precise results, at the cost of greater runtime. It all depends on how stringent you wish your results to be.

Regarding HBond detection, it's typically set up in Rosetta to be energy based. How the energy is computed is a complex statistics-based spline-fit polynomial. While you could theoretically change the polynomial such that the inter-atomic distances are longer, such a change would be non-trivial. Probably easier is to change the threshold energy at which one counts something as a hydrogen bond. That should allow for larger distances, but would also allow for less-ideal angles and orientations. Additionally at some point (I don't know which) the energy would drop to zero, and no more discrimination would be possible. -- If you're finding hydrogen bonds with the HBondSet class, I believe that any (backbone) pairing with non-zero energy is included in the set already, so you would simply adjust your exclusion criteria.

Thu, 2013-08-01 19:25
rmoretti

Actually the the reason for increasing the bond length for h-bond was to find some weak interactions/salt bridges for which the distance between 4-6 angs has been reported in literature. As you told that threshold energy can be changed , can you tell how can that be done ??

Thu, 2013-08-01 20:15
bharat_46010

I'm not quite sure how you're detecting hydrogen bonds currently. If you're using the HBondSet class, that should contain all the potential hydrogen bonds with non-zero energy already, so to change the threshold you would just change at what level you filter them as you iterate through the set.

If you're looking for weak salt bridge interactions, I might suggest forgoing the internal HBond code altogether, and making a custom scanner for the interactions. You can make a double loop through all the residues looking for potential interacting pairs, and then check the relevant atom geometries. This way you can tailor your criteria to whatever specifics you're looking for.

Fri, 2013-08-02 10:53
rmoretti

Hi,

Is it possible to calculate the packing fraction of main chain atoms of residues in a loop/turn. Does pyrosetta has this facility ??

Thu, 2013-08-08 19:43
bharat_46010

Bharat, can you go ahead and create a new topic for this? Thanks!

Fri, 2013-08-09 09:31
jadolfbr

I have posted that question as a new topic in the forum. Can you tell whether its possible to calculate the energy of a certain region of protein in it's native state, I mean the structure taken from PDB. For eg, if I have a peptide of 4 residues with a turn like conformation I am interested only in calculating the replusive/steric forces for the main chain atoms alone?? How can I do this with pyrosetta

Fri, 2013-08-09 17:44
bharat_46010

You can get the energy of a certain subset of residues with the get_sub_score(pose, residues) method of ScoreFunction. Here residues is a vector1 of bools the same length as the pose which indicates which residues to include in the scoring. (You can specify residues to exclude with get_sub_score_exclude_res() method.)

If you want specific scoreterms, probably the easiest way in PyRosetta is to simply create a new ScoreFunction with just those energy terms you want.

Unfortunately, if you just want mainchain atoms, things get a little harder. Rosetta is residue-based, so there isn't a simple way to pull out the sum. Instead you'd probably have to put together something yourself. One easy way is to make a copy of the pose and just mutate all the residues to glycine. This would leave you with just the backbone-backbone interactions. You could try to sum up just the backbone components of the score, but that'd probably be harder to do.

Fri, 2013-08-09 18:22
rmoretti

Dear Sir,

The electrostatic energy between two atoms (N&O) of two different 4 residue peptides is -1.3 and -0.7 (REU). I want to whether the difference between the two energy values can be considered significant or not ??

Sun, 2013-08-18 18:59
bharat_46010

It somewhat depends on the size of the complex and how much noise you expect from the modeling. As a gauge, 0.7 to 1.3 REU would be about what one would expect for a single hydrogen bond.

It's at the level where it's worth paying attention to, but there's certainly no guarantee that a difference in that range would necessarily translate to an experimentally discernible difference. I'd be especially wary if it's an atom-atom energy. The Rosetta energy function is optimized more on a whole residue level (or whole system level), rather than a per atom residue, so you may get large favorable interaction between certain atoms which are balanced by disfavorable interactions between closely associated atoms. Also, atomwise electrostatics in Rosetta aren't as well benchmarked as other scoreterms, so caution would be advised.

One possible way to get a handle on things is to look at sets of randomly generated interactions of a similar type. Most of those won't be good interactions, so you can get a better handle on what the noise level in your particular case is.

Mon, 2013-08-19 19:59
rmoretti

Thanks for the advice, Sir. Actually I am looking into atom-atom energies only, mainly electrostatic interaction and repulsion. I am trying to calculate these energies on two different conformations in the pdb structures. As you told to look into randomly generated interactions of a similar type, in that case my data-set seems to be apt for that ... Am I right here ??

Also I want to know that how can I check for a constraint during minimization.. For eg right now I am using movemap and minmover method of minization and while doing that I have to restrict the bond length between two c- alpha atoms to be around 7 angs. The reason I need to do this is to check that the minimized structure should not deviate from an ideal geometry. How can this be achieved ??

In addition to that, is it possible to sample the side chain rotamers of a residue and look for the lowest energy one amongst them?

I found that I can use packer_task for side chain packing and scoring. As I was not able to find a way to put constraints during minimization , I used MC simulation to perturb the backbone by changing phi psi angles and used the lowest scoring to repack the sidechain.

Here's the script

scorefxn = ScoreFunction()
scorefxn.set_weight(fa_rep,1.0)
scorefxn.set_weight(fa_atr,1.0)

#Initial structure with defined phi psi values for a 4 residue turn
p1 = pose_from_pdb("2fi9_a.pdb")
p1.set_psi(1,135)
p1.set_phi(4,-135)
p1.set_phi(2,60)
p1.set_psi(2,30)
p1.set_phi(3,90)
p1.set_psi(3,0)
p1.dump_pdb("new.pdb")

p = pose_from_pdb("new.pdb")
mc = MonteCarlo(p, scorefxn, 1.0)

packer_task = standard_packer_task(p)
packer_task.restrict_to_repacking() #Only repacking no design
packer_task.temporarily_fix_everything() # Fix all residues
packer_task.temporarily_set_pack_residue(2,True) # allow side chain packing for residue 2
packer=PackRotamersMover(scorefxn,packer_task)

# Monte Carlo approach to vary the phi psi angles within 30 degree

for i in range(1,60000):
mc.reset(p)
p.set_phi(2, p.phi(2)-20+random.random()*10)
p.set_psi(2, p.psi(2)-20+random.random()*10)
p.set_phi(3, p.phi(3)-20+random.random()*10)
p.set_psi(3, p.psi(3)-20+random.random()*10)
mc.boltzmann(p)

if (i%1000 == 0):
mc.recover_low(p)
print mc.lowest_score()
print p.phi(2),p.psi(2),"\t",p.phi(3),p.psi(3)
# Pack the side chain of lowest scoring structure
packer.apply(p)
p.dump_pdb('packed.pdb')

scorefxn(p)
scorefxn.show(p)

# Calculation interactions between side chain and main chain of Asn(residue 2) of original pdb
r1 = p.residue(2)
r2 = p.residue(2)
a1 = r1.atom("O")
a2 = r2.atom("ND2")
print etable_atom_pair_energies(a1, a2, scorefxn)

pp = pose_from_pdb("packed.pdb")
scorefxn(pp)
scorefxn.show(pp)

# Calculation interactions between side chain and main chain of Asn(residue 2) of modified pdb

r1 = pp.residue(2)
r2 = pp.residue(2)
a1 = r1.atom("O")
a2 = r2.atom("ND2")
print etable_atom_pair_energies(a1, a2, scorefxn)

print pp.psi(1)
print pp.phi(2),pp.psi(2)
print pp.phi(3),pp.psi(3)
print pp.phi(4)

Here, I want to know whether is there any option to simultaneously increase and decrease an angle, during MC simulation. For eg, in one of the MC tutorials it is written that :-

"Perturbations are made to the structure by randomly selecting a residue and then perturbing its φ and ψ by a random magnitude, from -25° to 25° " and then they have used the following code

pose.set_phi(resnum, pose.phi(resnum)-25+random()*50)
pose.set_psi(resnum, pose.psi(resnum)-25+random()*50)

In my case I need to perturb the phi psi angles +/30 deg. How can I do this ??

Still, I am looking for a method to constraint distance between two atoms during minimization.

---------
BHARAT

Wed, 2013-08-21 22:19
bharat_46010

To apply a constraint during minimization, you'd need to add a constraint object to the pose and make sure the appropriate score term is on in the scorefunction you're using. In your case it would be a core.scoring.constraints.AtomPairConstraint, added with the pose.add_constraint() method. You also need to make sure the atom_pair_constraint scoreterm has a non-zero weight in the scorefunction you're using during minimization.

If you want to scan through rotamers looking for the best one, I'd look at the RotamerTrialsMover (in protocols.simple_moves). This will randomly pick a variable position in the associated packer task, then scan through all the possible rotamers in the current context, replace the current rotamer with the best one, and then move on to the next position. If you give a packer task with only a single residue allowed to repack, it will optimize that one. You can also use the RotamerTrialsMinMover to minimize each rotamer in the current context prior to scoring.

Regarding +/-30, the random() function returns a number between 0 and 1. random()*50 would thus be between 0 and 50, and -25+random()*50 would be between -25 and 25. You then add that to the current phi/psi value, and reset the angle to that value. To get a different range, you'd adjust the multiplier and the offset to get a random number in the range you want. E.g. for +/-30, you'd use -30+random()*60

Thu, 2013-08-22 10:52
rmoretti

Dear Sir,

Thanks for your response. I used the constraints in the following way :

ca1 = AtomID(2,12)
ca2 = AtomID(2,15)
GF = constraints.GaussianFunc( 5.0, 1.0 )
apc = constraints.AtomPairConstraint( ca1,ca2,GF )
p.add_constraint(apc)

I want to know whether is it possible to keep two strands of a beta-hairpins static and just mutate and change the properties of residues within the turn region of beta-hairpin ?? Does pyrosetta has a secondary structure constraint or something of that sort ??

Thu, 2013-08-22 22:56
bharat_46010

If you're looking to modify the backbone coordinates of the loop while keeping the strands fixed, you're going to need to reorder the foldtree.

Rosetta typically uses a foldtree where the effect of backbone changes propagate to the next residue. What you want to do instead is to set up a foldtree where the two sheets are connected to each other at the root of the foldtree, and the loop is represented by two "tails" which go out to meet each other.

The way this is typically done is with loop modeling. Set up a protocols.loops.Loop object defining the loop - the constructor takes the pose-numbered start residue, end residue and cutpoint residue for the loop. Then add it to a protocols.loops.Loops object with the add_loop() method. Then you can use the protocols.loops.FoldTreeFromLoops mover. Give it the Loops object with the loops() method, and then call the apply() function on your pose. That should set up a foldtree where changes to the loop residues don't propagate to non-loop residues.

The other caveat you need to be aware of is that changing the backbone coordinates of the loop will likely result in chainbreaks (non-natural bond lengths, angles and dihedrals at the cutpoint. You'll want to turn on the chainbreak scoreterm in the scorefunction to tell you how bad it is, and then you can try various loop modeling techniques to attempt to close the chainbreak. (Regular minimization, CCD closure, kinematic closure, cartesian minimization, etc.)

Fri, 2013-08-23 10:50
rmoretti

I applied to fold tree method in the following manner :-

p = pose_from_pdb("VNGH.pdb")
print p.sequence(), p.total_residue()

# Loop consists of two residues only i.e. 13,14

i = 10
j = 17
cutpoint = 13
loop1 = Loop(i,j,cutpoint)
loop_list=Loops()
loop_list.add_loop(loop1)

ft = FoldTree()
ft.add_edge(1,10,-1)
ft.add_edge(10,13,-1)
ft.add_edge(10,17,1)
ft.add_edge(17,14,-1)
ft.add_edge(17,23,-1)

print ft
ft.check_fold_tree()

# change the turn type by changing dihedral angles
p.set_phi(13,60)
p.set_psi(13,-120)
p.set_phi(14,-80)
p.set_psi(14,0)

p.fold_tree(ft)
p.dump_pdb('test.pdb')

# High Resolution Loop Refinement
scorefxn = create_score_function('standard')
loop_refine = LoopMover_Refine_CCD(loop_list,scorefxn)
loop_refine.apply(p)

p.dump_pdb('test1.pdb')

After checking the minimized loop structure, I found that the hairpin structure is entirely disrupted i.e. the strand are far apart from each other. In my case there are only two residues in the loop(res 13 and 14). So, how can I keep the hairpin structure intact...

Fri, 2013-08-23 18:34
bharat_46010

First off, when you print the foldtree, what does it print out?

Secondly, you need to apply the reordered foldtree to the pose *before* you do your backbone manipulation. As you have it now, you change the dihedrals with the regular linear foldtree, and then impose your loop foldtree to the now-separated pose. Pose manipulation happens with the foldtree that's currently in the pose.

Sat, 2013-08-24 15:02
rmoretti

Dear Sir,

I am using the following script to first convert the loop residues to straight chain and then change dihedral angles + loop refinement (fold back the loop to strands) and finally minimization of entire beta-hairpin.

===========================
array1 = ['G','A','V','L','I','C','M','W','F','P','S','T','D','N','E','Q','K','R','Y','H']

fh = open("H:/pyrosetta/turn/scores_hp.txt","a+")
scorefxn = ScoreFunction()
scorefxn.set_weight(fa_rep,1.0)
scorefxn.set_weight(fa_atr,1.0)

for k in range(0,len(array1)):
p = pose_from_pdb("EGDT.pdb")
print p.sequence(), p.total_residue()

i = 9
j = 16
cutpoint = 13
loop1 = Loop(i,j,cutpoint)
loop_list=Loops()
loop_list.add_loop(loop1)

ft = FoldTree()
ft.add_edge(1,9,-1)
ft.add_edge(9,13,-1)
ft.add_edge(9,16,1)
ft.add_edge(16,14,-1)
ft.add_edge(16,24,-1)

print ft
ft.check_fold_tree()
p.fold_tree(ft)

p.set_phi(12,180)
p.set_psi(12,-180)
p.set_phi(13,180)
p.set_psi(13,-180)

mutate_residue(p,11,'G')
mutate_residue(p,12,'%s'%array1[k])
mutate_residue(p,13,'G')
mutate_residue(p,14,'G')
p.dump_pdb("check.pdb")

p.set_phi(12,60)
p.set_psi(12,30)
p.set_phi(13,90)
p.set_psi(13,0)

loop_refine = LoopMover_Refine_CCD(loop_list,scorefxn)
loop_refine.apply(p)

p.dump_pdb('typeI_%s.pdb'%array1[k])

pose = pose_from_pdb("typeI_%s.pdb"%array1[k])
sc = scorefxn(pose)
scorefxn.show(pose)
fh.write("%s\t%s\t"%(array1[k],sc))

for i in range (11,14):
print pose.energies().show(i)

mm1 = MoveMap()
mm1.set_bb(True)
mm1.set_bb_true_range(1,24)
mm1.set_chi_true_range(1,24)
minmover = MinMover(mm1, scorefxn, 'dfpmin', 0.001, True)
minmover.apply(pose)

pose.dump_pdb("aftermin_type-I%s.pdb"%array1[k])
pp1 = pose_from_pdb("aftermin_type-I%s.pdb"%array1[k])
print pp1.residue(12).name(),"\t",pp1.phi(12),"\t",pp1.psi(12)
print pp1.residue(13).name(),"\t",pp1.phi(13),"\t",pp1.psi(13)
print pp1.sequence(), pp1.total_residue()

sc1 = scorefxn(pp1)
phi1 = pp1.phi(12); psi1 = pp1.psi(12)
phi2 = pp1.phi(13); psi2 = pp1.psi(13)

fh.write("%s |\t%s %s %s %s\n"%(sc1,phi1,psi1,phi2,psi2))
fh.close()
==================

I am getting almost similar energies after loop refinement for all the the mutants of beta-hairpin . However, the energies after complete beta-hairpin minimization do not match when I run the same program another time.Why is it so ??... How should I know that when to consider the best score for different trail runs?? Another question is regarding the minimization protocol. I am confused about how to select tolerance value or can I apply the same tolerance value for same amino acid mutation but with different dihedral angles ??

I also checked that, when I change the tolerance value from 0.01 to 0.001 I use to get a folded beta-hairpin for one mutation and if I apply the same tolerance value for another mutant I don't get a folded beta-hairpin. Why is it so ??

If I want to report difference in energies after changing the structure of loop, whether the energy after minimization of entire beta-hairpin structure is correct or energy after loop refinement is correct??

Is there any command to superimpose two structures in pyrosetta ?? Also, I want to know that after applying foldtree + some changes in loop backbone + minimization of the entire beta-hairpin, I use to run for some 50 cycles of minimization and select the structure that has least RMSD from the original structure. In such case the energy of that structure is slightly higher than the one with the least energy (say a difference of 1-2 REU). Is my selection of the model correct as the energies between the best model and the one that I select is quite similar. During the minimization cycle I keep the count at which I get my desired model. Out the two conformations that I am studying, I use to find the one of conformation more early during minimization (eg. round 5-10) whereas, the second one takes time (around 25-30). Can I say that energetically the first conformation is more favorable and easy to obtain compared to the second one ??

Pls find the attached scores file (with strategy used in program in schematics) .

Hope my questions are clear ...

Wed, 2013-08-28 23:51
bharat_46010

The difference in energies after minimization isn't too surprising - the different residue identities likely support the entire hairpin differently, a change that might not be captured in the CCD refinement step. If you take a look at the structures, you might see changes like the sheets of the hairpin coming apart, or other such issues. To a limited extent this can be starting structure dependent - if you repeated the CCD closure multiple times you might get different values (as I believe that CCD is stochastic.) I'd recommend trying exact replicates and see the spread of result energies. I would go with the after-full-minimization energy as indicative of the appropriateness of the mutation, but with the knowledge that minimization can cause the structure to fall apart if the starting structure is not quite right.

Minimization tolerance is typically specified as a relative value. A value of 0.01 means that the minimization algorithm thinks that the final energy value is within 1% of the true value of the local minimum. How stringent you want to set this depends on how precise you want the results. Lower values will give more minimized structures, at the expense of longer runtimes.

Regarding aligning structures, there are a number of functions in the core.scoring namespace. Take a look at calpha_superimpose_pose(modified_pose, reference_pose), which assumes that the two poses are the same length, and superimpose_pose (mod_pose, ref_pose, atom_map) which allows you to specify which atoms to superimpose over.

I wouldn't say that when you get a structure with respect to rounds of minimization tells you anything. It's too sensitive to subtle changes in starting conformation.

Thu, 2013-08-29 13:49
rmoretti

I am using calpha_superimpose_pose command to superimpose two structures and the RMSD comes around 151.9 .I checked the superposition in pymol and there I am got RMSD = 0.056 Angs. Why the value is so large in case of pyrosetta ?

Sat, 2013-08-31 17:28
bharat_46010

Hmm, best I can tell, calpha_superimpose_pose returns the rmsd prior to superposition. Try inserting calls to the non-modifying function core::scoring::CA_rmsd(pose1,pose2) before and after superposition to check.

Mon, 2013-09-02 14:24
rmoretti

Thanks... using CA_rmsd gives the correct value after superposition.

Mon, 2013-09-02 16:49
bharat_46010

Hi,

I have some doubts regarding minimization using Pyrosetta. I am current using Movemap() method of minimization and I am confused about choosing appropriate value for dfp_min. Actually I am studying two different sets of conformation and when I reduce the value for dfp_min lower for one conformer than the other I get better (i.e. lower) energy score for that conformation. If, I apply the same dfp_min to other conformer I get higher energy values. So, is it ok if I use two different values of dfp_min for different sets of conformation and if that's the case, how shall I account for the results I obtain ?? Hope my query is clear ?

Thu, 2013-09-05 16:57
bharat_46010

What do you mean by "reducing the value" for dfp_min? The dfpmin settings are for setting the type of minimizer - there's really not any value that can be increased/reduced there.

I'm wondering if what you mean is that you're changing the tolerance associated with the "dfpmin" minimizer. If that's the case, I'm not quite sure what's going on in your situation. In general, starting with the exact same pose, minimizing with a smaller tolerance shouldn't give you higher energies - it should generally be the same or lower (though internal implementation details may make it such that there's a very slight anomalous difference - especially if you're using the nonmonotone variants). For different poses, though, I don't think there's necessarily any consistency. I wouldn't be surprised if two very similar structures showed an energy difference at one tolerance setting that disappeared or inverted when the tolerance was lowered.

I think you may be putting too much though/effort into finding an "optimal" tolerance setting. My general recommendation is to set the tolerance such that the associated delta is in the noise level for your particular application. (e.g. for the default relative tolerance, a setting of 0.001 on a -300 REU protein will have an associated nominal delta of ~0.3 REU; if you're using the "atol" versions of the minimizer types, then the set tolerance is the delta value.) If you're really looking to find the lowest energy, other techniques are likely more worth your time, such as using the repulsive ramping of the relax protocol (I believe you can use a movemap to do a fixed backbone relax), doing loop remodeling, or using things like cartesian or flexible bond length/angle minimization.

Fri, 2013-09-06 12:42
rmoretti