Re: Using a database of loop conformations together with de novo folding protocol

[Anonymous]

Hi,

I want to know two things about a design experiment involving loop modeling:-
1. I have a database of loop conformations, that I have extracted from PDB. I want to model those conformations onto a beta-hairpin structure. How can I do that ??
2. Is it possible to apply de-novo folding protocol only for the loop-sequence, keeping the hairpin sequence unaltered.

---

Bharat

[Anonymous]

1. I don’t know of a way to ask Rosetta to thread particular loop conformations onto a structure. It may be possible to do this by converting your database to either a traditional fragments database or a LoopHash database. In both cases you’d need to modify the C++ code significantly to get them to sample exhaustively and once-through instead of randomly.

I would instead probably write a script (probably in PyRosetta) to generate semi-closed structures of each of your loops. Take your starting structure and copy the internal coordinates and sequence of your database loop onto the structure, write the pose, put the next loop on, write the pose, … Depending on your exact needs I can put together some pseudocode for you. After creating these starting structures, feed them to loop modeling.

Since you know what PDBs your loop fragments are from, you can do fragment picking (https://www.rosettacommons.org/manuals/archive/rosetta3.4_user_guide/dc/d10/app_fragment_picker.html) using only those PDBs as input, then use fragment-based CCD loop closure (https://www.rosettacommons.org/manuals/archive/rosetta3.4_user_guide/de/d6f/ccd_loop_modeling.html) to stay close to the inserts.

2. I’m not sure what you mean. I think you mean, “If I have a loop of 10 residues, can I keep 4 in the middle that form a hairpin fixed?” The answer is yes, but it’s kind of hacky; the AnchoredDesign interface-design module modifies loops in this fashion. Is that what you meant; if so I can explain further?

[Anonymous]

Like Steven suggested, I think PyRosetta is the way to go, if you can’t get the fragment picker setup and are up for some scripting. There is now code in Rosetta to make grafting easier, most of which is made up of functions he wrote for AnchoredDesign.

Lewis SM, Kuhlman BA. Anchored design of protein-protein interfaces. PLoS One. 2011;6(6):e20872. Epub 2011 Jun 17.

You would import the namespace like this: import rosetta.protocols.grafting as graft
Once you have the newest PyRosetta binaries from http://www.pyrosetta.org , take a look at this file where grafting is employed in the GUI for an example: /GUIs/pyrosetta_toolkit/protocols/GraftingProtocols.py

The function that actually runs the graft in this file is run_graft(), with the grafting class being: graft.AnchoredGraftMover(start, end)
There are utility functions you can use as well, use iPython to get an idea of what you can do: graft.delete_region(pose, start, stop); etc.

Here is the description of the class from the C++ file:

///@brief Grafting class adapted from Steven Lewis’ pose_into_pose algorithm. Basic, and quick, but with many options.
///
/// example:
/// mover = AnchoredGraftMover(start, end)
/// mover.set_piece(piece, cter_overhang, nter_overhang)
/// mover.apply(pose)
///
/// see also: grafting/util.hh.
///
///@details Uses a single loop and a single arm to close the loop by default.
/// ****Nter_loop_start—->Piece—-> | Cter_loop_end****
/// Default movemap keeps insert Frozen in dihedral angle space, But think of the insert as part of a giant arm.
/// Default flexibility on Nter and Cter is only two residues (–> part of diagram).
/// Will delete any residues between start and end, and any overhang residues from the insert.
///
/// Algorithm originally from pose_into_pose:
/// The insert will be left unchanged in internal-coordinate space except for the phi on the first residue, and the psi/omega on the last residue, and atoms whose bonding partners change as a result of the insertion.
/// Internally, apply performs the insertion, idealizes the loop residues (omegas to 180, peptide bonds idealized) and the newly made polymer connections at the insert point, and then attempts to close the loop(s).
/// It is intended, but not guaranteed, to produce a loop with good rama, omega, and chainbreak/peptide_bond scores. It does NOT attempt to give a loop with good sidechains (it does not repack at all) or worry overmuch about van der Waals
///

Just for more information that could be useful:
These are the simplest functions of the class to control which residues are used to sample and do the insert (the class also accepts a movemap if you know what that is):

///@brief Sets scaffold flexiblity on either end of scaffold
set_scaffold_flexibility(Size const Nter_scaffold_flexibility, Size const Cter_scaffold_flexibility);

///@brief Sets insert flexibility on either end of insert
set_insert_flexibility(Size const Nter_insert_flexibility, Size const Cter_insert_flexibility);

The example should help you further, or we could both help you with the script.

-J

[Anonymous]

Note that the part 1 grafting’s relationship to AnchoredDesign has no bearing on the appearance of AnchoredDesign as my answer to part 2 of your question – it’s coincidental that I’ve worked on both…

[Anonymous]

Thank you J for ur detailed advice. I have some questions regarding the script :-

1. The default flexibility mode needs/makes two residues of either side of loop as flexible. In my case I have only 2 residue in the turn and other two are present at the end of beta-hairpin strand. What I want is to allow the flexibility for these four residues only. How can that be done ??

2. mc = MonteCarlo(scaffold, scorefxn, 1.0); #We will use the montecarlo object to get the lowest graft found. This will also set the lowest score at this starting model. If you don’t want this as a bias, construct the MonteCarlo after you delete the region. I didn’t understand the meaning of the comment the comment .

3. Another thing that I want to know is that , can I specifically look at the reason for high energy structures during loop grafting. For eg change in the dihedral angles or some constrained bonds or clashing atoms. As my study is not related to design. It is more of a theoretical study to see why the energy increases when one type of loop is grafted to another type.

I hope I made my points clear.

Thanks

---

BHARAT

[Anonymous]

Must be a glitch – I just checked all the links and all of them work except for PyRosetta v1.1 which could be downloaded by using this direct link: http://graylab.jhu.edu/pyrosetta/downloads/release/ but as Jared said: v1.1 is really outdated and will be of really little practical use.

Bottom line: please try to download again and if you still got an error post it here with full URL that does not work.

Thanks!

[Anonymous]

For the twist sign, you don’t need to use the whole acos/degree thing – you’re just looking for the sign (whether the two vectors are in the same direction or not).

triple_product = cross_prdt.dot(ca_vec)
if triple_product < 0:
….print “Left Handed”
….twist_sign = -1
else:
….print “Right Handed”
….twist_sign = +1

Other than that, I can’t see any obvious issues. (Again, it’d be worth testing the procedure on a known case to double check.)

[Anonymous]

For Koga et al. when they refer to the chirality of the beta hairpins, they’re talking about which direction the loop of the hairpin is turning with respect to the orientation of the residue immediately prior to the loop of the hairpin – or equivalently, which side of the N-terminal beta strand the C-terminal beta strand is on.

To my understanding, the difference between type I’ and type II’ beta hairpins is the orientation of the peptide bond (the orientation of the carbonyl) between the two residues of the two-residue turn (see http://swissmodel.expasy.org/course/text/chapter2.htm for a diagram). This sense of handedness is different from the handedness that Koga et al. talk about. For example, both of the hairpins in the linked figure have the same Koga-handedness, but one is type I’ and one is type II’. — This goes back to the point I was trying to make before. There are multiple senses of handedness at play here, and they don’t all measure the same thing, so they’re not going to agree for all structures. You need to carefully consider what it is you actually want to use the results of the chirality calculation for, in order to pick the correct calculation and get interpretable results.

[Anonymous]

Thank you both for your detailed suggestions. Lewis, for my question, pyrosetta code would be suitable for me.Since I have not done any programming in pyrosetta or python, I think your pseudo-code would provide me some help. As far as the second question is concerned, I don’t want to fix the loop . I want to fix the surrounding strand regions, thereby allowing loop residue to move. What I want to look here is how the dihedral angles of loop residues affect the energy of the structure and also, how steric forces come into play to either stabilize or destabilize the turn(in terms of energy). I hope my explanations are clear to you …

[Anonymous]

Jared – do you already have a working python insertion example? I don’t want to spend time writing poor pseudocode when you already solved this.

[Anonymous]

Sorry, yes, it’s in the file I pointed to.

For clarity, that code is splicing a region from one protein (from_pose) to another (pose) using the grafting functions/class.

Here is some more detail, as how to use the class and functions are within the header files in trunk:

The region includes itself plus a few residues on either side which are used for a superposition (graftmover.superimpose_overhangs_heavy), then deleted once the graft is started. Superposition is really only used if either of these settings are called: graftmover.set_use_double_loop_double_CCD_arms(True) or graftmover.set_use_double_loop_quad_CCD_arms(True). These will keep your loop frozen in cartesian space and use residues on either side to close the loop. This is less aggressive then the default. If your loops don’t fit well into your scaffold, use the default.

Here is some code close to what your script will look like (I should really add rosettascript support at some point too..):
#I’m assuming you will have your loops in separate PDB files in a directory, and a scaffold protein. You should generate a text file that lists the path to each loop. I’m also assuming you don’t have any overhang residues and are using the default method, default flexibility (2 residues on either side in the scaffold).

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
from rosetta import *
import rosetta.protocols.graft as graft

rosetta.init()

scaffold = Pose()
pose_from_pdb(scaffold, “path_to_scaffold.pdb”)

scorefxn = create_score_function_ws_patch(“standard”, “score12”)
mc = MonteCarlo(scaffold, scorefxn, 1.0); #We will use the montecarlo object to get the lowest graft found. This will also set the lowest score at this starting model. If you don’t want this as a bias, construct the MonteCarlo after you delete the region.

#Here we delete the region:
start_residue = scaffold.pdb_info.pdb2pose(“A”, 23)
stop_residue = scaffold.pdb_info.pdb2pose(“A”, 27)
graft.delete_region(scaffold, start_residue, stop_residue); #Will delete from start->stop including start+stop residues

graftmover = AnchoredGraftMover(start_residue-1, stop_residue+1)

#Settings
graftmover.set_cycles(500); #May need to experiment with this. It depends on how hard it will be to insert your loop.
graftmover.set_use_smooth_centroid_settings(True); #This seems to improve results a bit. Uses a different centroid scorefunction

#Try out each loop (Forum eats the indentation)
LIST = open(“path_to_pdb_list.txt”, ‘r’)

//////////Start For Loop////////////
for pdb_path in LIST:
pdb_path = pdb_path.strip()
temp_pose = Pose()
loop_pose = Pose()
pose_from_pdb(loop_pose, pdb_path)

#Copy the pose
temp_pose.assign(scaffold); #This is because the numbering will change each time we graft. And instead of updating this works.

#Initialize grafter
graftmover.set_piece(loop_pose, 0, 0)

#Run + see if it is lower in energy. If it is, keep it.
graftmover.apply(temp_pose)

##Here you can do some sampling. Remodeling the loop, relaxing the loop, repacking, etc. Whatever you want to do. I’m going to repack the loop and residues used to close the loop.
graftmover.repack_connection_and_residues_in_movemap_and_piece(temp_pose, scorefxn)

score = scorefxn(temp_pose)
if score < mc.lowest_score():
mc.set_lowest_score_pose(temp_pose)
//////////End For Loop//////////////

LIST.close()
final_pose = mc.lowest_score_pose()
final_pose.dump_pdb(“Lowest_pose_encountered.pdb”)

[Anonymous]

1) This flexibility is the default – 2 residues on either end. Though grafting may be overkill in this regard. Its up to you, but you may want to just copy the residues into your pose instead of doing all the grafting. If you have 2 residues in the turn in your starting pose, this can be done using the function graft.replace_residues(from_pose, to_pose, from_pose_start_residue, to_pose_start_residue, insertion_length) as Steven suggested originally.

2) Hard to explain without first explaining what the monte carlo object is. Basically, upon construction, it will set its lowest_energy_pose variable to the pose you gave it. If you were doing design, it may favor that native pose you pass to it. If you’re interested in all this, I suggest spending a day going over the PyRosetta tutorials. You may find them useful for more projects other then this. http://www.pyrosetta.org/tutorials Specifically, the montecarlo tutorials are Workshops 4 + 5.

3) Interesting! What you would want to do is print the energy information for the pose scorefxn.show(pose), and dump the pose after each graft/replacement (instead of only at the end) to look at it. There are ways of getting specific energy term information per residue or pair of residues as well. To do per residue or a pair of residues, you will want to look into the EmapVector object. For an example of this, checkout /GUIs/pyrosetta_toolkit/modules/RegionalScoring.py. Something else that may be useful to you here is the PyMOLMover object. This is covered in the tutorials. Basically it will allow you to send your pose to PyMOL and color it by a specific energy term (or total_score) in PyMOL after each graft. It can also send H-Bonds, etc. What’s nice is that it can send them as different frames, and you can compare them to the energies that are printed out.

I know all of this is a bit too much information, but I guess I like to be thorough. Have a look at PyRosetta and the tutorials, and things will seem a bit more clear. Definately a very cool project!

-J

[Anonymous]

I misunderstood 1. You have 2 residues in the turn you want flexible. You also have two residues on one end, or one residue on each end of each strand? Regardless, this can be done, but for graft, you would not want to allow flexibility in the entire turn. This will completely mess up the structure of the loop, and it may even crash PyRosetta.

What you will want to do is graft with those other 2 residues, or at least try it since you only have a 2 residue insert. Like I said in the other post, replace_residues may work better. After that, you will want to use the MoveMap to specify which residues you want flexible, and choose how you want to optimize all 4 residues. You could repack sidechains. Do loop modeling, or simply minimize/FastRelax the 4 residues. I would recommend this, as you won’t mess up your turn too much. I can help with how to do this, once you get here. Movemap is covered in the tutorial, and I can show you how to use FastRelax as it’s not really anywhere in the tutorial or guide.

[Anonymous]

Thanks J. As per you suggestion, I think I should spend some time learning pyrosetta in order to understand clearly what you are saying. As far as the loop refinement is concerned, simple minimization would be better. The reason being , I want to look at the behavior of backbone and not the sidechains.

Right now I am using Fedora 15(32-bit). There is no download version of Pyrosetta for fedora/32 bit linux system. What about the ubuntu -32 bit ver, would it work in fedora ??

[Anonymous]

“Right now I am using Fedora 15(32-bit). There is no download version of Pyrosetta for fedora/32 bit linux system. What about the ubuntu -32 bit ver, would it work in fedora ??”

Probably. Try it and let us know.

[Anonymous]

I tried downloading the windows version, it seems that the link fails to open …

[Anonymous]

I tried and also got a weird result. I let the sysadmin know to take a look at it.

[Anonymous]

The windows version is extremely old and buggy. Were talking V1.1, which only had a fraction of the available namespaces bound. Many more were buggy…

I take it the Fedora version did not work? Have you thought about switching to Ubuntu?

[Anonymous]

My weird result (file does not exist) was with the Ubuntu version – I suspect there may just be a server problem.

[Anonymous]

V2.012, Sorry. Still a few years old, but if Sergey is able to fix the link, should be fairly useful.

[Anonymous]

Hi,

I tried downloading the windows version 2.012 (http://graylab.jhu.edu/pyrosetta/downloads/release/download.release.Windows.32Bit.html) when I was in windows and the link didn’t work. But when I accessed the same link in linux it worked. I don’t know the reason why it happened like that .

I have installed the PyRosetta in / and when I run the command pose_from_pdb(“TEST”) , an error from windows is flashed and then the ipython prompt gets closed.

[Anonymous]

Ok. So, you are first running SetPyRosettaEnvironment.sh? Then, in IPython calling

from rosetta import *
rosetta.init()

?

Then, what you need is: p = Pose()
pose_from_pdb(p, “TEST.pdb”)

TEST.pdb should be the full path if it’s not in the directory you started IPython from…

-J

[Anonymous]

Hi,

I have started working on your script, but I am getting an error at graft.delete_region command. I checked the pyrosetta folder and I didn’t get any protocol related to grafting .. Also, instead of grafting, if I want to replace the the central two residues with ones in my local db, how could this be achieved.

[Anonymous]

First, I’m glad that you were able to get PyRosetta working. It can be tricky sometimes.

For the grafting (and the GUI): It’s only in the newest PyRosetta releases, not the old ones (which the windows version is).

Since the utility functions are fairly small, I might be able to rewrite them in python as a module. Which you would then import. If you can’t use Ubuntu 12.04 for some reason (VirtualBox works great in this instance and is easy to setup), I can write this up. You would need to test it, as even the newer windows version is fairly old. Let me know what you think.

[Anonymous]

I will try setting up virtual box for ubuntu and will contact you after that..

Thanks

---

BHARAT

[Anonymous]

Hi,

I have installed virtual machine with ubuntu 12 (32 bit). I installed pyrosetta and ran a test script, which give the following error :

Traceback (most recent call last):
File “test.py”, line 1, in
from rosetta import *
File “/home/bharat/Downloads/PyRosetta.Ubuntu-12.04LTS-r54814.64Bit/rosetta/__init__.py”, line 27, in
import utility
File “/home/bharat/Downloads/PyRosetta.Ubuntu-12.04LTS-r54814.64Bit/rosetta/utility/__init__.py”, line 1, in
from __utility_all_at_once_ import *
ImportError: /home/bharat/Downloads/PyRosetta.Ubuntu-12.04LTS-r54814.64Bit/rosetta/utility/__utility_all_at_once_.so: wrong ELF class: ELFCLASS64

[Anonymous]

Your running 64 bit PyRosetta on a 32 bit ubuntu. Is your system architecture 64 bit? Try installing 64 bit ubuntu, and you should be good to go.

[Anonymous]

I tried that .. but my system’s architecture is 32 bit and VM gives error while installing 64 bit Ubuntu

[Anonymous]

The last hora is to compile from source, which is a bit of a pain. It can also take quite a long time. I recommend this, but to use the grafting stuff you would need the newest developers version… It looks like you would need to either upgrade your system, or try another route for the grafting…

[Anonymous]

Hi, I have succesfully installed pyrosetta(64 bit) on centos 6.0. I am getting the error while running the scipt :
name ‘AnchoredGraftMover’ is not defined

[Anonymous]

My fault – instead of

graftmover = AnchoredGraftMover(start_residue-1, stop_residue+1)

do

graftmover = graft.AnchoredGraftMover(start_residue-1, stop_residue+1)

[Anonymous]

Hi,

I want to calculate the distance vector for the mid-point of C-N of one peptide(i) bond and midpoint C-N of successive peptide bond(i+1). How can I do that?? (Please see the attached file)

[Anonymous]

There’s no compiled application which can do it, but if you can use PyRosetta, or don’t mind coding C++, it’s relatively easy.

You can get a particular residue n (pose numbered) from the pose with:
residue = pose.residue(n)

You can get the atomic coordinates for a particular atom in the residue with:
xyzvector = residue.xyz(atomname)

numeric::xyzVector has all of the common vector manipulation routines, so you can get the midpoint between two atoms by averaging:
midpoint = (atom1_vector + atom2_vector)/2

Then you can get the vector between them by subtraction, or the distance using:
dist = xyzvector1.distance(xyzvector2)

[Anonymous]

Dear Sir,

I tested a code for calculating the angles between two vectors. One vector is the midpoint between the adjacent C-N(i,i+1 residue). Another vector is the midpoint between the adjacent C-N(i,i+1). Both vectors form one strand of anti-parallel beta-hairpin.In order to calculate the angle between these to vectors, I used the following code :

C1_xyz = p.residue(18).xyz(“C”)

N1_xyz = p.residue(19).xyz(“N”)

C2_xyz = p.residue(19).xyz(“C”)

N2_xyz = p.residue(20).xyz(“N”)

C3_xyz = p.residue(25).xyz(“C”)

N3_xyz = p.residue(26).xyz(“N”)

C4_xyz = p.residue(26).xyz(“C”)

N4_xyz = p.residue(27).xyz(“N”)

midpoint1 = (C1_xyz + N1_xyz)/2

midpoint2 = (C2_xyz + N2_xyz)/2

midpoint3 = (C3_xyz + N3_xyz)/2

midpoint4 = (C4_xyz + N4_xyz)/2

#vector between midpoints 1&2 and between midpoints 3&4

vec1 = midpoint2 – midpoint1
vec2 = midpoint4 – midpoint3

#To calculate the angle between the two vectors
vec3 = vec1.dot(vec2)

#print vec3
-12.03716225

(Please see the attached image for clear explanation of my query )

Can u please tell whether my assumptions and calculations are correct or not ??

[Anonymous]

The dot product of two vectors is equal to the cosine of the angle, times the two magnitudes, rather than the dot product itself.

There’s actually a utility function to compute the angle between two vectors, which will take care of all that:

angle_of(vec1, vec2)

This gives the value in radians (zero to pi), but you can use the numeric::conversions::degrees() function to do the conversion to degrees, if that’s what you need.

[Anonymous]

Actually, I followed or misunderstood the dot product notation from Workshop tutorial 2 of Pyrosetta [..find the N–Cα–C bond angle using the vector dot product function,v3 = v1.dot(v2)..].
I checked for angle_of() in pyrosetta shell, which gives me an error : NameError: name ‘angle_of’ is not defined. How can I use it then ??

Sir, my final objective is to calculate the twist-angle and direction of twist for anti-parallel beta-hairpins. For this I am following the JMB paper (Twist and Shear in beta-sheets and beta-ribbons). As per the paper –

The sheet twist is measure of twisting of beta-sheet. It is defined as the angle between the backbone vectors of the two residues in the inter-strand pair.

Backbone vector for a residue is defined as the vector joining the point equidistant from N and C atoms of the two peptide units that form the bacobone of residues. Sheet twist for a hydrogen bonded and non-hydrogen bonded pair as the acute angle between the backbone vectors, b1 and b2, of a pair of inter-strand neighboring residues(See Figure 1). If the interstrand pair is right-hand twisted, then the cross product b1 X b2 will be pointing roughly in direction of d21, the vector from Calpha1 to Calpha2. Otherwise, the interstrand pair is left-hand twisted. The sign of sheet twist is defined as the sign of scalar product between b1Xb2 and d21 i.e. ((b1Xb2).d21).
Please refer figure 1 for the above description.

In my case I have slight changed the diagram for my calculations and understanding(See figure 2). I used the pyrosetta code for calculating the twist angle and sign of twist-angle ??. Here’s the code

In [252]: C1 = p.residue(15).xyz(“C”)

In [253]: N1 = p.residue(16).xyz(“N”)

In [254]: C2 = p.residue(16).xyz(“C”)

In [255]: N2 = p.residue(17).xyz(“N”)

In [256]: mp1 =(C1+N1)/2

In [257]: mp2 =(C2+N2)/2

In [258]: vec2 = mp2 – mp1

In [259]: C3 = p.residue(18).xyz(“C”)

In [260]: N3 = p.resi
p.residue p.residue_type

In [260]: N3 = p.residue(19).xyz(“N”)

In [261]: C4 = p.residue(19).xyz(“C”)

In [262]: N4 = p.residue(20).xyz(“N”)

In [263]: mp3 = (C3+N3)/2

In [264]: mp4 = (C4+N4)/2

In [265]: vec1 = mp3-mp4

In [266]: twist_angle = vec2.dot(vec1)

In [267]: print twist_angle
10.69854225

In [270]: ca2 = p.residue(16).xyz(“CA”)

In [271]: ca1 = p.residue(19).xyz(“CA”)

In [272]: ca_vec= ca2 – ca1

In [273]: twist_sign = (vec1.cross(vec2)).dot(ca_vec)

In [274]: print twist_sign
10.644989098

I kindly request you to check whether my calculations are correct or not .

Regards

---

BHARAT

[Anonymous]

I’m not too experienced with PyRosetta, so it could be that the function isn’t exposed to PyRosetta – it’s a templated inlined friend function of xyzVector, so the Python wrapping might not be picking it up.

You can get the same effect directly:

#indents as leading dots, so the forum doesn’t eat them.
import math

mag = vec1.length() * vec2.length()
if mag == 0:
….twist_angle = 0 # Shouldn’t happen, if the input vectors are reasonable
else:
….twist_angle = math.degrees(math.acos( vec1.dot(vec2)/mag ))

Note for anti-parallel strands, the angle calculated will be obtuse, not acute, unless you flip the direction of one of the vectors. But if you’re looking also at the cross product, I wouldn’t suggest arbitrarily flipping the vector direction – doing so would reverse the sign of the cross product. Instead, just subtract from 180 degrees to get the acute angle:

if twist_angle > 90.:
….twist_angle = 180.0 – twist_angle

Regarding handedness, the key thing is the sign of the triple product. If the triple product is greater than zero, it’s one handedness, and if it’s below zero it’s the other. Which is which depends on how the directions of the vectors and handedness is defined. I’d recommend closely examining the examples you want to be consistent with, and checking that the way you’re computing things comes out with the answer that it should be.

[Anonymous]

Dear Sir,

I have calculated the twist_angle and twist_sign for one of the turns. The value for twist is around 30 deg. and twist sign in positive.
Here’s my calculation. Please see the attached figure for vector definition and calculation. This time I didn’t reverse the direction of the first vector. The twist_angle came around 149 deg. Therefore, the actual twist calculated is 30.9 deg (180-149.1). For calculating the twist direction I used the convention followed in the JMB paper that I am referring. As per the paper , if the interstrand pair is right-hand twisted, then the cross product b1 X b2 will be pointing roughly in direction of d21, the vector from Calpha1 to Calpha2. Otherwise, the interstrand pair is left-hand twisted. The sign of sheet twist is defined as the sign of scalar product between b1Xb2 and d21 i.e. ((b1Xb2).d21). Here’s the code

In [310]: C1 = q.residue(21).xyz(“C”)

In [311]: N1 = q.residue(22).xyz(“N”)

In [312]: C2 = q.residue(22).xyz(“C”)

In [313]: N2 = q.residue(23).xyz(“N”)

In [314]: mp1 = (C1+N1)/2

In [315]: mp2 = (C2+N2)/2

In [316]: vec1 = mp1 – mp2

In [317]: C3 =q.residue(24).xyz(“C”)

In [318]: N3 = r.residue(25).xyz(“N”)

In [319]: N3 = q.residue(25).xyz(“N”)

In [320]: C4 = q.residue(25).xyz(“C”)

In [321]: N4 = q.residue(26).xyz(“N”)

In [322]: mp3 = (C3+N3)/2

In [323]: mp4 = (C4+N4)/2

In [324]: vec2 = mp3-mp4

In [325]: mag = vec1.length * vec2.length

In [326]: print mag
9.85137085158

In [327]: twist_angle = math.degrees(math.acos(vec1.dot(vec2)/mag))

In [328]: print twist_angle
149.120449995

In [329]: print 180 – twist_angle
30.8795500048

In [330]: ca1 = q.residue(22).xyz(“CA”)

In [331]: ca2 = q.residue(25).xyz(“CA”)

In [332]: ca_vec = ca1 – ca2

In [343]: cross_prdt = vec1.cross(vec2)

In [344]: twist_sign = math.degrees(math.acos(cross_prdt.dot(ca_vec)/mag2)

In [345]: print twist_sign
159.579574946

Now, is my calculation correct as the figure(and its nomenclature)

[Anonymous]

Dear Sir,

I understood the calculation for twist angles. I still have a problem about sign/direction of twist. I understood the description and when I apply the same to code, the sign comes negative . This means that the hairpin in left-handed. However, the turn connecting the hairpin is right-handed which means that the hairpin should also be right handed. In calculation when I change the vectors in cross product, the sign also changes to positive. The question that I have is then how to decide the order of two vectors in cross product. Here’s the output

In [511]: twist_axis = vec1.cross(vec2)

In [512]: print twist_axis.dot(ca_vec)
-33.7640500425

In [513]: twist_axis = vec2.cross(vec1)

In [514]: print twist_axis.dot(ca_vec)
33.7640500425

Also, I have changed the nomenclature according to my understanding . If you see the figure 1 (from jmb paper, figure 1) and compare it with mine (figure 2) , you can find that both the figures are opposite.

[Anonymous]

There’s two handednesses going on here. The first is the handedness of the hairpin itself – the direction in which the loop connecting them is turning with respect to the orientation of the sidechains of the strands. Ho and Curmi are talking about something else when they discuss handedness – the overall curvature of the beta sheet (like the curvature which makes the TIM barrel a barrel). These aren’t necessarily going to be the same handedness.

Regarding the order in which to do the cross product, just be very careful to follow the conventions used in the source paper. Luckily, Ho & Curmi note that the way they define the curvature is invariant to swapping the labels on the two vectors. So the sign of the handedness is given by ((b1 × b2)·d21), regardless of which strand you pick as 1. The trick is to make sure that between-strand d21 vector is going from strand 2 to 1, rather than the other way around.

[Anonymous]

Sir,

I am interested only in finding the handedness of the beta-hairpin alone. So, if you refer to figure 3 from the previous post, I have defined the following for twist_sign:
ca_vec = ca1-ca2
twist_axis = vec2.cross(vec1)
twist_sign = twist_axis.dot(ca_vec)

If the sign is positive the twist is right handed else if the sign is negative the twist is left handed. Is this correct ??

[Anonymous]

If you’re interested in the handedness of the strand-loop-strand hairpin – that is, the direction in which the second strand comes off of the first, rather than the cupping of the two strands, the method of Ho & Curmi isn’t going to help. You’ll want a reference that discusses the handedness of the hairpin specifically. Koga et al. “Principles for designing ideal protein structures.” http://dx.doi.org/10.1038/nature11600 would be a better bet.

In it they discuss the handedness of a beta-beta hairpin in terms of the N-C vector of the residue just preceding the loop (vecNC), the Calpha-Cbeta vector (perpendicular to the plane of the sheet) of that residue (vecAB), and the Calpha-Calpha vector from the residue just preceding the loop to the one just following the loop (vecPF).

In this case, a right handed hairpin has vecAB.cross(vecNC) in the same direction as vecPF, and left handed hairpins have it pointed in the opposite direction. Or in other words, if vecAB.cross(vecNC).dot(vecPF) is positive, the hairpin is right handed, and if it’s negative, it’s left handed. (In other words, if you place your thumb of your hand along the Calpha-Cbeta vector of the residue preceding the loop, the loop curls in the same direction the fingers do of the corresponding hand.)

[Anonymous]

Thanks for the reference , Sir.

I was wondering whether there is any method to find out hydrogen bonded pairs and non-hydrogen bonded pairs for an anti-parallel beta hairpin. Actually I tried to figure it out from the dssp output. For eg.,

# RESIDUE AA STRUCTURE BP1 BP2 ACC N-H–>O O–>H-N N-H–>O O–>H-N
83 85 A L E -G 77 0A 37 20,-0.3 20,-2.7 -6,-0.2 2,-0.5
84 86 A W E -GH 76 102A 2 -8,-2.6 -8,-2.5 -2,-0.4 2,-0.3
85 87 A Q E -GH 75 101A 59 16,-2.9 16,-2.1 -2,-0.5 2,-1.0
86 88 A I E + H 0 100A 16 -12,-2.1 14,-0.2 -2,-0.3 3,-0.1
87 89 A Q E – 0 0A 78 12,-3.3 2,-0.3 -2,-1.0 -1,-0.2
88 90 A H E – H 0 99A 60 11,-1.1 11,-2.3 -3,-0.2 2,-0.3
89 91 A H E + H 0 98A 24 -2,-0.3 2,-0.3 9,-0.2 -84,-0.2
90 92 A L E – H 0 97A 43 7,-1.9 7,-3.3 -2,-0.3 2,-0.6
91 93 A M E – H 0 96A 81 -2,-0.3 2,-0.5 5,-0.2 5,-0.2
92 94 A V E > S- H 0 95A 35 3,-3.1 3,-1.9 -2,-0.6 -2,-0.0
93 95 A R T 3 S- 0 0 239 -2,-0.5 -1,-0.1 1,-0.3 3,-0.1
94 96 A G T 3 S+ 0 0 66 1,-0.2 2,-0.4 0, 0.0 -1,-0.3

if the NH-CO and CO-NH is equal for a particular residue then it is hydrogen bonded to that residue no., which is equal to the NH-CO bonded residue. For eg residue 87 is hydrogen bonded to 87+16 =103. How to look for non-hydrogen bonded pairs then , is my query ?? As of now, I am considering the pair next to hydrogen bonded pair as non-hydrogen bonded one, but that’s not the case with every beta-hairpin.

In addition to that regarding the calculation of twist of beta-sheet twist, I have a query regarding the length of beta-strands to consider. As I have separately calculated the twist of beta-turn as pseudo-dihedral angles between four C-alpha atoms. Therefore in the calculation of twist of beta-strands connecting the loop should I remove those end-residues that I included in calculating the twist of beta-turn ??

Also, there is another query regarding adding the vectors. I am trying to append the set of vectors from one strand to a list and I print the list I get an output as : [, 0AE068F0>, ]

[Anonymous]

Regarding determining hydrogen bonds, you can get Rosetta to detect them with the core.scoring.hbonds.HBondSet class, initialized with HBondSet(pose). You can then access hydrogen bonds with hbondset.hbond( index ), with indexes from 1 to hbondset.nhbonds(). You can then use the returned core.scoring.hbonds.HBond object to look at the residues and atoms involved with the hydrogen bond. To find if a particular pair is not in a hydrogen bond, the simplest way is probably to look though the list of hydrogen bonds and keep a flag if you see it or not.

I don’t entirely understand your second to last paragraph, though my guess to the simple answer is to think carefully about what you actually want to know about. Slightly different things will call for slightly different calculations.

Regarding the printing, while we try to get PyRosetta to print usable things for Rosetta objects (and more recent version should be much better at this), not all of the Rosetta objects have “useful” representations when printed at the Python prompt. This may be the case here, where the objects in your list don’t give you much in the way of useful representations. You may need to manually iterate through the list and print off relevant information yourself.

[Anonymous]

Dear Sir,

Whenever I use this code for storing the info., pyrosetta crashes every time.

p=pose_from_pdb(“C:/cygwin/home/blast-dataset/pdb-files/clean-pdb/%s.pdb”%pdb)
dssp = secstr.Dssp(p)
dssp.insert_ss_into_pose(p)
sec_str = p.secstruct()

pose_hbonds = hb.HBondSet() # imported rosetta.core.scoring.hbonds as hb
hb.fill_hbond_set(p,False,pose_hbonds)
for hbond in range(1,pose_hbonds.nhbonds()+1):
hbond=pose_hbonds.hbond(hbond)
acceptor_residue = hbond.acc_res()
donor_residue = hbond.don_res()
print donor_residue,acceptor_residue

[Anonymous]

When you say crashes, what exactly happens? What message do you get? Which line of the code is it crashing on?

One possibility is the repeated use of the variable “hbonds” – you use it for the loop index as well as the hbond object. I don’t think that would necessarily give you a problem with Python, but I’d recommend changing it just to be safe.

[Anonymous]

actually it stops working and windows reports an error that pyrosetta is not responding, But when I use the command – p.update_residue_neighbors(); it works fine.

---

BHARAT

[Anonymous]

Dear Sir,

I have calculated the handedness of beta-hairpin as per your advice. Here’ the calculation (please see attached fig.)

In [35]: vecN = p.residue(115).xyz(“N”)

In [36]: vecC = p.residue(115).xyz(“C”)

In [37]: vecNC = vecN-vecC

In [38]: vecCA = p.residue(115).xyz(“CA”)

In [39]: vecCB = p.residue(115).xyz(“CB”)

In [40]: vecAB = vecCA – vecCB

In [42]: vecP = p.residue(115).xyz(“CA”)

In [43]: vecF = p.residue(118).xyz(“CA”)

In [44]: vecPF = vecP-vecF

In [45]: print vecAB.cross(vecNC).dot(vecPF)
20.452565332

The structure used for this analysis was 1GFL, residue 22-25.Here, does the positive sign means right-handed twist?? Can u please verify if this calculation matched with the figure ??

That paper (Koga et.al.) also talks about the relation between the chirality of loop and beta-hairpin. I want to know how can I find the torsional energy of the loops using pyrosetta ??

I have a doubt about one of the results mentioned in that paper(Koga et.al.). It says that for the Beta-loop-beta motif, where loop length is 2 prefers the L chirality. Does it mean that the loop twist in direction different from twist of beta-strands ??? I am confused about this point. I have read in many papers that as far the two residue loops are concerned , they mostly have the right-handed twist compatible with the twist of beta-strands. Can u please clarify this point ??

Also, if you look at the attached figure, can you tell whether the chirality of this motif is L or R?? (The figure contains a beta-hairpin , residue 12-36 of pdb-id 1GFL)

[Anonymous]

Dear Sir,

I have been waiting for a long time for your reply … Please find time to kindly respond to my queries ???

[Anonymous]

From what I understand (I’m no expert), the hairpin chirality of the the posted loop would be right handed. (Thumb pointed Calpha to Cbeta on the residue just before the loop, the loop curls in the same direction as the fingers of your right hand.) This does indeed correspond to a positive value in your calculation.

Regarding discrepencies between Koga et al and other papers, I presume the issue may be differing definitions of chirality (The chirality definitions are slightly arbitrary, and other people may be using other definitions), or the more likely case is that they’re talking about the chirality of something different. Again, there’s multiple chiralities at play here, and the chirality of the loop in a hairpin may not be the same as the chirality of the twist of the hairpin or the chirality of the curvature of the hairpin.

Regarding loop torsional energy, your best bet would be to look at the “rama” scores of the residues in the loop, which is a statistical energy based on the Ramachandran phi/psi bin the residue is occupying. You may also want to look at “p_aa_pp” (probability of amino acid given phi and psi), which is a measure of how much a particular amino acid identity likes that phi/psi bin.

[Anonymous]

Dear Sir,

Sorry for late response to your reply. So, it means that my calculation and understanding is correct. Then I still have a doubt, about the word chirality used in the paper. As per the description in Supplementary material of the paper you mentioned (Principles for designing ideal protein structures), it says that :-
“the chirality of a ββ-unit was considered using the vector ! along the axis of
the first strand, and the vector ! perpendicular to ! between the centers of the two strands. Since twisted
strands lead to inaccurate assignments of the chirality, however, we used atom coordinates close to the
loop between the strands for the definition of ! and !: ! is a vector from the N (backbone amide nitrogen)
to C (backbone carbonyl carbon) atoms of the strand residue preceding the connecting loop and ! is a
vector from the Cα atom of the strand residue preceding the loop to the Cα atom of the strand residue
following the loop”.

Using this rule can we say whether the beta-hairpin is right or left hand twisted as the calculation is based entirely on the residue around the turn ??

[Anonymous]

You can for the definition of chirality as used in the paper. It won’t if you’re thinking of chirality in a different sense than the paper is using.

There’s multiple “handednesses” going on in a beta hairpin, and even for a given handedness sense which one is left and which one is right can sometimes depend on which convention (which definition) you’re using (in the cases where one definition hasn’t become a widely accepted standard). If you’re trying to match the handedness sense used by a particular paper, you want to make sure you’re using the definition and doing the calculations the same way the paper does. If you’re not necessarily trying to match a paper, you need to think carefully about why you want to classify things as right handed/left handed – what you’re achieving by splitting them into those two groups – and which handedness and which definition of chirality will help you achieve that goal.

[Anonymous]

Dear Sir,

Thank you for your reply…I guess I am not able to explain my queries clearly, sorry for that ..

My sole objective is to find the direction of the twist(right/left) of the beta-strands. I have already calculated the direction of twist(whether right /left) for the turns. Now, I want to check whether the twist direction of beta-hairpin’s strands match with that of it’s beta-turn.

Now, as per the paper, it says that “The chirality of β-hairpins is determined by the length of the loop between the two strands “. What is the meaning of the word chirality here ?? Does it mean that the word “chirality” refers to right/left handedness of beta-hairpin strands alone. Also, what is the meaning of right/left handed beta strands then, I guess the answer to this should clarify my doubt (ref ??). Just by looking at the strands can we say whether they are right or left hand twisted ??

In the supplementary materials its mentioned that , “the chirality of a ββ-unit was considered using the vector ! along the axis of
the first strand, and the vector ! perpendicular to ! between the centers of the two strands. Since twisted strands lead to inaccurate assignments of the chirality, however, we used atom coordinates close to the loop between the strands for the definition of vec u and vec v . u is a vector from the N (backbone amide nitrogen) to C (backbone carbonyl carbon) atoms of the strand residue preceding the connecting loop and v is a vector from the Cα atom of the strand residue preceding the loop to the Cα atom of the strand residue following the loop”, this is similar to the explanation that you provided me, when I asked for the calculation of handedness of beta-hairpins.

In order to confirm my calculation, I calculated the twist handedness of the beta-hairpins having two different types of turns – type I'(right-handed) and type II(left-handed)??. For both I got the handedness of hairpin (calculated using the above mentioned method)
as right-handed.

For type-I’ eg I used pdb-id 1a1x (pdb residues 94-97) and for type-I eg I used pdb-id 1dqg (pdb-residues 20-23). They both give the angle positive, which means that they are the strands in both cases are right hand twisted.

Also, further check my calculation , I did the same calculation for on of the NMR models from the same paper. I used the pdb id 2LN3 (residues 64-67). In this case also, the handedness comes in right hand direction ?? Since the pdb contained NMR models , I therefore used the first model details for calculation, Here’s the result for that calculation :-

In [10]: vecN = p.residue(64).xyz(“N”)

In [11]: vecC = p.residue(64).xyz(“C”)

In [12]: vecNC = vecN-vecC

In [13]: vecCA = p.residue(64).xyz(“CA”)

In [14]: vecCB = p.residue(64).xyz(“CB”)

In [15]: vecAB = vecCA-vecCB

In [16]: vecP = p.residue(64).xyz(“CA”)

In [17]: vecF = p.residue(67).xyz(“CA”)

In [18]: vecPF = vecP-vecF

In [19]: print vecAB.cross(vecNC).dot(vecPF)
19.063025436

.. Is the calculation correct or not ??. As the paper mentioned that all in all the simulations with the 2 residue turns the chirality was L ?? Pls find the attached figure containing all the beta-hairpins used in the calculation.

Hope I have made my points clear this time … looking forward for your reply

Regards

---

BHARAT

[Anonymous]

Yes Sir, I understand the difference between the type I’&II’ turns. I know this sense of handedness is different from the one described in the paper.

What I want to know that you just looked at the structure at told whether it’s chirality is L/R ??

Sorry for repeatedly asking the same thing , but I am not able to understand that why my calculation’s result doesn’t match with that of the paper. What about my calculation for the first hairpin from the paper, is it correct ?? Also, I don’t understand the meaning of this sentence “which side of the N-terminal beta strand the C-terminal beta strand is on” ??

You also said that “For example, both of the hairpins in the linked figure have the same Koga-handedness, but one is type I’ and one is type II’. ” It means that my calculation doesnot seem to be right. So, here’s the calculation for the second hairpin :-

In [13]: vecN = p.residue(92).xyz(“N”)

In [14]: vecC = p.residue(92).xyz(“C”)

In [15]: vecNC = vecN-vecC

In [16]: vecCA = p.residue(92).xyz(“CA”)

In [17]: vecCB = p.residue(92).xyz(“CB”)

In [18]: vecAB = vecCA-vecCB

In [19]: vecP = p.residue(95).xyz(“CA”)

In [20]: vecP = p.residue(92).xyz(“CA”)

In [21]: vecF = p.residue(95).xyz(“CA”)

In [22]: vecPF = vecP-vecF

In [23]: print vecAB.cross(vecNC).dot(vecPF)
18.693263828

In response to your query “You need to carefully consider what it is you actually want to use the results of the chirality calculation for, in order to pick the correct calculation and get interpretable results.” , I am rephrasing my question again . Can I use the chirality to tell/calculate whether a hairpin is right/left hand twisted ??.. Is there any other reference that you can provide (if possible) in this regard ??

[Anonymous]

Regarding “which side of the N-terminal beta strand the C-terminal beta strand is on”: Imagine lying down on the N-terminal beta strand of the hairpin, with your head at the c-terminal end, your feet toward the N-terminal end, and looking up at the Cbeta atom of the last residue. The loop of the hairpin comes up off of the top of your head, curves around and twists down to lie by your side. Which side is it? The left or the right? That’s what the chirality from the Kogas’ paper is measuring.

I’m at a bit of a loss to provide further help, as I’m not clear as to what you mean by “whether a hairpin is right/left hand twisted”. Twisted how? In what sense? What’s the geometry of the situation? In some sense a chirality is simply a rotation and a direction. What’s being rotated with respect to what? – You say you’re comparing chiralities, but what do you hope to achieve by this comparison? What insight do you intend to obtain? Perhaps by being clear on the big-picture goals, the smaller details will also become clearer.

[Anonymous]

Dear Sir,

I am really thankful to you for helping me clearing my doubts. After such a long discussion on this topic, I am now totally clear about what I need and what actually chirality means. Your last post cleared all my doubts. I once again thank you for lending your time in clearing my concepts.

Regards

---

BHARAT

[Author]

Posts