PRS directional response

prody/dynamics/perturb.py

@@ -45,10 +45,27 @@ def calcPerturbResponse(model, **kwargs):
         ATPase subdomain IA is a mediator of interdomain allostery in Hsp70
         molecular chaperones. *PLoS Comput. Biol.* **2014** 10:e1003624.
 
-    If *turbo* is **True** (default), then PRS is approximated by the limit of 
-    large numbers of forces and no perturbation forces are explicitly applied. 
-    If set to **False**, then each residue/node is perturbed *repeats* times (default 100) 
+    If *turbo* is **True** (default), then PRS is approximated by the limit of
+    large numbers of forces and no perturbation forces are explicitly applied.
+    If set to **False**, then each residue/node is perturbed *repeats* times (default 100)
     with a random unit force vector as in ProDy v1.8 and earlier.
+
+    If *return_vectors* is **True**, the full directional responses are kept
+    instead of being squared and averaged into magnitudes. This overrides
+    *turbo*: explicit perturbing forces are always applied (as when *turbo* is
+    **False**), so the *turbo* setting is ignored. The response displacement
+    vectors ``cov . f`` are returned as a :class:`.ModeEnsemble` with one modeset
+    per perturbed node, each holding the *repeats* response vectors to random
+    unit forces, with eigenvalues set to the squared response magnitudes. A
+    specific direction (or
+    set of directions) may be supplied via *force* (a 3-vector or an ``(m, 3)``
+    array) to override the random forces; *perturb_node* (an index or list of
+    indices) restricts the scan to selected nodes. The resulting ensemble plugs
+    into the SignDy tools (:func:`.calcOverlap`, :meth:`.ModeEnsemble.match`,
+    :func:`.saveModeEnsemble`), so each response can be overlapped with a target
+    conformational change (e.g. from :func:`.calcDeformVector`) to find which
+    perturbations trigger a given transition and map a functional cycle. Requires
+    a 3d *model* (ANM or PCA); GNM responses have no direction.
     """
 
     if not isinstance(model, (NMA, ModeSet, Mode)):
@@ -76,6 +93,60 @@ def calcPerturbResponse(model, **kwargs):
 
     cov = model.getCovariance()
 
+    return_vectors = kwargs.get('return_vectors', False)
+    if return_vectors:
+        # Directional PRS: keep the response displacement vectors cov.f rather
+        # than squaring/averaging them into magnitudes, so their direction can be
+        # overlapped with a target conformational change. Returns a ModeEnsemble
+        # (one modeset per perturbed node, one "mode" per force).
+        from .signature import ModeEnsemble
+
+        if not model.is3d():
+            raise ValueError('return_vectors=True requires a 3d model (ANM or '
+                             'PCA); GNM responses have no direction')
+
+        repeats = kwargs.pop('repeats', 100)
+        force = kwargs.get('force', None)
+        perturb_node = kwargs.get('perturb_node', None)
+
+        if perturb_node is None:
+            nodes = list(range(n_atoms))
+        elif np.isscalar(perturb_node):
+            nodes = [int(perturb_node)]
+        else:
+            nodes = [int(n) for n in perturb_node]
+
+        if force is not None:
+            forces = np.atleast_2d(np.asarray(force, dtype=float))
+            if forces.shape[1] != 3:
+                raise ValueError('force must be a 3-vector or an (m, 3) array')
+
+        resnums = atoms.getResnums() if atoms is not None else None
+        ens = ModeEnsemble('{0} PRS response'.format(model.getTitle()))
+
+        LOGGER.progress('Calculating directional perturbation response',
+                        len(nodes), '_prody_prs_vec')
+        for k, i in enumerate(nodes):
+            i3 = i * 3
+            if force is None:
+                f = np.random.rand(repeats * 3).reshape((repeats, 3))
+                f /= ((f**2).sum(1)**0.5).reshape((repeats, 1))
+            else:
+                f = forces
+            # response vectors cov.f, shape (3*n_atoms, n_forces)
+            resp = np.dot(cov[:, i3:i3 + 3], f.T)
+            vals = (resp**2).sum(0)  # squared response magnitude per force
+            label = 'perturb_{0}'.format(i if resnums is None else resnums[i])
+            nma_i = NMA(label)
+            nma_i.setEigens(resp, vals)
+            ens.addModeSet(nma_i, label=label)
+            LOGGER.update(k, label='_prody_prs_vec')
+        LOGGER.clear()
+
+        if atoms is not None:
+            ens.setAtoms(atoms)
+        return ens
+
     turbo = kwargs.get('turbo', True)
     if turbo:
         if not model.is3d():

prody/tests/dynamics/test_perturb.py

@@ -0,0 +1,117 @@
+"""This module contains unit tests for :func:`.calcPerturbResponse`, in
+particular the directional ``return_vectors`` response ensemble, exercised on
+the adenylate kinase open->closed transition (4ake chain A -> 1ake chain A)."""
+
+import numpy as np
+from numpy.testing import assert_allclose
+
+from prody import (ANM, GNM, matchChains, superpose,
+                   calcDeformVector, calcOverlap, calcPerturbResponse, LOGGER)
+from prody.dynamics.signature import ModeEnsemble
+from prody.tests import unittest
+from prody.tests.datafiles import parseDatafile
+
+LOGGER.verbosity = 'none'
+
+
+class TestCalcPerturbResponse(unittest.TestCase):
+
+    @classmethod
+    def setUpClass(cls):
+        # open (4ake, chain A) and closed (1ake) adenylate kinase; bundled data
+        op = parseDatafile('pdb4ake_fixed').select('calpha')
+        cl = parseDatafile('pdb1ake').select('calpha')
+        matches = matchChains(op, cl, pwalign=True, seqid=90, overlap=90)
+        op_ca, cl_ca = matches[0][0], matches[0][1]
+
+        # closure target: open->closed internal deformation (after superposition)
+        cl_fit = cl_ca.copy()
+        superpose(cl_fit, op_ca)
+        cls.target = calcDeformVector(op_ca, cl_fit)
+
+        cls.n_atoms = op_ca.numAtoms()
+        cls.anm = ANM('4ake_A')
+        cls.anm.buildHessian(op_ca)
+        cls.anm.calcModes(n_modes=None, zeros=False)
+        cls.gnm = GNM('4ake_A')
+        cls.gnm.buildKirchhoff(op_ca)
+        cls.gnm.calcModes(n_modes=None, zeros=False)
+
+    def testLegacyReturnsTriple(self):
+        """Default call is unchanged: (matrix, effectiveness, sensitivity)."""
+        prs, eff, sens = calcPerturbResponse(self.anm)
+        self.assertEqual(prs.shape, (self.n_atoms, self.n_atoms))
+        self.assertEqual(eff.shape, (self.n_atoms,))
+        self.assertEqual(sens.shape, (self.n_atoms,))
+
+    def testReturnVectorsEnsembleShape(self):
+        """return_vectors gives a ModeEnsemble: one modeset/node, repeats modes."""
+        repeats = 7
+        ens = calcPerturbResponse(self.anm, return_vectors=True, repeats=repeats)
+        self.assertIsInstance(ens, ModeEnsemble)
+        self.assertEqual(ens.numModeSets(), self.n_atoms)
+        self.assertEqual(ens.numModes(), repeats)
+        self.assertEqual(ens.numAtoms(), self.n_atoms)
+
+    def testResponseVectorEqualsCovDotForce(self):
+        """The kept response vector is exactly cov . f (no square/average)."""
+        cov = self.anm.getCovariance()
+        node = 5
+        force = np.array([1.0, -2.0, 0.5])
+        ens = calcPerturbResponse(self.anm, return_vectors=True,
+                                  perturb_node=node, force=force)
+        self.assertEqual(ens.numModeSets(), 1)
+        self.assertEqual(ens.numModes(), 1)
+        resp = np.asarray(ens[0].getEigvecs()).flatten()
+        expected = np.dot(cov[:, node * 3:node * 3 + 3], force)
+        assert_allclose(resp, expected, rtol=1e-6, atol=1e-8)
+
+    def testEigvalsAreSquaredResponseMagnitude(self):
+        """Eigenvalues store the squared response magnitude of each force."""
+        cov = self.anm.getCovariance()
+        node = 3
+        force = np.array([0.3, 0.4, 0.0])
+        ens = calcPerturbResponse(self.anm, return_vectors=True,
+                                  perturb_node=node, force=force)
+        expected = np.dot(cov[:, node * 3:node * 3 + 3], force)
+        eigval = np.atleast_1d(ens[0].getEigvals())[0]
+        assert_allclose(eigval, (expected ** 2).sum(), rtol=1e-6, atol=1e-8)
+
+    def testPerturbNodeList(self):
+        """perturb_node restricts the scan to the requested nodes."""
+        nodes = [0, 4, 9]
+        ens = calcPerturbResponse(self.anm, return_vectors=True,
+                                  repeats=3, perturb_node=nodes)
+        self.assertEqual(ens.numModeSets(), len(nodes))
+
+    def testGNMRaisesForVectors(self):
+        """Directional responses require a 3d model; GNM has no direction."""
+        self.assertRaises(ValueError, calcPerturbResponse, self.gnm,
+                          return_vectors=True)
+
+    def testDirectionalPRSRecoversClosure(self):
+        """A directed force at some node drives the open->closed closure: the
+        best response overlaps the transition better than any single ANM mode."""
+        np.random.seed(42)
+        ens = calcPerturbResponse(self.anm, return_vectors=True, repeats=100)
+        tv = self.target.getArray()
+        tv = tv / np.linalg.norm(tv)
+
+        best_per_node = []
+        for i in range(ens.numModeSets()):
+            R = np.asarray(ens[i].getEigvecs())
+            if R.shape[0] != tv.shape[0]:
+                R = R.T
+            R = R / np.linalg.norm(R, axis=0, keepdims=True)
+            best_per_node.append(np.abs(tv @ R).max())
+        best_directional = float(np.max(best_per_node))
+
+        best_single_mode = float(np.abs(calcOverlap(self.anm[:20], self.target)).max())
+
+        # closure is genuinely recovered, and beats the best single normal mode
+        self.assertGreater(best_directional, 0.6)
+        self.assertGreater(best_directional, best_single_mode)
+
+
+if __name__ == '__main__':
+    unittest.main()