Appendix: Additional Examples in Python

These are full listings of programming examples that are offered for download, elsewhere in this chapter. For an guide to programming examples, see the list elsewhere in this chapter.

Listing 1: Generate Szmap points and probe orientations at ligand atoms.

#!/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.

# generate szmap points and probe orientations at ligand atoms

import sys
from openeye import oechem
from openeye import oeszmap

InterfaceData = """
!BRIEF SzmapBestOrientations.py [-prob #.#] [-du] <designunitfile> [-o] <molfile>
!PARAMETER -du
  !TYPE string
  !BRIEF Input designunit
  !REQUIRED true
  !KEYLESS 1
!END
!PARAMETER -o
  !TYPE string
  !BRIEF Output file for points and probes molecules
  !REQUIRED true
  !KEYLESS 2
!END
!PARAMETER -prob
  !TYPE double
  !DEFAULT 0.5
  !BRIEF Cutoff for cumulative probability of probes
  !REQUIRED false
!END
"""


def GenerateSzmapProbes(oms, cumulativeProb, lig, prot):
    """
    generate multiconf probes and data-rich points at ligand coords
    @rtype : None
    @param oms: output mol stream for points and probes
    @param cumulativeProb: cumulative probability for cutoff of point set
    @param lig: mol defining coordinates for szmap calcs
    @param prot: context mol for szmap calcs (must have charges and radii)
    """

    sz = oeszmap.OESzmapEngine(prot)

    rslt = oeszmap.OESzmapResults()

    points = oechem.OEGraphMol()
    points.SetTitle("points %s" % lig.GetTitle())
    probes = oechem.OEMol()

    coord = oechem.OEFloatArray(3)

    for i, atom in enumerate(lig.GetAtoms()):
        lig.GetCoords(atom, coord)

        if not oeszmap.OEIsClashing(sz, coord):
            oeszmap.OECalcSzmapResults(rslt, sz, coord)

            rslt.PlaceNewAtom(points)

            clear = False
            rslt.PlaceProbeSet(probes, cumulativeProb, clear)

            name = oeszmap.OEGetEnsembleName(oeszmap.OEEnsemble_NeutralDiffDeltaG)
            nddG = rslt.GetEnsembleValue(oeszmap.OEEnsemble_NeutralDiffDeltaG)
            print("%2d (%7.3f, %7.3f, %7.3f): %s = %.3f"
                  % (i, coord[0], coord[1], coord[2], name, nddG))
        else:
            print("%2d (%7.3f, %7.3f, %7.3f): CLASH"
                  % (i, coord[0], coord[1], coord[2]))

    oechem.OEWriteMolecule(oms, points)
    oechem.OEWriteMolecule(oms, probes)


def main(argv=(__name__)):
    """
    the protein should have charges and radii but the ligand need not
    """
    itf = oechem.OEInterface()

    if not oechem.OEConfigure(itf, InterfaceData):
        oechem.OEThrow.Fatal("Problem configuring OEInterface!")

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

    duFile = itf.GetString("-du")
    du = oechem.OEDesignUnit()
    if not oechem.OEReadDesignUnit(duFile, du):
        oechem.OEThrow.Fatal("Unable to open %s for reading" % duFile)

    if not du.HasLigand():
        oechem.OEThrow.Fatal("Input designunit %s does not have a ligand bound" % du.GetTitle())

    lig = oechem.OEGraphMol()
    du.GetLigand(lig)
    prot = oechem.OEGraphMol()
    du.GetComponents(prot, oechem.OEDesignUnitComponents_Protein)

    outputFile = itf.GetString("-o")
    if not oechem.OEIsWriteable(outputFile):
        oechem.OEThrow.Fatal("Invalid file extension for output %s" % outputFile)

    oms = oechem.oemolostream()
    if not oms.open(outputFile):
        oechem.OEThrow.Fatal("Unable to open %s for writing" % outputFile)

    GenerateSzmapProbes(oms, itf.GetDouble("-prob"), lig, prot)


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

Listing 2: Example of Szmap calculations.

#!/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.

import sys
from openeye import oechem
from openeye import oeszmap

InterfaceData = """
!BRIEF SzmapCalcs.py [-prob #.#] [-p] <molfile> [-l] <molfile> [-o] <molfile>
!PARAMETER -p
  !TYPE string
  !BRIEF Input protein (or other) context mol
  !REQUIRED true
  !KEYLESS 1
!END
!PARAMETER -l
  !TYPE string
  !BRIEF Input ligand coordinates for calculations
  !REQUIRED true
  !KEYLESS 2
!END
!PARAMETER -o
  !TYPE string
  !BRIEF Output file for points and probes molecules
  !REQUIRED true
  !KEYLESS 3
!END
!PARAMETER -prob
  !TYPE double
  !DEFAULT 0.5
  !BRIEF Cutoff for cumulative probability of probes
  !REQUIRED false
!END
"""


def RunSzmap(lig, prot):
    """
    run szmap at ligand coordinates in the protein context
    @rtype : None
    @param lig: mol defining coordinates for szmap calcs
    @param prot: context mol for szmap calcs (must have charges and radii)
    """

    opt = oeszmap.OESzmapEngineOptions()

    mol = oechem.OEGraphMol()
    opt.GetProbeMol(mol)
    print("probe mol smiles: %s" % oechem.OEMolToSmiles(mol))

    opt.SetProbe(360)
    print("norient = %d" % opt.NumOrientations())

    opt.SetProbe(24).SetMaskCutoff(999.99)
    print("norient = %d" % opt.NumOrientations())
    print("cutoff  = %.3f" % opt.GetMaskCutoff())

    opt.SetMaskCutoff(0.0)  # reset cutoff
    print("cutoff  = %.3f" % opt.GetMaskCutoff())

    p = opt.GetProbe()
    print("probe nconf = %d" % p.NumConfs())

    sz = oeszmap.OESzmapEngine(prot, opt)

    opt = sz.GetOptions()
    print("norient= %d" % opt.NumOrientations())
    print("name   = %s" % opt.GetProbeName())
    print("cutoff = %.3f" % opt.GetMaskCutoff())

    mol = sz.GetContext()
    print("context mol: %s" % mol.GetTitle())

    print("probe smiles: %s" % oechem.OEMolToSmiles(opt.GetProbe()))

    coord = oechem.OEFloatArray(3)

    rslt = oeszmap.OESzmapResults()

    for atom in lig.GetAtoms():
        lig.GetCoords(atom, coord)

        name = oeszmap.OEGetEnsembleName(oeszmap.OEEnsemble_NeutralDiffDeltaG)
        print("ensemble name: %s" % name)

        longName = True
        name = oeszmap.OEGetEnsembleName(oeszmap.OEEnsemble_NeutralDiffDeltaG, longName)
        nddg = oeszmap.OECalcSzmapValue(sz, coord)
        print("ensemble name: %s, value: %.3f" % (name, nddg))

        name = oeszmap.OEGetComponentName(oeszmap.OEComponent_Interaction)
        print("component name: %s" % name)

        if not oeszmap.OEIsClashing(sz, coord):
            oeszmap.OECalcSzmapResults(rslt, sz, coord)
            nddg = rslt.GetEnsembleValue(oeszmap.OEEnsemble_NeutralDiffDeltaG)

        if not oeszmap.OEIsClashing(sz, coord):
            oeszmap.OECalcSzmapResults(rslt, sz, coord)
            # ...

            print("result norient = %d" % rslt.NumOrientations())

            nddg = rslt.GetEnsembleValue(oeszmap.OEEnsemble_NeutralDiffDeltaG)
            print("nddG = %.3f" % nddg)
            print("vdw  = %.3f" % rslt.GetEnsembleValue(oeszmap.OEEnsemble_VDW))

            coulomb = rslt.GetComponent(oeszmap.OEComponent_Interaction)
            print("interaction:")
            print(" ".join("%.3f" % c for c in coulomb))

            order = rslt.GetProbabilityOrder()
            print("conf with greatest prob = %d" % order[0])

            prob = rslt.GetProbabilities()
            print("greatest prob = %.3f" % prob[order[0]])

            vdw = oechem.OEDoubleArray(rslt.NumOrientations())
            rslt.GetComponent(vdw, oeszmap.OEComponent_VDW)
            print("vdw greatest prob = %.3f" % vdw[order[0]])

            point = oechem.OEFloatArray(3)
            rslt.GetCoords(point)
            print("calc point: (%.3f, %.3f, %.3f)" % (point[0], point[1], point[2]))

            mcmol = oechem.OEMol()
            rslt.PlaceProbeSet(mcmol)

            probCutoff = 0.5
            rslt.PlaceProbeSet(mcmol, probCutoff)
            print("nconf to yield 50pct = %d" % mcmol.NumConfs())

            clear = False
            cumulativeProb = rslt.PlaceProbeSet(mcmol, 10, clear)
            print("best 10 cumulative prob = %.3f" % cumulativeProb)

            amol = oechem.OEGraphMol()
            atom = rslt.PlaceNewAtom(amol)

            print("vdw = %s" % atom.GetData("vdw"))

            pmol = oechem.OEGraphMol()
            rslt.PlaceProbeMol(pmol, order[0])


def GetSzmapEnergies(lig, prot):
    """
    run szmap at ligand coordinates in the protein context
    @rtype : None
    @param lig: mol defining coordinates for szmap calcs
    @param prot: context mol for szmap calcs (must have charges and radii)
    """

    print("num\tatom\t%s\t%s\t%s\t%s\t%s"
          % (oeszmap.OEGetEnsembleName(oeszmap.OEEnsemble_NeutralDiffDeltaG),
             oeszmap.OEGetEnsembleName(oeszmap.OEEnsemble_PSolv),
             oeszmap.OEGetEnsembleName(oeszmap.OEEnsemble_WSolv),
             oeszmap.OEGetEnsembleName(oeszmap.OEEnsemble_VDW),
             oeszmap.OEGetEnsembleName(oeszmap.OEEnsemble_OrderParam)))

    coord = oechem.OEFloatArray(3)

    sz = oeszmap.OESzmapEngine(prot)
    rslt = oeszmap.OESzmapResults()

    for i, atom in enumerate(lig.GetAtoms()):
        lig.GetCoords(atom, coord)

        if not oeszmap.OEIsClashing(sz, coord):
            oeszmap.OECalcSzmapResults(rslt, sz, coord)

            print("%2d\t%s\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f"
                  % (i, atom.GetName(),
                     rslt.GetEnsembleValue(oeszmap.OEEnsemble_NeutralDiffDeltaG),
                     rslt.GetEnsembleValue(oeszmap.OEEnsemble_PSolv),
                     rslt.GetEnsembleValue(oeszmap.OEEnsemble_WSolv),
                     rslt.GetEnsembleValue(oeszmap.OEEnsemble_VDW),
                     rslt.GetEnsembleValue(oeszmap.OEEnsemble_OrderParam)))
        else:
            print("%2d\t%s CLASH" % (i, atom.GetName()))


def GenerateSzmapProbes(oms, cumulativeProb, lig, prot):
    """
    generate multiconf probes and data-rich points at ligand coords
    @rtype : None
    @param oms: output mol stream for points and probes
    @param cumulativeProb: cumulative probability for cutoff of point set
    @param lig: mol defining coordinates for szmap calcs
    @param prot: context mol for szmap calcs (must have charges and radii)
    """

    coord = oechem.OEFloatArray(3)

    sz = oeszmap.OESzmapEngine(prot)
    rslt = oeszmap.OESzmapResults()

    points = oechem.OEGraphMol()
    points.SetTitle("points %s" % lig.GetTitle())

    probes = oechem.OEMol()

    for i, atom in enumerate(lig.GetAtoms()):
        lig.GetCoords(atom, coord)

        if not oeszmap.OEIsClashing(sz, coord):
            oeszmap.OECalcSzmapResults(rslt, sz, coord)

            rslt.PlaceNewAtom(points)

            clear = False
            rslt.PlaceProbeSet(probes, cumulativeProb, clear)

    oechem.OEWriteMolecule(oms, points)
    oechem.OEWriteMolecule(oms, probes)


def main(argv=(__name__)):
    """
    the protein should have charges and radii but the ligand need not
    """
    itf = oechem.OEInterface()

    if not oechem.OEConfigure(itf, InterfaceData):
        oechem.OEThrow.Fatal("Problem configuring OEInterface!")

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

    ligFile = itf.GetString("-l")
    ims = oechem.oemolistream()
    if not ims.open(ligFile):
        oechem.OEThrow.Fatal("Unable to open %s for reading" % ligFile)
    lig = oechem.OEGraphMol()
    oechem.OEReadMolecule(ims, lig)

    contextFile = itf.GetString("-p")
    ims = oechem.oemolistream()
    if not ims.open(contextFile):
        oechem.OEThrow.Fatal("Unable to open %s for reading" % contextFile)
    prot = oechem.OEGraphMol()
    oechem.OEReadMolecule(ims, prot)

    RunSzmap(lig, prot)

    GetSzmapEnergies(lig, prot)

    outputFile = itf.GetString("-o")
    oms = oechem.oemolostream()
    if not oms.open(outputFile):
        oechem.OEThrow.Fatal("Unable to open %s for writing" % outputFile)

    GenerateSzmapProbes(oms, itf.GetDouble("-prob"), lig, prot)


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

Listing 3: Analyze binding site energies.

#!/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.

# analyze binding site energies with szmaptk

import sys
from openeye import oechem
from openeye import oeszmap

InterfaceData = """
!BRIEF SzmapEnergies.py [-high_res] [-du] <designunitfile>
!PARAMETER -du
  !TYPE string
  !BRIEF Input designunit
  !REQUIRED true
  !KEYLESS 1
!END
!PARAMETER -high_res
  !TYPE bool
  !DEFAULT false
  !BRIEF If true, increase the number of rotations to 360
  !REQUIRED false
!END
"""


def GetSzmapEnergies(lig, prot, highRes):
    """
    run szmap at ligand coordinates in the protein context
    @rtype : None
    @param lig: mol defining coordinates for szmap calcs
    @param prot: context mol for szmap calcs (must have charges and radii)
    @param highRes: if true, use 360 rotations rather than the default of 60
    """

    opt = oeszmap.OESzmapEngineOptions()
    if highRes:
        opt.SetProbe(360)

    sz = oeszmap.OESzmapEngine(prot, opt)
    coord = oechem.OEFloatArray(3)
    rslt = oeszmap.OESzmapResults()

    print("num\tatom\t%s\t%s\t%s\t%s\t%s"
          % (oeszmap.OEGetEnsembleName(oeszmap.OEEnsemble_NeutralDiffDeltaG),
             oeszmap.OEGetEnsembleName(oeszmap.OEEnsemble_PSolv),
             oeszmap.OEGetEnsembleName(oeszmap.OEEnsemble_WSolv),
             oeszmap.OEGetEnsembleName(oeszmap.OEEnsemble_VDW),
             oeszmap.OEGetEnsembleName(oeszmap.OEEnsemble_OrderParam)))

    for i, atom in enumerate(lig.GetAtoms()):
        lig.GetCoords(atom, coord)

        if not oeszmap.OEIsClashing(sz, coord):
            oeszmap.OECalcSzmapResults(rslt, sz, coord)

            print("%2d\t%s\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f"
                  % (i, atom.GetName(),
                     rslt.GetEnsembleValue(oeszmap.OEEnsemble_NeutralDiffDeltaG),
                     rslt.GetEnsembleValue(oeszmap.OEEnsemble_PSolv),
                     rslt.GetEnsembleValue(oeszmap.OEEnsemble_WSolv),
                     rslt.GetEnsembleValue(oeszmap.OEEnsemble_VDW),
                     rslt.GetEnsembleValue(oeszmap.OEEnsemble_OrderParam)))
        else:
            print("%2d\t%s CLASH" % (i, atom.GetName()))


def main(argv=(__name__)):
    """
    the protein should have charges and radii but the ligand need not
    """
    itf = oechem.OEInterface()

    if not oechem.OEConfigure(itf, InterfaceData):
        oechem.OEThrow.Fatal("Problem configuring OEInterface!")

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

    duFile = itf.GetString("-du")
    du = oechem.OEDesignUnit()
    if not oechem.OEReadDesignUnit(duFile, du):
        oechem.OEThrow.Fatal("Unable to open %s for reading" % duFile)

    if not du.HasLigand():
        oechem.OEThrow.Fatal("Input designunit %s does not have a ligand bound" % du.GetTitle())

    lig = oechem.OEGraphMol()
    du.GetLigand(lig)
    prot = oechem.OEGraphMol()
    du.GetComponents(prot, oechem.OEDesignUnitComponents_Protein)

    highRes = itf.GetBool("-high_res")
    GetSzmapEnergies(lig, prot, highRes)


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