(by claude) Fix PDB sequence clusters: migrate off deprecated RCSB bc-*.out files

prody/proteins/pdbclusters.py

@@ -41,6 +41,9 @@ def loadPDBClusters(sqid=None):
                                 'bc-{0}.out.gz'.format(sqid))
         if not os.path.isfile(filename):
             fetchPDBClusters(sqid)
+        if not os.path.isfile(filename):
+            raise IOError('PDB sequence clusters for {0}% sequence identity '
+                          'could not be downloaded.'.format(sqid))
 
         if PDB_CLUSTERS_UPDATE_WARNING:
             import time
@@ -57,7 +60,12 @@ def loadPDBClusters(sqid=None):
         for cluster_str in clusters_str.split('\n'):
             cluster_str = cluster_str.strip()
             if len(cluster_str):
-                cluster = [tuple(item.split('_')) for item in cluster_str.split()]
+                # split on the last underscore only: identifiers may themselves
+                # contain underscores (e.g. computed structure models such as
+                # ``AF_AFP12345F1``), so each member is an (identifier, entity)
+                # pair where *entity* is the trailing entity number.
+                cluster = [tuple(item.rsplit('_', 1))
+                           for item in cluster_str.split()]
                 clusters.append(cluster)
 
         PDB_CLUSTERS[sqid] = clusters
@@ -69,22 +77,21 @@ def loadPDBClusters(sqid=None):
         return clusters
 
 
-def listPDBCluster(pdb, ch, sqid=95):
-    """Returns the PDB sequence cluster that contains chain *ch* in structure
-    *pdb* for sequence identity level *sqid*.  PDB sequence cluster will be
-    returned in as a list of tuples, e.g. ``[('1XXX', 'A'), ]``.  Note that
-    PDB clusters individual chains, so the same PDB identifier may appear
-    twice in the same cluster if the corresponding chain is present in the
-    structure twice.
+def listPDBCluster(pdb, entity=1, sqid=95):
+    """Returns the PDB sequence cluster that contains polymer *entity* in
+    structure *pdb* for sequence identity level *sqid*.  The cluster is
+    returned as a list of tuples, e.g. ``[('1XXX', '1'), ]``.  Note that RCSB
+    now clusters by polymer *entity* (not by individual chain), so *entity* is
+    an entity identifier such as ``1`` and the same PDB identifier may appear
+    more than once across clusters if the structure has multiple entities.
 
     Before this function is used, :func:`fetchPDBClusters` needs to be called.
     This function will load the PDB sequence clusters for *sqid* automatically
     using :func:`loadPDBClusters`."""
 
     assert isinstance(pdb, str) and len(pdb) == 4, \
         'pdb must be 4 char long string'
-    assert isinstance(ch, str) and len(ch) == 1, \
-        'ch must be a one char long string'
+    entity = str(entity)
     try:
         sqid = int(sqid)
     except TypeError:
@@ -96,23 +103,24 @@ def listPDBCluster(pdb, ch, sqid=95):
     if clusters is None:
         loadPDBClusters(sqid)
         clusters = PDB_CLUSTERS[sqid]
-    pdb_ch = (pdb.upper(), ch)
+    pdb_entity = (pdb.upper(), entity)
 
     for cluster in clusters:
-        if pdb_ch in cluster:
+        if pdb_entity in cluster:
             return cluster
-    return 
+    return
 
 def fetchPDBClusters(sqid=None):
     """Retrieve PDB sequence clusters.  PDB sequence clusters are results of
-    the weekly clustering of protein chains in the PDB generated by blastclust.
-    They are available at FTP site: https://cdn.rcsb.org/resources/sequence/clusters/
+    the regular clustering of polymer entities in the PDB.  They are available
+    at: https://cdn.rcsb.org/resources/sequence/clusters/ as
+    :file:`clusters-by-entity-{sqid}.txt` (the legacy :file:`bc-{sqid}.out`
+    files have been deprecated by RCSB).
 
-    This function will download about 10 Mb of data and save it after
-    compressing in your home directory in :file:`.prody/pdbclusters`.
-    Compressed files will be less than 4 Mb in size.  Cluster data can
-    be loaded using :func:`loadPDBClusters` function and be accessed
-    using :func:`listPDBCluster`."""
+    This function will download the cluster data and save it after compressing
+    in your home directory in :file:`.prody/pdbclusters`.  Cluster data can be
+    loaded using :func:`loadPDBClusters` function and be accessed using
+    :func:`listPDBCluster`."""
 
     if sqid is not None:
         if isListLike(sqid):
@@ -134,18 +142,29 @@ def fetchPDBClusters(sqid=None):
                     '_prody_fetchPDBClusters')
     count = 0
     for i, x in enumerate(keys):
-        filename = 'bc-{0}.out'.format(x)
-        url = ('https://cdn.rcsb.org/resources/sequence/clusters/' + filename)
+        # RCSB deprecated the legacy bc-*.out files; the current data is served
+        # as clusters-by-entity-*.txt (clustered by polymer entity, not chain).
+        urlname = 'clusters-by-entity-{0}.txt'.format(x)
+        url = ('https://cdn.rcsb.org/resources/sequence/clusters/' + urlname)
         try:
             inp = openURL(url)
         except IOError:
             LOGGER.warning('Clusters at {0}% sequence identity level could '
                            'not be downloaded.'.format(x))
             continue
         else:
-            out = openFile(filename+'.gz', 'w', folder=PDB_CLUSTERS_PATH)
-            out.write(inp.read())
+            data = inp.read()
             inp.close()
+            # Guard against empty bodies or HTML error pages (e.g. a 404 that
+            # is returned with a 200 status) being saved as valid cluster data
+            # and reported as a successful download.
+            if not data or pystr(data[:200]).lstrip().startswith('<'):
+                LOGGER.warning('Clusters at {0}% sequence identity level could '
+                               'not be downloaded.'.format(x))
+                continue
+            out = openFile('bc-{0}.out.gz'.format(x), 'w',
+                           folder=PDB_CLUSTERS_PATH)
+            out.write(data)
             out.close()
             count += 1
         LOGGER.update(i, label='_prody_fetchPDBClusters')

prody/tests/proteins/test_pdbclusters.py

@@ -0,0 +1,107 @@
+"""This module contains unit tests for :mod:`~prody.proteins.pdbclusters`."""
+
+import os
+
+from numpy.testing import *
+
+from prody.utilities import importDec
+dec = importDec()
+
+from prody import fetchPDBClusters, loadPDBClusters, listPDBCluster
+from prody import LOGGER, SETTINGS, getPackagePath
+from prody.proteins import pdbclusters
+from prody.tests import unittest, TEMPDIR
+
+LOGGER.verbosity = 'none'
+
+SQID = 40
+
+# Cluster files are stored under getPackagePath(); when the package path is not
+# set/writable (e.g. a fresh CI environment) prody prompts for it via input(),
+# which fails under captured output ("reading from stdin while output is
+# captured!").  Point the package path at the writable test TEMPDIR so the tests
+# never block on a prompt; the original value is restored afterwards.
+_ORIGINAL_PACKAGE_PATH = None
+
+
+def setUpModule():
+    global _ORIGINAL_PACKAGE_PATH
+    _ORIGINAL_PACKAGE_PATH = SETTINGS.get('package_path', None)
+    SETTINGS['package_path'] = TEMPDIR
+
+
+def tearDownModule():
+    SETTINGS['package_path'] = _ORIGINAL_PACKAGE_PATH
+
+
+class TestFetchPDBClusters(unittest.TestCase):
+    """Test downloading and querying the RCSB sequence clusters.  RCSB serves
+    these as :file:`clusters-by-entity-{sqid}.txt` (the legacy
+    :file:`bc-{sqid}.out` files were deprecated)."""
+
+    @dec.slow
+    def testFetchSavesUsableFile(self):
+
+        fetchPDBClusters(SQID)
+        filename = os.path.join(getPackagePath(), 'pdbclusters',
+                                'bc-{0}.out.gz'.format(SQID))
+        self.assertTrue(os.path.isfile(filename),
+                        'fetchPDBClusters did not save a cluster file')
+        self.assertGreater(os.path.getsize(filename), 100000,
+                           'saved cluster file is implausibly small '
+                           '(likely an error page, not cluster data)')
+
+    @dec.slow
+    def testListPDBClusterEntityFormat(self):
+
+        fetchPDBClusters(SQID)
+        # 101M is myoglobin; its single polymer entity (1) belongs to a large
+        # cluster at low sequence identity.
+        cluster = listPDBCluster('101M', 1, sqid=SQID)
+        self.assertIsNotNone(cluster,
+                             'listPDBCluster found no cluster for 101M entity 1')
+        self.assertIn(('101M', '1'), cluster,
+                      'query entity is missing from its own cluster')
+        self.assertGreater(len(cluster), 1,
+                           'cluster unexpectedly contains a single member')
+        # Every member must be an (identifier, entity) pair, even for computed
+        # structure models whose identifiers contain underscores (AF_..., MA_...).
+        self.assertTrue(all(len(member) == 2 for member in cluster),
+                        'cluster members are not (identifier, entity) pairs')
+        # New data clusters by polymer entity, so the trailing field is an
+        # entity number (e.g. "1"), not a chain letter (e.g. "A").
+        self.assertTrue(all(ent.isdigit() for _, ent in cluster),
+                        'cluster members are not in IDENTIFIER_ENTITY format')
+
+
+class TestLoadPDBClustersFailure(unittest.TestCase):
+    """A failed download must raise a clear error rather than crashing later in
+    :func:`os.path.getmtime` with a :exc:`FileNotFoundError`."""
+
+    def testFailedDownloadRaisesIOError(self):
+
+        sqid = 50
+        filename = os.path.join(getPackagePath(), 'pdbclusters',
+                                'bc-{0}.out.gz'.format(sqid))
+        orig_openURL = pdbclusters.openURL
+        orig_clusters = pdbclusters.PDB_CLUSTERS[sqid]
+        existed = os.path.isfile(filename)
+        if existed:
+            os.rename(filename, filename + '.bak')
+
+        def _fail(*args, **kwargs):
+            raise IOError('forced download failure')
+
+        pdbclusters.openURL = _fail
+        pdbclusters.PDB_CLUSTERS[sqid] = None
+        try:
+            self.assertRaises(IOError, loadPDBClusters, sqid)
+        finally:
+            pdbclusters.openURL = orig_openURL
+            pdbclusters.PDB_CLUSTERS[sqid] = orig_clusters
+            if existed:
+                os.rename(filename + '.bak', filename)
+
+
+if __name__ == '__main__':
+    unittest.main()