🔄 Visualizing Protein-Ligand Contacts

Problem

You want to generate an interactive image (in svg file format) that depicts the active site i.e. the ligand and the nearby residues. The interactive image allows revealing or hiding the atoms that interact with a given residue. See example in Figure 1.

Click on any residue to reveal which part of the ligand it interacts with

Figure 1. Example of depicting the contact map of 1A1B complex

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

contactmap2img.py

See also the Usage subsection.

Source Code

contactmap2img

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

"""Depicts active site with residues interacting with the ligand."""

import argparse
import os
import sys
from pathlib import Path

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

__SCRIPT_NAME__ = Path(__file__).absolute().stem
__SCRIPT_DESC__ = "Depict contact map of an active site."
__SCRIPT_TOOLKITS__ = ["oechem", "oedepict", "oegrapheme"]
__SCRIPT_CATEGORIES__ = ["visualization", "protein-ligand 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]",
    )

    # 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=True,
        metavar="IMAGE-FILE",
        help="output image file (SVG)",
    )
    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)",
    )

    parser.add_argument("--help-image", action=HelpPreviewAction)
    return parser.parse_args()


def main() -> int:
    """Depict the contact map 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!")

    image = oedepict.OEImage(args.width, args.height)
    opts = oegrapheme.OE2DActiveSiteDisplayOptions(args.width, args.height)

    depict_contact_map(image, protein, ligand, opts)

    icon_scale = 0.5
    oedepict.OEAddInteractiveIcon(image, oedepict.OEIconLocation_TopRight, icon_scale)
    oedepict.OEDrawCurvedBorder(image, oedepict.OELightGreyPen, 10.0)
    oedepict.OEWriteImage(args.image, image)

    return os.EX_OK


def depict_contact_map(
    image: oedepict.OEImageBase,
    protein: oechem.OEMolBase,
    ligand: oechem.OEMolBase,
    depict_options: oegrapheme.OE2DActiveSiteDisplayOptions,
) -> None:
    """Depict protein-ligand contact 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.OERenderContactMap(image, active_site_disp)


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 Path(args.image).suffix[1:].upper() != "SVG":
        oechem.OEThrow.Fatal(
            "This script only accepts SVG as the output image file format!"
        )

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


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

1. First the 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 OERenderContactMap function generates an image in which each residue cycle behaves as a toggle button. By clicking on a residue cycle, the ligand atoms with which the given residue interacts are revealed or hidden. The generated image must be written out in svg image file format, which supports interactive elements.

def depict_contact_map(
    image: oedepict.OEImageBase,
    protein: oechem.OEMolBase,
    ligand: oechem.OEMolBase,
    depict_options: oegrapheme.OE2DActiveSiteDisplayOptions,
) -> None:
    """Depict protein-ligand contact 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.OERenderContactMap(image, active_site_disp)

Usage

See Download section to download the script.

> contactmap2img --help

Visualizing 1YWR_DU_0.oedu design unit of 1YWR.

> contactmap2img --design-unit 1YWR_DU_0.oedu --image image.svg

Discussion

Currently the OEPerceiveInteractionHints function perceives the following interaction types:

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

All of these interaction types are considered when positioning the residues around the ligand.

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.

Related Python Cookbook Recipes

• Visualizing Protein-Ligand B-factor

• Visualizing Protein-Ligand B-factor Map

• 🔄 Visualizing Protein-Ligand Interactions

• Visualizing Protein-Ligand Unpaired Interactions

• Visualizing Protein-Ligand Maps

See also in OEChem TK manual

Theory

• Biopolymers chapter

• Protein Preparation chapter

API

• OEDesignUnit class

• OEInteractionHintContainer class

• OEPerceiveInteractionHints function

See also in Spruce TK manual

• Simple Macromolecular Preparation Workflow chapter

See also in OEDepict TK manual

Theory

• Molecule Depiction chapter

• Generating Interactive SVG images chapter

API

• OEImage class

See also in Grapheme^TM TK manual

API

• OE2DActiveSiteDisplay class

• OEPrepareActiveSiteDepiction function

• OERenderContactMap function