Cavity examples

Compare cavities in proteins using histograms

This example allows you to detect and compare cavities from PDB files stored in a user-defined directory. The cavity histograms comparison uses feature distance histograms and the similarity scores are written out to a .csv file. Note that cavity detection is the rate-determining step and therefore it is best practice to store the detected cavities as XML files or in a cavity database for future reuse and analysis.

#!/usr/bin/env python
#
# This script can be used for any purpose without limitation subject to the
# conditions at http://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx
#
# This permission notice and the following statement of attribution must be
# included in all copies or substantial portions of this script.
#
# 2015-02-25: created by the Cambridge Crystallographic Data Centre
#
'''
    cavity_compare.py - compare the cavities in all proteins (PDB files) in a certain directory.

    Usage: python cavity_compare.py pdb_directory
'''
#####################################################################

import sys
import argparse
import os

from ccdc.cavity import Cavity
from glob import glob

#####################################################################


def parse_args():
    usage_txt = ('Detect cavities from PDB files stored in a user defined directory'
                 ' and perform an all-by-all comparison of the detected cavities.')
    parser = argparse.ArgumentParser(
        description=usage_txt,
        formatter_class=argparse.ArgumentDefaultsHelpFormatter)

    parser.add_argument('-d', '--directory',
                        help='Directory of PDB files to find cavities in.',
                        required=True)
    parser.add_argument('-m', '--maximum-incomplete-residues-per-chain',type=int,
                        help='Number of incomplete residues that are permissible when generating cavities',
                        default=100) # High default - presumably we want all cavities included unless they are really poor otherwise they would not be in the directory

    return parser.parse_args()


def main():
    # Get command line argument:
    args = parse_args()
    pdb_dir = args.directory

    # Collect PDB files
    pdbs = glob(os.path.join(pdb_dir, '*.pdb'))

    # At least two PDB files are needed to perform comparisons
    if len(pdbs) < 2:
        print('Error: At least two PDB files are required. Not enough found in %s' % pdb_dir)
        sys.exit(1)

    print(' -> Detecting the cavities for all PDB structures in %s' % pdb_dir)
    cavs = []
    for pdb in pdbs:
        cavs += Cavity.from_pdb_file(pdb, maximum_incomplete_residues_per_chain = args.maximum_incomplete_residues_per_chain)

    # Generate CSV output
    with open('similarities.csv', 'w') as f:
        # header:
        cav_ids = [''] + [cav.identifier for cav in cavs]
        f.write(','.join(cav_ids) + '\n')

        # rows:
        for i in range(0, len(cavs)):
            cells = [cavs[i].identifier]
            for j in range(0, len(cavs)):
                score = cavs[i].compare(
                    cavs[j],
                    comparison_method=Cavity.ComparisonMethod.CAVITY_HISTOGRAMS_COMPARISON)
                cells.append(str(round(score, 3)))
            f.write(','.join(cells) + '\n')

        print(' -> Similarity matrix written to similarities.csv')

#####################################################################


if __name__ == '__main__':
    main()

Overlay pairs of proteins based on their contained cavities

This example allows you to detect and compare cavities from pairs of PDB files. The script reads and detects cavities in 2 protein structures then overlays pairs of cavities. The resultant transformation matrices are applied to the input proteins and new output files are written that contain the superimposed protein with the transformation applied.

#!/usr/bin/env python
#
# This script can be used for any purpose without limitation subject to the
# conditions at http://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx
#
# This permission notice and the following statement of attribution must be
# included in all copies or substantial portions of this script.
#
# 2015-02-25: created by the Cambridge Crystallographic Data Centre
#
'''
    cavity_pair_view.py - generate cavity overlays for all cavities in a  pair of protein structures and then generate mol2 files of the proteins in the overlaid configuration. Also writes a pairwise similarity matrix for the cavities

    Usage: python cavity_pair_view.py <pdb file 1> <pdb file 2>
'''
#####################################################################

import sys
import argparse
import os

from ccdc.cavity import Cavity
from ccdc.protein import Protein
from ccdc.io import MoleculeWriter

from glob import glob

#####################################################################


def parse_args():
    usage_txt = ('Detect cavities from a pair of PDB files and then overlays them based on the best match')
    parser = argparse.ArgumentParser(
        description=usage_txt,
        formatter_class=argparse.ArgumentDefaultsHelpFormatter)

    parser.add_argument('filenames', metavar='<filename>',type=str, nargs='+',
                        help='Filenames of PDB files to find cavities in.')
    return parser.parse_args()


def main():
    # Get command line argument:
    args = parse_args()


    # At least two PDB files are needed to perform comparisons
    if len(args.filenames) != 2:
        print('Error: You must specify exactly two PDB filenames')
        sys.exit(1)

    print(' -> Detecting the cavities for all PDB structures in %s and %s' % ( args.filenames[0],args.filenames[1] ))


    file_lookup = {}

    first_cavs = Cavity.from_pdb_file(args.filenames[0])
    second_cavs =  Cavity.from_pdb_file(args.filenames[1])

    for cav in first_cavs:
        file_lookup[cav.identifier] = args.filenames[0]

    for cav in second_cavs:
        file_lookup[cav.identifier] = args.filenames[1]

    # Generate CSV output
    with open('similarities.csv', 'w') as f:
        # header:
        cav_ids = [''] + [cav.identifier for cav in second_cavs]
        f.write(','.join(cav_ids) + '\n')

        for first_cav in first_cavs:
            cells = [first_cav.identifier]
            for second_cav in second_cavs:

                comparison_object = first_cav.compare(
                    second_cav,
                    comparison_method=Cavity.ComparisonMethod.CAVITY_GRAPH_COMPARISON)

                score = comparison_object.score
                cells.append(str(round(score, 3)))
                filename = first_cav.identifier + "_" + second_cav.identifier + ".mol2"

                m = Protein.from_file(file_lookup[second_cav.identifier])
                m.transform(comparison_object.transformation_matrix)
                print(' -> Writing %s - a file containing PDB entry %s - superimposition using cavity %s in %s onto cavity %s' % ( filename,m.identifier,second_cav.identifier,m.identifier, first_cav.identifier))

                er = MoleculeWriter(filename)
                er.write(m)

            f.write(','.join(cells) + '\n')

        print(' -> Similarity matrix written to similarities.csv')

#####################################################################


if __name__ == '__main__':
    main()

Compare XML cavity files

This script allows you to compare cavity information stored in XML cavity files. You can switch between a cavity histograms comparison, a fast cavity graph comparison and a cavity graph comparison, using the appropriate parameter in the command line argument. The script furthermore provides a more elaborate example on how to make use of Python’s argument parser.

#!/usr/bin/env python
#
# This script can be used for any purpose without limitation subject to the
# conditions at http://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx
#
# This permission notice and the following statement of attribution must be
# included in all copies or substantial portions of this script.
#
# 2015-02-25: created by the Cambridge Crystallographic Data Centre
#
'''
    cavity_compare_args.py - compare a cavity against a set of other
        cavities using command line arguments.

    Usage: python cavity_compare_args.py -d C:/path/to/xml_cavities
                                         -q query_cavity
                                         -c compare_cavities
                                         -m method
'''
#####################################################################

import argparse
from ccdc.cavity import Cavity
import os

#####################################################################


def parse_args():
    usage_txt = 'Compare a cavity against a set of other cavities using command line arguments.'
    parser = argparse.ArgumentParser(description=usage_txt,
                                     formatter_class=argparse.ArgumentDefaultsHelpFormatter)
    parser.add_argument('-m', '--method', default='GRAPH',
                        choices=['HISTOGRAMS', 'FAST_GRAPH', 'GRAPH'],
                        help='Comparison method to use.')
    parser.add_argument('-q', '--query',
                        help='Query cavity name (e.g. "pdb3jw3.1" for cavity_pdb3jw3.1.xml).',
                        required=True)
    parser.add_argument('-c', '--comparisons', nargs='+', required=True,
                        help='Cavity names with which to compare, separated by spaces.')
    parser.add_argument('-d', '--directory',
                        help='Directory of XML cavity files.')
    return parser.parse_args()


def main():
    # Get command line arguments:
    args = parse_args()
    xml_dir = args.directory
    query = args.query
    method = args.method
    comps = args.comparisons

    # Load the query cavity:
    query_cav = Cavity.from_xml_file(os.path.join(xml_dir, 'cavity_%s.xml' % query))

    # Do the comparisons:
    if method == 'HISTOGRAMS':
        # Create histograms for cavity histograms comparison
        _ = query_cav.cavity_distance_histograms()
        print('=== Cavity histogram comparison ===')
        print('Query\tComparison\tScore')
        for comp in comps:
            try:
                comp_cav = Cavity.from_xml_file(os.path.join(xml_dir, 'cavity_%s.xml' % comp))
                score = query_cav.compare(
                    comp_cav,
                    comparison_method=Cavity.ComparisonMethod.CAVITY_HISTOGRAMS_COMPARISON)
                print('\t'.join([query, comp, str(score)]))
            except RuntimeError as err:
                print('\t'.join([query, comp, 'Error: %s' % err]))
    elif method == 'FAST_GRAPH':
        print('=== Fast graph comparison ===')
        print('Query\tComparison\tScore')
        for comp in comps:
            try:
                comp_cav = Cavity.from_xml_file(os.path.join(xml_dir, 'cavity_%s.xml' % comp))
                score = query_cav.compare(
                    comp_cav,
                    comparison_method=Cavity.ComparisonMethod.FAST_CAVITY_GRAPH_COMPARISON).score
                print('\t'.join([query, comp, str(score)]))
            except RuntimeError as err:
                print('\t'.join([query, comp, 'Error: %s' % err]))
    elif method == 'GRAPH':
        print('=== Graph comparison ===')
        print('Query\tComparison\tScore')
        for comp in comps:
            try:
                comp_cav = Cavity.from_xml_file(os.path.join(xml_dir, 'cavity_%s.xml' % comp))
                score = query_cav.compare(
                    comp_cav,
                    comparison_method=Cavity.ComparisonMethod.CAVITY_GRAPH_COMPARISON).score
                print('\t'.join([query, comp, str(score)]))
            except RuntimeError as err:
                print('\t'.join([query, comp, 'Error: %s' % err]))

#####################################################################


if __name__ == '__main__':
    main()

Compare cavities using multiple cores

This example shows how to use multiple cores to speed up the comparison of large numbers of cavities.

#!/usr/bin/env python
#
# This script can be used for any purpose without limitation subject to the
# conditions at http://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx
#
# This permission notice and the following statement of attribution must be
# included in all copies or substantial portions of this script.
#
# 2018-02-02: created by the Cambridge Crystallographic Data Centre
#
'''
    cavity_compare_db_multiple_cores.py - Pairwise comparison of cavities
        taken from a database using multiple cores.

    Usage: python cavity_compare_db_multiple_cores.py -n n_processes
                                                      -q query_cavity
                                                      -v volume_min
                                                      -w volume_max
'''
#####################################################################

import argparse
from multiprocessing import Pool, get_context
import os
from ccdc.cavity import *

#####################################################################

process_db = None


def initializer(initargs):
    """  Create a database connection for each process. """
    global process_db
    process_db = CavityDatabase(initargs)


def comparison(cavities):
    """ Compare two cavities in a tuple and output the scoring. """
    query_cav = process_db.get_cavity_by_name(cavities[0])
    comp_cav = process_db.get_cavity_by_name(cavities[1])
    try:
        result = query_cav.compare(
            comp_cav,
            comparison_method=Cavity.ComparisonMethod.CAVITY_GRAPH_COMPARISON)
        print('%s\t%s' % ('\t'.join(cavities[:2]), result.score))
    except RuntimeError as exc:
        print('%s\tComparison unsuccessful: %s' % ('\t'.join(cavities[:2]), str(exc)))


def mp_handler(comparisons_list, n_processes, db_path):
    """ Create a pool of worker processes and distribute the cavity pairs between them. """
    with get_context('spawn').Pool(n_processes, initializer, initargs=(db_path,)) as p:
        p.map(comparison, comparisons_list)


def parse_args():
    usage_txt = 'Pairwise comparison of cavities taken from a database using multiple cores.'
    parser = argparse.ArgumentParser(
        description=usage_txt,
        formatter_class=argparse.ArgumentDefaultsHelpFormatter
    )
    parser.add_argument(
        '-n', '--n_processes',
        help='Number of parallel processes',
        default=2
    )
    parser.add_argument(
        '-q', '--query',
        help='Query cavity name',
        required=True
    )
    parser.add_argument(
        '-v', '--volume_min',
        help='Minimum volume for target cavities search',
        required=True
    )
    parser.add_argument(
        '-w', '--volume_max',
        help='Maximum volume for target cavities search',
        required=True
    )
    parser.add_argument(
        '-d', '--db_path',
        help='Cavity database to use as source data.',
        default=os.path.join(CavityDatabase.drugbank_database_dir(), 'drugbank_cavities.sqlite')
    )
    return parser.parse_args()


def main():
    args = parse_args()
    n_processes = int(args.n_processes)
    query = args.query
    volume_min = float(args.volume_min)
    volume_max = float(args.volume_max)
    db_path = args.db_path

    if not os.path.exists(db_path):
        print('Cavity database %s does not exist!' % db_path)
        exit(1)

    print('Comparing cavity %s with a volume range of %s-%s in database %s using %s processes.' %
          (query, volume_min, volume_max, db_path, n_processes))

    # Create database connection for initial search
    db = CavityDatabase(db_path)

    query_cav = db.get_cavity_by_name(query)
    if query_cav is None:
        raise RuntimeError('Query cavity {0} not found in database'.format(query))

    # Define database search settings
    settings = CavityDatabase.Settings()
    settings.volume_range = [volume_min, volume_max]

    # Search using settings
    comparisons = db.search(settings=settings)

    # Make list of comparison cavities fulfilling the search criteria
    comparisons_list = []
    for c in comparisons:
        comparisons_list.append((query, c.identifier))

    print('Number of comparisons to be performed: %s' % len(comparisons_list))
    print('=== Cavity Graph Comparison ===')
    print('Query\tComparison\tScore')

    mp_handler(comparisons_list, n_processes, db_path=db_path)

#####################################################################


if __name__ == '__main__':
    import ccdc.utilities
    ccdc.utilities._fix_multiprocessing_on_macos()
    main()

Compare cavities stored in an SQLite database

Cavity comparison using an SQLite database is faster than processing individual XML files. Refer to this example to set up a comparison using SQLite input.

#!/usr/bin/env python
#
# This script can be used for any purpose without limitation subject to the
# conditions at http://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx
#
# This permission notice and the following statement of attribution must be
# included in all copies or substantial portions of this script.
#
# 2015-02-25: created by the Cambridge Crystallographic Data Centre
#
'''
    cavity_compare_db_pairwise.py - Pairwise comparison of cavities taken from a database.

    Usage: python cavity_compare_db_pairwise.py -q query_cavity
                                                -c comparison1 comparison2 comparison3
'''
#####################################################################

import argparse
from ccdc.cavity import *

#####################################################################


def parse_args():
    usage_txt = 'Pairwise comparison of cavities taken from a database.'
    parser = argparse.ArgumentParser(description=usage_txt,
                                     formatter_class=argparse.ArgumentDefaultsHelpFormatter)
    parser.add_argument('-c', '--comparisons', nargs='+',
                        default=['pdb3dar.1', 'pdb3f8s.1', 'pdb2hha.1', 'pdb21bi.2'],
                        help='Names of cavities in the database to compare against.')
    parser.add_argument('-q', '--query', default='pdb21bi.1',
                        help='The name of the cavity to use as a reference.')
    parser.add_argument('-d', '--database',
                        default=os.path.join(CavityDatabase.drugbank_database_dir(),
                                             'drugbank_cavities.sqlite'),
                        help='Database file to use a source for cavities.')
    return parser.parse_args()


def main():
    args = parse_args()
    db = CavityDatabase(args.database)

    # Load the query cavity:
    query_cav = db.get_cavity_by_name(args.query)

    print('=== Cavity comparison ===')
    print('Query\tComparison\tScore')

    for comp in args.comparisons:
        comp_cav = db.get_cavity_by_name(comp)
        try:
            result = query_cav.compare(
                comp_cav,
                comparison_method=Cavity.ComparisonMethod.CAVITY_GRAPH_COMPARISON)
            # Output the comparison result:
            print('\t'.join([args.query, comp, str(result.score)]))
        except RuntimeError as e:
            print('\t'.join([args.query, comp, 'Failed: %s' % e.message]))

#####################################################################


if __name__ == '__main__':
    main()

A script to compile an SQLite database from your own data is provided as ccdc/utilities/create_cavity_database/create_cavity_database.py.

Download and extract the PDB header for a cavity

This example shows how to download a PDB file from the wwPDB for a cavity and extract header information using Biopython.

#!/usr/bin/env python
#
# This script can be used for any purpose without limitation subject to the
# conditions at http://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx
#
# This permission notice and the following statement of attribution must be
# included in all copies or substantial portions of this script.
#
# 2018-03-01: created by the Cambridge Crystallographic Data Centre
#
'''
    cavity_fetch_pdb_header.py - Retrieve pdb file header information for a cavity
                                 from wwwPDB using Biopython

    Usage: python cavity_fetch_pdb_header.py --query [cavity]
                                             --database [database_path]
'''
#####################################################################

import argparse
from Bio import PDB
from ccdc.cavity import *
from ccdc.utilities import output_file, to_utf8
import os
import requests

import tempfile

#####################################################################


def parse_args():
    """Parse command-line arguments."""
    usage = 'Get pdb file header information for a cavity from the wwwPDB using Biopython.'
    parser = argparse.ArgumentParser(description=usage,
                                     formatter_class=argparse.ArgumentDefaultsHelpFormatter)
    parser.add_argument('-q', '--query', default='pdb3k2h.8',
                        help='The name of the cavity to use as a reference.')
    parser.add_argument('-d', '--database',
                        default=os.path.join(CavityDatabase.drugbank_database_dir(),
                                             'drugbank_cavities.sqlite'),
                        help='Database file to use a source for cavities.')
    return parser.parse_args()


def fetch_from_pdb(pdb_id, filename):
    """Try to retrieve the given structure from the RCSB web database.

    :param pdb_id: (:obj:`str`) The PDB structure identifier to retrieve.
    :param filename: (:obj:`str`) The file name to write the structure to.

    :raises: requests.HTTPError if an error occurs trying to download the structure.
    """

    # Try to retrieve the structure from the online database.
    pdb_req = requests.get('https://files.rcsb.org/download/%s.pdb' % pdb_id)
    # Raise a requests.HTTPError if any HTTP errors have occurred.
    pdb_req.raise_for_status()

    # Write the structure returned from the PDB to the given file name
    with output_file(filename) as pdb_file:
        pdb_file.write(to_utf8(pdb_req.text))


def extract_pdb_id(cavity_id):
    """Convert Relibase chain id to CSD Python API chain id"""
    pdb_id = cavity_id.split('.')[0]
    pdb_id = pdb_id.replace('pdb', '')
    return pdb_id


def safe_key_value(a_dict, key):
    """Return a value corresponding to the specified header key."""
    value = a_dict.get(key, '')
    return value.strip() if isinstance(value, str) else value


class Header(object):
    def __init__(self, structure):
        self.header = structure.header
        self.structure_id = structure.id.upper()
        self.resolution = safe_key_value(self.header, 'resolution')
        self.structure_method = safe_key_value(self.header, 'structure_method')
        self.deposition_date = safe_key_value(self.header, 'deposition_date')
        self.pdb_title = safe_key_value(self.header, 'name').capitalize()
        self.pdb_class = safe_key_value(self.header, 'head').capitalize()
        self.compound_dict = safe_key_value(self.header, 'compound')
        self.source_dict = safe_key_value(self.header, 'source')
        try:
            self.pdb_chain_list = list({chain.id.upper() for chain in structure.get_chains()})
        except Exception:
            self.pdb_chain_list = []

    def __getitem__(self, key):
        return self.header[key]


def main():
    args = parse_args()
    db = CavityDatabase(args.database)

    # Load the query cavity
    query_cavity = db.get_cavity_by_name(args.query)
    if query_cavity is None:
        raise Exception('Failed to find {} in database {}!'.format(args.query, args.database))
    pdb_id = extract_pdb_id(query_cavity.identifier)

    # load the pdb file for the protein from wwwPDB
    output_dir = tempfile.mkdtemp()
    pdb_file = os.path.join(output_dir, '{}.pdb'.format(pdb_id))
    fetch_from_pdb(pdb_id, pdb_file)

    # strict parser, will assign altloc
    struct_parser = PDB.PDBParser(PERMISSIVE=False, get_header=True)
    pdb_structure = struct_parser.get_structure(pdb_id, pdb_file)
    pdb_header = Header(pdb_structure)

    print("Identifier: {}".format(pdb_id))
    print("Title: {}".format(pdb_header.pdb_title))
    print("Class: {}".format(pdb_header.pdb_class))
    print("Resolution: {}".format(pdb_header.resolution))
    print("Structure method: {}".format(pdb_header.structure_method))
    print("Deposition date: {}".format(pdb_header.deposition_date))
    print("Compound dictionary: {}".format(pdb_header.compound_dict))
    print("Source dictionary: {}".format(pdb_header.source_dict))
    print("Chain list: {}".format(pdb_header.pdb_chain_list))


#####################################################################

if __name__ == '__main__':
    main()

Find the most buried cavity area using k-means clustering and use it for subcavity comparisons

This script allows you to locate the most buried area of a cavity using a k-means clustering algorithm and compare it against other cavities. The value of k giving the number of clusters and the minimum size of the clusters can be adjusted within the script.

#!/usr/bin/env python
#
# This script can be used for any purpose without limitation subject to the
# conditions at http://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx
#
# This permission notice and the following statement of attribution must be
# included in all copies or substantial portions of this script.
#
# 2018-01-15: created by the Cambridge Crystallographic Data Centre
#
'''
    cavity_most_buried.py - use pseudocenter depth values and k-means clustering
    to detect the most buried area of a cavity.

    Usage: python cavity_most_buried.py -x c:/path/to/input/xml/directory
                                        -i c:/path/to/input/directory/reference_file.xml
                                        [-v c:/path/to/visualisation/output/directory]
'''
#####################################################################

import argparse
from ccdc.cavity import *
import glob
from math import sqrt
import os
import random

#####################################################################


def parse_args():

    usage_txt = ('Use pseudocenter depth values and k-means clustering to detect',
                 ' the most buried area of a cavity.')
    parser = argparse.ArgumentParser(description=usage_txt,
                                     formatter_class=argparse.ArgumentDefaultsHelpFormatter)

    parser.add_argument(
        '-x', '--xml_directory',
        help='Directory for input XML files.',
        required=True
    )
    parser.add_argument(
        '-i', '--input_file',
        help='Reference cavity XML file.',
        required=True)
    parser.add_argument(
        '-v', '--visualisation_directory',
        help='Output directory for PyMol visualisation file.',
        default=None
    )
    parser.add_argument(
        '-s', '--min_subcavity_size',
        help='The minimum size of the buried subcavity (number of pseudocenters).',
        type=int, default=15
    )
    parser.add_argument(
        '-d', '--min_depth',
        help=('Minimum depth value of the pseudocenters to be used for subcavity detection.',
              ' This is in the interval [0,7] where 7 means most buried.'),
        type=float, default=7.0)
    parser.add_argument(
        '-k', '--k_value',
        help='The number of clusters to be detected.',
        type=int, default=3)
    return parser.parse_args()


def k_means_clustering(features, min_depth, num_clusters):

    def euclidean_dist(p, q):
        return sqrt(sum([pow(p[i] - q[i], 2) for i in range(3)]))

    def get_buried_indices(features, min_depth):
        indices = []
        for i in range(len(features)):
            if features[i].burial >= min_depth:
                indices.append(i)
        return indices

    def get_centroid(cluster):
        x = y = z = 0.0
        for feature in cluster:
            x += feature.coordinates[0]
            y += feature.coordinates[1]
            z += feature.coordinates[2]
        n = len(cluster)
        return (x / n, y / n, z / n)

    buried_indices = get_buried_indices(features, min_depth)

    # Return set of empty clusters if not enough points with the specified
    # depth were found:
    if len(buried_indices) < num_clusters:
        return [[] for i in range(num_clusters)]

    clusters = []

    # Initialisation: generate 'num_clusters' clusters
    # each of which contains one random feature.
    used_positions = []
    for i in range(num_clusters):
        new_cluster = []
        while True:
            pos = random.randint(0, len(buried_indices) - 1)
            if pos not in used_positions:
                new_cluster.append(features[buried_indices[pos]])
                used_positions.append(pos)
                break
        clusters.append(new_cluster)

    # Iteration: re-generate clusters until convergence.
    while True:
        cluster_centers = [get_centroid(c) for c in clusters]
        cluster_sizes = [len(c) for c in clusters]
        clusters = [[] for i in range(num_clusters)]   # clear previous clusters

        for bi in buried_indices:
            closest_cluster_center_index = -1
            closest_distance = 10000.0
            for ci in range(len(cluster_centers)):
                d = euclidean_dist(features[bi].coordinates, cluster_centers[ci])
                if d < closest_distance:
                    closest_cluster_center_index = ci
                    closest_distance = d
            clusters[closest_cluster_center_index].append(features[bi])

        new_cluster_sizes = [len(c) for c in clusters]
        if new_cluster_sizes == cluster_sizes:
            break

    return clusters


def visualise_points(cluster, scale=0.50, color=(255, 120, 50)):
    serial = str(random.randint(0, 10000))
    spheres = []
    for f in cluster:
        p = f.coordinates
        spheres.append('7.0, %f, %f, %f, %f' % (p[0], p[1], p[2], scale))
    out = []
    out.append('points_' + serial + ' = [' + ', '.join(spheres) + ']')
    out.append('cmd.load_cgo(points_' + serial + ', \'points_' + serial + '\', 1)')
    out.append('cmd.set_color(\'color_' + serial + '\', [%i, %i, %i])' % color)
    out.append('cmd.color(\'color_' + serial + '\', \'points_' + serial + '\')')
    return os.linesep.join(out)


def main():
    args = parse_args()

    # Load a previously detected cavity:
    cav = Cavity.from_xml_file(args.input_file)

    min_subcavity_size = args.min_subcavity_size

    # Define the initial minimum depth value of the pseudocenters that should be used for
    # the subcavity detection. This is in the interval [0,7] where 7 means most buried.
    min_depth = args.min_depth

    # Define the number of clusters to be detected
    k_value = args.k_value

    # Perform a k-means clustering with a decreasing depth threshold until a cluster was
    # found that is large enough:
    clusters = []
    while min_depth >= 1.0:
        print('Running k-means clustering with min_depth %s' % min_depth)

        clusters = k_means_clustering(cav.features, min_depth, k_value)

        print('Cluster sizes: %s' % [len(c) for c in clusters])
        if max([len(c) for c in clusters]) >= min_subcavity_size:
            break
        min_depth -= 0.5
    print('Detected subcavity has size %s' % max([len(c) for c in clusters]))

    # Write a PyMol visualisation file for the detected clusters if an output
    # path is given
    if args.visualisation_directory is not None:
        out_file_name = os.path.join(args.visualisation_directory, 'cav_clusters.py')
        with open(out_file_name, 'w') as out_file:
            out_file.write('from pymol.cgo import *%s' % os.linesep)
            out_file.write('from pymol import cmd%s' % os.linesep)
            for i in range(len(clusters)):
                color = (80 + 80 * i, 50, 50)
                out_file.write(visualise_points(clusters[i], 0.50, color))
                out_file.write(os.linesep)

    # Find the index of the largest detected cluster
    max_cluster_index = max_cluster_size = -1
    for i in range(len(clusters)):
        if len(clusters[i]) > max_cluster_size:
            max_cluster_index = i
            max_cluster_size = len(clusters[i])

    # Create a subcavity from the largest cluster
    subcav = cav.subcavity(clusters[max_cluster_index])
    # Cavities to be compared with the deeply buried subcavity
    comp_cavs = glob.glob(os.path.join(args.xml_directory, '*.xml'))
    print('%s==========%s' % (os.linesep, os.linesep))
    print('Comparison\tScore')
    for cav_xml in comp_cavs:
        c = Cavity.from_xml_file(cav_xml)
        try:
            comp = subcav.compare(
                c, comparison_method=Cavity.ComparisonMethod.CAVITY_GRAPH_COMPARISON)
            score = comp.score
        except RuntimeError as exc:
            score = 'Failed: %s' % str(exc)
        print('\t'.join([c.identifier, str(score)]))

#####################################################################


if __name__ == '__main__':
    main()

Evaluate enrichment in cavity comparison

Screen a cavity against an entire database and obtain a list of all comparison cavities ordered by their degree of similarity. In the example script the cavity pdb3k2h.8, which hosts the ligand LYA (Pemetrexed), is used as query and all cavities that accommodate the same ligand are being defined as true positives. Also, the ROC AUC of the enrichment rates is calculated.

#!/usr/bin/env python
#
# This script can be used for any purpose without limitation subject to the
# conditions at http://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx
#
# This permission notice and the following statement of attribution must be
# included in all copies or substantial portions of this script.
#
# 2015-02-25: created by the Cambridge Crystallographic Data Centre
#
'''
    cavity_retrieval.py - Screen a cavity against an entire database and
    calculate the ROC AUC of the enrichment rates of cavities binding the same
    ligand.

    Usage: python cavity_retrieval.py
'''
#####################################################################

import argparse
from ccdc.cavity import *
import os

#####################################################################


def parse_args():
    usage_txt = ('Screen a cavity against an entire database and calculate',
                 ' the ROC AUC of the enrichment rates of cavities',
                 ' binding the same ligand.')
    parser = argparse.ArgumentParser(description=usage_txt,
                                     formatter_class=argparse.ArgumentDefaultsHelpFormatter)
    parser.add_argument('-q', '--query', default='pdb3k22.2',
                        help='The name of the cavity to use as a reference.')
    parser.add_argument('-l', '--ligand', default='LYA',
                        help='The identifier of the ligand within the cavity to define the pocket.')
    parser.add_argument('-d', '--database',
                        default=os.path.join(CavityDatabase.drugbank_database_dir(),
                                             'drugbank_cavities.sqlite'),
                        help='Database file to use a source for cavities.')
    return parser.parse_args()


def get_roc_coords(hits_positions, num_results):
    """ Calculates the coordinates of a ROC plot
        using the obtained positions of true positives and the total
        number of results in the list.
    """
    coords = []
    pos = 0
    for i in range(len(hits_positions)):
        while pos < hits_positions[i]:
            pos += 1
        x = float(pos - i) / (num_results - len(hits_positions))
        y = float(i) / len(hits_positions)
        coords.append((x, y))
    coords.append((1.0, 1.0))
    return coords


def get_auc(coords):
    """ Calculates the area under curve (AUC) for a given set of ROC coordinates """
    auc = 0.0
    for i in range(1, len(coords)):
        auc += (coords[i][0] - coords[i - 1][0]) * coords[i][1]
    return auc


def main():
    args = parse_args()

    # Load the cavity database and query cavity
    db = CavityDatabase(args.database)
    query_cav = db.get_cavity_by_name(args.query)

    if query_cav is None:
        raise RuntimeError(f'Cavity {args.query} not found in database {args.database}')

    # Collect all cavities from the database that bind the target ligand
    hits = db.get_cavity_identifiers_by_ligand(args.ligand)

    # Compare the query cavity against all entries of the database using
    # the cavity histograms comparison method:
    results = db.search(query_cav, Cavity.ComparisonMethod.CAVITY_HISTOGRAMS_COMPARISON)

    # Collect and print the positions of the true positives in the result list:
    print('Hit\tPos\tDistance')
    hits_positions = []

    for i in range(len(results) - 1):
        cav_id = results[i][1]

        # Don't output the comparison of the query cavity with itself
        if cav_id == args.query:
            pass

        # Output any other results
        elif cav_id in hits:
            hits_positions.append(i)
            print('\t'.join([cav_id, str(i), str(results[i][0])]))

    # Calculate the ROC coordinates and the corresponding AUC:
    roc = get_roc_coords(hits_positions, len(results))
    auc = get_auc(roc)
    print('%sAUC: %s' % (os.linesep, auc))

#####################################################################


if __name__ == '__main__':
    main()

Create a subcavity and use it for comparisons

This example shows how to generate a subcavity on the basis of a regular cavity and use it for a database screening afterwards. We are using a Thrombin structure (PDB code: 5af9) as the starting point and are interested in creating a new cavity that only comprises the S1 pocket of that enzyme. In order to identify the features that make up the subcavity, we employ the methoxy-benzamide fragment of the bound ligand, as we know that it points directly into the S1 subpocket. The newly created subcavity is then compared against all the other cavities in a database by making use the cavity graph comparison method. The cavity histograms comparison method is not recommended for subcavity comparisons.

#!/usr/bin/env python
#
# This script can be used for any purpose without limitation subject to the
# conditions at http://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx
#
# This permission notice and the following statement of attribution must be
# included in all copies or substantial portions of this script.
#
# 2016-12-13: created by the Cambridge Crystallographic Data Centre
#
'''
    cavity_subcavity.py - Create a subcavity based on a ligand fragment and
    compare it against a database.

    Usage: python cavity_subcavity.py
'''
#####################################################################

import argparse
from ccdc.cavity import *
import os
import sys

#####################################################################


def parse_args():
    usage_txt = ('Create a subcavity based on a ligand fragment and ',
                 'compare it against a database')
    parser = argparse.ArgumentParser(description=usage_txt,
                                     formatter_class=argparse.ArgumentDefaultsHelpFormatter)
    parser.add_argument('-p', '--pdb_file', default='5af9.pdb',
                        help='The pdb file to work from.')
    parser.add_argument('-l', '--ligand', default='SJR',
                        help='The identifier of the ligand to define the subcavity from.')
    parser.add_argument('-f', '--fragment_labels', nargs='+',
                        default=['C', 'O', 'C1', 'C12', 'C11', 'C2', 'C3', 'C4'],
                        help='The atoms defining the subcavity fragment within the ligand.')
    parser.add_argument('-r', '--subcavity_radius', type=float, default=6.0,
                        help='The radius around the fragment atoms defining the subcavity volume.')
    parser.add_argument('-s', '--max_hit_structures', type=int, default=50,
                        help='The maximum number of structures to compare against.')
    parser.add_argument('-m', '--max_results', type=int, default=10,
                        help='The maximum number of results to output.')
    parser.add_argument('-d', '--database',
                        default=os.path.join(CavityDatabase.drugbank_database_dir(),
                                             'drugbank_cavities.sqlite'),
                        help='Database file to use a source for cavities.')
    return parser.parse_args()


def main():
    args = parse_args()
    cavs = Cavity.from_pdb_file(args.pdb_file)
    pdb = Protein.from_file(args.pdb_file)
    db = CavityDatabase(args.database)

    # Get the ligand matching the given identifier
    ligand = None
    for lig in pdb.ligands:
        if args.ligand in lig.identifier:
            ligand = lig
            break

    # Get the cavity that hosts the ligand:
    query_cav = None
    for cav in cavs:
        for lig in cav.ligand_identifiers:
            if args.ligand in lig:
                query_cav = cav
                break

    if ligand is None:
        print('Could not find the ligand %s in the PDB structure.' % args.ligand)
        sys.exit(1)
    if query_cav is None:
        print('Could not find the cavity hosting the ligand %s.' % args.ligand)
        sys.exit(1)

    # Get the ligand atoms that make up the defined fragment
    atoms = [a for a in ligand.atoms if a.label in args.fragment_labels]

    # Create the subcavity as the volume within the given radius around the fragment atoms
    subcav = query_cav.subcavity(query_cav.features_by_atom_distance(atoms, args.subcavity_radius))

    # Apply settings that limits the maximum number of cavities found in the search
    settings = CavityDatabase.Settings()
    settings.max_hit_structures = args.max_hit_structures
    settings.verbose = True

    # Find matches for the subcavity using the fast cavity graph comparison method
    results = db.search(subcav, settings=settings,
                        comparison_method=Cavity.ComparisonMethod.FAST_CAVITY_GRAPH_COMPARISON)

    # Print the 10 most similar cavities
    print('%sCavity\tScore' % os.linesep)
    for result in results[:args.max_results]:
        print('%s\t%s' % (result[1], result[0]))

#####################################################################


if __name__ == '__main__':
    main()

Perform subcavity comparisons and filter out results from structures homologous to the query

This examples demonstrates the creation of a subcavity using a specified radius around a fragment of a ligand. This is then used in fast cavity graph comparisons against a cavity database and the results of this search are filtered on the basis of a FASTA sequence search. All cavities from proteins with a chain that has statistically significant sequence similarity to the chain associated with the query ligand, as judged by an expectation value below a certain threshold (default 0.001), are removed from the results. Cavity identifiers and scores are written out to a CSV file.

#!/usr/bin/env python
#
# This script can be used for any purpose without limitation subject to the
# conditions at http://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx
#
# This permission notice and the following statement of attribution must be
# included in all copies or substantial portions of this script.
#
# 2018-03-21: created by the Cambridge Crystallographic Data Centre
#
'''
    cavity_subcavity_compare_filter.py - generate a subcavity of a certain
    radius around a specified ligand and perform comparisons against a
    database, filtering out results for proteins within the same homologous
    family using FASTA and Biopython. Expectation (E()) values are used in filtering
    the results, which depend on the size of the database searched. For protein
    searches, sequences with E()-values < 0.001 for searches of a 10,000 entry database
    are almost always homologous. Here, expectation values are calculated for each chain
    separately, and if any of these fall below the specified expectation threshold, all
    cavities associated with this PDB will be excluded from the results. Expectation
    values printed in the output files are the lowest found for the corresponding PDB
    structure, rather than being particularly associated with that cavity. The printed
    sequence similarity is that associated with the chain that led to the lowest expectation
    value.
'''
#####################################################################

import os
import argparse
import requests
import tempfile
import shutil
import errno

from subprocess import call
from collections import Counter
from time import sleep

from Bio import SeqIO, PDB, SearchIO

from ccdc.utilities import _temporary_copy, output_file, to_utf8
from ccdc.protein import Protein
from ccdc.cavity import Cavity, CavityDatabase


class Runner(argparse.ArgumentParser):
    '''Define and parse arguments, and run the comparisons according to the arguments.'''

    def __init__(self):
        '''Defines and parses the arguments.'''
        super(self.__class__, self).__init__(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)

        cavity_dir = CavityDatabase.drugbank_database_dir()
        if cavity_dir is not None:
            default_db = os.path.join(cavity_dir, 'drugbank_cavities.sqlite')
            have_default_db = os.path.exists(default_db)
        else:
            have_default_db = False
        if have_default_db:
            self.add_argument(
                '-d', '--database', default=default_db,
                help='The cavity database to search [%s]' % default_db
            )
        else:
            self.add_argument(
                '-d', '--database', required=True,
                help='The cavity database to search. See create_cavity_database.py if you need to build one'
            )
        self.add_argument(
            '-q',
            '--query_pdb_id',
            help='PDB identifier for the query protein, e.g. 3a7h',
            required=True
        )
        self.add_argument(
            '-l',
            '--ligand_id',
            help='Ligand identifier with residue number for the cavity of interest, e.g. ATP400',
            required=True
        )
        self.add_argument(
            '-c',
            '--ligand_chain_id',
            help='Chain identifier for ligand in the cavity of interest, e.g. B',
            default=None
        )
        self.add_argument(
            '-a',
            '--ligand_atoms', nargs='+',
            help='Ligand atoms used in defining the binding site, '
                 'e.g. N3 C2 N1 C6 C5 N7 C8 N9 C4',
            default=None
        )
        self.add_argument(
            '-r',
            '--radius_around_ligand',
            help='Radius around ligand used in defining the binding site [6.0]',
            type=float,
            default=6.0
        )
        self.add_argument(
            '-m',
            '--max_hit_structures',
            help='Maximum number of structures returned in the database search',
            type=int,
            default=-1
        )
        group = self.add_mutually_exclusive_group()
        group.add_argument(
            '--pdb_mirror', nargs='?',
            help='PDB mirror file path template, e.g. /local/mirror/pub/pdb/data/structures/all/pdb/pdb{}s.ent.gz',
            default=None
        )
        group.add_argument(
            '--pdb_url', nargs='?',
            help='PDB URL file template',
            default='https://files.rcsb.org/download/{}.pdb'
        )
        if 'FASTA_SEQUENCE_SEARCH_TOOL' in os.environ:
            self.add_argument(
                '-f',
                '--fasta_location',
                help='Path to FASTA executable',
                default=os.environ['FASTA_SEQUENCE_SEARCH_TOOL']
            )
        else:
            self.add_argument(
                '-f',
                '--fasta_location',
                help='Path to FASTA executable',
                required=True
            )
        self.add_argument(
            '-e',
            '--expectation_threshold',
            help='Threshold for the FASTA expectation value, below which structures will be considered homologous [0.001]',
            type=float,
            default=0.001
        )
        self.add_argument(
            '-s',
            '--score-filter',
            help='Threshold below which comparison scores will not be written to the output file',
            type=float,
            default=0.0
        )
        self.add_argument(
            '-o',
            '--results_output_file',
            help='Name of output file for results',
            default='output.csv'
        )
        self.add_argument(
            '-x',
            '--excluded_results_output_file',
            help='Name of output file for excluded results',
            default='excluded_output.csv'
        )

        self.args = self.parse_args()
        self.temp_dir = tempfile.mkdtemp()

        self.query_sequence_file = 'query.fasta'
        open(self.query_sequence_file, 'w').close()
        self.targets_sequence_file = 'targets.fasta'
        open(self.targets_sequence_file, 'w').close()
        self.fasta_output_file = 'output.fasta'

        if not os.path.exists(self.args.database):
            raise ValueError('The specified database {} does not exist'.format(self.args.database))

        self.db = CavityDatabase(self.args.database)
        self.settings = CavityDatabase.Settings()
        self.query_pdb_file = 'pdb{}.pdb'.format(self.args.query_pdb_id)
        self.query_protein = self.load_query_protein()
        self.query_cavities = Cavity.from_pdb_file(self.query_pdb_file)
        self.ligand_chain_id = self.args.ligand_chain_id

    @staticmethod
    def safe_key_value(a_dict, key):
        '''Return a value corresponding to the specified header key.'''
        value = a_dict.get(key, '')
        return value.strip() if isinstance(value, str) else value

    def get_pdb_class(self, pdb_id, pdb_file):
        struct_parser = PDB.PDBParser(PERMISSIVE=False, get_header=True)  # strict parser, will assign altloc
        pdb_structure = struct_parser.get_structure(pdb_id, pdb_file)
        classification = self.safe_key_value(pdb_structure.header, 'head').capitalize()
        return classification

    def load_query_protein(self):
        self.fetch_pdb(self.args.query_pdb_id, self.query_pdb_file)
        if not os.path.exists(self.query_pdb_file):
            raise RuntimeError('Unable to fetch file {} for query PDB {}'.format(self.query_pdb_file, self.args.query_pdb_id))
        protein = Protein.from_file(self.query_pdb_file)
        return protein

    def fetch_from_pdb(self, pdb_id, filename):
        """Try to retrieve the given structure from the RCSB web database.

        :param pdb_id: (:obj:`str`) The PDB structure identifier to retrieve.
        :param filename: (:obj:`str`) The file name to write the structure to.

        :raises: requests.HTTPError if an error occurs trying to download the structure.
        """

        # Try to retrieve the structure from the online database.
        pdb_req = requests.get(self.args.pdb_url.format(pdb_id))
        # Raise a requests.HTTPError if any HTTP errors have occurred.
        pdb_req.raise_for_status()

        # Write the structure returned from the PDB to the given file name
        with output_file(filename) as pdb_file:
            pdb_file.write(to_utf8(pdb_req.text))

    def fetch_pdb_from_mirror(self, pdb_id, filename):
        full_path_input_file = self.args.pdb_mirror.format(pdb_id)
        with _temporary_copy(full_path_input_file) as pdb_file:
            shutil.copy2(pdb_file, filename)

    def fetch_pdb(self, pdb_id, filename):
        if self.args.pdb_mirror is not None:
            self.fetch_pdb_from_mirror(pdb_id, filename)
        else:
            for attempt in range(15):
                try:
                    self.fetch_from_pdb(pdb_id, filename)
                except:
                    sleep(1)
                else:
                    break
            else:
                print('Unable to fetch PDB {} from URL {}'.format(pdb_id, self.args.pdb_url.format(pdb_id)))

    def create_subcavity_around_ligand(self):
        '''Return the subcavity of a specified radius around the query ligand.'''
        # Get the ligand that contains the ligand identifier
        query_ligand = None
        for ligand in self.query_protein.ligands:
            chain_id_name = ligand.identifier.split(':')
            if self.args.ligand_id == chain_id_name[-1]:
                if len(chain_id_name) == 1 or self.ligand_chain_id is None or self.ligand_chain_id == chain_id_name[0]:
                    query_ligand = ligand
                    # Set chain ID for selected ligand if none was specified
                    if len(chain_id_name) == 2 and self.ligand_chain_id is None:
                        self.ligand_chain_id = chain_id_name[0]
                    break

        if query_ligand is None:
            raise IndexError('Ligand not found for ligand identifier {}'.format(self.args.ligand_id))
        print('Selected ligand ' + query_ligand.identifier)

        # Get the cavity that hosts the ligand
        query_cav = None
        for cav in self.query_cavities:
            for lig in cav.ligand_identifiers:
                chain_id_name = lig.split(':')
                if chain_id_name[-1] in self.args.ligand_id:
                    if len(chain_id_name) == 1 or self.ligand_chain_id == chain_id_name[0]:
                        query_cav = cav
                        break

        if query_cav is None:
            raise RuntimeError('Could not find the cavity hosting the ligand.')

        print('Selected cavity ' + query_cav.identifier + ' containing ligand ')
        [print(cav) for cav in query_cav.ligand_identifiers]

        # Get the ligand atoms that make up the desired fragment of the ligand
        # around which the subcavity will be defined
        if not self.args.ligand_atoms:
            fragment_labels = [str(a.label) for a in query_ligand.atoms]
        else:
            fragment_labels = self.args.ligand_atoms
        atoms = [a for a in query_ligand.atoms if a.label in fragment_labels]

        # Create the subcavity by using the pseudocenters within a specified radius of the fragment atoms
        subcavity = query_cav.subcavity(query_cav.features_by_atom_distance(atoms, self.args.radius_around_ligand))
        return subcavity

    @staticmethod
    def get_subcavity_chain(subcavity):
        '''Return the ID for the chain associated with the majority of subcavity pseudocenters.'''
        feature_chain_ids = []
        for feature in subcavity.features:
            feature_chain_ids.append(feature.chain.identifier)
        count = Counter(feature_chain_ids)
        chain_id = str(count.most_common(1)[0][0])
        return chain_id

    @staticmethod
    def write_sequence_file(pdb_file, sequence_file, chain_id=None):
        '''Write sequence file for FASTA.'''
        sequences = []
        records = list(SeqIO.parse(pdb_file, 'pdb-seqres'))
        for record in records:
            if not record.seq:
                continue
            if ':' in record.id:
                record_chain_id = record.id.split(':')[1]
            else:
                # Handle the case where Biopython cannot determine the PDB ID by getting this from the file name
                record_chain_id = record.id
                record.id = pdb_file[-8:-4].upper() + ':' + record_chain_id
            if chain_id is not None:
                if record_chain_id == chain_id:
                    sequences.append(record)
                    break
            else:
                sequences.append(record)

        if not sequences:
            raise IndexError('No chain sequences found.')

        with open(sequence_file, 'a') as output_handle:
            for sequence in sequences:
                SeqIO.write(sequence, output_handle, 'fasta')

    def find_expectation_values(self):
        '''Run FASTA and return a dictionary associating a PDB identifier with a tuple of the lowest
        expectation value found during FASTA sequence alignment and the associated sequence similarity.'''
        max_evalue = 0
        for record in SeqIO.parse("query.fasta", "fasta"):
            max_evalue += 1
        for record in SeqIO.parse("targets.fasta", "fasta"):
            max_evalue += 1

        args = ['-q', '-m', '10', '-E', str(max_evalue),
                self.query_sequence_file,
                self.targets_sequence_file]
        print('Running FASTA with arguments ' + str(args) + ' and writing the results to ' + self.fasta_output_file)
        with open(self.fasta_output_file, 'w') as out:
            program = self.args.fasta_location
            call([program] + args, stdout=out)

        pdb_dictionary = {}
        results = SearchIO.read('output.fasta', "fasta-m10")
        for target in results:
            target_pdb_id = target.id.split(':')[0]
            if target_pdb_id.lower() in pdb_dictionary:
                if target[0].evalue < pdb_dictionary[target_pdb_id.lower()][0]:
                    pdb_dictionary[target_pdb_id.lower()] = (target[0].evalue, target[0].pos_pct)
            else:
                pdb_dictionary[target_pdb_id.lower()] = (target[0].evalue, target[0].pos_pct)

        return pdb_dictionary

    def _extract_pdb_code(self, entry_id):
        pdb_code = entry_id[3:7] if entry_id[:3] == 'pdb' else entry_id[:4]
        return pdb_code

    def write_targets_sequence_file(self, results, temp_dir):
        '''Write FASTA sequence file for the PDB structures returned from the database search
        and return the classification of each PDB.'''
        header_classes = []
        target_pdbs = []
        for i in range(len(results)):
            target_pdb_id = self._extract_pdb_code(results[i][1])
            target_pdb_file = os.path.join(temp_dir, 'pdb{}.pdb'.format(target_pdb_id))
            self.fetch_pdb(target_pdb_id, target_pdb_file)
            if not os.path.exists(target_pdb_file):
                continue
            try:
                header_classes.append(self.get_pdb_class(target_pdb_id, target_pdb_file))
            except PDB.PDBExceptions.PDBConstructionException as e:
                print('Unable to determine PDB classification: {}'.format(e))
                header_classes.append('')
            if target_pdb_id.lower() not in target_pdbs:
                self.write_sequence_file(target_pdb_file, self.targets_sequence_file)
                target_pdbs.append(target_pdb_id.lower())
            os.remove(target_pdb_file)
        return header_classes

    def filter_homologous_structures(self, header_classes, results):
        '''Return two lists of results (non-homologous and homologous with query structure).'''
        pdb_dictionary = self.find_expectation_values()
        filtered_results = []
        homologous_results = []
        for result, header_class in zip(results, header_classes):
            expectation_value, similarity = pdb_dictionary[self._extract_pdb_code(result[1]).lower()]
            if expectation_value >= self.args.expectation_threshold:
                filtered_results.append((result[1], result[0], header_class, expectation_value, similarity))
            else:
                homologous_results.append((result[1], result[0], header_class, expectation_value, similarity))
        return filtered_results, homologous_results

    def write_results_file(self, subcavity, results, output_file):
        with open(output_file, 'w') as output:
            output.write('Identifier, Comparison score, Classification, No. pseudocenters, Lowest expectation, Associated similarity\n')
            # Results for subcavity self-self comparison
            output.write(', '.join((str(subcavity.identifier) + ' subcavity', str(subcavity.compare(subcavity).score),
                                    str(self.get_pdb_class(self.args.query_pdb_id, self.query_pdb_file)),
                                    str(len(subcavity.features)), '-1', '-1')) + '\n')
            for i in range(len(results)):
                if results[i][1] >= self.args.score_filter:
                    # Filtered database comparison results
                    output.write(', '.join((results[i][0], str(results[i][1]), results[i][2],
                                            str(len(self.db.get_cavity_by_name(results[i][0]).features)),
                                            str(results[i][3]), str(results[i][4]))) + '\n')

    def run(self):
        '''Run comparisons of the subcavity against the database and filter out structures homologous to the query.'''
        atoms = ''
        if self.args.ligand_atoms is not None:
            atoms = 'atoms {} of '.format(self.args.ligand_atoms)
        chain = ''
        if self.ligand_chain_id is not None:
            chain = self.ligand_chain_id + ':'
        print(u'Creating subcavity of radius {} Angstroms around {}'
              'the {}{} ligand for PDB structure {}'
              .format(self.args.radius_around_ligand,
                      atoms,
                      chain,
                      self.args.ligand_id,
                      self.query_protein.identifier).encode('utf-8'))

        subcavity = self.create_subcavity_around_ligand()
        print('The query subcavity has ' + str(len(subcavity.features)) + ' pseudocenters')

        chain_id = self.get_subcavity_chain(subcavity)
        self.write_sequence_file(self.query_pdb_file, self.query_sequence_file, chain_id)

        self.settings.max_hit_structures = self.args.max_hit_structures
        self.settings.verbose = True
        # Compare the subcavity against the database using the fast cavity graph comparison method
        results = self.db.search(subcavity,
                                 settings=self.settings,
                                 comparison_method=Cavity.ComparisonMethod.FAST_CAVITY_GRAPH_COMPARISON)

        # Make temporary directory for fetching PDB targets for writing sequence file
        try:
            temp_dir = tempfile.mkdtemp()
            header_classes = self.write_targets_sequence_file(results, temp_dir)
        finally:
            try:
                shutil.rmtree(temp_dir)
            except OSError as exc:
                if exc.errno != errno.ENOENT:
                    raise

        filtered_results, homologous_results = self.filter_homologous_structures(header_classes, results)

        # Write cavity identifiers, comparison scores and PDB classifications to a file (sorted by score)
        print('Writing filtered results to {}'.format(self.args.results_output_file))
        self.write_results_file(subcavity, filtered_results, self.args.results_output_file)
        print('Writing excluded results to {}'.format(self.args.excluded_results_output_file))
        self.write_results_file(subcavity, homologous_results, self.args.excluded_results_output_file)

########################################################################################

if __name__ == '__main__':
    Runner().run()

Perform a k-leave-one-out cross-validation

This example shows how to perform a leave-one-out cross-validation (LOOCV) on two classes of binding sites using a k-nearest-neighbor classifier. In the example script, the LOOCV is carried out for k = 1, 3, 5, 7, 9. The two classes of binding sites are Imatinib binders (also known as Gleecev/Glivec) and Pemetrexed binders (also known as Alimta). Both are cytostatics used in cancer treatment. However, Imatinib is a tyrosine kinase inhibitor while Pemetrexed inhibits the folate-dependent enzymes thymidylate synthase, dihydrofolate reductase and glycinamide ribonucleotide formyltransferase. The LOOCV is a method to evaluate how well two classes of binding sites can be separated by making use of the results of a cavity comparison method.

#!/usr/bin/env python
#
# This script can be used for any purpose without limitation subject to the
# conditions at http://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx
#
# This permission notice and the following statement of attribution must be
# included in all copies or substantial portions of this script.
#
# 2015-02-25: created by the Cambridge Crystallographic Data Centre
#
'''
    k-leave-one-out_cross-validation.py - Generate a matrix of distances by
    using two classes of binding sites and run a k leave-one-out cross-validation
    to determine classification success. First, an all-against-all similarity
    matrix is generated. This matrix is subsequently subjected to a
    leave-one-out cross-validation by using a k-nearest-neighbor method as
    classifier and the success rates are printed for five different values of
    k (1, 3, 5, 7 ,9).

    Usage: python k-leave-one-out_cross-validation.py
'''
#####################################################################

import sys
import argparse

from ccdc.cavity import *

#####################################################################

# Perform a k-leave-one-out cross-validation using a similarity matrix, a value
# for k, and an index that defines the boundary between the first and the
# second class of cavities.
def kloocv(scores, boundary, k):
    def get_lowest_positions(values, n, exclude):
        positions = []
        while len(positions) < n:
            lowest = sys.maxsize
            lowest_pos = -1
            for i in range(len(values)):
                if values[i] < lowest and i not in positions and i != exclude:
                    lowest = values[i]
                    lowest_pos = i
            positions.append(lowest_pos)
        return positions

    total = 0
    correct = 0
    for i in range(len(scores)):
        line = list(scores[i])
        lowests = get_lowest_positions(line, k, i)
        votes_for_class_b = sum(v > boundary for v in lowests)
        votes_for_class_a = (k - votes_for_class_b)
        total += 1
        if i > boundary:
            if votes_for_class_b > votes_for_class_a:
                correct += 1
        else:
            if votes_for_class_a > votes_for_class_b:
                correct += 1
    return (100.0 * correct / total)

def main(cavity_database):
    db = CavityDatabase(cavity_database)
    # Collect the cavities of the two classes to be used:
    sti_cavities = db.get_cavities_by_ligand("STI")  # cavities binding Imatinib (Gleevec/Glivec)
    lya_cavities = db.get_cavities_by_ligand("LYA")  # cavities binding Pemetrexed

    print("Found", len(sti_cavities), "cavities binding STI.")
    print("Found", len(lya_cavities), "cavities binding LYA.")

    cavities = sti_cavities + lya_cavities

    # Build scoring matrix
    scores = []
    for i in range(len(cavities)):
        scores.append([])
        for j in range(len(cavities)):
            scores[i].append(0.0)

    # Insert similarity scores
    print("Performing the comparisons")
    c = 0
    start = time.time()
    for i in range(len(cavities)):
        for j in range(i + 1, len(cavities)):
            d = cavities[i].compare(cavities[j], comparison_method=Cavity.ComparisonMethod.CAVITY_HISTOGRAMS_COMPARISON)
            scores[i][j] = d
            scores[j][i] = d
            c += 1
    elapsed = round(time.time() - start, 3)
    print("  -> ran", c, "comparisons in", elapsed, "sec.")

    # Run k leave-one-out cross-validation
    print("---")
    print("k leave-one-out cross-validation results:")
    boundary = len(sti_cavities) - 1
    for k in range(1, 11, 2):
        success = kloocv(scores, boundary, k)
        print("Success rate for k =", k, ":", round(success, 1), "%")

#####################################################################

if __name__ == '__main__':

    drugbank_cavity_database = \
        os.path.join(CavityDatabase.drugbank_database_dir(), 'drugbank_cavities.sqlite')
    parser = argparse.ArgumentParser(description='Cavity K Leave One Out Cross Validation.')
    parser.add_argument(
        '-d',
        '--cavity_database',
        help='Cavity sqlite database for testing.',
        default=drugbank_cavity_database)
    args = parser.parse_args()
    main(args.cavity_database)