OEBio Examples

Extracting the backbone of a protein

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2003-2015 OpenEye Scientific Software, Inc.
#############################################################################
# This program demonstrates how to extract the backbone of a protein.
#############################################################################
import os, sys
from openeye.oechem import *

def BackBone(ifs, ofs):
    adjustHCount = True
    mol = OEGraphMol()
    backboneMol = OEGraphMol()

    while OEReadMolecule(ifs, mol):
        if not OEHasResidues(mol):
            OEPerceiveResidues(mol, OEPreserveResInfo_All)
        aiter = mol.GetAtoms(OEIsBackboneAtom())
        member = OEIsAtomMember(aiter)

        OESubsetMol(backboneMol, mol, member, adjustHCount)
        OEWriteMolecule(ofs, backboneMol)

def main(argv=[__name__]):
    if len(sys.argv) != 3:
        OEThrow.Usage("%s <infile> <outfile>" % argv[0])

    ifs = oemolistream()
    if not ifs.open(argv[1]):
        OEThrow.Fatal("Unable to open %s for reading" % argv[1])

    ofs = oemolostream()
    if not ofs.open(argv[2]):
       OEThrow.Fatal("Unable to open %s for writing" % argv[2])

    BackBone(ifs, ofs)

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

See also

Splitting a Macro-molecular Complex

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2015 OpenEye Scientific Software, Inc.
#############################################################################
# split a mol complex (a PDB structure, for example) into basic categories
#############################################################################

import sys
from openeye.oechem import *

def main(argv=[__name__]):

    itf = OEInterface(InterfaceData)
    OEConfigureSplitMolComplexOptions(itf)

    if not OEParseCommandLine(itf, argv):
        OEThrow.Fatal("Unable to interpret command line!")

    if itf.GetBool("-verbose"):
        OEThrow.SetLevel(OEErrorLevel_Verbose)

    iname = itf.GetString("-in")
    oname = itf.GetString("-out")

    ims = oemolistream()
    if not itf.GetUnsignedInt("-modelnum") == 1:
        ims.SetFlavor(OEFormat_PDB,
                      OEGetDefaultIFlavor(OEFormat_PDB) & ~OEIFlavor_PDB_ENDM)
    if not ims.open(iname):
        OEThrow.Fatal("Cannot open input file!")

    oms = oemolostream()
    if not oms.open(oname):
        OEThrow.Fatal("Cannot open output file!")

    inmol = OEGraphMol()
    if not OEReadMolecule(ims, inmol):
        OEThrow.Fatal("Cannot read input file!")

    opts = OESplitMolComplexOptions()
    OESetupSplitMolComplexOptions(opts, itf)

    if itf.GetBool("-verbose"):
        # don't bother counting sites unless we're going to print them
        numSites = OECountMolComplexSites(inmol, opts)
        OEThrow.SetLevel(OEErrorLevel_Verbose)
        OEThrow.Verbose("sites %d" % numSites)

    lig   = OEGraphMol()
    prot  = OEGraphMol()
    wat   = OEGraphMol()
    other = OEGraphMol()

    if not OESplitMolComplex(lig, prot, wat, other, inmol, opts):
        OEThrow.Fatal("Unable to split mol complex from %s" % iname)

    if not lig.NumAtoms() == 0:
        OEThrow.Verbose("  lig %s" % lig.GetTitle())
        OEWriteMolecule(oms, lig)

    if not prot.NumAtoms() == 0:
        OEThrow.Verbose(" prot %s" % prot.GetTitle())
        OEWriteMolecule(oms, prot)

    if not wat.NumAtoms() == 0:
        OEThrow.Verbose("  wat %s" % wat.GetTitle())
        OEWriteMolecule(oms, wat)

    if not other.NumAtoms() == 0:
        OEThrow.Verbose("other %s" % other.GetTitle())
        OEWriteMolecule(oms, other)

    oms.close()

#############################################################################
# INTERFACE
#############################################################################

InterfaceData = '''
!BRIEF splitmolcomplex.py <inmol> [<outmol>]

!CATEGORY "input/output options :"
   !PARAMETER -in
      !ALIAS -i
      !TYPE string
      !BRIEF Input molecule (usually a pdb file)
      !VISIBILITY simple
      !REQUIRED true
      !KEYLESS 1
   !END

   !PARAMETER -out
      !ALIAS -o
      !TYPE string
      !DEFAULT splitmolcomplex.oeb.gz
      !BRIEF Output molecule (usually an oeb)
      !VISIBILITY simple
      !REQUIRED false
      !KEYLESS 2
   !END
!END

!CATEGORY "Display options :"
   !PARAMETER -verbose
      !ALIAS -v
      !TYPE bool
      !DEFAULT false
      !BRIEF If true, show molecule titles and number of binding sites
      !VISIBILITY simple
      !REQUIRED false
   !END
!END
'''

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

Splitting a Macro-molecular Complex Efficiently and Flexibly

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2016 OpenEye Scientific Software, Inc.
#############################################################################
# split a mol complex (a PDB structure, for example) using low-level funcs
#############################################################################

import sys
from openeye.oechem import *

def main(argv=[__name__]):

    itf = OEInterface(InterfaceData)
    OEConfigureSplitMolComplexOptions(itf)

    if not OEParseCommandLine(itf, argv):
        OEThrow.Fatal("Unable to interpret command line!")

    if itf.GetBool("-verbose"):
        OEThrow.SetLevel(OEErrorLevel_Verbose)

    iname = itf.GetString("-in")
    oname = itf.GetString("-out")

    ims = oemolistream()
    if not itf.GetUnsignedInt("-modelnum") == 1:
        ims.SetFlavor(OEFormat_PDB,
                      OEGetDefaultIFlavor(OEFormat_PDB) & ~OEIFlavor_PDB_ENDM)
    if not ims.open(iname):
        OEThrow.Fatal("Cannot open input file!")

    oms = oemolostream()
    if not oms.open(oname):
        OEThrow.Fatal("Cannot open output file!")

    inmol = OEGraphMol()
    if not OEReadMolecule(ims, inmol):
        OEThrow.Fatal("Cannot read input file!")

    opts = OESplitMolComplexOptions()
    OESetupSplitMolComplexOptions(opts, itf)

    frags = OEAtomBondSetVector()
    if not OEGetMolComplexFragments(frags, inmol, opts):
        OEThrow.Fatal("Unable to split mol complex from %s" % iname)

    numSites = OECountMolComplexSites(frags)

    if itf.GetBool("-verbose"):
        OEThrow.SetLevel(OEErrorLevel_Verbose)
        OEThrow.Verbose("sites %d" % numSites)

    lig   = OEGraphMol()
    prot  = OEGraphMol()
    wat   = OEGraphMol()
    other = OEGraphMol()

    if not OECombineMolComplexFragments(lig, frags, opts, opts.GetLigandFilter()):
        OEThrow.Fatal("Unable to split ligand from %s" % iname)

    if not OECombineMolComplexFragments(prot, frags, opts, opts.GetProteinFilter()):
        OEThrow.Fatal("Unable to split protein complex from %s" % iname)

    if not OECombineMolComplexFragments(wat, frags, opts, opts.GetWaterFilter()):
        OEThrow.Fatal("Unable to split waters from %s" % iname)

    if not OECombineMolComplexFragments(other, frags, opts, opts.GetOtherFilter()):
        OEThrow.Fatal("Unable to split other mols from %s" % iname)

    if not lig.NumAtoms() == 0:
        OEThrow.Verbose("  lig %s" % lig.GetTitle())
        OEWriteMolecule(oms, lig)

    if not prot.NumAtoms() == 0:
        OEThrow.Verbose(" prot %s" % prot.GetTitle())
        OEWriteMolecule(oms, prot)

    if not wat.NumAtoms() == 0:
        OEThrow.Verbose("  wat %s" % wat.GetTitle())
        OEWriteMolecule(oms, wat)

    if not other.NumAtoms() == 0:
        OEThrow.Verbose("other %s" % other.GetTitle())
        OEWriteMolecule(oms, other)

    oms.close()

#############################################################################
# INTERFACE
#############################################################################

InterfaceData = '''
!BRIEF splitmolcomplexlowlevel.py <inmol> [<outmol>]

!CATEGORY "input/output options :"
   !PARAMETER -in
      !ALIAS -i
      !TYPE string
      !BRIEF Input molecule (usually a pdb file)
      !VISIBILITY simple
      !REQUIRED true
      !KEYLESS 1
   !END

   !PARAMETER -out
      !ALIAS -o
      !TYPE string
      !DEFAULT splitmolcomplex.oeb.gz
      !BRIEF Output molecule (usually an oeb)
      !VISIBILITY simple
      !REQUIRED false
      !KEYLESS 2
   !END
!END

!CATEGORY "Display options :"
   !PARAMETER -verbose
      !ALIAS -v
      !TYPE bool
      !DEFAULT false
      !BRIEF If true, show molecule titles and number of binding sites
      !VISIBILITY simple
      !REQUIRED false
   !END
!END
'''

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

Splitting a Macro-molecular Complex Into Fragments

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2015 OpenEye Scientific Software, Inc.
#############################################################################
# split a mol complex (a PDB structure, for example) into fragments
#############################################################################

import sys
from openeye.oechem import *

def main(argv=[__name__]):

    itf = OEInterface(InterfaceData)
    OEConfigureSplitMolComplexOptions(itf)

    if not OEParseCommandLine(itf, argv):
        OEThrow.Fatal("Unable to interpret command line!")

    iname = itf.GetString("-in")
    oname = itf.GetString("-out")

    ims = oemolistream()
    if not itf.GetUnsignedInt("-modelnum") == 1:
        ims.SetFlavor(OEFormat_PDB,
                      OEGetDefaultIFlavor(OEFormat_PDB) & ~OEIFlavor_PDB_ENDM)
    if not ims.open(iname):
        OEThrow.Fatal("Cannot open input file!")

    oms = oemolostream()
    if not oms.open(oname):
        OEThrow.Fatal("Cannot open output file!")

    inmol = OEGraphMol()
    if not OEReadMolecule(ims, inmol):
        OEThrow.Fatal("Cannot read input file!")

    opts = OESplitMolComplexOptions()
    OESetupSplitMolComplexOptions(opts, itf)

    if itf.GetBool("-verbose"):
        # don't bother counting sites unless we're going to print them
        numSites = OECountMolComplexSites(inmol, opts)
        OEThrow.SetLevel(OEErrorLevel_Verbose)
        OEThrow.Verbose("sites %d" % numSites)

    for frag in OEGetMolComplexComponents(inmol, opts):
        OEThrow.Verbose("frag %s" % frag.GetTitle())
        OEWriteMolecule(oms, frag)

    oms.close()

#############################################################################
# INTERFACE
#############################################################################

InterfaceData = '''
!BRIEF splitmolcomplexfrags.py <inmol> [<outmol>]

!CATEGORY "input/output options :"
   !PARAMETER -in
      !ALIAS -i
      !TYPE string
      !BRIEF Input molecule (usually a pdb file)
      !VISIBILITY simple
      !REQUIRED true
      !KEYLESS 1
   !END

   !PARAMETER -out
      !ALIAS -o
      !TYPE string
      !DEFAULT splitmolcomplex.oeb.gz
      !BRIEF Output molecule (usually an oeb)
      !VISIBILITY simple
      !REQUIRED false
      !KEYLESS 2
   !END
!END

!CATEGORY "Display options :"
   !PARAMETER -verbose
      !ALIAS -v
      !TYPE bool
      !DEFAULT false
      !BRIEF If true, show fragment titles and number of binding sites
      !VISIBILITY simple
      !REQUIRED false
   !END
!END
'''

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

Preparing a Protein

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2016 OpenEye Scientific Software, Inc.
#############################################################################
# prepare a molecule: alts, hydrogens, split ligand
#############################################################################
import sys
from openeye.oechem import *

def WaterProcess(processName):
    if processName == "fullsearch":
        return OEPlaceHydrogensWaterProcessing_FullSearch
    elif processName == "focused":
        return OEPlaceHydrogensWaterProcessing_Focused
    return OEPlaceHydrogensWaterProcessing_Ignore

def main(argv=[__name__]):

    itf = OEInterface(InterfaceData)
    OEConfigureSplitMolComplexOptions(itf,
                                      OESplitMolComplexSetup_All              & \
                                  ~ ( OESplitMolComplexSetup_CovBondTreatment | \
                                      OESplitMolComplexSetup_CovCofactor ) )

    if not OEParseCommandLine(itf, argv):
        OEThrow.Fatal("Unable to interpret command line!")

    verbose = itf.GetBool("-verbose")
    if verbose:
        OEThrow.SetLevel(OEErrorLevel_Verbose)

    altProcess = itf.GetString("-alts")
    keepAlts     = (altProcess != "a")
    highestOcc   = (altProcess == "occupancy")
    compareAlts  = (altProcess == "compare")

    siteNum    = itf.GetUnsignedInt("-bindingsitenum")
    allSites   = (siteNum == 0)
    otherModel = (itf.GetUnsignedInt("-modelnum") != 1)

    placeHyd = itf.GetBool("-placehydrogens")
    splitlig = itf.HasString("-ligout")

    watProcessName = itf.GetString("-waterprocessing")
    waterProcess = WaterProcess(watProcessName)

    standardize = itf.GetBool("-standardizehyd")
    badclash  = itf.GetDouble("-clashcutoff")
    flipbias  = itf.GetDouble("-flipbias")
    maxStates = itf.GetDouble("-maxsubstates")

    flavor = OEIFlavor_PDB_Default | OEIFlavor_PDB_DATA
    if keepAlts:
        flavor = flavor | OEIFlavor_PDB_ALTLOC
    if otherModel:
        flavor = flavor & ~ OEIFlavor_PDB_ENDM

    ims = oemolistream()
    ims.SetFlavor(OEFormat_PDB, flavor)

    inputFile = itf.GetString("-in")
    if not ims.open(inputFile):
        OEThrow.Fatal("Unable to open %s for reading." % inputFile)

    if not OEIs3DFormat(ims.GetFormat()):
        OEThrow.Fatal("%s is not in a 3D format." % inputFile)

    inftype = OEGetFileType(OEGetFileExtension(inputFile))
    if (inftype == OEFormat_PDB) and not keepAlts:
        OEThrow.Verbose("Default processing of alt locations (keep just 'A' and ' ').")

    sopt = OESplitMolComplexOptions()
    OESetupSplitMolComplexOptions(sopt, itf)

    inmol = OEGraphMol()
    if not OEReadMolecule(ims, inmol):
        OEThrow.Fatal("Unable to read %s." % inputFile)

    ims.close()

    if (inmol.NumAtoms() == 0):
        OEThrow.Fatal("Input molecule %s contains no atoms." % inputFile)

    if inmol.GetTitle() == "":
      inmol.SetTitle("input mol")

    OEThrow.Verbose("Processing %s." % inmol.GetTitle())

    if not OEHasResidues(inmol):
        OEPerceiveResidues(inmol, OEPreserveResInfo_All)

    if highestOcc or compareAlts:
        alf = OEAltLocationFactory(inmol)
        if alf.GetGroupCount() != 0:
            if highestOcc:
                OEThrow.Verbose("Dropping alternate locations from input.")
                alf.MakePrimaryAltMol(inmol)
            elif compareAlts:
                OEThrow.Verbose("Fixing alternate location issues.")
                inmol = alf.GetSourceMol()

    outmol = OEGraphMol()
    if allSites:
        outmol = inmol
    else:
        OEThrow.Verbose("Splitting out selected complex.")

        soptSiteSel = OESplitMolComplexOptions(sopt)
        soptSiteSel.SetSplitCovalent(False) # do any cov lig splitting later

        frags = OEAtomBondSetVector()
        if not OEGetMolComplexFragments(frags, inmol, soptSiteSel):
            OEThrow.Fatal("Unable to fragment %s." % inmol.GetTitle())

        howManySites = OECountMolComplexSites(frags)
        if howManySites < siteNum:
              OEThrow.Warning(("Binding site count (%d) " + \
                               "less than requested site (%d) in %s.") % \
                              (howManySites, siteNum, inmol.GetTitle()))
              exit(0)

        if not OECombineMolComplexFragments(outmol, frags, soptSiteSel):
            OEThrow.Fatal("Unable to collect fragments from %s." % inmol.GetTitle())

        if (outmol.NumAtoms() == 0):
            OEThrow.Fatal("No fragments selected from %s." % inmol.GetTitle())

    if placeHyd:
        OEThrow.Verbose("Adding hydrogens to complex.")

        hopt = OEPlaceHydrogensOptions()
        hopt.SetAltsMustBeCompatible(compareAlts)
        hopt.SetStandardizeBondLen(standardize)
        hopt.SetWaterProcessing(waterProcess)
        hopt.SetBadClashOverlapDistance(badclash)
        hopt.SetFlipBiasScale(flipbias)
        hopt.SetMaxSubstateCutoff(maxStates)

        if verbose:
            details = OEPlaceHydrogensDetails()
            if not OEPlaceHydrogens(outmol, details, hopt):
                OEThrow.Fatal("Unable to place hydrogens and get details on %s." % inmol.GetTitle())
            OEThrow.Verbose(details.Describe())
        else:
            if not OEPlaceHydrogens(outmol, hopt):
                OEThrow.Fatal("Unable to place hydrogens on %s." % inmol.GetTitle())

    oms1 = oemolostream()
    cplxFile = itf.GetString("-cplxout")
    if not oms1.open(cplxFile):
        OEThrow.Fatal("Unable to open %s for writing." % cplxFile)

    if splitlig:
        OEThrow.Verbose("Splitting ligand from complex.")

        frags = OEAtomBondSetVector()
        if not OEGetMolComplexFragments(frags, outmol, sopt):
            OEThrow.Fatal("Unable to fragment complex from %s." % inmol.GetTitle())

        lfilter = sopt.GetLigandFilter()

        protComplex = OEGraphMol()
        if not OECombineMolComplexFragments(protComplex,
                                            frags,
                                            sopt,
                                            OENotRoleSet(lfilter)):
            OEThrow.Fatal("Unable to collect complex from %s." % inmol.GetTitle())

        if (protComplex.NumAtoms() == 0):
            OEThrow.Warning("No complex identified in %s." % inmol.GetTitle())
        else:
            OEWriteMolecule(oms1, protComplex)

        lig = OEGraphMol()
        if not OECombineMolComplexFragments(lig, frags, sopt, lfilter):
            OEThrow.Fatal("Unable to collect ligand from %s." % inmol.GetTitle())

        if (lig.NumAtoms() == 0):
            OEThrow.Warning("No ligand identified in %s." % inmol.GetTitle())
        else:
            oms2 = oemolostream()
            if splitlig:
                ligFile = itf.GetString("-ligout")
                if not oms2.open(ligFile):
                    OEThrow.Fatal("Unable to open %s for writing." % ligFile)

            OEThrow.Verbose("Ligand: %s" % lig.GetTitle())
            OEWriteMolecule(oms2, lig)
            oms2.close()
    else:
        OEWriteMolecule(oms1, outmol)

    oms1.close()

#############################################################################
# INTERFACE
#############################################################################

InterfaceData = '''
!BRIEF proteinprep.py [-options] <inmol> [<outcplx> [<outlig>]]

!CATEGORY "input/output options :" 1
   !PARAMETER -in 1
      !ALIAS -i
      !TYPE string
      !BRIEF Input molecule filename (must have 3D coordinates)
      !SIMPLE true
      !REQUIRED true
      !KEYLESS 1
   !END

   !PARAMETER -cplxout 2
      !ALIAS -p
      !TYPE string
      !DEFAULT proteinprep.oeb.gz
      !BRIEF Output complex filename
      !SIMPLE true
      !REQUIRED false
      !KEYLESS 2
   !END

   !PARAMETER -ligout 3
      !ALIAS -l
      !TYPE string
      !BRIEF Output ligand filename
      !SIMPLE true
      !REQUIRED false
      !KEYLESS 3
   !END
!END

!CATEGORY "Calculation options :" 2
    !PARAMETER -alts 1
       !TYPE string
       !LEGAL_VALUE occupancy
       !LEGAL_VALUE a
       !LEGAL_VALUE ignore
       !LEGAL_VALUE compare
       !DEFAULT occupancy
       !BRIEF Alternate location atom handling (affects atom:atom interactions)
       !SIMPLE true
       !REQUIRED false
       !DETAIL
         occupancy - keep just the highest average occupancy for each alt group
         a - keep only loc code A (and blank)
         ignore - assume alts already selected appropriately
         compare - keep all alts but only interact if same loc code (or blank)
    !END

    !PARAMETER -placehydrogens 2
       !TYPE bool
       !DEFAULT true
       !BRIEF If false, hydrogens will not be added
       !SIMPLE true
       !REQUIRED false
    !END

    !PARAMETER -waterprocessing 3
       !TYPE string
       !LEGAL_VALUE ignore
       !LEGAL_VALUE focused
       !LEGAL_VALUE fullsearch
       !DEFAULT fullsearch
       !BRIEF How waters are processed
       !SIMPLE true
       !REQUIRED false
       !DETAIL
         ignore - leave water hydrogens in a random orientation
         focused - search orientations based on neighboring polar groups
         fullsearch - do an extensive search of water orientations
    !END

    !PARAMETER -standardizehyd 4
       !ALIAS -stdhyd
       !TYPE bool
       !DEFAULT true
       !BRIEF If false, bonds for hydrogens are not adjusted to standard lengths
       !SIMPLE false
       !REQUIRED false
    !END

    !PARAMETER -clashcutoff 5
        !TYPE double
        !DEFAULT 0.4
        !BRIEF Van der Waals overlap (in Angstroms) defined to be a bad clash
        !SIMPLE false
        !REQUIRED false
    !END

    !PARAMETER -flipbias 6
        !TYPE double
        !DEFAULT 1.0
        !BRIEF Scale factor for the bias against flipping sidechains such as HIS
        !SIMPLE false
        !REQUIRED false
    !END

    !PARAMETER -maxsubstates 7
        !TYPE double
        !DEFAULT 1.0e8
        !BRIEF Maximum number of substates in a single step of hydrogen placement optimization
        !SIMPLE false
        !REQUIRED false
    !END
!END

!CATEGORY "Display options :" 9
   !PARAMETER -verbose 1
      !ALIAS -v
      !TYPE bool
      !DEFAULT false
      !BRIEF Display more information about the process
      !SIMPLE true
      !REQUIRED false
   !END
!END
'''

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

Perceive and Print Protein-Ligand Interactions

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2015-2016 OpenEye Scientific Software, Inc.
#############################################################################
# Print protein-ligand interactions
#############################################################################
from __future__ import print_function
import sys

from openeye.oechem import *


def main(argv=[__name__]):

    itf = OEInterface()
    OEConfigure(itf, InterfaceData)
    OEConfigureSplitMolComplexOptions(itf, OESplitMolComplexSetup_LigName)

    if not OEParseCommandLine(itf, argv):
        return 1

    iname = itf.GetString("-complex")

    ifs = oemolistream()
    if not ifs.open(iname):
        OEThrow.Fatal("Unable to open %s for reading" % iname)

    complexmol = OEGraphMol()
    if not OEReadMolecule(ifs, complexmol):
        OEThrow.Fatal("Unable to read molecule from %s" % iname)

    if not OEHasResidues(complexmol):
        OEPerceiveResidues(complexmol, OEPreserveResInfo_All)

    # Separate ligand and protein

    sopts = OESplitMolComplexOptions()
    OESetupSplitMolComplexOptions(sopts, itf)

    ligand = OEGraphMol()
    protein = OEGraphMol()
    water = OEGraphMol()
    other = OEGraphMol()

    pfilter = sopts.GetProteinFilter()
    wfilter = sopts.GetWaterFilter()
    sopts.SetProteinFilter(OEOrRoleSet(pfilter, wfilter))
    sopts.SetWaterFilter(OEMolComplexFilterFactory(OEMolComplexFilterCategory_Nothing))

    OESplitMolComplex(ligand, protein, water, other, complexmol, sopts)

    if ligand.NumAtoms() == 0:
        OEThrow.Fatal("Cannot separate complex!")

    # Perceive interactions

    asite = OEInteractionHintContainer(protein, ligand)
    if not OEIsValidActiveSite(asite):
        OEThrow.Fatal("Cannot initialize active site!")

    OEPerceiveInteractionHints(asite)

    print ("Number of interactions:", asite.NumInteractions())
    for itype in OEGetActiveSiteInteractionHintTypes():
        numinters = asite.NumInteractions(OEHasInteractionHintType(itype))
        if numinters == 0:
            continue
        print("%d %s :" % (numinters, itype.GetName()))

        inters = [s for s in  asite.GetInteractions(OEHasInteractionHintType(itype))]
        print ("\n".join(sorted(GetInteractionString(s) for s in inters)))


    print ("\nResidue interactions:")
    for res in OEGetResidues(asite.GetMolecule(OEProteinInteractionHintComponent())):
        PrintResidueInteractions(asite, res)

    print ("\nLigand atom interactions:")
    for atom in asite.GetMolecule(OELigandInteractionHintComponent()).GetAtoms():
        PrintLigandAtomInteractions(asite, atom)


def GetResidueName(residue):
    return "%3s %4d %s" % (residue.GetName(), residue.GetResidueNumber(), residue.GetChainID())


def GetInteractionString(inter):

    fragstrs = []
    for frag in [inter.GetBgnFragment(), inter.GetEndFragment()]:
        if frag is None:
            continue
        fragtype = frag.GetComponentType()
        if fragtype == OELigandInteractionHintComponent():
            fragstrs.append("ligand:" + " ".join(sorted(str(a) for a in frag.GetAtoms())))
        if fragtype == OEProteinInteractionHintComponent():
            fragstrs.append("protein:" + " ".join(sorted(str(a) for a in frag.GetAtoms())))

    return  " ".join(sorted(f for f in fragstrs))


def PrintResidueInteractions(asite, residue):

    ligatomnames = set()
    for inter in asite.GetInteractions(OEHasResidueInteractionHint(residue)):
        ligfrag = inter.GetFragment(OELigandInteractionHintComponent())
        if ligfrag is None:
            continue
        for latom in ligfrag.GetAtoms():
            ligatomnames.add(str(latom))

    if len(ligatomnames) != 0:
        print(GetResidueName(residue), ": ", " ".join(sorted(a for a in ligatomnames)))


def PrintLigandAtomInteractions(asite, atom):

    resnames = set()
    for inter in asite.GetInteractions(OEHasInteractionHint(atom)):
        profrag = inter.GetFragment(OEProteinInteractionHintComponent())
        if profrag is None:
            continue
        for patom in profrag.GetAtoms():
            residue = OEAtomGetResidue(patom)
            resnames.add(GetResidueName(residue))

    if len(resnames) != 0:
        print(atom, ":", " ".join(sorted(r for r in resnames)))


InterfaceData = '''
!BRIEF printinteractions [-complex] <input>

!CATEGORY "input/output options :"

  !PARAMETER -complex
    !ALIAS -c
    !TYPE string
    !KEYLESS 1
    !REQUIRED true
    !VISIBILITY simple
    !BRIEF Input filename of the protein complex
  !END

!END
'''

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