Finding pockets and their properties (like, Surface Area) from a PDB file

This example shows how to find pockets in an input PDB file using the OEFindPockets function. OEDesignUnit file (OEDesignUnit) can also be used as input to find the pockets.

Command Line Interface

This example uses an input PDB file, and will print surface area and list of residues for each found pocket in given PDB file.

prompt> findpockets <input protein molecule or DesignUnit file>

Code

Download code

findpockets.py and the 7mhf.pdb (the input PDB file),

#!/usr/bin/env python
# (C) 2022 Cadence Design Systems, Inc. (Cadence) 
# All rights reserved.
# TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
# provided to current licensees or subscribers of Cadence products or
# SaaS offerings (each a "Customer").
# Customer is hereby permitted to use, copy, and modify the Sample Code,
# subject to these terms. Cadence claims no rights to Customer's
# modifications. Modification of Sample Code is at Customer's sole and
# exclusive risk. Sample Code may require Customer to have a then
# current license or subscription to the applicable Cadence offering.
# THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
# NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
# PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
# liable for any damages or liability in connection with the Sample Code
# or its use.

#############################################################################################################
# This program demonstrates how to find pockets and some of its properties from a protein or DesignUnit file.
#############################################################################################################
import sys
from openeye import oechem
from openeye import oespruce


def readProteinMol(ifilename):
    if oechem.OEGetFileExtension(ifilename) == 'oedu':
        du = oechem.OEDesignUnit()
        if not oechem.OEReadDesignUnit(ifilename, du):
            oechem.OEThrow.Fatal("Unable to open %s for reading OEDesignUnit" % ifilename)
        return du
    else:
        ifs = oechem.oemolistream()
        if not ifs.open(ifilename):
            oechem.OEThrow.Fatal("Unable to open %s for reading" % ifilename)
        mol = oechem.OEGraphMol()
        oechem.OEReadMolecule(ifs, mol)
        return mol


def main(argv=None):
    if argv is None:
        argv = [__name__]
    if len(sys.argv) != 2:
        oechem.OEThrow.Usage("%s <protein or DesignUnit input file>" % argv[0])

    mol = readProteinMol(sys.argv[1])

    pockets = oespruce.OEFindPockets(mol)
    print("pockets count: %s" % len(list(pockets)))
    pockets.ToFirst()
    pocket_cntr = 0
    for pocket in pockets:
        pocket_cntr += 1
        pocket_residues = pocket.GetResidues()
        print("pocket_%s Residues count: " % pocket_cntr, len(list(pocket_residues)))
        pocket_residues.ToFirst()
        print("pocket_%s Residues: " % pocket_cntr)
        for res in pocket_residues:
            print(res)
        print("pocket_%s Surface Area: " % pocket_cntr,  pocket.GetSurfaceArea())


if __name__ == "__main__":
    sys.exit(main(sys.argv))

Example

prompt> findpockets.py 7mhf.pdb

will generate the following output:

pockets count: 1
pocket_1 Residues count:  22
pocket_1 Residues:
HOH 623   A 2
GLY 143   A 1
ASN 142   A 1
HOH 502   A 2
DMS 404   A 2
MET 165   A 1 A
HOH 624   A 2
ASP 187   A 1
HOH 671   A 2
HOH 655   A 2
HOH 589   A 2
MET 49   A 1
HOH 782   A 2
ARG 188   A 1
THR 25   A 1
HIS 41   A 1
GLN 189   A 1
CYS 44   A 1
HOH 750   A 2
HOH 548   A 2
THR 45   A 1
SER 46   A 1
pocket_1 Surface Area:  211.0912322998047

See also

• OEPocketOptions class