Merge branch 'master' of github.com:prody/ProDy into clustenm_app

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)')
 

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

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

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