Merge branch 'prody-cavifinder' of github.com:karolamik13/ProDy into …

…prody-cavifinder

prody/dynamics/clustenm.py

@@ -1076,7 +1076,7 @@ def run(self, cutoff=15., n_modes=3, gamma=1., n_confs=50, rmsd=1.0,
         self._sparse = kwargs.get('sparse', False)
         self._kdtree = kwargs.get('kdtree', False)
         self._turbo = kwargs.get('turbo', False)
-        self._nproc = kwargs.pop('nproc', None)
+        self._nproc = kwargs.pop('nproc', 0)
         if kwargs.get('zeros', False):
             LOGGER.warn('ClustENM cannot use zero modes so ignoring this kwarg')
 

prody/proteins/cifheader.py

@@ -745,8 +745,11 @@ def _getReference(lines):
             except:
                 continue
             if what == 'AUTH':
-                surname, initials = value.split(',')
-                author = initials+surname
+                try:
+                    surname, initials = value.split(',')
+                    author = initials+surname
+                except ValueError:
+                    author = value
                 authors.append(author.strip().upper())
 
         ref['authors'] = authors

prody/proteins/interactions.py

@@ -711,9 +711,11 @@ def calcChHydrogenBonds(atoms, **kwargs):
 
         ChainsHBs = [ i for i in HBS_calculations if str(i[2]) != str(i[5]) ]
         if not ChainsHBs:
-            ligand_name = list(set(atoms.select('all not protein and not ion').getResnames()))[0]
-            ChainsHBs = [ ii for ii in HBS_calculations if ii[0][:3] == ligand_name or ii[3][:3] == ligand_name ]
-
+            ligand_sel = atoms.select('all not protein and not ion')
+            if ligand_sel:
+                ligand_name = list(set(ligand_sel.getResnames()))[0]
+                ChainsHBs = [ ii for ii in HBS_calculations if ii[0][:3] == ligand_name or ii[3][:3] == ligand_name ]
+
         return ChainsHBs 
 
 
@@ -3559,8 +3561,43 @@ def saveInteractionsPDB(self, **kwargs):
             writePDB('filename', atoms, **kw)
             LOGGER.info('PDB file saved.')
 
-
-    def getFrequentInteractors(self, contacts_min=3):
+
+    def getInteractors(self, residue_name):
+        """ Provide information about interactions for a particular residue
+        
+        :arg residue_name: name of a resiude
+                            example: LEU234A, where A is a chain name
+        :type residue_name: str
+        """
+
+        atoms = self._atoms   
+        interactions = self._interactions
+
+        InteractionsMap = np.empty([atoms.select('name CA').numAtoms(),atoms.select('name CA').numAtoms()], dtype=object)
+        resIDs = list(atoms.select('name CA').getResnums())
+        resChIDs = list(atoms.select('name CA').getChids())
+        resIDs_with_resChIDs = list(zip(resIDs, resChIDs))
+        interaction_type = ['hb','sb','rb','ps','pc','hp','dibs']
+        ListOfInteractions = []
+
+        for nr,i in enumerate(interactions):
+            if i != []:
+                for ii in i: 
+                    m1 = resIDs_with_resChIDs.index((int(ii[0][3:]),ii[2]))
+                    m2 = resIDs_with_resChIDs.index((int(ii[3][3:]),ii[5]))
+                    ListOfInteractions.append(interaction_type[nr]+':'+ii[0]+ii[2]+'-'+ii[3]+ii[5])
+
+        aa_ListOfInteractions = []
+        for i in ListOfInteractions:
+            inter = i.split(":")[1:][0]
+            if inter.split('-')[0] == residue_name or inter.split('-')[1] == residue_name:
+                LOGGER.info(i)
+                aa_ListOfInteractions.append(i)
+
+        return aa_ListOfInteractions
+
+
+    def getFrequentInteractors(self, contacts_min=2):
         """Provide a list of residues with the most frequent interactions based 
         on the following interactions:
             (1) Hydrogen bonds (hb)
@@ -3572,7 +3609,7 @@ def getFrequentInteractors(self, contacts_min=3):
             (7) Disulfide bonds (disb)
         
         :arg contacts_min: Minimal number of contacts which residue may form with other residues, 
-                           by default 3.
+                           by default 2.
         :type contacts_min: int  """
 
         atoms = self._atoms   
@@ -3589,36 +3626,59 @@ def getFrequentInteractors(self, contacts_min=3):
                 for ii in i: 
                     m1 = resIDs_with_resChIDs.index((int(ii[0][3:]),ii[2]))
                     m2 = resIDs_with_resChIDs.index((int(ii[3][3:]),ii[5]))
-                    InteractionsMap[m1][m2] = interaction_type[nr]+':'+ii[0]+ii[2]+'-'+ii[3]+ii[5]
+
+                    if InteractionsMap[m1][m2] is None:
+                        InteractionsMap[m1][m2] = []
+
+                    InteractionsMap[m1][m2].append(interaction_type[nr] + ':' + ii[0] + ii[2] + '-' + ii[3] + ii[5])
 
-        ListOfInteractions = [ list(filter(None, InteractionsMap[:,j])) for j in range(len(interactions[0])) ]
-        ListOfInteractions = list(filter(lambda x : x != [], ListOfInteractions))
-        ListOfInteractions = [k for k in ListOfInteractions if len(k) >= contacts_min ]
-        ListOfInteractions_list = [ (i[0].split('-')[-1], [ j.split('-')[0] for j in i]) for i in ListOfInteractions ]
-        LOGGER.info('The most frequent interactions between:')
-        for res in ListOfInteractions_list:
+        ListOfInteractions = [list(filter(None, [row[j] for row in InteractionsMap])) for j in range(len(InteractionsMap[0]))]
+        ListOfInteractions = list(filter(lambda x: x != [], ListOfInteractions))
+        ListOfInteractions_flattened = [j for sublist in ListOfInteractions for j in sublist]
+
+        swapped_ListOfInteractions_list = []
+        for interaction_group in ListOfInteractions_flattened:
+            swapped_group = []
+            for interaction in interaction_group:
+                interaction_type, pair = interaction.split(':')
+                swapped_pair = '-'.join(pair.split('-')[::-1])
+                swapped_group.append("{}:{}".format(interaction_type, swapped_pair))
+            swapped_ListOfInteractions_list.append(swapped_group)
+
+        doubleListOfInteractions_list = ListOfInteractions_flattened+swapped_ListOfInteractions_list
+        ListOfInteractions_list = [(i[0].split('-')[-1], [j.split('-')[0] for j in i]) for i in doubleListOfInteractions_list]
+
+        merged_dict = {}
+        for aa, ii in ListOfInteractions_list:
+            if aa in merged_dict:
+                merged_dict[aa].extend(ii)
+            else:
+                merged_dict[aa] = ii
+
+        ListOfInteractions_list = [(key, value) for key, value in merged_dict.items()] 
+        ListOfInteractions_list2 = [k for k in ListOfInteractions_list if len(k[-1]) >= contacts_min]
+
+        for res in ListOfInteractions_list2:
             LOGGER.info('{0}  <--->  {1}'.format(res[0], '  '.join(res[1])))
 
-        LOGGER.info('Legend: hb-hydrogen bond, sb-salt bridge, rb-repulsive ionic bond, ps-Pi stacking interaction,'
-                             'pc-Cation-Pi interaction, hp-hydrophobic interaction, dibs-disulfide bonds')
+        LOGGER.info('\nLegend: hb-hydrogen bond, sb-salt bridge, rb-repulsive ionic bond, ps-Pi stacking interaction,'
+                             '\npc-Cation-Pi interaction, hp-hydrophobic interaction, dibs-disulfide bonds')
 
         try:
             from toolz.curried import count
         except ImportError:
             LOGGER.warn('This function requires the module toolz')
             return
 
-        LOGGER.info('The biggest number of interactions: {}'.format(max(map(count, ListOfInteractions))))
-
-        return ListOfInteractions_list
+        return ListOfInteractions_list2
 
 
-    def showFrequentInteractors(self, cutoff=5, **kwargs):
+    def showFrequentInteractors(self, cutoff=4, **kwargs):
         """Plots regions with the most frequent interactions.
         
         :arg cutoff: minimal score per residue which will be displayed.
                      If cutoff value is to big, top 30% with the higest values will be returned.
-                     Default is 5.
+                     Default is 4.
         :type cutoff: int, float
 
         Nonstandard resiudes can be updated in a following way:

prody/proteins/localpdb.py

@@ -212,12 +212,15 @@ def fetchPDB(*pdb, **kwargs):
     if len(pdb) == 1 and isinstance(pdb[0], list):
         pdb = pdb[0]
 
-    identifiers = checkIdentifiers(*pdb)
-
     folder = kwargs.get('folder', '.')
     compressed = kwargs.get('compressed')
     format_ = kwargs.get('format', 'pdb')
 
+    if format_ != 'emd':
+        identifiers = checkIdentifiers(*pdb)
+    else:
+        identifiers = pdb
+
     # check *folder* specified by the user, usually pwd ('.')
     filedict = findPDBFiles(folder, compressed=compressed, 
                             format=format_)

prody/proteins/pdbfile.py

@@ -214,13 +214,15 @@ def _parsePDB(pdb, **kwargs):
         if filename is None:
             try:
                 LOGGER.warn("Trying to parse mmCIF file instead")
+                chain = kwargs.pop('chain', chain)
                 return parseMMCIF(pdb+chain, **kwargs)
-            except:
+            except OSError:
                 try:
                     LOGGER.warn("Trying to parse EMD file instead")
+                    chain = kwargs.pop('chain', chain)
                     return parseEMD(pdb+chain, **kwargs)
                 except:
-                    raise IOError('PDB file for {0} could not be downloaded.'
+                    raise IOError('PDB file for {0} could not be parsed.'
                                 .format(pdb))
         pdb = filename
     if title is None:

prody/proteins/wwpdb.py

@@ -254,7 +254,22 @@ def fetchPDBviaHTTP(*pdb, **kwargs):
     output_folder = kwargs.pop('folder', None)
     compressed = bool(kwargs.pop('compressed', True))
 
-    extension = '.pdb'
+    format = kwargs.get('format', 'pdb')
+    noatom = bool(kwargs.pop('noatom', False))
+    if format == 'pdb':
+        extension = '.pdb'
+    elif format == 'xml':
+        if noatom:
+            extension = '-noatom.xml'
+        else:
+            extension = '.xml'
+    elif format == 'cif':
+        extension = '.cif'
+    elif format == 'emd' or format == 'map':
+        extension = '.map'
+    else:
+        raise ValueError(repr(format) + ' is not valid format')
+
     local_folder = pathPDBFolder()
     if local_folder:
         local_folder, is_divided = local_folder
@@ -294,7 +309,10 @@ def fetchPDBviaHTTP(*pdb, **kwargs):
             filenames.append(None)
             continue
         try:
-            handle = openURL(getURL(pdb))
+            url = getURL(pdb)
+            if kwargs.get('format', 'pdb') != 'pdb':
+                url = url.replace('.pdb', extension)
+            handle = openURL(url)
         except Exception as err:
             LOGGER.warn('{0} download failed ({1}).'.format(pdb, str(err)))
             failure += 1

prody/tests/proteins/test_wwpdb.py

@@ -15,15 +15,15 @@
 LOGGER.verbosity = 'none'
 
 
-class TestFTP(unittest.TestCase):
+class TestHTTP(unittest.TestCase):
 
     def setUp(self):
 
         self.pdb = ['1ubi', '1aar', 'arg', 1234]
         self.fns = []
         self.len = [683, 1218, None, None]
-        self.fetch = fetchPDBviaFTP
-        self.protocol = 'FTP'
+        self.fetch = fetchPDBviaHTTP
+        self.protocol = 'HTTP'
 
 
     @dec.slow
@@ -69,12 +69,3 @@ def tearDown(self):
                 pass
             except:
                 pass
-
-class TestHTTP(TestFTP):
-
-    def setUp(self):
-
-        TestFTP.setUp(self)
-        self.fetch = fetchPDBviaHTTP
-        self.protocol = 'HTTP'
-

prody/utilities/catchall.py

@@ -739,6 +739,8 @@ def showMatrix(matrix, x_array=None, y_array=None, **kwargs):
         im = ax3.imshow(matrix, aspect=aspect, vmin=vmin, vmax=vmax,
                         norm=norm, cmap=cmap, origin=origin, **kwargs)
     else:
+        zeros_matrix = np.zeros(matrix.shape)
+        im = ax3.imshow(zeros_matrix, vmin=-1, vmax=1, cmap='seismic')
         plot_list = []
         for rows,cols in zip(np.where(matrix!=0)[0],np.where(matrix!=0)[1]):
             plot_list.append([cols,rows,matrix[rows,cols]])