eval_ci_2b() with hydrogen bonds

[Anonymous]

I am trying to compute the residue-residue pairwise energy for several residues inside a protein structure. I found the method scorefxn.eval_ci_2b() can compute the pairwise energy between residues for a given scorefxn. This works for the LJ, Solvation and electrostatic potentials. Also as epxected, it ignores intra-residue interactions and backbone parameters when computing with eval_ci_2b(). However, I was wondering why the Hydrogen bonded potential is also not included? Is there an easy way to include the hydrogen bond energy in the pairwise energy?

Not sure if this is helpful to even see, but I have a minimum example of what I am trying to do below. Not sure if it would run as it was cut/pasted from various sections of my code, but it should get at the gist of what I was trying to do.

import numpy as np

import pyrosetta as pyr

from pyrosetta import Pose

from pyrosetta.toolbox import get_hbonds

from pyrosetta import teaching as pyrt



pose = pyr.pose_from_pdb(my_pdb_file)



stuff = "(fa_atr 1) (fa_rep 0.55) (fa_sol 1) (fa_intra_rep 0.005) (fa_intra_sol_xover4 1) (lk_ball_wtd 1) (fa_elec 1) (pro_close 1.25) (hbond_sr_bb 1) (hbond_lr_bb 1) (hbond_bb_sc 1) (hbond_sc 1) (dslf_fa13 1.25) (omega 0.4) (fa_dun 0.7) (p_aa_pp 0.6) (yhh_planarity 0.625) (ref 1) (rama_prepro 0.45)"



things = stuff.strip().split(")")

weights = []

order = []

for i_thing_count,thing in enumerate(things):

if len(thing) > 0:

if i_thing_count == 0:

useful_stuff = thing[1:].split()

else:

useful_stuff = thing[2:].split()

order.append(eval("pyrt." + useful_stuff[0]))

weights.append(float(useful_stuff[1]))



assert len(weights) == len(order)

scorefxn = pyrt.ScoreFunction()

for thing, wt in zip(order,weights):

scorefxn.set_weight(thing, wt)



hbond_set = get_hbonds(pose)

pair_E = np.zeros((nresidues,nresidues))

for contact in these_contacts:



idx = contact[0]

jdx = contact[1]

emap = pyrt.EMapVector()

scorefxn.eval_ci_2b(pose.residue(idx), pose.residue(jdx), pose, emap)

this_E = 0.

for thing,wt in zip(order,weights):

if idx == 8 and jdx == 20:

print thing

print emap[thing]

this_E += emap[thing] * wt

pair_E[idx-1, jdx-1] = this_E

 

[Anonymous]

Hydrogen bonding is a context-dependent two body energy, so it wouldn’t be seen with scorefxn.eval_ci_2b()  (the “ci” means “context independent”) — However I think it should be availible with scorefxn.eval_cd_2b().

If you’re playing around with the energy function on a lower-level like you are, I’d highly encourage you to read through Leaver-Fay et al. https://doi.org/10.1016/B978-0-12-381270-4.00019-6 to get an overview of the structure, if you haven’t already.

[Anonymous]

Thanks for the reply. I just tried using calc_cd_2b(), but it unfortunately did not work either. I have an example script in a .txt file attached as well as the pdb file I’ve been testing with.

System Specs:

OS: Ubuntu

Python: Python2.7.12

PyRosetta: Downloaded binary from the website http://www.pyrosetta.org/dow .

Also tested this script on a RedHat system where PyRosetta was compiled from source with similar results. Any advice you can give on how to proceed will be appreciated. Many Thanks.

[Author]

Posts