You are here

HBondEnergy

16 posts / 0 new
Last post
HBondEnergy
#1

Hi,

I want to know that how can I access the HBondEnergy of residue residue pair of the HBondSet(). I checked that there is a function called hbond.energy(), but don't know how to use that ?? But when I use it for the residue in pdb to store information using :
for residue in range(1,p.total_residue()+1):

for hbond in range (1,p_hb.nhbonds()+1):
hbond = p_hb.hbond(hbond)
donor = hbond.don_res()
acceptor = hbond.acc_res()
if residue == donor or residue == acceptor:
donor_residues.append(donor); acceptor_residues.append(acceptor)
donor_array.append(p.residue(donor).atom_name(hbond.don_hatm())),acceptor_array.append(p.residue(acceptor).atom_name(hbond.acc_atm())),energy.append(hbond.energy())

The information gets repeated, if a residue is found . Hope I have made myself clear here.

I also want to calculate the dihedral energy of a residue, how can I calculate that ??

I am working on peptide model , for which I have find h-bond. It works well with full-atom mode but when I apply centroid mode it doesn't show any information. Basically, I need to find the h-bond when the peptide contains only backbone atoms. how can I do that .. currently I am doing this :-

In [30]: q = pose_from_pdb("VNGH.pdb")
core.pack.task: Packer task: initialize from command line()

In [31]: q_hb = hb.HBondSet()

In [32]: q.update_residue_neighbors()

In [33]: hb.fill_hbond_set(q,False,q_hb)

In [34]: q_hb.show(q)
#Dch Dn Dres Da Ach An Ares Aa length AHDang BAHang BAtor weight energy
#A 1 VAL N A 4 HIS O 1.71 134.6 119.0 -115.0 1.000 -0.578
#A 4 HIS N A 1 VAL O 1.91 150.1 103.7 99.1 1.000 -1.104

In [35]: centroid_mode = SwitchResidueTypeSetMover('centroid')

In [36]: centroid_mode.apply(q)
core.chemical.ResidueTypeSet: Finished initializing centroid residue type set. Created 1980 residue types

In [37]: q_hb = hb.HBondSet()

In [38]: hb.fill_hbond_set(q,False,q_hb)

In [39]: q_hb.show(q)
"No output for this command"

-----
BHARAT

Post Situation: 
Wed, 2013-05-15 19:03
bharat_46010

Bharat,

First, I would change the name of the first use of hbond to hbond_index, then use this to get the real Hbond object: p_hb.hbond(hbond_index)
For the energies, you will want to use the hbond.energy() like you found and multiply that by hbond.weight(). This is because different types of hbonds (bb-bb vs bb-sc, ect.) have different weights in the scorefunction.

For dihedral energy, you need to be more specific. Are you looking for a torsional potential? If so, you can create a scorefunction using the mm_std.wts, or just make an empty scorefunction with the mm_twist scoreterm I believe. For more information on this: http://www.plosone.org/article/info:doi/10.1371/journal.pone.0032637

The 'dihedral score' in the standard score12 scorefunction is a statistical potential based on the Ramachandran plot. This is the scoreterm rama.

For individual residue energies, take a look at the Emap. Its mentioned here: http://www.pyrosetta.org/documentation#TOC-9.-ScoreTypes and also in the PyRosetta appendix.

For centroid hbonds, It looks like there are 3 centroid scoreterms for hbonds: hbond_sr_bb and hbond_lr_bb and cen_hb. However, since they are in centroid mode where you don't have true atoms and the true geometry of the hbonds, it probably doesnt not use the HbondSet object.
For this, you will want to construct a scorefunction with hbond_sr_bb and hbond_lr_bb, and use it and the emap to get the info. Checkout GUIs/pyrosetta_toolkit/modules/RegionalScoring.py There may also be a PyRosetta script or demo using the emaps. Your going to be using either score.eval_cd_2b or score.eval_ci_2b function to get the energy between the two residues. I'm not sure if the hbond energies are context dependant or context independent.

Thu, 2013-05-16 09:20
jadolfbr

What about the energies that p_hb.show(p) gives?? Are those energy values weighted or not ?? I have another doubt regarding the name of donor and acceptor atoms. I am attaching my script, if you can have a look and rectify according to my need. What I am doing is to find the hbond between first and last residue of turn and their energies and distance between the NH and CO.

Thu, 2013-05-16 09:51
bharat_46010

HBond.show() should show the *unweighted energy* along with the associated weight as the second to last and last item on the line. Which one is which depends on which function that call is wrapping. My guess is the energy is second-to-last and the weight is last. It should be easy to tell apart. The weight will probably be a positive number between zero and one (although could be greater), and the energy can be any number, but will likely be negative.

The donor atom is the one that's bonded to the hydrogen in the hydrogen bond, and the acceptor is the atom that has the lone pairs for the hydrogen bond. For backbone hydrogen bonds, the donor atom name will be "N" and the acceptor atom name will be "O".

Thu, 2013-05-16 11:09
rmoretti

Hi J,

If possible can u pls brief about scorefunction using the mm_std.wts for calculating the dihedral/torsional energy ??

Sun, 2013-05-19 20:09
bharat_46010

Hey Bharat,

I don't know much about it besides what I have read in the paper. However, I read it to get a general idea on the the scorefunction and not specifics, so you will have to read it and see for yourself if mm_twist is what you need.

For the actual code, you can use the emap or look at the classes MMTorsionScore and MMTorsionLibrary in rosetta.core.scoring.mm. I have emailed Doug Renfrew to see if he can help more. The MMTorsionScore takes a Key4Tuple in order to get the score of one torsion - I'm not sure if this is available in PyRosetta, so I have also emailed Sergey to get his input.

-J

Sun, 2013-05-19 20:47
jadolfbr

Okay and thank you for your prompt response. I will wait for their reply ...

Sun, 2013-05-19 22:09
bharat_46010

I also want to calculate the dihedral energy of a residue, how can I calculate that ??

Like Jared said, you will need to be MUCH more specific about what you mean when you say "dihedral energy" as it is completely unclear what you mean.
Do you mean the torsional component of the orientation dependent hydrogen bonding term?
Do you mean the knowledge-based back bone torsion energy (we usually call this the Ramachandran term)?
Do you mean the knowledge-based side chain torsion energy (we usually call this the Dunbrack term)?
Do you mean the CHARMM-based torsion term called mm_twist used in the mm_std scoring function?
Do you mean something else entirely?

I am working on peptide model , for which I have find h-bond. It works well with full-atom mode but when I apply centroid mode it doesn't show any information. Basically, I need to find the h-bond when the peptide contains only backbone atoms. how can I do that ..

First of all, centroid mode doesn't contain only backbone atoms, there is an additional "atom" (sometimes called a bead) that represents the interaction sphere of the side chain. If your hydrogen bond involves the side chain, then those atoms will not be present. Secondly, if you are interested in specific hydrogen bonds, why are you even thinking of working in centroid mode? The energy function usually used in centroid mode has completely different terms to model the interactions within and between residues.

If possible can u pls brief about scorefunction using the mm_std.wts for calculating the dihedral/torsional energy ??
The mm_std scoring function was created to score models with non-canonical amino acids. It removes the knowledge-based terms that are evaluated based on residue-type and replaces them with some molecular mechanics terms, notably the CHARMM27 torsion and LJ terms. It is described in the link Jared provided. The mm_twist (as it it called in Rosetta) term is evaluated on all sets of 4 bonded atoms. This term is not usually used in Rosetta so I doubt it is what you want.

Mon, 2013-05-20 07:57
renfrew

Dear Sir,

Referring to this paper "http://www.nature.com/nature/journal/v491/n7423/full/nature11600.html#f1", I asked about how to calculate the torsional energy as I want to check the similar torsional energies for my data-set of two residue loops. Also, if possible it would be interesting to find the torsional component of orientation dependent h-bonding term...can u provide some ref. for this ??

Tue, 2013-05-21 01:40
bharat_46010

Bharat,

The paper mentions Rama in the supplementary information (pg 7) when talking about torsional energies:

The torsion energies for L- and R- hairpins were calculated for each loop length using Rosetta with the
rama score term 1.0. There are no R-hairpins at the loop length 2. For 2 and 3 residue loops, L-hairpins
have lower energy, and for 5 residue loop, R-hairpins have lower energy.

To use this: create a scorefunction + set the rama weight to 1.0
scorefxn = ScoreFunction()
scorefxn.set_weight(rama, 1.0)

Score the pose, and use the emap to get individual residue energies your interested in.

One way to do that is this (after you score the pose)
In [24]: emap = p.energies().residue_total_energies(10)

In [25]: emap[rama]
Out[25]: -0.03849164103910052

I am moving some functions from HBond.show method to their own methods to allow access to this kind of information. Until then: You can calculate it like this:
base_atomno = (pose.residue(hbond.acc_res()).atom_base(acc_atm()));
bbase_atomno =pose.residue(hbond.acc_res()).atom_base(base_atomno));
Bxyz = pose.residue(hbond.acc_res()).xyz(pose.residue(base_atomno)));
BBxyz = pose.residue(hbond.acc_res()).xyz(bbase_atomno));
Hxyz = pose.residue(hbond.don_res()).xyz(hbond.don_hatm()));
Axyz = pose.residue(hbond.acc_res()).xyz(hbond.acc_atm()));

Each of these are xyz vectors. To calculate the dihedral there is a function rosetta.numeric.dihedral_degrees ( dihedral_degrees(BBxyz, Bxyz, Axyz, Hxyz) )that you could use, but my PyRosetta doesn't seem to have it. It may be a newer function. If you can't access this, just calculate the dihedral by hand. I don't have the math on hand for this, so you will need to figure this one out on your own.

Tue, 2013-05-21 08:13
jadolfbr

Thank you, Sir.

Could you please comment on analyzing the rama score in relation with p_aa_pp. As far rama score is concerned, the lower the better , if I am correct ??

Tue, 2013-05-21 19:55
bharat_46010

Lower the better for Rama,
p_aa_pp is basically the statistical opposite of rama. It is the probability of amino acid given phi/psi; while rama is the probability of phi/psi given amino acid. To get the energy, both of these are -log(prob). This paper: http://www.biostat.wisc.edu/bmi576/papers/rohl.pdf covers most of the default scorefunction.

Our lab doesn't think p_aa_pp is necessary in the scorefunction other than for design - but from the general knowledge of the community, it seems that if you remove it from regular modelling runs, everything gets worse. So, for now at least, it is in the default scorefunction.

Wed, 2013-05-22 07:21
jadolfbr

Hi,

How does rosetta identify the hydrogen bonded pair in beta-hairpins. Is it based on the DSSP rule that a NH-CO and CO-NH bonds between donor and acceptor residue should be present ?

Also, I want to know how can I find hydrogen-bonded and non-hydrogen bonded pairs of a beta-hairpin . I wrote a script to find HB and NHB pairs but it works only when the bhp is symmetric.

Pls find the attached figure in order to understand my query

Sun, 2013-06-16 04:25
bharat_46010

Rosetta identifies hydrogen bonds by atom geometry. If the donor, hydrogen and acceptor atoms are in the correct geometry to make a hydrogen bond, Rosetta counts it as a hydrogen bond. As Rosetta has an energy term for hydrogen bonds, usually the calling if something is "really" a hydrogen bond (rather than just something where the atoms are in vaguely the correct region) is usually made by an energy cutoff. If the core.scoring.hbonds.HBond object is initialized properly, you can get the energy associated with each hydrogen bond from it.

I don't believe there's any special consideration for the hydrogen bonds in beta-hairpins. (At least for full atom mode - there might be something in centroid mode, but I get the sense there isn't.)

Probably the easiest way to find hydrogen bonds in a hairpin is probably to initialize an core.scoring.hbonds.HBondSet object with the pose, then iterate through the HBonds in the set, pulling out the ones which are between residues in the hairpin and more specifically the ones between backbone atoms ( there's don_hatm_is_protein_backbone() and acc_atm_is_protein_backbone() methods in the HBond object for that purpose(). You can further refine this search by whatever additional energy or geometry criteria you wish to use.

Sun, 2013-06-16 12:41
rmoretti

Yes, Sir I know these methods. But query is related to find the pair of non-hydrogen bonded pair also. I have posted the figure of the case when the beta-hairpin incorporates a bulge and this is where I don't understand how to count for non-hydrogen bonded pairs.

Sorry to ask this question here, but I have read many papers about pairing preferences in anti-parallel beta hairpins. Almost all of them mention that they use DSSP conventions find the HB and NHB pairs. I tried to understand it pairing preferences from DSSP output , but I was not able to find out how the HB and NHB pairings were identified. For eg, I have attached a sample Dssp output of a bhp. Can anybody help me out in this regard ??

Sun, 2013-06-16 19:12
bharat_46010

I'm not familiar with the DSSP conventions. Unless someone else on the board is, your best bet is probably to track down the original DSSP references and see how they describe looking for "non-hydrogen bonded pairs".

From the bhp-problem.png figure, it looks like the hydrogen bonded pair might require that the N and O hydrogen bonds be going to the same residue. For the non-hydrogen bonded pair on the left, there's no hydrogen bonds, and for the bulge, there's hydrogen bonds, but the N and O hydrogen bonds are going to different residues.

But if you want to match a particular convention, you'll need to look at the original reference to see how they define things.

Mon, 2013-06-17 10:40
rmoretti