diff --git a/prody/apps/prody_apps/prody_anm.py b/prody/apps/prody_apps/prody_anm.py index 17d1ccd9d..63072ecc3 100644 --- a/prody/apps/prody_apps/prody_anm.py +++ b/prody/apps/prody_apps/prody_anm.py @@ -13,7 +13,7 @@ for key, txt, val in [ ('model', 'index of model that will be used in the calculations', 1), ('altloc', 'alternative location identifiers for residues used in the calculations', "A"), - ('cutoff', 'cutoff distance (A)', 15.), + ('cutoff', 'cutoff distance (A)', '15.'), ('gamma', 'spring constant', '1.'), ('sparse', 'use sparse matrices', False), ('kdtree', 'use kdtree for Hessian', False), @@ -57,7 +57,6 @@ def prody_anm(pdb, **kwargs): selstr = kwargs.get('select') prefix = kwargs.get('prefix') - cutoff = kwargs.get('cutoff') sparse = kwargs.get('sparse') kdtree = kwargs.get('kdtree') nmodes = kwargs.get('nmodes') @@ -90,6 +89,19 @@ def prody_anm(pdb, **kwargs): raise NameError("Please provide gamma as a float or ProDy Gamma class") except TypeError: raise TypeError("Please provide gamma as a float or ProDy Gamma class") + + try: + cutoff = float(kwargs.get('cutoff')) + LOGGER.info("Using cutoff {0}".format(cutoff)) + except ValueError: + try: + import math + cutoff = eval(kwargs.get('cutoff')) + LOGGER.info("Using cutoff {0}".format(cutoff)) + except NameError: + raise NameError("Please provide cutoff as a float or equation using math") + except TypeError: + raise TypeError("Please provide cutoff as a float or equation using math") anm = prody.ANM(pdb.getTitle()) @@ -275,7 +287,7 @@ def addCommand(commands): group = addNMAParameters(subparser) - group.add_argument('-c', '--cutoff', dest='cutoff', type=float, + group.add_argument('-c', '--cutoff', dest='cutoff', type=str, default=DEFAULTS['cutoff'], metavar='FLOAT', help=HELPTEXT['cutoff'] + ' (default: %(default)s)') diff --git a/prody/dynamics/editing.py b/prody/dynamics/editing.py index 9ce1d47fe..bccd5db58 100644 --- a/prody/dynamics/editing.py +++ b/prody/dynamics/editing.py @@ -5,8 +5,9 @@ from prody.atomic import Atomic, AtomGroup, AtomMap, AtomSubset from prody.atomic import Selection, SELECT, sliceAtoms, extendAtoms -from prody.utilities import importLA, isListLike -from prody import _PY3K +from prody.measure import calcDistance +from prody.utilities import importLA, isListLike, getCoords +from prody import _PY3K, LOGGER from .nma import NMA from .mode import VectorBase, Mode, Vector @@ -19,7 +20,8 @@ __all__ = ['extendModel', 'extendMode', 'extendVector', 'sliceMode', 'sliceModel', 'sliceModelByMask', 'sliceVector', - 'reduceModel', 'reduceModelByMask', 'trimModel', 'trimModelByMask'] + 'reduceModel', 'reduceModelByMask', 'trimModel', 'trimModelByMask', + 'interpolateModel'] def extendModel(model, nodes, atoms, norm=False): @@ -122,10 +124,14 @@ def sliceVector(vector, atoms, select): which, sel = sliceAtoms(atoms, select) - vec = Vector(vector.getArrayNx3()[ - which, :].flatten(), - '{0} slice {1}'.format(str(vector), select), - vector.is3d()) + if vector.is3d(): + vec = Vector(vector.getArrayNx3()[which, :].flatten(), + '{0} slice {1}'.format(str(vector), select), + vector.is3d()) + else: + vec = Vector(vector.getArray()[which].flatten(), + '{0} slice {1}'.format(str(vector), select), + vector.is3d()) return (vec, sel) def trimModel(model, atoms, select): @@ -260,9 +266,14 @@ def sliceMode(mode, atoms, select): which, sel = sliceAtoms(atoms, select) - vec = Vector(mode.getArrayNx3()[which, :].flatten() * - mode.getVariance()**0.5, - '{0} slice {1}'.format(str(mode), select), mode.is3d()) + if mode.is3d(): + vec = Vector(mode.getArrayNx3()[which, :].flatten() * + mode.getVariance()**0.5, + '{0} slice {1}'.format(str(mode), select), mode.is3d()) + else: + vec = Vector(mode.getArray()[which].flatten() * + mode.getVariance()**0.5, + '{0} slice {1}'.format(str(mode), select), mode.is3d()) return (vec, sel) @@ -485,3 +496,151 @@ def _reduceModel(matrix, system): matrix = ss return matrix + + +def interpolateModel(model, nodes, coords, norm=False, **kwargs): + """Interpolate a coarse grained *model* built for *nodes* to *coords*. + + *model* may be :class:`.ANM`, :class:`.PCA`, or :class:`.NMA` + instance + + :arg nodes: the coordinate set or object with :meth:`getCoords` method + that corresponds to the model + :type nodes: :class:`.Atomic`, :class:`~numpy.ndarray` + + :arg coords: a coordinate set or an object with :meth:`getCoords` method + onto which the model should be interpolated + :type coords: :class:`.Atomic`, :class:`~numpy.ndarray` + + This function will take the part of the normal modes for each node + (i.e. Cα atoms) and extend it to nearby atoms. + + If *norm* is **True**, extended modes are normalized. + + Adapted from ModeHunter as described in [JS09]_. + + .. [JS09] Stember JN, Wriggers W. Bend-twist-stretch model for coarse + elastic network simulation of biomolecular motion. *J Chem Phys* **2009** 131:074112. + + # Legal notice: + # + # This software is copyrighted, (c) 2009-10, by Joseph N. Stember and Willy Wriggers + # under the following terms: + # + # The authors hereby grant permission to use, copy, modify, and re-distribute this + # software and its documentation for any purpose, provided that existing copyright + # notices are retained in all copies and that this notice is included verbatim in + # any distributions. No written agreement, license, or royalty fee is required for + # any of the authorized uses. + # + # IN NO EVENT SHALL THE AUTHORS OR DISTRIBUTORS BE LIABLE TO ANY PARTY FOR DIRECT, + # INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE + # OF THIS SOFTWARE, ITS DOCUMENTATION, OR ANY DERIVATIVES THEREOF, EVEN IF THE + # AUTHORS HAVE BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + # + # THE AUTHORS AND DISTRIBUTORS SPECIFICALLY DISCLAIM ANY WARRANTIES, INCLUDING, + # BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A + # PARTICULAR PURPOSE, AND NON-INFRINGEMENT. THIS SOFTWARE IS PROVIDED ON AN "AS + # IS" BASIS, AND THE AUTHORS AND DISTRIBUTORS HAVE NO OBLIGATION TO PROVIDE + # MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. + # + # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + """ + + quiet = kwargs.pop('quiet', False) + + try: + eigvecs = model._getArray() + eigvals = model.getEigvals() + except AttributeError: + raise TypeError('model must be an NMA instance') + + if model.numAtoms() != nodes.numAtoms(): + raise ValueError('atom numbers must be the same') + + if not model.is3d(): + raise ValueError('model must be 3D') + + cg_coords = nodes.getCoords() + len_cg = nodes.numAtoms() + + fg_coords = getCoords(coords) + len_fg = fg_coords.shape[0] + + n_modes = model.numModes() + + eigvecs_fg = np.zeros((3*len_fg, n_modes)) + + if not quiet: + LOGGER.progress('Interpolating {0} coarse modes to fine modes...'.format(n_modes), + n_modes, '_prody_iterp') + + for k in np.arange(n_modes): # loop through cg eigenvecs + coarse_grained_nms = eigvecs[:, k] + coarse_grained_nms.shape = (-1,3) + + lmat = np.zeros((len_cg+4,len_cg+4)) + for i in np.arange(len_cg): + for j in np.arange(len_cg): + lmat[i,j] = np.linalg.norm(cg_coords[i]-cg_coords[j]) + + for i in np.arange(len_cg): + lmat[i,len_cg] = 1.0 + lmat[i,len_cg+1] = cg_coords[i,0] + lmat[i,len_cg+2] = cg_coords[i,1] + lmat[i,len_cg+3] = cg_coords[i,2] + lmat[len_cg,i] = 1.0 + lmat[len_cg+1,i] = cg_coords[i,0] + lmat[len_cg+2,i] = cg_coords[i,1] + lmat[len_cg+3,i] = cg_coords[i,2] + + lmat_inv = np.linalg.inv(lmat) + + vxprime = coarse_grained_nms[:,0] + for i in np.arange(4): + vxprime = np.append(vxprime,0.0) + + vyprime = coarse_grained_nms[:,1] + for i in np.arange(4): + vyprime = np.append(vyprime,0.0) + + vzprime = coarse_grained_nms[:,2] + for i in np.arange(4): + vzprime = np.append(vzprime,0.0) + + vx = np.dot(lmat_inv,vxprime) + vy = np.dot(lmat_inv,vyprime) + vz = np.dot(lmat_inv,vzprime) + + nmx=np.zeros(len_fg) + nmy=np.zeros(len_fg) + nmz=np.zeros(len_fg) + + for j in np.arange(len_fg): + nmx[j] += vx[len_cg] + fg_coords[j,0]*vx[len_cg+1] + fg_coords[j,1]*vx[len_cg+2] + fg_coords[j,2]*vx[len_cg+3] + nmy[j] += vy[len_cg] + fg_coords[j,0]*vy[len_cg+1] + fg_coords[j,1]*vy[len_cg+2] + fg_coords[j,2]*vy[len_cg+3] + nmz[j] += vz[len_cg] + fg_coords[j,0]*vz[len_cg+1] + fg_coords[j,1]*vz[len_cg+2] + fg_coords[j,2]*vz[len_cg+3] + + dist_j = calcDistance(cg_coords, fg_coords[j]) + for i in np.arange(len_cg): + nmx[j] += vx[i]*dist_j[i] + nmy[j] += vy[i]*dist_j[i] + nmz[j] += vz[i]*dist_j[i] + + eigvecs_fg[::3, k] = nmx + eigvecs_fg[1::3, k] = nmy + eigvecs_fg[2::3, k] = nmz + + if not quiet: + LOGGER.update(k+1, label='_prody_iterp') + + if not quiet: + LOGGER.finish() + + if norm: + for k in np.arange(np.size(cg_coords)): + eigvecs_fg[k] = eigvecs_fg[k]/np.linalg.norm(eigvecs_fg[k]) + + interpolated = NMA('Interpolated ' + str(model)) + interpolated.setEigens(eigvecs_fg, eigvals) + return interpolated, coords diff --git a/prody/dynamics/gamma.py b/prody/dynamics/gamma.py index 7566aa869..19e64f673 100644 --- a/prody/dynamics/gamma.py +++ b/prody/dynamics/gamma.py @@ -6,7 +6,7 @@ from prody.atomic import Atomic -__all__ = ['Gamma', 'GammaStructureBased', 'GammaVariableCutoff', 'GammaGOdMD'] +__all__ = ['Gamma', 'GammaStructureBased', 'GammaVariableCutoff', 'GammaED', 'GammaGOdMD'] class Gamma(object): @@ -37,7 +37,7 @@ class GammaStructureBased(Gamma): """Facilitate setting the spring constant based on the secondary structure and connectivity of the residues. - A recent systematic study [LT10]_ of a large set of NMR-structures analyzed + A systematic study [LT10]_ of a large set of NMR-structures analyzed using a method based on entropy maximization showed that taking into consideration properties such as sequential separation between contacting residues and the secondary structure types of the interacting @@ -301,15 +301,37 @@ def gamma(self, dist2, i, j): return gamma -class GammaGOdMD(Gamma): +class GammaED(Gamma): """Facilitate setting the spring constant based on - sequence distance and spatial distance as in GOdMD. + sequence distance and spatial distance as in ed-ENM [OL10]_. + This ENM is refined based on comparison to essential dynamics + from MD simulations and can reproduce flexibility in NMR ensembles. + + It has also been implemented in FlexServ [CJ09]_ and + used in MDdMD [SP12]_ and GOdMD [SP13]_. The sequence distance-dependent term is Cseq/(S**2) for S=abs(i,j) <= Slim The structure distance-dependent term is (Ccart/dist)**Ex + .. [OL10] Orellana L, Rueda M, Ferrer-Costa C, Lopez-Blanco JR, Chacón P, Orozco M. + Approaching Elastic Network Models to Molecular Dynamics Flexibility. + *J Chem Theory Comput* **2010** 6(9):2910-23. + + .. [CJ09] Camps J, Carrillo O, Emperador A, Orellana L, Hospital A, Rueda M, + Cicin-Sain D, D'Abramo M, Gelpí JL, Orozco M. + FlexServ: an integrated tool for the analysis of protein flexibility. + *Bioinformatics* **2009** 25(13):1709-10. + + .. [SP12] Sfriso P, Emperador A, Orellana L, Hospital A, Gelpí JL, Orozco M. + Finding Conformational Transition Pathways from Discrete Molecular Dynamics Simulations. + *J Chem Theory Comput* **2012** 8(11):4707-18. + + .. [SP13] Sfriso P, Hospital A, Emperador A, Orozco M. + Exploration of conformational transition pathways from coarse-grained simulations. + *Bioinformatics* **2013** 29(16):1980-6. + **Example**: Let's parse coordinates from a PDB file. @@ -396,3 +418,5 @@ def gamma(self, dist2, i, j): else: dist = dist2**0.5 return (Ccart/dist)**Ex + +GammaGOdMD = GammaED diff --git a/prody/proteins/compare.py b/prody/proteins/compare.py index 6c4171e71..9c624560f 100644 --- a/prody/proteins/compare.py +++ b/prody/proteins/compare.py @@ -923,12 +923,12 @@ def mapChainOntoChain(mobile, target, **kwargs): :type overlap: float :keyword mapping: what method will be used if the trivial mapping based on residue numbers - fails. If ``"ce"`` or ``"cealign"``, then the CE algorithm [IS98]_ will be + fails. If ``"ce"`` or ``"cealign"``, then the CE structural alignment [IS98]_ will be performed. It can also be a list of prealigned sequences, a :class:`.MSA` instance, or a dict of indices such as that derived from a :class:`.DaliRecord`. If set to **True** then the sequence alignment from Biopython will be used. If set to **False**, only the trivial mapping will be performed. - Default is **"auto"** + Default is **"auto"**, which means try sequence alignment then CE. :type mapping: list, str, bool :keyword pwalign: if **True**, then pairwise sequence alignment will diff --git a/prody/tests/dynamics/test_editing.py b/prody/tests/dynamics/test_editing.py index f86185363..c4f7f75f5 100644 --- a/prody/tests/dynamics/test_editing.py +++ b/prody/tests/dynamics/test_editing.py @@ -1,11 +1,13 @@ """This module contains unit tests for :mod:`~prody.KDTree` module.""" -from numpy.testing import assert_array_equal +from numpy.testing import assert_array_equal, assert_equal -from numpy import concatenate +from numpy import concatenate, array from prody.dynamics import calcGNM, calcANM -from prody.dynamics import extendModel, extendMode, extendVector +from prody.dynamics import (extendModel, extendMode, extendVector, + sliceModel, sliceMode, sliceVector, + interpolateModel) from prody.tests import unittest from prody.tests.datafiles import parseDatafile @@ -28,6 +30,14 @@ EXT3D = concatenate([list(range(i*3, (i+1)*3)) * len(res) for i, res in enumerate(ATOMS.iterResidues())]) +SLC1D = concatenate([[i] for i in ATOMS.ca.getIndices()]) + +SLC3D = concatenate([list(range(i*3, (i+1)*3)) + for i in ATOMS.ca.getIndices()]) + +ANM_AA, NODES_AA = calcANM(ATOMS, selstr='all', cutoff=8) +GNM_AA = calcGNM(ATOMS, selstr='all', cutoff=4)[0] + class TestExtending(unittest.TestCase): @@ -67,3 +77,50 @@ def testVector1d(self): ext = extendVector(vector, NODES, ATOMS)[0] assert_array_equal(ext._getArray(), vector._getArray()[EXT1D]) + +class TestSlicing(unittest.TestCase): + + def testModel3d(self): + + slc = sliceModel(ANM_AA, NODES_AA, NODES)[0] + assert_array_equal(slc._getArray(), ANM_AA._getArray()[SLC3D, :]) + + def testModel1d(self): + + slc = sliceModel(GNM_AA, NODES_AA, NODES)[0] + assert_array_equal(slc._getArray(), GNM_AA._getArray()[SLC1D, :]) + + def testMode3d(self): + + mode = ANM_AA[0] + slc = sliceMode(mode, NODES_AA, NODES)[0] + assert_array_equal(slc._getArray(), mode._getArray()[SLC3D] * + mode.getVariance()**0.5) + + def testMode1d(self): + + mode = GNM_AA[0] + slc = sliceMode(mode, NODES_AA, NODES)[0] + assert_array_equal(slc._getArray(), mode._getArray()[SLC1D] * + mode.getVariance()**0.5) + + def testVector3d(self): + + vector = ANM_AA[0] * 1 + slc = sliceVector(vector, NODES_AA, NODES)[0] + assert_array_equal(slc._getArray(), vector._getArray()[SLC3D]) + + def testVector1d(self): + + vector = GNM_AA[0] * 1 + slc = sliceVector(vector, NODES_AA, NODES)[0] + assert_array_equal(slc._getArray(), vector._getArray()[SLC1D]) + + +class TestInterpolating(unittest.TestCase): + + def testModel3d(self): + + interp, atoms = interpolateModel(ANM, NODES, ATOMS) + assert_equal(interp._getArray().shape, ANM._getArray()[EXT3D, :].shape) + assert_equal(atoms.numAtoms(), NODES_AA.numAtoms())