make sure we don't collect contact atoms of the query residue

deeprank/models/atom.py

@@ -18,5 +18,5 @@ def __lt__(self, other):
         return self.id < other.id
 
     def __repr__(self):
-        return "Atom {} ({} {} {}) at {}".format(self.id, self.chain_id, self.residue, self.name, self.position)
+        return "Atom {} ({}) from {} {} at {}".format(self.id, self.chain_id, self.residue, self.name, self.position)
 

deeprank/models/residue.py

@@ -1,10 +1,10 @@
 
 
 class Residue:
-    def __init__(self, number, name, atoms=[]):
+    def __init__(self, number, name):
         self.number = number
         self.name = name
-        self.atoms = atoms
+        self.atoms = []
 
     @property
     def chain_id(self):

deeprank/operate/pdb.py

@@ -8,18 +8,6 @@
 
 _log = logging.getLogger(__name__)
 
-def get_squared_distance(position1, position2):
-    """Get squared euclidean distance, this is computationally cheaper that the euclidean distance
-
-        Args:
-            position1 (numpy vector): first position
-            position2 (numpy vector): second position
-
-        Returns (float): the squared distance
-    """
-
-    return numpy.sum(numpy.square(position1 - position2))
-
 
 def get_distance(position1, position2):
     """ Get euclidean distance between two positions in space.
@@ -31,7 +19,7 @@ def get_distance(position1, position2):
         Returns (float): the distance
     """
 
-    return numpy.sqrt(get_squared_distance(position1, position2))
+    return numpy.sqrt(numpy.sum(numpy.square(position1 - position2)))
 
 
 def get_atoms(pdb2sql):
@@ -81,9 +69,6 @@ def get_residue_contact_atom_pairs(pdb2sql, chain_id, residue_number, max_intera
         Returns ([Pair(int, int)]): pairs of atom numbers that contact each other
     """
 
-    # Square the max distance, so that we can compare it to the squared euclidean distance between each atom pair.
-    squared_max_interatomic_distance = numpy.square(max_interatomic_distance)
-
     # get all the atoms in the pdb file:
     atoms = get_atoms(pdb2sql)
 
@@ -97,12 +82,15 @@ def get_residue_contact_atom_pairs(pdb2sql, chain_id, residue_number, max_intera
     contact_atom_pairs = set([])
     for atom in atoms:
 
+        # Check that the atom is not one of the residue's own atoms:
+        if atom.chain_id == chain_id and atom.residue.number == residue_number:
+            continue
+
         # Within the atom iteration, iterate over the atoms in the residue:
         for residue_atom in residue_atoms:
 
-            # Check that the two atom numbers are not the same and check their distance:
-            if atom != residue_atom and \
-                    get_squared_distance(atom.position, residue_atom.position) < squared_max_interatomic_distance:
+            # Check that the atom is close enough:
+            if get_distance(atom.position, residue_atom.position) < max_interatomic_distance:
 
                 contact_atom_pairs.add(Pair(residue_atom, atom))
 

test/operate/test_pdb.py

@@ -35,23 +35,41 @@ def _find_atom(atoms, chain_id, residue_number, name):
     return matching_atoms[0]
 
 
+def _find_residue(atoms, chain_id, residue_number):
+    matching_atoms = [atom for atom in atoms if atom.chain_id == chain_id and
+                                                atom.residue.number == residue_number]
+
+    assert len(matching_atoms) > 0, "Expected at least one matching atom, got zero"
+
+    return matching_atoms[0].residue
+
+
 def test_residue_contact_atoms():
 
     pdb_path = os.path.join(pkg_resources.resource_filename(__name__, ''),
                             "../1AK4/native/1AK4.pdb")
+    chain_id = 'D'
+    residue_number = 145
 
     try:
         pdb = pdb2sql(pdb_path)
 
         atoms = get_atoms(pdb)
+        query_residue = _find_residue(atoms, chain_id, residue_number)
 
-        contact_atom_pairs = get_residue_contact_atom_pairs(pdb, 'D', 145, 8.5)
+        contact_atom_pairs = get_residue_contact_atom_pairs(pdb, chain_id, residue_number, 8.5)
 
         # List all the atoms in the pairs that we found:
         contact_atoms = set([])
         for atom1, atom2 in contact_atom_pairs:
             contact_atoms.add(atom1)
             contact_atoms.add(atom2)
+
+            # Be sure that this atom pair does not pair two atoms of the same residue:
+            ok_(atom1.residue != atom2.residue)
+
+            # Be sure that one of the atoms in the pair is from the query residue:
+            ok_(atom1.residue == query_residue or atom2.residue == query_residue)
     finally:
         pdb._close()
 
@@ -62,7 +80,7 @@ def test_residue_contact_atoms():
     distant = _find_atom(atoms, 'C', 134, 'OE2')  # this residue is very far away
 
     # Check that the close residue is present in the list and that the distant residue is absent.
+    # Also, the function should not pair a residue with itself.
 
     ok_(neighbour in contact_atoms)
-
     ok_(distant not in contact_atoms)