Grapheme Examples

Depicting BFactor

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2013-2015 OpenEye Scientific Software, Inc.
#############################################################################
# Depicts the BFactor of a ligand and its environment
#############################################################################

import sys
from openeye.oechem import *
from openeye.oedepict import *
from openeye.oegrapheme import *


def main(argv=[__name__]):

    itf = OEInterface()
    OEConfigure(itf, InterfaceData)
    OEConfigureImageWidth(itf, 600.0)
    OEConfigureImageHeight(itf, 600.0)
    OEConfigure2DMolDisplayOptions(itf, OE2DMolDisplaySetup_AromaticStyle)
    OEConfigureSplitMolComplexOptions(itf, OESplitMolComplexSetup_LigName)

    if not OEParseCommandLine(itf, argv):
        return 1

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

    ifs = oemolistream()
    if not ifs.open(iname):
        OEThrow.Fatal("Cannot open input file!")

    ext = OEGetFileExtension(oname)
    if not OEIsRegisteredImageFile(ext):
        OEThrow.Fatal("Unknown image type!")

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

    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()

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

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

    # Calculate average BFactor of the whole complex

    avgbfactor = GetAverageBFactor(complexmol)

    # Calculate minimum and maximum BFactor of the ligand and its environment

    minbfactor, maxbfactor = GetMinAndMaxBFactor(ligand, protein)

    # Attach to each ligand atom the average BFactor of the nearby protein atoms

    stag = "avg residue BFfactor"
    itag = OEGetTag(stag)
    SetAverageBFactorOfNearbyProteinAtoms(ligand, protein, itag)

    OEThrow.Info("Average BFactor of the complex = %+.3f" % avgbfactor)
    OEThrow.Info("Minimum BFactor of the ligand and its environment = %+.3f" % minbfactor)
    OEThrow.Info("Maximum BFactor of the ligand and its environment = %+.3f" % maxbfactor)

    # Create image

    imagewidth, imageheight = OEGetImageWidth(itf), OEGetImageHeight(itf)
    image = OEImage(imagewidth, imageheight)

    mframe = OEImageFrame(image, imagewidth, imageheight * 0.90, OE2DPoint(0.0, 0.0))
    lframe = OEImageFrame(image, imagewidth, imageheight * 0.10, OE2DPoint(0.0, imageheight * 0.90))

    opts = OE2DMolDisplayOptions(mframe.GetWidth(), mframe.GetHeight(), OEScale_AutoScale)
    OESetup2DMolDisplayOptions(opts, itf)
    opts.SetAtomColorStyle(OEAtomColorStyle_WhiteMonochrome)

    # Create BFactor color gradient

    colorg = OELinearColorGradient()
    colorg.AddStop(OEColorStop(0.0, OEDarkBlue))
    colorg.AddStop(OEColorStop(10.0, OELightBlue))
    colorg.AddStop(OEColorStop(25.0, OEYellowTint))
    colorg.AddStop(OEColorStop(50.0, OERed))
    colorg.AddStop(OEColorStop(100.0, OEDarkRose))

    # Prepare ligand for depiction

    OEPrepareDepictionFrom3D(ligand)
    arcfxn = BFactorArcFxn(colorg, itag)
    for atom in ligand.GetAtoms():
        OESetSurfaceArcFxn(ligand, atom, arcfxn)
    opts.SetScale(OEGetMoleculeSurfaceScale(ligand, opts))

    # Render ligand and visualize BFactor

    disp = OE2DMolDisplay(ligand, opts)
    colorbfactor = ColorLigandAtomByBFactor(colorg)
    OEAddGlyph(disp, colorbfactor, OEIsTrueAtom())
    OEDraw2DSurface(disp)
    OERenderMolecule(mframe, disp)

    # Draw color gradient

    opts = OEColorGradientDisplayOptions()
    opts.SetColorStopPrecision(1)
    opts.AddMarkedValue(avgbfactor)
    opts.SetBoxRange(minbfactor, maxbfactor)

    OEDrawColorGradient(lframe, colorg, opts)

    OEWriteImage(oname, image)

    return 0

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


def GetAverageBFactor(mol):

    sumbfactor = 0.0
    for atom in mol.GetAtoms():
        res = OEAtomGetResidue(atom)
        sumbfactor += res.GetBFactor()
    avgbfactor = sumbfactor / mol.NumAtoms()

    return avgbfactor


def ConsiderResidueAtom(atom, res):
    if atom.GetAtomicNum() == OEElemNo_H:
        return False
    if res.GetName() == "HOH":
        return False
    return True


def GetMinAndMaxBFactor(ligand, protein, maxdistance=4.0):

    minbfactor = float("inf")
    maxbfactor = float("-inf")

    # Ligand atoms
    for latom in ligand.GetAtoms(OEIsHeavy()):
        res = OEAtomGetResidue(latom)
        minbfactor = min(minbfactor, res.GetBFactor())
        maxbfactor = max(maxbfactor, res.GetBFactor())

    # Protein atoms close to ligand atoms
    nn = OENearestNbrs(protein, maxdistance)
    for latom in ligand.GetAtoms(OEIsHeavy()):
        for neigh in nn.GetNbrs(latom):
            ratom = neigh.GetBgn()
            res = OEAtomGetResidue(ratom)
            if ConsiderResidueAtom(ratom, res):
                minbfactor = min(minbfactor, res.GetBFactor())
                maxbfactor = max(maxbfactor, res.GetBFactor())

    return minbfactor, maxbfactor


def SetAverageBFactorOfNearbyProteinAtoms(ligand, protein, itag, maxdistance=4.0):

    nn = OENearestNbrs(protein, maxdistance)
    for latom in ligand.GetAtoms(OEIsHeavy()):
        sumbfactor = 0.0
        neighs = []
        for neigh in nn.GetNbrs(latom):
            ratom = neigh.GetBgn()
            res = OEAtomGetResidue(ratom)
            if ConsiderResidueAtom(ratom, res):
                sumbfactor += res.GetBFactor()
                neighs.append(ratom)

        avgbfactor = 0.0
        if len(neighs) > 0:
            avgbfactor = sumbfactor / len(neighs)
        latom.SetDoubleData(itag, avgbfactor)

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


class BFactorArcFxn(OESurfaceArcFxnBase):
    def __init__(self, colorg, itag):
        OESurfaceArcFxnBase.__init__(self)
        self.colorg = colorg
        self.itag = itag

    def __call__(self, image, arc):
        adisp = arc.GetAtomDisplay()
        if adisp is None or not adisp.IsVisible():
            return False

        atom = adisp.GetAtom()
        if atom is None:
            return False

        avgresiduebfactor = atom.GetDoubleData(self.itag)
        if avgresiduebfactor == 0.0:
            return True
        color = self.colorg.GetColorAt(avgresiduebfactor)

        pen = OEPen(color, color, OEFill_Off, 5.0)

        center = arc.GetCenter()
        bAngle = arc.GetBgnAngle()
        eAngle = arc.GetEndAngle()
        radius = arc.GetRadius()

        OEDrawDefaultSurfaceArc(image, center, bAngle, eAngle, radius, pen)

        return True

    def CreateCopy(self):
        return BFactorArcFxn(self.colorg, self.itag).__disown__()

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


class ColorLigandAtomByBFactor(OEAtomGlyphBase):
    def __init__(self, colorg):
        OEAtomGlyphBase.__init__(self)
        self.colorg = colorg

    def RenderGlyph(self, disp, atom):
        adisp = disp.GetAtomDisplay(atom)
        if adisp is None or not adisp.IsVisible():
            return False

        res = OEAtomGetResidue(atom)
        bfactor = res.GetBFactor()
        color = self.colorg.GetColorAt(bfactor)

        pen = OEPen(color, color, OEFill_On, 1.0)
        radius = disp.GetScale() / 3.0

        layer = disp.GetLayer(OELayerPosition_Below)
        OEDrawCircle(layer, OECircleStyle_Default, adisp.GetCoords(), radius, pen)
        return True

    def CreateCopy(self):
        return ColorLigandAtomByBFactor(self.colorg).__disown__()

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

InterfaceData = '''
!BRIEF [-complex] <input> [-out] <output pdf>

!CATEGORY "input/output options :"

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

  !PARAMETER -out
    !ALIAS -o
    !TYPE string
    !REQUIRED true
    !KEYLESS 2
    !VISIBILITY simple
    !BRIEF Output filename
  !END

!END
'''

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

See also

Depicting Active Site Interactions

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2014-2016 OpenEye Scientific Software, Inc.
#############################################################################
# Depicts the interactions of an active site
#############################################################################

import sys
from openeye.oechem import *
from openeye.oedepict import *
from openeye.oegrapheme import *


def main(argv=[__name__]):

    itf = OEInterface()
    OEConfigure(itf, InterfaceData)
    OEConfigureImageWidth(itf, 900.0)
    OEConfigureImageHeight(itf, 600.0)
    OEConfigure2DMolDisplayOptions(itf, OE2DMolDisplaySetup_AromaticStyle)
    OEConfigureSplitMolComplexOptions(itf, OESplitMolComplexSetup_LigName)

    if not OEParseCommandLine(itf, argv):
        return 1

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

    ifs = oemolistream()
    if not ifs.open(iname):
        OEThrow.Fatal("Cannot open input file!")

    ext = OEGetFileExtension(oname)
    if not OEIsRegisteredImageFile(ext):
        OEThrow.Fatal("Unknown image type!")

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

    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 asite.IsValid():
        OEThrow.Fatal("Cannot initialize active site!")
    asite.SetTitle(ligand.GetTitle())

    OEPerceiveInteractionHints(asite)

    OEPrepareActiveSiteDepiction(asite)

    # Depict active site with interactions

    width, height = OEGetImageWidth(itf), OEGetImageHeight(itf)
    image = OEImage(width, height)

    cframe = OEImageFrame(image, width * 0.80, height, OE2DPoint(0.0, 0.0))
    lframe = OEImageFrame(image, width * 0.20, height, OE2DPoint(width * 0.80, 0.0))

    opts = OE2DActiveSiteDisplayOptions(cframe.GetWidth(), cframe.GetHeight())
    OESetup2DMolDisplayOptions(opts, itf)

    adisp = OE2DActiveSiteDisplay(asite, opts)
    OERenderActiveSite(cframe, adisp)

    lopts = OE2DActiveSiteLegendDisplayOptions(10, 1)
    OEDrawActiveSiteLegend(lframe, adisp, lopts)

    OEWriteImage(oname, image)

    return 0


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

InterfaceData = '''
!BRIEF [-complex] <input> [-out] <output image>

!CATEGORY "input/output options :"

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

  !PARAMETER -out
    !ALIAS -o
    !TYPE string
    !REQUIRED true
    !KEYLESS 2
    !VISIBILITY simple
    !BRIEF Output filename
  !END

!END
'''

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

Depicting Shape and Color Atom Overlaps

Note

OEShape TK license is not required to run this example.

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2013-2015 OpenEye Scientific Software, Inc.
#############################################################################
# Depicts shape and color overlap between a 3D reference structure and
# a sets of 3D fit molecules. The molecules have to be pre-aligned.
# The first molecule is expected be the reference.
#############################################################################

import sys
from openeye.oechem import *
from openeye.oedepict import *
from openeye.oegrapheme import *
from openeye.oeshape import *


def main(argv=[__name__]):

    itf = OEInterface(InterfaceData)

    if not OEParseCommandLine(itf, argv):
        return 1

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

    ext = OEGetFileExtension(oname)
    if not OEIsRegisteredMultiPageImageFile(ext):
        OEThrow.Fatal("Unknown multipage image type!")

    ifs = oemolistream()
    if not ifs.open(iname):
        OEThrow.Fatal("Cannot open input molecule file!")

    refmol = OEGraphMol()
    if not OEReadMolecule(ifs, refmol):
        OEThrow.Fatal("Cannot read reference molecule!")

    ropts = OEReportOptions(3, 1)
    ropts.SetHeaderHeight(40.0)
    ropts.SetFooterHeight(20.0)
    report = OEReport(ropts)

    cff = OEColorForceField()
    cff.Init(OEColorFFType_ImplicitMillsDean)
    cffdisplay = OEColorForceFieldDisplay(cff)

    qopts = GetShapeQueryDisplayOptions()
    sopts = GetShapeOverlapDisplayOptions()
    copts = GetColorOverlapDisplayOptions()

    refdisp = OEShapeQueryDisplay(refmol, cff, qopts)

    dots = OEDots(100, 10, "shape overlaps")

    for fitmol in ifs.GetOEGraphMols():

        if maxhits > 0 and dots.GetCounts() >= maxhits:
            break
        dots.Update()

        maincell = report.NewCell()

        grid = OEImageGrid(maincell, 1, 3)
        grid.SetMargins(5.0)
        grid.SetCellGap(5.0)

        # TITLE + SCORE GRAPH + QUERY
        cell = grid.GetCell(1, 1)
        cellw, cellh = cell.GetWidth(), cell.GetHeight()

        font = OEFont(OEFontFamily_Default, OEFontStyle_Bold, 10, OEAlignment_Left, OEBlack)
        pos = OE2DPoint(10.0, 10.0)
        cell.DrawText(pos, fitmol.GetTitle(), font, cell.GetWidth())

        rframe = OEImageFrame(cell, cellw, cellh * 0.35, OE2DPoint(0.0, cellh * 0.10))
        mframe = OEImageFrame(cell, cellw, cellh * 0.50, OE2DPoint(0.0, cellh * 0.50))

        RenderScoreRadial(rframe, fitmol)
        OERenderShapeQuery(mframe, refdisp)

        font = OEFont(OEFontFamily_Default, OEFontStyle_Bold, 8, OEAlignment_Center, OEGrey)
        pos = OE2DPoint(20.0, 10.0)
        mframe.DrawText(pos, "query", font)
        OEDrawCurvedBorder(mframe, OELightGreyPen, 10.0)

        odisp = OEShapeOverlapDisplay(refdisp, fitmol, sopts, copts)

        # SHAPE OVERLAP
        cell = grid.GetCell(1, 2)
        OERenderShapeOverlap(cell, odisp)
        RenderScore(cell, fitmol, "ROCS_ShapeTanimoto", "Shape Tanimoto")

        # COLOR OVERLAP
        cell = grid.GetCell(1, 3)
        OERenderColorOverlap(cell, odisp)
        RenderScore(cell, fitmol, "ROCS_ColorTanimoto", "Color Tanimoto")

        OEDrawCurvedBorder(maincell, OELightGreyPen, 10.0)

    dots.Total()

    cffopts = OEColorForceFieldLegendDisplayOptions(1, 6)
    for header in report.GetHeaders():
        OEDrawColorForceFieldLegend(header, cffdisplay, cffopts)
        OEDrawCurvedBorder(header, OELightGreyPen, 10.0)

    font = OEFont(OEFontFamily_Default, OEFontStyle_Default, 12, OEAlignment_Center, OEBlack)
    for idx, footer in enumerate(report.GetFooters()):
        OEDrawTextToCenter(footer, "-" + str(idx + 1) + "-", font)

    OEWriteReport(oname, report)

    return 0


def AddCommonDisplayOptions(opts):
    opts.SetTitleLocation(OETitleLocation_Hidden)
    opts.SetAtomLabelFontScale(1.5)
    pen = OEPen(OEBlack, OEBlack, OEFill_Off, 1.5)
    opts.SetDefaultBondPen(pen)


def GetShapeQueryDisplayOptions():

    qopts = OEShapeQueryDisplayOptions()
    AddCommonDisplayOptions(qopts)
    arcpen = OEPen(OELightGreyPen)
    qopts.SetSurfaceArcFxn(OEDefaultArcFxn(arcpen))
    qopts.SetDepictOrientation(OEDepictOrientation_Square)
    return qopts


def GetShapeOverlapDisplayOptions():

    sopts = OEShapeOverlapDisplayOptions()
    AddCommonDisplayOptions(sopts)
    arcpen = OEPen(OEGrey, OEGrey, OEFill_Off, 1.0, 0x1111)
    sopts.SetQuerySurfaceArcFxn(OEDefaultArcFxn(arcpen))
    sopts.SetOverlapColor(OEColor(110, 110, 190))

    return sopts


def GetColorOverlapDisplayOptions():

    copts = OEColorOverlapDisplayOptions()
    AddCommonDisplayOptions(copts)
    arcpen = OEPen(OEGrey, OEGrey, OEFill_Off, 1.0, 0x1111)
    copts.SetQuerySurfaceArcFxn(OEDefaultArcFxn(arcpen))

    return copts


def GetScore(mol, sdtag):

    if OEHasSDData(mol, sdtag):
        return float(OEGetSDData(mol, sdtag))
    return 0.0


def RenderScoreRadial(image, mol):

    sscore = max(min(GetScore(mol, "ROCS_ShapeTanimoto"), 1.0), 0.00)
    cscore = max(min(GetScore(mol, "ROCS_ColorTanimoto"), 1.0), 0.00)
    if sscore > 0.0 or cscore > 0.0:
        scores = OEDoubleVector([sscore, cscore])
        OEDrawROCSScores(image, scores)


def RenderScore(image, mol, sdtag, label):

    score = GetScore(mol, sdtag)
    if score == 0.0:
        return

    w, h = image.GetWidth(), image.GetHeight()
    frame = OEImageFrame(image, w, h * 0.10, OE2DPoint(0.0, h * 0.90))
    font = OEFont(OEFontFamily_Default, OEFontStyle_Default, 9, OEAlignment_Center, OEBlack)
    OEDrawTextToCenter(frame, label + " = " + str(score), font)


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

InterfaceData = '''
!BRIEF [-in] <input> [-out] <output pdf>

!CATEGORY "input/output options:" 0

  !PARAMETER -in
    !ALIAS -i
    !TYPE string
    !REQUIRED true
    !KEYLESS 1
    !VISIBILITY simple
    !BRIEF Input molecule filename
    !DETAIL
         The first molecule in the file is expected to be the reference
         molecule
  !END

  !PARAMETER -out
    !ALIAS -o
    !TYPE string
    !REQUIRED true
    !KEYLESS 2
    !VISIBILITY simple
    !BRIEF Output image filename
  !END

!END

!CATEGORY "general options:" 1

  !PARAMETER -maxhits
    !ALIAS -mhits
    !TYPE int
    !REQUIRED false
    !DEFAULT 0
    !LEGAL_RANGE 0 500
    !VISIBILITY simple
    !BRIEF Maximum number of hits depicted
    !DETAIL
            The default of 0 means there is no limit.
  !END

!END
'''

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