Exclude backbone oxygens when fitting missing side chains

pdbfixer/pdbfixer.py

@@ -216,13 +216,13 @@ def _overlayPoints(points1, points2):
 
     Parameters
     ----------
-    points1 (numpy array of openmm.unit.Quantity with units compatible with distance) - reference set of coordinates
-    points2 (numpy array of openmm.unit.Quantity with units compatible with distance) - set of coordinates to be rotated
+    points1 (list of Vec3, in nm) - reference set of coordinates
+    points2 (list of Vec3, in nm) - set of coordinates to be rotated
 
     Returns
     -------
     translate2 - vector to translate points2 by in order to center it
-    rotate - rotation matrix to apply to centered points2 to map it on to points1
+    rotate - proper rotation matrix (det = +1) to apply to centered points2 to map it on to points1
     center1 - center of points1
 
     Notes
@@ -234,7 +234,7 @@ def _overlayPoints(points1, points2):
     if len(points1) == 0:
         return (mm.Vec3(0, 0, 0), np.identity(3), mm.Vec3(0, 0, 0))
     if len(points1) == 1:
-        return (points1[0], np.identity(3), -1*points2[0])
+        return (-1*points2[0], np.identity(3), points1[0])
 
     # Compute centroids.
 
@@ -251,10 +251,15 @@ def _overlayPoints(points1, points2):
             for j in range(3):
                 R[i][j] += y[i]*x[j]
 
-    # Use an SVD to compute the rotation matrix.
+    # Use an SVD to compute the rotation matrix.  If the best fit is a reflection, flip the
+    # direction with the smallest singular value (the last one) so the result is a proper
+    # rotation; otherwise nearly coplanar points can map a template onto its mirror image.
 
-    (u, s, v) = lin.svd(R)
-    return (-1*center2, np.dot(u, v).transpose(), center1)
+    (u, _, vh) = lin.svd(R)
+    d = np.sign(lin.det(u)*lin.det(vh))
+    rotate = np.dot(np.dot(u, np.diag([1, 1, d])), vh).transpose()
+
+    return (-1*center2, rotate, center1)
 
 def _dihedralRotation(points, angle):
     """Given four points that form a dihedral, compute the matrix that rotates the last point around the axis to
@@ -647,10 +652,13 @@ def _addAtomsToTopology(self, heavyAtomsOnly, omitUnknownMolecules):
 
                     template = self._getTemplate(residue.name)
                     atomPositions = dict((atom.name, self.positions[atom.index]) for atom in residue.atoms())
+                    # Backbone oxygens depend on psi and can tilt the fit enough to put a
+                    # missing side chain on the wrong side of the N-CA-C plane.
+                    skip = ('O', 'OXT') if all(name in atomPositions for name in ('N', 'CA', 'C')) else ()
                     points1 = []
                     points2 = []
                     for atom in template.topology.atoms():
-                        if atom.name in atomPositions:
+                        if atom.name in atomPositions and atom.name not in skip:
                             points1.append(atomPositions[atom.name].value_in_unit(unit.nanometer))
                             points2.append(template.positions[atom.index].value_in_unit(unit.nanometer))
 

pdbfixer/tests/test_find_positions.py

@@ -1,26 +1,99 @@
 from pathlib import Path
 
-from openmm import app, unit
+import numpy as np
+import pytest
+from openmm import Vec3, app, unit
 
 import pdbfixer
+from pdbfixer.pdbfixer import _overlayPoints
 
 
-def test_findOXTPosition():
-    """Test that OXT is added at the correct position."""
-    fixer = pdbfixer.PDBFixer(filename=(Path(__file__).parent / "data" / "1BHL.pdb").as_posix())
+def _load1BHL():
+    """Load the 1BHL test structure."""
+    return pdbfixer.PDBFixer(filename=(Path(__file__).parent / "data" / "1BHL.pdb").as_posix())
 
-    # Record the original position of OXT, then delete it.
 
-    originalPos = [pos for pos, atom in zip(fixer.positions, fixer.topology.atoms()) if atom.name == "OXT"]
+def _deleteAndRebuild(fixer, shouldDelete):
+    """Delete the atoms selected by shouldDelete, then have PDBFixer add them back without building missing residues."""
     modeller = app.Modeller(fixer.topology, fixer.positions)
-    modeller.delete([atom for atom in modeller.topology.atoms() if atom.name == "OXT"])
+    modeller.delete([atom for atom in modeller.topology.atoms() if shouldDelete(atom)])
     fixer.topology = modeller.topology
     fixer.positions = modeller.positions
-
-    # Have PDBFixer add it back and make sure it is sufficiently close to the original position.
-
     fixer.missingResidues = {}
     fixer.findMissingAtoms()
     fixer.addMissingAtoms()
-    newPos = [pos for pos, atom in zip(fixer.positions, fixer.topology.atoms()) if atom.name == "OXT"]
-    assert unit.norm(newPos[0] - originalPos[0]) < 0.01 * unit.nanometer
+
+
+def _oxtPositions(fixer):
+    """Return the positions of all OXT atoms."""
+    return [pos for pos, atom in zip(fixer.positions, fixer.topology.atoms()) if atom.name == "OXT"]
+
+
+def _residueHandedness(fixer, shouldInclude):
+    """Return the handedness (+1 or -1) of the selected alpha carbons, in residue order."""
+    positions = fixer.positions.value_in_unit(unit.nanometer)
+    handedness = []
+    for residue in fixer.topology.residues():
+        if shouldInclude(residue):
+            pos = {atom.name: positions[atom.index] for atom in residue.atoms()}
+            toN = pos["N"]-pos["CA"]
+            toC = pos["C"]-pos["CA"]
+            toCB = pos["CB"]-pos["CA"]
+            handedness.append(np.sign(np.dot(toN, np.cross(toC, toCB))))
+    return handedness
+
+
+def test_overlayPointsNeverReflects():
+    """Test that _overlayPoints returns a proper rotation even when a reflection would fit better."""
+    points2 = [Vec3(0, 0, 0), Vec3(1, 0, 0), Vec3(0, 1, 0), Vec3(1, 1, 0.1)]
+    points1 = [Vec3(p[0], p[1], -p[2]) for p in points2]
+    (_, rotate, _) = _overlayPoints(points1, points2)
+    assert np.isclose(np.linalg.det(rotate), 1)
+
+
+def test_findOXTPosition():
+    """Test that OXT is added at the correct position."""
+    fixer = _load1BHL()
+
+    # Record the original position of OXT, then delete it and have PDBFixer add it back.
+
+    originalPos = _oxtPositions(fixer)
+    _deleteAndRebuild(fixer, lambda atom: atom.name == "OXT")
+
+    # Make sure it is sufficiently close to the original position.
+
+    assert unit.norm(_oxtPositions(fixer)[0] - originalPos[0]) < 0.01 * unit.nanometer
+
+
+def test_addedCBChirality():
+    """Test that missing CB atoms are added with the correct chirality."""
+    fixer = _load1BHL()
+
+    # Record the handedness of every alanine, then delete its CB and have PDBFixer add it back.
+
+    originalHandedness = _residueHandedness(fixer, lambda residue: residue.name == "ALA")
+    assert len(originalHandedness) > 0
+    _deleteAndRebuild(fixer, lambda atom: atom.residue.name == "ALA" and atom.name == "CB")
+
+    # Make sure none of them was built as its mirror image.  Adding atoms keeps the residues
+    # in the same order, so the lists can be compared directly.
+
+    assert _residueHandedness(fixer, lambda residue: residue.name == "ALA") == originalHandedness
+
+
+@pytest.mark.parametrize("residueId", ["57", "72", "112"])
+def test_addedSideChainChirality(residueId):
+    """Test that the backbone oxygen does not put a rebuilt side chain on the wrong side."""
+    fixer = _load1BHL()
+
+    def isSelected(residue):
+        return residue.chain.id == "A" and residue.id == residueId
+
+    # SER57, VAL72, and THR112 have backbone oxygen positions that tilt the template
+    # fit enough to put CB on the wrong side, even when the fit is a proper rotation.
+
+    originalHandedness = _residueHandedness(fixer, isSelected)
+    assert len(originalHandedness) == 1
+    _deleteAndRebuild(fixer, lambda atom: isSelected(atom.residue) and atom.name not in ("N", "CA", "C", "O", "OXT"))
+
+    assert _residueHandedness(fixer, isSelected) == originalHandedness