ANMD and ClustENM fixes and test fixes

environment.yml

@@ -20,7 +20,8 @@ dependencies:
   - mdtraj
   - openmm
   - clustalw
-
+  - clustalo
+
   # Pip packages
   - pip
   - pip:

prody/dynamics/anmd.py

@@ -25,22 +25,20 @@
 __credits__ = ['James Krieger']
 __email__ = ['anupam.banerjee@stonybrook.edu', 'jamesmkrieger@gmail.com']
 
-from numbers import Number
 import os
+from numbers import Number
 
 from prody import LOGGER
 from prody.atomic.atomic import Atomic
-from prody.ensemble.ensemble import Ensemble
-from prody.proteins.pdbfile import parsePDB, writePDB
-
 from prody.dynamics.anm import ANM
 from prody.dynamics.clustenm import ClustENM
 from prody.dynamics.editing import extendModel
 from prody.dynamics.modeset import ModeSet
 from prody.dynamics.nma import NMA
 from prody.dynamics.pca import PCA
 from prody.dynamics.sampling import traverseMode
-
+from prody.ensemble.ensemble import Ensemble
+from prody.proteins.pdbfile import parsePDB, writePDB
 
 __all__ = ['runANMD']
 
@@ -89,15 +87,6 @@ def runANMD(atoms, num_modes=2, max_rmsd=2., num_steps=5, tolerance=10.0,
             A molecular assessment of the alterations in the spike-host protein 
             interactions. *iScience* **2022** 25(3):103939.
     """
-    try:
-        from simtk.openmm.app import PDBFile, ForceField, \
-            Simulation, HBonds, NoCutoff
-        from simtk.openmm import LangevinIntegrator
-        from simtk.unit import nanometer, kelvin, picosecond, picoseconds, \
-            angstrom, kilojoule, mole
-    except ImportError:
-        raise ImportError('Please install PDBFixer and OpenMM to use ANMD')
-
     if not isinstance(atoms, Atomic):
         raise TypeError('atoms should be an Atomic object')
 
@@ -112,7 +101,6 @@ def runANMD(atoms, num_modes=2, max_rmsd=2., num_steps=5, tolerance=10.0,
 
     if not isinstance(tolerance, Number):
         raise TypeError('tolerance should be a float')
-    tolerance = tolerance * kilojoule/mole/nanometer
 
     pos = kwargs.get('pos', True)
     if not isinstance(pos, bool):
@@ -134,6 +122,44 @@ def runANMD(atoms, num_modes=2, max_rmsd=2., num_steps=5, tolerance=10.0,
     if not isinstance(anm, (type(None), NMA, ModeSet)):
         raise TypeError('anm should be an NMA or ModeSet object')
 
+    # OpenMM 7.6 moved these out of the legacy simtk namespace, which is
+    # still shipped as a deprecated shim, so try the modern names first.
+    try:
+        try:
+            from openmm import LangevinIntegrator
+            from openmm.app import ForceField, HBonds, NoCutoff, PDBFile, Simulation
+            from openmm.unit import (
+                angstrom,
+                kelvin,
+                kilojoule,
+                mole,
+                nanometer,
+                picosecond,
+                picoseconds,
+            )
+        except ImportError:
+            from simtk.openmm import LangevinIntegrator
+            from simtk.openmm.app import (
+                ForceField,
+                HBonds,
+                NoCutoff,
+                PDBFile,
+                Simulation,
+            )
+            from simtk.unit import (
+                angstrom,
+                kelvin,
+                kilojoule,
+                mole,
+                nanometer,
+                picosecond,
+                picoseconds,
+            )
+    except ImportError:
+        raise ImportError('Please install PDBFixer and OpenMM to use ANMD')
+
+    tolerance = tolerance * kilojoule/mole/nanometer
+
     pdb_name=atoms.getTitle().replace(' ', '_')
 
     fix_name = pdb_name + '_fixed.pdb'
@@ -142,7 +168,7 @@ def runANMD(atoms, num_modes=2, max_rmsd=2., num_steps=5, tolerance=10.0,
     else:
         clustenm=ClustENM()
         clustenm.setAtoms(atoms)
-        clustenm.writePDBFixed()
+        fix_name = clustenm.writePDBFixed(fix_name)
     pdb_fix = PDBFile(fix_name)
 
     fixmin_name=pdb_name + '_fixedmin.pdb'
@@ -163,7 +189,8 @@ def runANMD(atoms, num_modes=2, max_rmsd=2., num_steps=5, tolerance=10.0,
         simulation.context.setPositions(pdb_fix.positions)
         simulation.minimizeEnergy(tolerance=tolerance)
         positions = simulation.context.getState(getPositions=True).getPositions()
-        PDBFile.writeFile(simulation.topology, positions, open(fixmin_name, 'w'))
+        with open(fixmin_name, 'w') as stream:
+            PDBFile.writeFile(simulation.topology, positions, stream)
         LOGGER.report('The fixed structure was minimised in %.2fs.\n',
                         label='_anmd_min')
 

prody/dynamics/clustenm.py

@@ -47,7 +47,7 @@
 from prody.measure import calcTransformation, applyTransformation, calcRMSD
 from prody.ensemble import Ensemble
 from prody.proteins import writePDB, parsePDB, writePDBStream, parsePDBStream
-from prody.utilities import createStringIO, importLA, mad
+from prody.utilities import checkBlocks, createStringIO, importLA, mad
 
 la = importLA()
 norm = la.norm
@@ -408,12 +408,14 @@ def _prep_sim(self, coords, external_forces=[]):
                                         0.002*picosecond)
 
         # precision could be mixed, but single is okay.
-        platform = self._platform if self._platform is None else Platform.getPlatformByName(self._platform)
+        platform = None if self._platform is None else Platform.getPlatformByName(self._platform)
+
+        # properties are platform specific, so they can only be given for a
+        # platform we chose. Without one OpenMM picks the fastest it has and
+        # uses its own defaults, and asking for properties as well raises.
         properties = None
 
-        if self._platform is None:
-            properties = {'Precision': 'single'}
-        elif self._platform in ['CUDA', 'OpenCL']:
+        if self._platform in ['CUDA', 'OpenCL']:
             properties = {'Precision': 'single'}
             # pin this Context to the device assigned to this parallel-sim worker (read at build time,
             # so it is immune to the CUDA_VISIBLE_DEVICES-vs-import timing race)
@@ -1106,19 +1108,27 @@ def _getCoordsets(self, indices=None, selected=True):
 
         return super(ClustENM, self)._getCoordsets(I, selected)
 
-    def writePDBFixed(self):
+    def writePDBFixed(self, filename=None, replace_spaces=False):
 
         'Write the fixed (initial) structure to a pdb file.'
 
+        if filename is None:
+            filename = self.getTitle()[:-8] + 'fixed.pdb'
+        if replace_spaces:
+            filename = filename.replace(" ", "_")
+
         try:
             from openmm.app import PDBFile
         except ImportError:
             raise ImportError('Please install PDBFixer and OpenMM 7.6 in order to use ClustENM.')
 
-        PDBFile.writeFile(self._topology,
-                          self._positions,
-                          open(self.getTitle()[:-8] + 'fixed.pdb', 'w'),
-                          keepIds=True)
+        with open(filename, 'w') as stream:
+            PDBFile.writeFile(self._topology,
+                            self._positions,
+                            stream,
+                            keepIds=True)
+
+        return filename
 
     def writePDB(self, filename=None, single=True, **kwargs):
 
@@ -1379,7 +1389,7 @@ def run(self, cutoff=15., n_modes=3, gamma=1., n_confs=50, rmsd=1.0,
 
             if len(self._maxclust) != self._n_gens + 1:
                 raise ValueError(
-                    'size mismatch: %d generations were set; %d maxclusts were given' % (
+                    'size mismatch: %d generations were set; %d maxclusts were given' % (  # noqa: UP031
                         self._n_gens + 1, len(self._maxclust)))
 
         if threshold is None:
@@ -1391,7 +1401,9 @@ def run(self, cutoff=15., n_modes=3, gamma=1., n_confs=50, rmsd=1.0,
                 self._threshold = (0,) + (threshold,) * n_gens
 
             if len(self._threshold) != self._n_gens + 1:
-                raise ValueError('size mismatch: %d generations were set; %d thresholds were given' % (self._n_gens + 1, self._threshold))
+                raise ValueError(
+                    'size mismatch: %d generations were set; %d thresholds were given' % (  # noqa: UP031
+                    self._n_gens + 1, len(self._threshold)))
 
         self._sol = solvent if self._nuc is None else 'exp'
         self._padding = kwargs.pop('padding', 1.0)
@@ -1418,6 +1430,17 @@ def run(self, cutoff=15., n_modes=3, gamma=1., n_confs=50, rmsd=1.0,
 
         self._cycle = 0
 
+        if self._atoms is None:
+            raise ValueError('atoms are not set, use `setAtoms`')
+
+        # the block based subclasses check that blocks are set before getting
+        # here, so this checks them against the nodes they have to describe,
+        # which fails the run before it minimises anything rather than later
+        # when its ANM is built
+        blocks = getattr(self, '_blocks', None)
+        if blocks is not None:
+            checkBlocks(blocks, self._n_cg, 'coarse grained nodes')
+
         # check for discontinuity in the structure
         gnm = GNM()
         gnm.buildKirchhoff(self._atoms[self._idx_cg], cutoff=self._cutoff)

prody/dynamics/rtb.py

@@ -5,7 +5,7 @@
 import numpy as np
 
 from prody import LOGGER
-from prody.utilities import checkCoords
+from prody.utilities import checkBlocks, checkCoords
 
 from .anm import ANMBase
 from .gnm import GNMBase
@@ -76,8 +76,7 @@ def buildHessian(self, coords, blocks, cutoff=15., gamma=1., **kwargs):
     def calcProjection(self, coords, blocks, **kwargs):
         natoms = self._n_atoms
 
-        if natoms != len(blocks):
-            raise ValueError('len(blocks) must match number of atoms')
+        checkBlocks(blocks, natoms)
 
         LOGGER.timeit('_rtb')
         from collections import defaultdict

prody/sequence/analysis.py

@@ -742,10 +742,15 @@ def buildMSA(sequences, title='Unknown', labels=None, **kwargs):
     :type align: bool
 
     :arg method: alignment method, one of either Biopython 'global',
-        Biopython 'local', 'clustalw', 'clustalw2', 'clustal' 
+        Biopython 'local', 'clustalw', 'clustalw2', 'clustal', 'clustalo'
         or another software in your path.
         Default is 'local'
     :type align: str
+
+    :arg outfilename: name of the file to write the alignment to, which
+        clustalw also writes beside as an ``.aln`` file.
+        Default is *title* with a ``.fasta`` extension
+    :type outfilename: str
     """
 
     align = kwargs.get('align', True)

prody/tests/dynamics/test_anmd.py

@@ -1,10 +1,16 @@
 """This module contains unit tests for :mod:`~prody.dynamics`."""
 
 import os
+import shutil
 import sys
+import tempfile
+import warnings
+
 import numpy as np
 from numpy.testing import *
+
 from prody.utilities import importDec
+
 dec = importDec()
 
 from prody import *
@@ -16,7 +22,7 @@
 try:
     from unittest.mock import MagicMock, patch
 except ImportError:
-    from mock import MagicMock, patch
+    from mock import MagicMock, patch  # noqa: UP026
 
 # Prevent threading hangs on remote servers
 os.environ['OMP_NUM_THREADS'] = '1'
@@ -29,6 +35,19 @@
 ENSEMBLE = PDBEnsemble(parseDatafile('anmd'))
 ENSEMBLE.setCoords(ENSEMBLE.getCoordsets()[2])
 
+# runANMD needs both of these to get as far as minimising anything. They are
+# optional dependencies, so the tests that need them are skipped rather than
+# failed when they are missing, with a warning so that the gap is not silent.
+OPENMM_SKIP_MSG = 'PDBFixer and OpenMM are needed to run ANMD'
+try:
+    import openmm
+    import pdbfixer
+    HAVE_OPENMM = True
+except ImportError:
+    HAVE_OPENMM = False
+    warnings.warn(OPENMM_SKIP_MSG + ', so the ANMD tests that run it are '
+                  'being skipped')
+
 class TestANMD(unittest.TestCase):
 
     def setUp(self):
@@ -141,5 +160,59 @@ def testResultsNumModes1(self):
                             rtol=1e-10, atol=0.25, # may not be so close
                             err_msg='runANMD with num_modes=1 failed to give expected RMSDs')
 
+class TestAnmdSpacedTitle(unittest.TestCase):
+    """runANMD hands the fixed structure to OpenMM through a PDB file whose
+    name it derives from the title of *atoms*, replacing spaces with
+    underscores. It used to look for that name while ClustENM wrote the
+    unreplaced one, so a title containing a space raised FileNotFoundError."""
+
+    def setUp(self):
+        if not prody.PY3K:
+            return
+
+        # runANMD writes its intermediate PDB files into the current directory
+        self.cwd = os.getcwd()
+        self.tmpdir = tempfile.mkdtemp()
+        os.chdir(self.tmpdir)
+
+        self.ATOMS = parseDatafile('1ubi')
+        self.ATOMS.setTitle('my prot')
+
+    def tearDown(self):
+        if not prody.PY3K:
+            return
+
+        os.chdir(self.cwd)
+        shutil.rmtree(self.tmpdir, ignore_errors=True)
+
+    @dec.slow
+    @unittest.skipUnless(HAVE_OPENMM, OPENMM_SKIP_MSG)
+    def testAnmdSpacedTitle(self):
+        """Test that a title containing a space is handled"""
+        if prody.PY3K:
+            # a loose tolerance keeps the minimisation short
+            ensembles = runANMD(self.ATOMS, num_modes=1, num_steps=1,
+                                tolerance=100.)
+
+            assert_equal(len(ensembles), 1,
+                         'runANMD with a spaced title failed to give 1 ensemble')
+            assert_equal(len(ensembles[0]), 3,
+                         'runANMD with a spaced title failed to give an '
+                         'ensemble with 3 conformers')
+
+    @dec.slow
+    @unittest.skipUnless(HAVE_OPENMM, OPENMM_SKIP_MSG)
+    def testAnmdSpacedTitleCleansUp(self):
+        """Test that the intermediate files of a spaced title are removed"""
+        if prody.PY3K:
+            runANMD(self.ATOMS, num_modes=1, num_steps=1, tolerance=100.)
+
+            leftover = [name for name in os.listdir(self.tmpdir)
+                        if name.endswith('.pdb')]
+            assert_equal(leftover, [],
+                         'runANMD with a spaced title left intermediate files '  # noqa: UP031
+                         'behind: %s' % leftover)
+
+
 if __name__ == '__main__':
     unittest.main()