Depicting Atom Contributions of XLogP

Problem

You want to depict the contribution of each atom to the total XLogP on your molecule diagram. See example in Figure 1.

Figure 1. Example of depicting the atom contributions of XLogP

Ingredients

OEChem TK - cheminformatics toolkit MolProp TK - molecular property calculation and filtering toolkit Quacpac TK - tautomers & charges toolkit OEDepict TK - molecule depiction toolkit Grapheme TK - molecule and property visualization toolkit

Difficulty Level

🌶️ 🌶️

Download

Download code

xlogp2img.py and xlogp2pdf.py

See also the Usage (xlogp2img) and Usage (xlogp2pdf) subsections.

Source Code

xlogp2img

#!/usr/bin/env python3
# (C) 2023 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.


"""Calculates XLogP and visualizes atom contributions using property map."""

import argparse
import io
import os
import sys
from pathlib import Path

from openeye import oechem, oedepict, oegrapheme, oemolprop, oequacpac
from PIL import Image
from rich_argparse import HelpPreviewAction, RichHelpFormatter

__SCRIPT_NAME__ = Path(__file__).absolute().stem
__SCRIPT_DESC__ = "Depict XLogP of molecule (atom based)."
__SCRIPT_TOOLKITS__ = ["oechem", "oedepict", "oegrapheme", "oemolprop", "oequacpac"]
__SCRIPT_CATEGORIES__ = ["depiction"]


def parse_options() -> argparse.Namespace:
    """Set up command line options."""
    parser = argparse.ArgumentParser(
        add_help=True,
        formatter_class=RichHelpFormatter,
        description="[yellow]" + __SCRIPT_DESC__ + "[/yellow]",
    )
    parser.add_argument("--help-image", action=HelpPreviewAction)

    input_group = parser.add_argument_group("Input options")
    input_group.add_argument(
        "--mol",
        type=str,
        required=True,
        metavar="MOL-FILE",
        help="input molecule file (oeb, sdf)",
    )

    image_group = parser.add_argument_group("Image options")
    image_group.add_argument(
        "--image",
        type=str,
        required=False,
        metavar="IMAGE-FILE",
        help="output image file (PNG, SVG) (required: %(required)s) -- if no output is provided the image will be displayed on the  screen",
    )
    image_group.add_argument(
        "--width",
        type=int,
        default=800,
        help="width of output image (default: %(default)s)",
    )
    image_group.add_argument(
        "--height",
        type=int,
        default=400,
        help="height of output image (default: %(default)s)",
    )
    return parser.parse_args()


def main() -> int:
    """Visualizes XLogP atom contributions."""
    args = parse_options()

    _check_image_file(args)
    mol = _get_molecule(args)

    image = oedepict.OEImage(args.width, args.height)

    opts = oedepict.OE2DMolDisplayOptions(
        args.width, args.height, oedepict.OEScale_AutoScale
    )
    opts.SetAtomColorStyle(oedepict.OEAtomColorStyle_WhiteMonochrome)

    # depict molecule with XLogP atom contributions

    oedepict.OEPrepareDepiction(mol)
    depict_molecule_xlogp(image, mol, opts)

    if args.image:
        oedepict.OEWriteImage(args.image, image)
    else:
        _img = Image.open(io.BytesIO(oedepict.OEWriteImageToBytes("png", image)))
        _img.show()

    return os.EX_OK


def set_atom_properties(mol: oechem.OEMolBase, tag: int) -> None:
    """Attach the XLogP atom contribution to each atom with the given tag."""
    oequacpac.OERemoveFormalCharge(mol)

    atom_values = oechem.OEFloatArray(mol.GetMaxAtomIdx())
    xlogp = oemolprop.OEGetXLogP(mol, atom_values)

    mol.SetTitle(f"{mol.GetTitle()} -- OEXLogP = {xlogp:.2f}")

    for atom in mol.GetAtoms():
        val = atom_values[atom.GetIdx()]
        atom.SetData(tag, val)


def depict_molecule_xlogp(
    image: oedepict.OEImageBase,
    mol: oechem.OEMolBase,
    opts: oedepict.OE2DMolDisplayOptions,
) -> None:
    """Generate an image of a molecule depicting the atom contribution of XLogP."""
    scale = oegrapheme.OEGetMoleculeSurfaceScale(mol, opts)
    opts.SetScale(scale)

    str_tag = "XLogP"
    int_tag = oechem.OEGetTag(str_tag)
    set_atom_properties(mol, int_tag)

    disp = oedepict.OE2DMolDisplay(mol, opts)

    prop_map = oegrapheme.OE2DPropMap(opts.GetBackgroundColor())
    prop_map.SetNegativeColor(oechem.OEDarkGreen)
    prop_map.SetPositiveColor(oechem.OEDarkPurple)
    prop_map.SetLegendLocation(oegrapheme.OELegendLocation_Left)
    prop_map.Render(disp, str_tag)

    oedepict.OERenderMolecule(image, disp)


def _check_image_file(args: argparse.Namespace) -> None:
    # script will terminate if there is some issues
    if not args.image:
        # image will be displayed on the screen
        return
    ext = Path(args.image).suffix[1:].upper()
    if not oedepict.OEIsRegisteredImageFile(ext):
        oechem.OEThrow.Fatal("Unknown image output type!")

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


def _get_molecule(args: argparse.Namespace) -> oechem.OEMolBase:
    ifs = oechem.oemolistream()
    if not ifs.open(args.mol):
        oechem.OEThrow.Fatal(f"Cannot open {args.mol} input file!")

    mol = oechem.OEGraphMol()
    if not oechem.OEReadMolecule(ifs, mol):
        oechem.OEThrow.Fatal(f"Cannot read molecule from {args.mol} input file!")

    return mol


setattr(main, "__SCRIPT_NAME__", __SCRIPT_NAME__)
setattr(main, "__SCRIPT_DESC__", __SCRIPT_DESC__)
setattr(main, "__SCRIPT_TOOLKITS__", __SCRIPT_TOOLKITS__)
setattr(main, "__SCRIPT_CATEGORIES__", __SCRIPT_CATEGORIES__)

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

xlogp2pdf

#!/usr/bin/env python3
# (C) 2023 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.


"""Calculates XLogP of set of molecules and visualizes the atom contributions using property map."""

import argparse
import os
import pathlib
import sys

import rich.console
from openeye import oechem, oedepict, oegrapheme, oemolprop, oequacpac
from rich_argparse import HelpPreviewAction, RichHelpFormatter

__SCRIPT_NAME__ = pathlib.Path(__file__).absolute().stem
__SCRIPT_DESC__ = "Depict XLogP of set of molecules."
__SCRIPT_TOOLKITS__ = ["oechem", "oedepict", "oegrapheme", "oemolprop", "oequacpac"]
__SCRIPT_CATEGORIES__ = ["depiction"]


def parse_options() -> argparse.Namespace:
    """Set up command line options."""
    parser = argparse.ArgumentParser(
        add_help=True,
        formatter_class=RichHelpFormatter,
        description="[yellow]" + __SCRIPT_DESC__ + "[/yellow]",
    )
    parser.add_argument("--help-image", action=HelpPreviewAction)

    input_group = parser.add_argument_group("Input options")
    input_group.add_argument(
        "--mol",
        type=str,
        required=True,
        metavar="MOL-FILE",
        help="input multi-conformer molecule file (oeb, sdf)",
    )

    report_group = parser.add_argument_group("Report options")
    report_group.add_argument(
        "--report",
        type=str,
        required=True,
        metavar="REPORT-FILE",
        help="output report file (PDF)",
    )
    report_group.add_argument(
        "--rows",
        type=int,
        default=3,
        choices=range(2, 6),
        metavar="N",
        help="number of rows per page (default: %(default)s)",
    )
    report_group.add_argument(
        "--cols",
        type=int,
        default=2,
        choices=range(1, 3),
        metavar="N",
        help="number of columns per page (default: %(default)s)",
    )
    report_group.add_argument(
        "--page-by-page",
        action="store_true",
        help="write pages of report to separate numbered image files",
    )
    return parser.parse_args()


def main() -> int:
    """Visualizes XLogP of set of molecules."""
    args = parse_options()

    _check_report_file(args)

    mols: list[oechem.OEMolBase] = _read_molecules(args.mol)
    console = rich.console.Console()
    console.print(f"Imported {len(mols)} molecules from {args.mol}")

    # initialize multi-page report

    report_opts = oedepict.OEReportOptions(args.rows, args.cols)
    report_opts.SetHeaderHeight(35)
    report_opts.SetFooterHeight(45)
    report_opts.SetPageMargins(10)
    report_opts.SetCellGap(5)
    report = oedepict.OEReport(report_opts)

    # setup depiction options

    width, height = report.GetCellWidth(), report.GetCellHeight()
    opts = oedepict.OE2DMolDisplayOptions(width, height, oedepict.OEScale_AutoScale)
    opts.SetAtomColorStyle(oedepict.OEAtomColorStyle_WhiteMonochrome)

    # depict molecule with XLogP atom contributions
    depict_molecules_xlogp(report, mols, opts)

    if args.page_by_page:
        oedepict.OEWriteReportPageByPage(args.report, report)
    else:
        oedepict.OEWriteReport(args.report, report)

    return os.EX_OK


def set_atom_properties(
    mol: oechem.OEMolBase, tag: int, min_value: float, max_value: float
) -> tuple[float, float]:
    """Attache the XLogP atom contribution to each atom with the given tag."""
    oequacpac.OERemoveFormalCharge(mol)

    atom_values = oechem.OEFloatArray(mol.GetMaxAtomIdx())
    xlogp = oemolprop.OEGetXLogP(mol, atom_values)

    mol.SetTitle(f"{mol.GetTitle()} -- OEXLogP = {xlogp:.2f}")

    for atom in mol.GetAtoms():
        val = atom_values[atom.GetIdx()]
        atom.SetData(tag, val)
        min_value = min(min_value, val)
        max_value = max(max_value, val)

    return min_value, max_value


def depict_molecules_xlogp(
    report: oedepict.OEReport,
    mols: list[oechem.OEMolBase],
    opts: oedepict.OE2DMolDisplayOptions,
) -> None:
    """Generate a report of molecules depicting the atom contribution of XLogP."""
    str_tag = "XLogP"
    int_tag = oechem.OEGetTag(str_tag)

    min_value, max_value = float("inf"), float("-inf")
    for mol in mols:
        min_value, max_value = set_atom_properties(mol, int_tag, min_value, max_value)

    mol_scale = float("inf")
    for mol in mols:
        oedepict.OEPrepareDepiction(mol)
        mol_scale = min(mol_scale, oegrapheme.OEGetMoleculeSurfaceScale(mol, opts))
    opts.SetScale(mol_scale)

    prop_map = oegrapheme.OE2DPropMap(opts.GetBackgroundColor())
    prop_map.SetNegativeColor(oechem.OEDarkGreen)
    prop_map.SetPositiveColor(oechem.OEDarkPurple)
    prop_map.SetLegendLocation(oegrapheme.OELegendLocation_Left)
    prop_map.SetMinValue(min_value)
    prop_map.SetMaxValue(max_value)

    for mol in mols:
        disp = oedepict.OE2DMolDisplay(mol, opts)
        prop_map.Render(disp, str_tag)
        cell = report.NewCell()
        oedepict.OERenderMolecule(cell, disp)


def _check_report_file(args: argparse.Namespace) -> bool:
    ext = pathlib.Path(args.report).suffix[1:]
    if not oedepict.OEIsRegisteredImageFile(ext):
        oechem.OEThrow.Fatal("Unknown image outout type!")

    if not args.page_by_page and not oedepict.OEIsRegisteredMultiPageImageFile(ext):
        oechem.OEThrow.Warning("Report will be generated into separate pages!")
        args.page_by_page = True

    return True


def _read_molecules(filename: str) -> list[oechem.OEMolBase]:
    ifs = oechem.oemolistream()
    if not ifs.open(filename):
        oechem.OEThrow.Fatal(f"Cannot open {filename} input file!")

    mols: list[oechem.OEMolBase] = [oechem.OEGraphMol(m) for m in ifs.GetOEGraphMols()]
    if not mols:
        oechem.OEThrow.Fatal(f"No molecules could be read from {filename}")
    return mols


setattr(main, "__SCRIPT_NAME__", __SCRIPT_NAME__)
setattr(main, "__SCRIPT_DESC__", __SCRIPT_DESC__)
setattr(main, "__SCRIPT_TOOLKITS__", __SCRIPT_TOOLKITS__)
setattr(main, "__SCRIPT_CATEGORIES__", __SCRIPT_CATEGORIES__)

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

Solution

The code snippet below shows how to calculate the total XLogP of a molecule along with the atom contributions by calling the OEGetXLogP function. Each atom contribution is then attached to the relevant atom as generic data with the given tag.

def set_atom_properties(mol: oechem.OEMolBase, tag: int) -> None:
    """Attach the XLogP atom contribution to each atom with the given tag."""
    oequacpac.OERemoveFormalCharge(mol)

    atom_values = oechem.OEFloatArray(mol.GetMaxAtomIdx())
    xlogp = oemolprop.OEGetXLogP(mol, atom_values)

    mol.SetTitle(f"{mol.GetTitle()} -- OEXLogP = {xlogp:.2f}")

    for atom in mol.GetAtoms():
        val = atom_values[atom.GetIdx()]
        atom.SetData(tag, val)

The depict_molecule_xlogp function below shows how to project the atom contributions of the total XLogP into a 2D molecular diagram using the OE2DPropMap class. After constructing the OE2DMolDisplay object to depict a molecule, an OE2DPropMap object is initialized and its properties are set that determine how the data is going to be visualized. For example in this case the negative color XLogP atom contributions are going to be represented by dark green, while the positive values are visualized by using dark purple.

When the OE2DPropMap.Render method is called with same tag that was used in the set_atom_properties function, the properties that were attached to the atoms as generic data are retrieved and a 2D grid is generated underneath the molecular diagram. The OELinearColorGradient object that is used to assign colors to the cells of the grid is also rendered based on the option set by the OE2DPropMap.SetLegendLocation method . You can see the result in Figure 1.

def depict_molecule_xlogp(
    image: oedepict.OEImageBase,
    mol: oechem.OEMolBase,
    opts: oedepict.OE2DMolDisplayOptions,
) -> None:
    """Generate an image of a molecule depicting the atom contribution of XLogP."""
    scale = oegrapheme.OEGetMoleculeSurfaceScale(mol, opts)
    opts.SetScale(scale)

    str_tag = "XLogP"
    int_tag = oechem.OEGetTag(str_tag)
    set_atom_properties(mol, int_tag)

    disp = oedepict.OE2DMolDisplay(mol, opts)

    prop_map = oegrapheme.OE2DPropMap(opts.GetBackgroundColor())
    prop_map.SetNegativeColor(oechem.OEDarkGreen)
    prop_map.SetPositiveColor(oechem.OEDarkPurple)
    prop_map.SetLegendLocation(oegrapheme.OELegendLocation_Left)
    prop_map.Render(disp, str_tag)

    oedepict.OERenderMolecule(image, disp)

Hint

You can easily adapt this example to visualize other atom properties by writing your own set_atom_properties function.

Usage (xlogp2img)

See Download section to download the script.

> xlogp2img --help

The following commands will generate the image shown in Figure 1 for input example.ism

> xlogp2img --mol example.ism --image image.svg

Discussion

The example above shows how to visualize the atom contributions for a single molecule, however you might want to visualize the XLogP data for a set of molecules. In this case the set_atom_properties function not only attaches the XLogP contributions to the relevant atom, but also calculates the minimum and maximum atom contributions for the whole molecule set.

def set_atom_properties(
    mol: oechem.OEMolBase, tag: int, min_value: float, max_value: float
) -> tuple[float, float]:
    """Attache the XLogP atom contribution to each atom with the given tag."""
    oequacpac.OERemoveFormalCharge(mol)

    atom_values = oechem.OEFloatArray(mol.GetMaxAtomIdx())
    xlogp = oemolprop.OEGetXLogP(mol, atom_values)

    mol.SetTitle(f"{mol.GetTitle()} -- OEXLogP = {xlogp:.2f}")

    for atom in mol.GetAtoms():
        val = atom_values[atom.GetIdx()]
        atom.SetData(tag, val)
        min_value = min(min_value, val)
        max_value = max(max_value, val)

    return min_value, max_value

These minimum and maximum values are used to initialize the value range of the linear color gradient of the property map depict_molecules_xlogp function below). Each molecule, along with its property map, is then rendered into a cell of an OEReport object. The OEReport class is a layout manager allowing generation of multi-page images in a convenient way. You can see the generated multi-page PDF in Table 1. While the value range of the color gradients depicted alongside the molecules represents the range for the whole set, the black box rendered on each color gradient represents the minimum and maximum XLogP atom contributions for the corresponding molecule.

def depict_molecules_xlogp(
    report: oedepict.OEReport,
    mols: list[oechem.OEMolBase],
    opts: oedepict.OE2DMolDisplayOptions,
) -> None:
    """Generate a report of molecules depicting the atom contribution of XLogP."""
    str_tag = "XLogP"
    int_tag = oechem.OEGetTag(str_tag)

    min_value, max_value = float("inf"), float("-inf")
    for mol in mols:
        min_value, max_value = set_atom_properties(mol, int_tag, min_value, max_value)

    mol_scale = float("inf")
    for mol in mols:
        oedepict.OEPrepareDepiction(mol)
        mol_scale = min(mol_scale, oegrapheme.OEGetMoleculeSurfaceScale(mol, opts))
    opts.SetScale(mol_scale)

    prop_map = oegrapheme.OE2DPropMap(opts.GetBackgroundColor())
    prop_map.SetNegativeColor(oechem.OEDarkGreen)
    prop_map.SetPositiveColor(oechem.OEDarkPurple)
    prop_map.SetLegendLocation(oegrapheme.OELegendLocation_Left)
    prop_map.SetMinValue(min_value)
    prop_map.SetMaxValue(max_value)

    for mol in mols:
        disp = oedepict.OE2DMolDisplay(mol, opts)
        prop_map.Render(disp, str_tag)
        cell = report.NewCell()
        oedepict.OERenderMolecule(cell, disp)

Usage (xlogp2pdf)

See Download section to download the script.

> xlogp2pdf --help

> xlogp2pdf --rows 2 --cols 1 --mol examples.ism --report report.pdf

The following commands will generate the image shown in Table 1 for input examples.ism

> xlogp2pdf --rows 2 --cols 1 --mol examples.ism --report report.pdf

Table 1. Example of depicting atom contributions of XLogP for a set of molecules (The pages are reduced here for visualization convenience)

page 1	page 2	page 3

Related Python Cookbook Recipes

• Depicting Fragment Contributions of XLogP

See also in OEChem TK manual

Theory

• Generic Data chapter

API

• OELinearColorGradient class

See also in MolProp TK manual

API

• OEGetXLogP function

See also in OEDepict TK manual

Theory

• Molecule Depiction chapter

• Multi Page Reports section

API

• OE2DMolDisplay class

• OE2DMolDisplayOptions class

• OEImage class

• OERenderMolecule function

• OEReport class

See also in Grapheme^TM TK manual

Theory

• Depicting Atom Property Maps chapter

API

• OE2DPropMap class