SequenceRelax on a loop

[pbradley]

Hi there,

Is there a way of doing SequenceRelax on a loop?

So far I have managed to do simple energy minimization with the attached script:

from rosetta import *
rosetta.init()
pose = Pose(“test_in.pdb”)
mm4060 = MoveMap()
mm4060.set_bb_true_range(15,24)
loop = Loop(15,24,20)
set_single_loop_fold_tree(pose, loop)
minmover = MinMover()
minmover.movemap(mm4060)
scorefxn = create_score_function(‘standard’)
minmover.score_function(scorefxn)
print scorefxn(pose)
minmover.apply(pose)
print scorefxn(pose)
minmover.apply(pose)
print scorefxn(pose)
minmover.apply(pose)
print scorefxn(pose)
minmover.apply(pose)
print scorefxn(pose)
pose.dump_pdb(“test_out.pdb”)

but I would really like to do something more like SequenceRelax to allow side chain repacking and gradient minimization (while ramping up ) if this is possible. I am also not sure how many steps of minimization MinMover does each time I apply it as it continues to find a better minimum each time I call it with apply().

Also it seems to me that the bond lengths and angles do not change with the my attached script except at the cut point but ideally I would like these to move as well if possible.

I hope these questions make sense and I am not asking anything too simple – I do not want to make large scale changes to my loop using CCD, etc just to find a local minimum for a loop conformation I already have.

thanks in advance for your help!

[Anonymous]

I can’t answer your actual question, but I can answer bits and pieces. Have you looked for a SequenceRelax Mover…?

“I am also not sure how many steps of minimization MinMover does each time I apply it as it continues to find a better minimum each time I call it with apply().”

There are several minimizer flavors. I know the “flavor” can be chosen as part of the constructor. linmin may be default; linmin is likely to behave in this fashion. The dfpmin and related algorithms are better in several senses. http://www.rosettacommons.org/manuals/archive/rosetta3.2.1_user_guide/minimization_overview.html

“Also it seems to me that the bond lengths and angles do not change with the my attached script except at the cut point but ideally I would like these to move as well if possible.”

Score12 has no terms to handle bond angle and length, only torsion. If the minimizer is allowed to change those DOFs, it will let them drift far outside the accepted range, because it doesn’t know any better. Thus, they are fixed by default. MM-based torsion, angle, and bond length terms exist, but I do not know if they can be accessed from that release, nor can they be trivially activated when using score12 otherwise.

If you want to add sidechain packing, it is easy to add to the code you describe. In pseudocode:

task = PackerTask(pose)
task.restrict_to_repacking()
vector = new vector of booleans, all set to false, length = pose.total_residue()
for i in range 15-24:
vector = true
task.restrict_to_residues(vector)
packer = PackRotamersMover(scorefunction?, packertask?)
#may need packer.task(task)
packer.apply(pose) # as desired


[pbradley]

Ah that makes sense yes I realized it was only doing linmin shortly after I posted this. I guess I’ll just have to ignore bond lengths/angles for this.

I can code something like SequenceRelax using the code above. There is a class called SequenceRelax but it doesn’t seem to take a MoveMap as input it is not obvious to me how to restrict it to move only a specified loop.

thanks again!

[pbradley]

one further question, sorry!

There must be some term in the score fxn that keeps the cut point closed when you minimize? Is it possible to define a loop without specifying a cut point as it seems to be redundant when you are doing gradient minimisation? Perhaps it is necessary for minimization too but I’m not sure why. I noticed in one of my minimization runs that it broke apart at the cut point.

[Anonymous]

The scorefunction term is chainbreak, or something similar to that (maybe cbreak?). Looking at your code, it seems likely it’s not turned on. If you can’t find it I’ll look for it in the C++.

set_single_loop_fold_tree(pose, loop)

This is the real heart of what the chain break does and why it’s important. The loop is kinematically isolated from the rest of the structure by routing protein folding AROUND the loop using a fixed jump between the loop’s termini. Because the fold tree must form a directed acyclic graph, and the termini are connected, a break is necessary. If you don’t isolate and break in this fashion, then when you change loop torsions, the entire downstream half of the protein moves with the loop.

Skipping the loop fold tree call is acceptable if you do only minimization, because the minimizer will penalize breaking the protein in two, but you will see some core shifts if you do this.

Read Chu’s paper about the fold tree for some more details (um, I think).

J Mol Biol. 2007 Oct 19;373(2):503-19. Epub 2007 Aug 2.
Protein-protein docking with backbone flexibility.

Wang C, Bradley P, Baker D.

[pbradley]

hi again,

sorry for all the questions.

if I add the line:

scorefxn.set_weight(chainbreak,1.0)

it gets added to the scorefxn but unfortunately it evaluates to 0 so doesn’t prevent a chain break and there doesn’t seem to be any documentation on how to use it. e.g. :

---

Scores Weight Raw Score Wghtd.Score

---

fa_atr 0.800 -1326.247 -1060.998
fa_rep 0.440 14937.280 6572.403
fa_sol 0.650 906.799 589.419
fa_intra_rep 0.004 731.007 2.924
pro_close 1.000 43.079 43.079
fa_pair 0.490 -52.043 -25.501
hbond_sr_bb 1.170 -48.035 -56.201
hbond_lr_bb 1.170 -84.287 -98.616
hbond_bb_sc 1.170 -25.401 -29.719
hbond_sc 1.100 -22.972 -25.269
dslf_ss_dst 1.000 -10.417 -10.417
dslf_cs_ang 1.000 -8.402 -8.402
dslf_ss_dih 1.000 -3.837 -3.837
dslf_ca_dih 1.000 -2.184 -2.184
fa_dun 0.560 438.271 245.432
p_aa_pp 0.640 -75.816 -48.522
ref 1.000 -38.230 -38.230
chainbreak 1.000 0.000 0.000

---

Total weighted score: 6045.362

[Anonymous]

If chainbreak isn’t triggering, probably the CUTPOINT_UPPER and CUTPOINT_LOWER patches haven’t been applied to the cut.

Try using the function add_cutpoint_variants(pose)?

[pbradley]

ah that did the trick. Thanks – you’ve been a great help!

[Anonymous]

I didn’t write any of the main loop modeling code, but I’ve re-written bits and pieces here and there within other protocols – in other words I’ve seen the bugs before. It also helps that you ask specific questions, which makes them easier to answer. In other words, you’re welcome!

[Anonymous]

Hello, I have a problem with the restrict_to_residues(vector) argument I have not been able to solve. As I understand, the argument needs to be a vector of booleans, so if I try to use a list of booleans as an argument I get the “C++ signature does not match” error. I’ve tried creating a vector with numpy.array but I still get the same error. Maybe this vector needs to be created using a rosetta library, but I have not been able to find a solution. Would someone be as kind as to point me in the right direction?

[Anonymous]

Hello, and thanks for your answers. I’ve been able to use rosetta.utility.vector1_bool() but only with limited success; I’m trying to, given a list of residues, create a vector of booleans with True for those residues and False for the rest. When I try:

to_pack = standard_packer_task(pose)
to_pack.restrict_to_repacking()
length = pose.total_residue()
vector = rosetta.utility.vector1_bool()
for i in range (1, length+1):
 vector.append(False)
to_pack.restrict_to_residues(vector)
packmover = PackRotamersMover(scorefxn, to_pack)
packmover.apply(pose)

Everything goes as planned and no sidechain is packed. However, when I try to set Falses and Trues according to a given list:

list1 = [1,2,3,4,10]
to_pack = standard_packer_task(pose)
vector = rosetta.utility.vector1_bool()
for i in range (1, length+1):
 if i in list1 == False:
  vector.append(False)
 elif i in list1 == True:
  vector.append(True)
to_pack.restrict_to_residues(vector)

I get an error on the lines of:

rosetta.utility.__utility_all_at_once_.vector1_bool object at 0x55322b8
Violación de segmento (`core’ generado)

I’m clueless as to what the reason of this error might be. Is .append the right method to use for this vector? I’ve tried with item assignment (vector = True) but it hasn’t worked.

[Anonymous]

“There is a class called SequenceRelax but it doesn’t seem to take a MoveMap as input”

This is a perennial irritation for me – I know partial abinitio and relax are possible but the authors refuse to make it easy…

[Anonymous]

Sure – its in the FAQ on the PyRosetta website:
http://www.pyrosetta.org/faq#TOC-2.-How-do-I-construct-Vector1-objects-

[Anonymous]

First off, please create a new thread for unrelated problems.

Your question happens to be covered in the FAQ (http://www.pyrosetta.org/faq#TOC-7.-Where-are-Vector-objects-)
Try seeing if rosetta.utility.vector1_bool works. If not, you can possibly make do with something like rosetta.utility.vector1_Size

[Anonymous]

I think it’s your python code:

Your getting an empty vector as your not appending.

-> Your if statements don’t do anything. Not exactly sure why, its a logical line, but it doesnt work.
To test, If you do this:
if 3 in list1 == True:
print “yes”

Nothing happens.

So…
Try this instead of if i in list1 == True/False:
for i in range(1, length+1):
if i in list1:
vector.append(True)
else:
vector.append(False)

In the future, try peppering print statements around your code when your debugging. Some things in Python SEEM like they should work, but don’t. I have run into this many times. Also, try iPython. Its great for debugging or writing python code. Just to see if something works or not.

-J

[Anonymous]

Thank you so much for your help, the problem is now solved.

[Author]

Posts