Visualizing Protein-Ligand Unpaired Interactions

Problem

You want to visualize protein-ligand unpaired and clash interactions. See example in Figure 1.

Figure 1. Example of visualizing protein-ligand unpaired and clash interactions (PDB: 1NQ2))

Ingredients

OEChem TK - cheminformatics toolkit (including OEBio TK) OEDepict TK - molecule depiction toolkit Grapheme TK - molecule and property visualization toolkit

Difficulty level

🌶️ 🌶️

Download

Download code

unpairedmap2img.py

See also Usage subsection.

Source Code

unpairedmap2img

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

"""Depict the unpaired and clash interactions of an active site."""

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

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

__SCRIPT_NAME__ = Path(__file__).absolute().stem
__SCRIPT_DESC__ = "Depict the unpaired and clash interactions of an active site."
__SCRIPT_TOOLKITS__ = ["oechem", "oedepict", "oegrapheme"]
__SCRIPT_CATEGORIES__ = ["visualization", "ligand-protein interactions"]


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 options
    input_group = parser.add_argument_group("Input ligand-protein complex")
    exclusive_input_group = input_group.add_mutually_exclusive_group(required=True)
    exclusive_input_group.add_argument(
        "--complex",
        type=str,
        required=False,
        metavar="PDB-FILE",
        help="input PDB file of the ligand-protein complex",
    )
    exclusive_input_group.add_argument(
        "--design-unit",
        "--du",
        type=str,
        metavar="DU-FILE",
        help="input design unit file",
    )
    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=900,
        help="width of output image (default: %(default)s)",
    )
    image_group.add_argument(
        "--height",
        type=int,
        default=600,
        help="height of output image (default: %(default)s)",
    )
    image_group.add_argument(
        "--interactive-legend",
        default=False,
        action="store_true",
        help="visualize legend on mouse hover (SVG-only feature) (default: %(default)s)",
    )
    return parser.parse_args()


def main() -> int:
    """Depict the unpaired and clash interactions of an active site."""
    args = parse_options()

    _check_image_file(args)

    if args.complex:
        protein, ligand = get_protein_and_ligand_from_pdb(args.complex)
    elif args.design_unit:
        protein, ligand = get_protein_and_ligand_from_design_unit(args.design_unit)
    else:
        oechem.OEThrow.Fatal("Invalid input option!")

    # depict unpaired interaction map
    image = oedepict.OEImage(args.width, args.height)

    cell_width, cell_height = args.width, args.height
    if not args.interactive_legend:
        cell_width = cell_width * 0.8

    opts = oegrapheme.OE2DActiveSiteDisplayOptions(cell_width, cell_height)
    opts.SetRenderInteractiveLegend(args.interactive_legend)

    if args.interactive_legend:
        depict_unpaired_map(image, protein, ligand, opts)
    else:
        main_frame = oedepict.OEImageFrame(
            image,
            args.width * 0.80,
            args.height,
            oedepict.OE2DPoint(args.width * 0.2, 0.0),
        )
        legend_frame = oedepict.OEImageFrame(
            image,
            args.width * 0.20,
            args.height,
            oedepict.OE2DPoint(args.width * 0.0, 0.0),
        )
        depict_unpaired_map(main_frame, protein, ligand, opts, legend_frame)

    if (
        args.image
        and Path(args.image).suffix[1:].lower() == "svg"
        and args.interactive_legend
    ):
        icon_scale = 0.5
        oedepict.OEAddInteractiveIcon(
            image, oedepict.OEIconLocation_TopRight, icon_scale
        )
    oedepict.OEDrawCurvedBorder(image, oedepict.OELightGreyPen, 10.0)

    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 depict_unpaired_map(
    image: oedepict.OEImageBase,
    protein: oechem.OEMolBase,
    ligand: oechem.OEMolBase,
    depict_options: oegrapheme.OE2DActiveSiteDisplayOptions,
    legend_frame: oedepict.OEImageBase | None = None,
) -> None:
    """Depict unpaired interaction map."""
    # perceive interactions
    active_site = oechem.OEInteractionHintContainer(protein, ligand)
    if not active_site.IsValid():
        oechem.OEThrow.Fatal("Cannot initialize active site!")
    active_site.SetTitle(ligand.GetTitle())

    oechem.OEPerceiveInteractionHints(active_site)

    # depiction

    oegrapheme.OEPrepareActiveSiteDepiction(active_site)
    active_site_disp = oegrapheme.OE2DActiveSiteDisplay(active_site, depict_options)
    oegrapheme.OERenderUnpairedInteractionMap(image, active_site_disp)

    if legend_frame:
        legend_options = oegrapheme.OE2DActiveSiteLegendDisplayOptions(12, 1)
        oegrapheme.OEDrawUnpairedInteractionMapLegend(
            legend_frame, active_site_disp, legend_options
        )


def get_protein_and_ligand_from_pdb(
    pdb_filename: str,
) -> tuple[oechem.OEMolBase, oechem.OEMolBase]:
    """Read protein and and ligand from from pdb/cif file."""
    ifs = oechem.oemolistream()
    if not ifs.open(pdb_filename):
        oechem.OEThrow.Fatal(f"Unable to open {pdb_filename} for reading")

    complex_mol = oechem.OEGraphMol()
    if not oechem.OEReadMolecule(ifs, complex_mol):
        oechem.OEThrow.Fatal(f"Unable to read complex from {pdb_filename}")

    if not oechem.OEHasResidues(complex_mol):
        oechem.OEPerceiveResidues(complex_mol, oechem.OEPreserveResInfo_All)

    # separate ligand and protein
    split_opts = oechem.OESplitMolComplexOptions()
    ligand = oechem.OEGraphMol()
    protein = oechem.OEGraphMol()
    water = oechem.OEGraphMol()
    other = oechem.OEGraphMol()

    split_opts.SetProteinFilter(
        oechem.OEOrRoleSet(split_opts.GetProteinFilter(), split_opts.GetWaterFilter())
    )
    split_opts.SetWaterFilter(
        oechem.OEMolComplexFilterFactory(oechem.OEMolComplexFilterCategory_Nothing)
    )

    oechem.OESplitMolComplex(ligand, protein, water, other, complex_mol, split_opts)

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

    return protein, ligand


def get_protein_and_ligand_from_design_unit(
    filename: str,
) -> tuple[oechem.OEMolBase, oechem.OEMolBase]:
    """Read protein and and ligand from from design unit file."""
    du = oechem.OEDesignUnit()
    if not oechem.OEIsReadableDesignUnit(filename) or not oechem.OEReadDesignUnit(
        filename, du
    ):
        oechem.OEThrow.Fatal("Cannot read design unit.")

    protein = oechem.OEGraphMol()
    if not du.GetComponents(protein, oechem.OEDesignUnitComponents_TargetComplex):
        oechem.OEThrow.Fatal("Could not extract protein from the design unit.")

    ligand = oechem.OEGraphMol()
    if not du.GetLigand(ligand):
        oechem.OEThrow.Fatal("Could not extract ligand from the design unit.")

    return (protein, ligand)


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!")


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 depict_unpairedmap illustrates how simple it is to generate these images.

1. OEInteractionHintContainer object is constructed that stores information about possible interactions between the ligand and the protein.

2. The interactions are perceived by calling the OEPerceiveInteractionHints function.

3. The active site is then prepared for 2D depiction by invoking the OEPrepareActiveSiteDepiction function.

4. When the OE2DActiveSiteDisplay object is constructed, residues are positioned around the ligand close to those atoms which they are interacting with.

5. The OERenderUnpairedInteractionMap function generates an image that displays the clash and unpaired interactions detected in the ligand and in nearby residues.

6. The legend associated with the unpaired map is rendered by invoking the OEDrawUnpairedInteractionMapLegend function.

def depict_unpaired_map(
    image: oedepict.OEImageBase,
    protein: oechem.OEMolBase,
    ligand: oechem.OEMolBase,
    depict_options: oegrapheme.OE2DActiveSiteDisplayOptions,
    legend_frame: oedepict.OEImageBase | None = None,
) -> None:
    """Depict unpaired interaction map."""
    # perceive interactions
    active_site = oechem.OEInteractionHintContainer(protein, ligand)
    if not active_site.IsValid():
        oechem.OEThrow.Fatal("Cannot initialize active site!")
    active_site.SetTitle(ligand.GetTitle())

    oechem.OEPerceiveInteractionHints(active_site)

    # depiction

    oegrapheme.OEPrepareActiveSiteDepiction(active_site)
    active_site_disp = oegrapheme.OE2DActiveSiteDisplay(active_site, depict_options)
    oegrapheme.OERenderUnpairedInteractionMap(image, active_site_disp)

    if legend_frame:
        legend_options = oegrapheme.OE2DActiveSiteLegendDisplayOptions(12, 1)
        oegrapheme.OEDrawUnpairedInteractionMapLegend(
            legend_frame, active_site_disp, legend_options
        )

Usage

See Download section to download the script.

> unpairedmap2img --help

Visualizing 1YWR_DU_0.oedu design unit of 1YWR.

> unpairedmap2img --design-unit 1YWR_DU_0.oedu --interactive-legend --image image.svg

Discussion

Interaction Perception

Currently the OEPerceiveInteractionHints function perceives the following interaction types:

Table 1. Interaction types currently available in OEChem TK

name	corresponding interaction class	corresponding interaction type namespace
cation-pi	OECationPiInteractionHint	OECationPiInteractionHintType
chelator	OEChelatorInteractionHint	OEChelatorInteractionHintType
clash	OEClashInteractionHint	None
contact	OEContactInteractionHint	None
covalent	OECovalentInteractionHint	None
halogen bond	OEHalogenBondInteractionHint	OEHalogenBondInteractionHintType
hydrogen bond	OEHBondInteractionHint	OEHBondInteractionHintType
salt-bridge	OESaltBridgeInteractionHint	OESaltBridgeInteractionHintType
stacking (T and Pi)	OEStackingInteractionHint	OEStackingInteractionHintType

The default geometric parameters used by the OEPerceiveInteractionHints function have been set based on literature data ([Kumar-2002], [Cavallo-2016], [Bissantz-2010], and [Marcou-2007] ). The interaction parameters can be customized by using the OEPerceiveInteractionOptions class.

Unpaired Interaction Depiction

The OERenderUnpairedInteractionMap function currently visualizes the following interactions detected by the OEPerceiveInteractionHints function.

• Atom clash interaction

  When visualizing protein-ligand atom clashes, a red outline of a residue circle indicates that there are one or more atoms in that residue which are too close to some ligand atom(s). Clashing ligand atoms are marked with a red arc that is directed towards the corresponding clashing residue. The red shading on the grey line representing the shape of the pocket is used to identify atom clashes easily.

  Table 2. Examples of visualizing atom clash(es)

  See also

  • OEClashInteractionHint class in the OEChem TK manual

• Unpaired types of the hydrogen bonding interaction

  An unpaired hydrogen bond interaction is detected:

  • if there is no ligand/protein acceptor atom that could interact with a protein/ligand donor atom

  • if there is no ligand/protein donor atom that could interact with a protein/ligand acceptor atom

  Different linker types are used to mark unpaired acceptor and donor hydrogen bond interactions. Please note that since these interactions are unpaired, the direction of the linkers has no real spatial meaning. Ligand linkers are directed away from the ligand, while protein linkers are directed towards the ligand.

  Table 3. Examples of visualizing unpaired hydrogen bond interactions

  unpaired ligand acceptor	unpaired protein acceptor	unpaired ligand donor	unpaired protein donor

  See also

  • OEHBondInteractionHint class in the OEChem TK manual

  • OEHBondInteractionHintType namespace in the OEChem TK manual

• Clash types of the hydrogen bonding interaction

  A hydrogen bond clash interaction is detected:

  • if an acceptor ligand atom interacts with an acceptor protein atom

  • if a donor ligand atom interacts with an donor protein atom

  Table 4. Examples of visualizing hydrogen bond clash interactions

  donor-donor clash	acceptor-acceptor clash

  See also

  • OEHBondInteractionHint class in the OEChem TK manual

  • OEHBondInteractionHintType namespace in the OEChem TK manual

• Unpaired types of the salt-bridge interaction

  An unpaired salt bridge interaction is detected if there is a positively / negatively charged functional group either in the ligand or in nearby protein without a matching negatively / positively charged functional group, respectively.

  Different linker types are used to mark unpaired positive and negative salt-bridge interactions. Please note that since these interactions are unpaired, the direction of the linkers has no real spatial meaning. Ligand linkers are directed away from the ligand, while protein linkers are directed towards the ligand.

  Table 5. Examples of visualizing unpaired salt-bridge interactions

  unpaired ligand positive	unpaired protein positive	unpaired ligand negative	unpaired protein negative

  See also

  • OESaltBridgeInteractionHint class in the OEChem TK manual

  • OESaltBridgeInteractionHintType namespace in the OEChem TK manual

Hydrogen Position Optimization

Since interaction perception depends on the position of hydrogens, it is highly recommended to optimize those positions prior to perceiving the interactions. The two images below reveal the effect of optimizing the hydrogen bond network in a protein-ligand complex: fewer atom clashes and fewer unpaired hydrogen bond interactions.

Table 6. Example of visualizing protein-ligand unpaired and clash interactions before and after optimizing hydrogen positions

original complex 1D3H.pdb downloaded from PDB Database	design unit 1D3H_DU_0.oedu generated with Spruce TK

See also

• OEPlaceHydrogens function in the OEChem TK manual

• Protein Preparation chapter in the OEChem TK manual

• Creating OEDesignUnits from a PDB file chapter in the Spruce TK manual

Unpaired Map vs Active Site Interaction Map

An unpaired interaction map provides a complementary view to the more common active site interaction map. While the interaction map (on the right) depicts interactions between the ligand and protein, the unpaired map (on the left) illustrates interactions that could contribute to binding but are not formed in the complex. Together, these two maps of the protein-ligand binding site provide insights into protein-ligand interactions and communicate complex 3D structural results to medicinal chemists in a directly actionable way.

Table 2. Example of visualizing unpaired and active site interaction maps

unpaired interaction map of 1D3H_DU_0.oedu	active site interaction map of 1D3H_DU_0.oedu

Related Python Cookbook Recipes

• Visualizing Protein-Ligand B-factor

• Visualizing Protein-Ligand B-factor Map

• 🔄 Visualizing Protein-Ligand Contacts

• 🔄 Visualizing Protein-Ligand Interactions

• Visualizing Protein-Ligand Maps

See also in OEChem TK manual

Theory

• Biopolymers chapter

• Protein Preparation chapter

API

• OEInteractionHintContainer class

• OEPerceiveInteractionHints function

• OEPlaceHydrogens function

See also in OEDepict TK manual

Theory

• Molecule Depiction chapter

API

• OEImage class

See also in Grapheme^TM TK manual

API

• OE2DActiveSiteDisplay class

• OEDrawUnpairedInteractionMapLegend function

• OEPrepareActiveSiteDepiction function

• OERenderUnpairedInteractionMap function