diff --git a/prody/proteins/pdbclusters.py b/prody/proteins/pdbclusters.py index b4452543e..be2f2fc0a 100644 --- a/prody/proteins/pdbclusters.py +++ b/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,13 +77,13 @@ 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 @@ -83,8 +91,7 @@ def listPDBCluster(pdb, ch, sqid=95): 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,8 +142,10 @@ 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: @@ -143,9 +153,18 @@ def fetchPDBClusters(sqid=None): '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') diff --git a/prody/tests/proteins/test_pdbclusters.py b/prody/tests/proteins/test_pdbclusters.py new file mode 100644 index 000000000..f2e75329b --- /dev/null +++ b/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()