Merge pull request #17 from weitzner/master

Bidentate placement and tests

.bumpversion.cfg

@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 0.2.3
+current_version = 0.3.0
 commit = True
 tag = True
 

.gitignore

@@ -60,3 +60,20 @@ target/
 
 # pyenv python configuration file
 .python-version
+
+# Jupyter Notebook
+.ipynb_checkpoints
+*/.ipynb_checkpoints/*
+
+# swap
+.sw[a-p]
+.*.sw[a-p]
+# session
+Session.vim
+# temporary
+.netrwhist
+*~
+tags
+
+local/
+

privileged_residues/__init__.py

@@ -2,7 +2,7 @@
 
 """Top-level package for Privileged Residues."""
 from .privileged_residues import *
-from . import hbond_ray_pairs, position_residue, process_networks
+from . import bidentify, hbond_ray_pairs, position_residue, process_bidentates, process_networks
 
 __author__ = """Brian D. Weitzner"""
 __email__ = '[email protected]'

privileged_residues/bidentify.py

@@ -0,0 +1,204 @@
+import os
+import numpy as np
+import itertools
+
+# The following packages are not pip-installable
+# The import calls are wrapped in a try/except block
+try:
+    import pyrosetta
+    from pyrosetta.rosetta.core.id import AtomID
+except ImportError:
+    print('Module "pyrosetta" not found in the current environment! '
+          'Go to http://www.pyrosetta.org to download it.')
+    pass
+
+
+from . import hbond_ray_pairs
+from . import position_residue as pr
+from . import process_networks as pn
+
+# the object of this submodule is to scan over all surface residues and 
+# identify pairs of rays to look up in each hash table
+def look_up_interactions(pairs_of_rays, ht, cart_resl, ori_resl, cart_bound):
+    """
+    """
+    hits = []
+    for r1, r2 in pairs_of_rays:
+        hashed_rays = int(hbond_ray_pairs.hash_rays(r1, r2))
+        try:
+            positioning_info = ht[hashed_rays]
+            print(positioning_info)
+        except KeyError:
+            continue
+
+        ray_frame = hbond_ray_pairs.get_frame_for_rays(r1, r2)
+        # now positioning_info is a set of possible values
+        fxnl_grps = list(pn.fxnl_groups.keys())
+
+        for pos_info in positioning_info:
+            xform = None
+            try:
+                xform = pr.get_ht_from_table(pos_info[1],
+                                             cart_resl,
+                                             ori_resl,
+                                             cart_bound)
+            except:
+                continue
+
+            rsd = fxnl_grps[pos_info[0]]
+            pos_grp = pn.fxnl_groups[rsd]
+
+            p = pyrosetta.Pose()
+            pyrosetta.make_pose_from_sequence(p, 'Z[{}]'.format(rsd), 'fa_standard')
+            coords = [np.array([*p.residues[1].xyz(atom)]) for atom in pos_grp.atoms]
+            c = np.stack(coords)
+
+            pos_frame = hbond_ray_pairs.get_frame_for_coords(c)
+            xf = np.dot(np.dot(ray_frame, xform), np.linalg.inv(pos_frame))
+            pr.transform_pose(p, xf)
+            hits.append(p)
+    return hits    
+
+
+def look_for_sc_bb_bidentates(p):
+    """
+    """
+
+    pairs_of_rays = []
+    for i in range(1, len(p.residues) + 1):
+        # look at the amino and carboxyl rays, swap the order and look ahead
+        # and back to see if the previous residue's carboxyl or the next 
+        # residue's amino ray can be used to place a bidentate interaction
+
+
+        rsd = p.residues[i]
+
+        # ray reference atoms in backbone
+        H = rsd.attached_H_begin(rsd.atom_index("N"))
+        N = rsd.atom_base(H)
+        O = rsd.atom_index("O")
+        C = rsd.atom_base(O)
+
+        n_ray = hbond_ray_pairs.create_ray(rsd.xyz(H), rsd.xyz(N))
+        c_ray = hbond_ray_pairs.create_ray(rsd.xyz(O), rsd.xyz(C))
+
+        pairs_of_rays += [(n_ray, c_ray), (c_ray, n_ray)]
+
+        if i != 1:
+            prev_rsd = p.residues[i - 1]
+
+            O = prev_rsd.atom_index("O")
+            C = prev_rsd.first_adjacent_heavy_atom(O)
+
+            prev_c = hbond_ray_pairs.create_ray(prev_rsd.xyz(O), 
+                                                prev_rsd.xyz(C))
+            pairs_of_rays += [(n_ray, prev_c), (prev_c, n_ray)]
+
+        if i != len(p.residues):
+            next_rsd = p.residues[i + 1]
+            H = next_rsd.attached_H_begin(rsd.atom_index("N"))
+            N = next_rsd.first_adjacent_heavy_atom(H)
+            next_n = hbond_ray_pairs.create_ray(next_rsd.xyz(H), 
+                                                next_rsd.xyz(N))
+            pairs_of_rays += [(next_n, c_ray), (c_ray, next_n)]
+    return pairs_of_rays
+
+def look_for_sc_scbb_bidentates(p):
+    """
+    """
+
+    pairs_of_rays = []
+    for rsd in p.residues:
+        H = rsd.attached_H_begin(rsd.atom_index("N"))
+        N = rsd.atom_base(H)
+        O = rsd.atom_index("O")
+        C = rsd.atom_base(O)
+
+        n_ray = hbond_ray_pairs.create_ray(rsd.xyz(H), rsd.xyz(N))
+        c_ray = hbond_ray_pairs.create_ray(rsd.xyz(O), rsd.xyz(C))
+
+        for i in rsd.Hpos_polar_sc():
+            ray = hbond_ray_pairs.create_ray(rsd.xyz(i), rsd.xyz(rsd.atom_base(i)))
+
+            pairs_of_rays += [(n_ray, ray), (ray, n_ray)]
+            pairs_of_rays += [(c_ray, ray), (ray, c_ray)]
+
+        for i in rsd.accpt_pos_sc():
+            ray = hbond_ray_pairs.create_ray(rsd.xyz(i), rsd.xyz(rsd.atom_base(i)))
+
+            pairs_of_rays += [(n_ray, ray), (ray, n_ray)]
+            pairs_of_rays += [(c_ray, ray), (ray, c_ray)]
+
+    return pairs_of_rays
+
+def look_for_sc_sc_bidentates(p):
+    """
+    """
+
+    pairs_of_rays = []
+    for rsd in p.residues:
+        for (i, j) in itertools.permutations(list(rsd.Hpos_polar_sc()) + list(rsd.accpt_pos_sc()), 2):
+            if i == j:
+                continue
+
+            first = hbond_ray_pairs.create_ray(rsd.xyz(i), rsd.xyz(rsd.atom_base(i)))
+            second = hbond_ray_pairs.create_ray(rsd.xyz(j), rsd.xyz(rsd.atom_base(j)))
+
+            pairs_of_rays += [(first, second)]
+
+    return pairs_of_rays
+
+def look_up_connected_network(p):
+    """
+    """
+
+    pairs_of_rays = []
+
+    hbondset = hbond_ray_pairs.identify_bonded_pairs(p, hbond_ray_pairs.find_hbonds(p))
+
+    for hbond in hbondset.hbonds():
+        don_res = p.residue(hbond.don_res())
+        acc_res = p.residue(hbond.acc_res())
+
+        don_list = list(don_res.Hpos_polar_sc()) + list(don_res.accpt_pos_sc())
+        acc_list = list(acc_res.Hpos_polar_sc()) + list(acc_res.accpt_pos_sc())
+
+        for (j, k) in itertools.product(don_list, acc_list):
+            if (j == hbond.don_hatm() or k == hbond.acc_atm()):
+                continue
+
+            print(*((y, x.atom_name(x.atom_base(y))) for (x, y) in zip([don_res, acc_res], [j, k])))
+
+            first = hbond_ray_pairs.create_ray(don_res.xyz(j), don_res.xyz(don_res.atom_base(j)))
+            second = hbond_ray_pairs.create_ray(acc_res.xyz(k), acc_res.xyz(acc_res.atom_base(k)))
+
+            pairs_of_rays += [(first, second)]
+            pairs_of_rays += [(second, first)]
+
+    return pairs_of_rays
+
+def find_hashable_positions(p, ht):
+    # look up each type of privileged interaction
+    #
+    hits = []
+    pairs_of_rays = []
+
+    # networks:
+    #     iterate over pairs of residues with side chains that are 
+    #     already hydrogen bonded to one another
+    # Sc_Sc bidentates:
+    #     iteratate over every residue, throw out non-bidentate capable 
+    #     residues
+    pairs_of_rays += look_for_sc_sc_bidentates(p)
+    # Sc_ScBb bidentates:
+    #     iterate over pairs of neighboring residues on the surface
+    pairs_of_rays += look_for_sc_scbb_bidentates(p)
+    # Sc_Bb bidentates:
+    #     iterate over pairs of neighboring residues on the surface
+    pairs_of_rays += look_for_sc_bb_bidentates(p)
+    # hits += look_up_interactions(pairs_of_rays, ht)
+    pairs_of_rays += look_up_connected_network(p)
+
+    # tmp
+    return pairs_of_rays
+

privileged_residues/hbond_ray_pairs.py

@@ -135,9 +135,9 @@ def create_ray(center, base):
         2. The point at which the ray is centered is at `center`.
 
     Args:
-        center (numpy.array): A (2, 3) array representing the
+        center (numpy.array): A (1, 3) array representing the
             coordinate at which to center the resulting ray.
-        base (numpy.array): A (2, 3) array representing the base used
+        base (numpy.array): A (1, 3) array representing the base used
             to determine the direction of the resulting ray.
 
     Returns:
@@ -147,12 +147,12 @@ def create_ray(center, base):
     direction = center - base
     direction /= np.linalg.norm(direction)
 
-    ray = np.empty((2, 4))
-    ray[0][:-1] = center
-    ray[0][-1] = 1.  # point!
+    ray = np.empty((4, 2))
+    ray[:-1, 0] = center
+    ray[-1, 0] = 1.  # point!
 
-    ray[1][:-1] = direction
-    ray[1][-1] = 0.  # direction!
+    ray[:-1, 1] = direction
+    ray[-1, 1] = 0.  # point!
 
     return ray
 
@@ -256,9 +256,9 @@ def hash_rays(r1, r2, resl=2, lever=10):
         `r1` and `r2` must be the same shape.
 
     Args:
-        r1 (numpy.array): An (N, 2, 4) array representing N rays in
+        r1 (numpy.array): An (N, 4, 2) array representing N rays in
             space each with a point and a unit direction.
-        r2 (numpy.array): An (N, 2, 4) array representing N rays in
+        r2 (numpy.array): An (N, 4, 2) array representing N rays in
             space each with a point and a unit direction.
 
     Returns:
@@ -298,7 +298,7 @@ def get_frame_for_coords(coords):
 
     assert(coords.shape == (3, 3))
 
-    frame = np.empty((4, 4))
+    frame = np.zeros((4, 4))
     frame[:3, 3] = coords[2]
     frame[3, 3] = 1.
 
@@ -309,7 +309,7 @@ def get_frame_for_coords(coords):
     y_hat = coords[0] - coords[1]
     y_hat -= np.dot(y_hat, z_hat) * z_hat
     y_hat /= np.linalg.norm(y_hat)
-    assert_allclose(np.dot(y_hat, z_hat), 0.)
+    assert_allclose(np.dot(y_hat, z_hat), 0., atol=1E-10)
 
     frame[:3, 1] = y_hat
 
@@ -332,20 +332,20 @@ def get_frame_for_rays(r1, r2):
         the point of r2 to lie in the xy-plane.
 
     Args:
-        r1 (numpy.array): A (2, 4) array representing a ray in space
+        r1 (numpy.array): A (4, 2) array representing a ray in space
             with a point and a unit direction.
-        r2 (numpy.array): A (2, 4) array representing a ray in space
+        r2 (numpy.array): A (4, 2) array representing a ray in space
             with a point and a unit direction.
 
     Returns:
         numpy.array: A (4, 4) array representing the coordinate frame
         as a homogenous transform.
     """
-    trans = r1[0]
-    x_hat = r1[1]
+    trans = r1[:, 0]
+    x_hat = r1[:, 1]
     assert_allclose(np.linalg.norm(x_hat), 1.)  # make sure this is unity
 
-    xy_pojnt = r2[0] - trans
+    xy_pojnt = r2[:, 0] - trans
     y_component = xy_pojnt - (np.dot(xy_pojnt, x_hat) * x_hat)
     y_hat = y_component / np.linalg.norm(y_component)
 

privileged_residues/position_residue.py

@@ -84,8 +84,8 @@ def get_ht_from_table(bin, cart_resl, ori_resl, cart_bound):
     return xh.get_center([bin])['raw'].squeeze()
 
 
-def tranform_pose(p, xform):
-    """Tranform the atomic coordinates of Pose by a specified
+def transform_pose(p, xform):
+    """Transform the atomic coordinates of Pose by a specified
     homogenous tranform.
 
     Args:
@@ -96,8 +96,7 @@ def tranform_pose(p, xform):
     coords = []
     for i in range(1, p.size() + 1):
         for j in range(1, p.residues[i].natoms() + 1):
-            c = np.ones(4)
-            c[:3] = np.array([*p.residues[i].xyz(j)])
+            c = np.array([*p.residues[i].xyz(j)] + [1])
             coords.append(c)
 
     new_coords = np.dot(xform, np.stack(coords).T)
@@ -109,5 +108,6 @@ def tranform_pose(p, xform):
             p.residues[i].atom(j).xyz(xyzVec(x, y, z))
             tot_atoms += 1
 
+
 if __name__ == '__main__':
     print('Don\'t execute me, bruh.')