Visualizing Protein-Ligand Maps

Problem

You would like to generate one image that contains all of the Grapheme TK maps that represent protein-ligand information. See example in Figure 1.

Figure 1. Example of depicting the various protein-ligand maps of 2A1B 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

activesitemaps2img.py

See also Usage subsection.

Source Code

activesitemaps2img

#!/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 active site maps."""

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 active site maps."
__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=True,
        metavar="IMAGE-FILE",
        help="output image file (SVG)",
    )
    image_group.add_argument(
        "--width",
        type=int,
        default=600,
        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:
    """Depict active site maps."""
    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)
    depict_active_site_maps(image, protein, ligand)

    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_active_site_maps(
    image: oedepict.OEImageBase, protein: oechem.OEMolBase, ligand: oechem.OEMolBase
) -> None:
    """Depict active site maps."""
    # 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)
    oegrapheme.OERenderActiveSiteMaps(image, active_site)


def get_protein_and_ligand_from_pdb(
    pdb_filename: str,
) -> tuple[oechem.OEMolBase, oechem.OEMolBase]:
    """Read protein and ligand 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 ligand 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 an issue
    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_activesite_maps 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 OERenderActiveSiteMaps function generates a multi-tab interactive image.

def depict_active_site_maps(
    image: oedepict.OEImageBase, protein: oechem.OEMolBase, ligand: oechem.OEMolBase
) -> None:
    """Depict active site maps."""
    # 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)
    oegrapheme.OERenderActiveSiteMaps(image, active_site)

Usage

See Download section to download the script.

> activesitemaps2img --help

Visualizing 1YWR_DU_0.oedu design unit of 1YWR.

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

Discussion

OERenderActiveSiteMaps currently visualizes protein-ligand information in four different maps:

• interaction map (using OERenderActiveSite)

• unpaired interaction map (using OERenderUnpairedInteractionMap)

• B-factor map (using OERenderBFactorMap)

• contact map (using OERenderContactMap)

Related Python Cookbook Recipes

• Visualizing Protein-Ligand B-factor

• Visualizing Protein-Ligand B-factor Map

• 🔄 Visualizing Protein-Ligand Interactions

• 🔄 Visualizing Protein-Ligand Contacts

• Visualizing Protein-Ligand Unpaired Interactions

See also in OEChem TK manual

Theory

• Biopolymers chapter

• Protein Preparation chapter

API

• OEInteractionHintContainer class

• OEPerceiveInteractionHints function

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

• OERenderActiveSiteMaps function