Visualizing Protein-Ligand B-factor Map

Problem

You want to visualize the B-factor of an active site in order to reveal regions with high flexibility. See example in Figure 1. Atoms with high B-factor (i.e. high flexibility) are colored red, while blue color indicates low B-factors.

hover the mouse over any residue circle to display its B-factor

Figure 1. Example of depicting the B-factor map of 2A1B complex

See also

• Visualizing Protein-Ligand B-factor recipe that shows an alternative visualization

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

bfactormap2img.py

See also Usage subsection.

Source Code

bfactormap2img

#!/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 B-factor of a ligand and its environment."""

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 B-factor of an active site."
__SCRIPT_TOOLKITS__ = ["oechem", "oedepict", "oegrapheme"]
__SCRIPT_CATEGORIES__ = ["visualization"]


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=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)",
    )
    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 B-factor."""
    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 active site with b-bfactor

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

    opts = oegrapheme.OE2DActiveSiteDisplayOptions(args.width, args.height)
    opts.SetRenderInteractiveLegend(args.interactive_legend)

    depict_bfactor_map(image, protein, ligand, opts)

    if args.image and Path(args.image).suffix[1:].lower() == "svg":
        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_bfactor_map(
    image: oedepict.OEImageBase,
    protein: oechem.OEMolBase,
    ligand: oechem.OEMolBase,
    opts: oegrapheme.OE2DActiveSiteDisplayOptions,
) -> None:
    """Depict B-factor map of active site."""
    # 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, opts)
    oegrapheme.OERenderBFactorMap(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 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_bfactormap 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. When generating the B-factor map the interactions are only used to position the nearby residues around the depicted ligand.

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 OERenderBFactorMap function generates the interactive image in which the B-factor is projected into the ligand using a color gradient (red indicates atoms with high B-factor while blue color indicates low B-factors)

   The residue circles in the image act as hover buttons that reveal the atomic representation of the residue, color-coded by the B-factor. The color of the residue circle indicates the average B-factor of the residue.

def depict_bfactor_map(
    image: oedepict.OEImageBase,
    protein: oechem.OEMolBase,
    ligand: oechem.OEMolBase,
    opts: oegrapheme.OE2DActiveSiteDisplayOptions,
) -> None:
    """Depict B-factor map of active site."""
    # 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, opts)
    oegrapheme.OERenderBFactorMap(image, active_site_disp)

Usage

See Download section to download the script.

> bfactormap2img --help

The following commands will generate the image shown in Figure 1 for the input file 1A1B.pdb

> bfactormap2img --complex 1A1B.pdb --interactive-legend --image image.svg

Discussion

See Discussion subsection of the Visualizing Protein-Ligand B-factor recipe.

Related Python Cookbook Recipes

• Visualizing Protein-Ligand B-factor

• Visualizing Protein-Ligand B-factor Heat Map

• 🔄 Visualizing Protein-Ligand Interactions

• 🔄 Visualizing Protein-Ligand Contacts

• Visualizing Protein-Ligand Unpaired Interactions

• Visualizing Protein-Ligand Maps

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

• OERenderBFactorMap function