Baker-Hubbard Hydrogen Bond IdentificationΒΆ

In [1]:
from __future__ import print_function
%matplotlib inline
import matplotlib.pyplot as plt
import itertools
import mdtraj as md
import mdtraj.testing

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-1-ce89838d5836> in <module>()
      4 import itertools
      5 import mdtraj as md
----> 6 import mdtraj.testing

/home/travis/miniconda3/envs/py27/lib/python2.7/site-packages/mdtraj/testing/__init__.py in <module>()
     23 
     24 from __future__ import print_function, division
---> 25 from mdtraj.testing.testing import *
     26 from mdtraj.testing.docstrings import *

/home/travis/miniconda3/envs/py27/lib/python2.7/site-packages/mdtraj/testing/testing.py in <module>()
     37   assert_raises, assert_string_equal, assert_warns)
     38 from numpy.testing.decorators import skipif, slow
---> 39 from nose.tools import ok_, eq_, raises
     40 from nose import SkipTest
     41 from pkg_resources import resource_filename

ImportError: No module named nose.tools

Load up some example data. This is a little 28 residue peptide

In [2]:
t = md.load_pdb('http://www.rcsb.org/pdb/files/2EQQ.pdb')
print(t)

<mdtraj.Trajectory with 20 frames, 423 atoms, 28 residues, without unitcells>

/home/travis/miniconda3/envs/py27/lib/python2.7/site-packages/mdtraj/formats/pdb/pdbfile.py:193: UserWarning: Unlikely unit cell vectors detected in PDB file likely resulting from a dummy CRYST1 record. Discarding unit cell vectors.
  warnings.warn('Unlikely unit cell vectors detected in PDB file likely '

md.baker_hubbard idenfies hydrogen bonds baced on cutoffs for the Donor-H...Acceptor distance and angle. The criterion employed is $\theta > 120$ and $r_\text{H...Acceptor} < 2.5 A$ in at least 10% of the trajectory. The return value is a list of the indices of the atoms (donor, h, acceptor) that satisfy this criteria.

In [3]:
hbonds = md.baker_hubbard(t, periodic=False)
label = lambda hbond : '%s -- %s' % (t.topology.atom(hbond[0]), t.topology.atom(hbond[2]))
for hbond in hbonds:
    print(label(hbond))

GLU1-N -- GLU1-OE2
GLU1-N -- GLU1-OE1
GLY6-N -- SER4-O
CYS7-N -- GLY5-O
TYR11-N -- VAL8-O
MET12-N -- LYS20-O
ARG13-NH1 -- TYR11-O
THR14-N -- ARG18-O
ASP16-N -- ASP16-OD1
GLY17-N -- THR14-O
ARG18-N -- THR14-OG1
ARG18-NE -- ASP16-OD2
LYS20-N -- MET12-O
THR22-N -- GLY10-O
THR14-OG1 -- ASP16-OD1
THR28-OG1 -- ILE27-O

Let's compute the actual distances between the donors and acceptors

In [4]:
da_distances = md.compute_distances(t, hbonds[:, [0,2]], periodic=False)

And plot a histogram for a few of them

In [5]:
color = itertools.cycle(['r', 'b', 'gold'])
for i in [2, 3, 4]:
    plt.hist(da_distances[:, i], color=next(color), label=label(hbonds[i]), alpha=0.5)
plt.legend()
plt.ylabel('Freq');
plt.xlabel('Donor-acceptor distance [nm]')
Out[5]:
<matplotlib.text.Text at 0x2b0ac82205d0>
In [6]:

(hbonds.ipynb; hbonds_evaluated.ipynb; hbonds.py)

Versions