Visualizing Protein-Ligand B-factor

Problem

You want to visualize the B-factor of a ligand and its environment in order to reveal regions with high flexibility at the active site. See example in Figure 1.

hover the mouse over any ligand atom to display its B-factor

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

Ligand atoms with high B-factor (i.e. high flexibility) are colored red, while blue color indicates low B-factors.

The 2D molecule surface (arcs around the ligand) represents the environment of the ligand in 3D. Each arc is associated with a ligand atom that is the center of the arc. A red arc indicates that the average B-factor of the nearby protein atoms is high i.e. the environment is more flexible. When this average B-factor is calculated, only protein atoms that are closer than 4.0 Ångströms (user-adjustable parameter) to the given ligand atom are considered. Blue arcs represent low average B-factor values of the nearby protein atoms. The lack of arc indicates that there is no protein heavy atom closer than 4.0 Ångströms to the corresponding ligand atom.

The linear color gradient, used for coloring B-factors, is drawn at the bottom of the image. The number (that is 18.00) marked on the color gradient is the average B-factor of the whole complex. The black box shown on the color gradient indicates the range of B-factor values close to the active site i.e. the minimum and maximum B-factor values of any protein atom that is closer than 4.0 Ångströms to any ligand atoms.

See also

• Visualizing Protein-Ligand B-factor Map and Visualizing Protein-Ligand B-factor Heat Map that show alternative visualizations

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

bfactor2img.py

See also the Usage subsection.

Source Code

bfactor2img

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

import rich.console
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 a ligand and its environment."
__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]",
    )

    # 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 options
    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)",
    )
    viz_group = parser.add_argument_group("Visualization options")
    viz_group.add_argument(
        "--max-dist",
        type=float,
        default=4.0,
        help="maximum distance of receptor atoms to be considered (default: %(default)s)",
    )
    parser.add_argument("--help-image", action=HelpPreviewAction)

    return parser.parse_args()


def main() -> int:
    """Depict B-factor."""
    args = parse_options()

    _check_image_file(args)
    is_svg: bool = args.image and Path(args.image).suffix[1:].upper() == "SVG"

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

    # calculate average BFactor of the whole complex
    avg_bfactor = get_average_bfactor(protein, ligand)

    console = rich.console.Console()
    console.print(f"Average B-factor in complex = {avg_bfactor:.2f}")

    # calculate minimum and maximum BFactor of the ligand and its environment
    min_bfactor, max_bfactor = get_min_and_max_bfactor(protein, ligand, args.max_dist)
    console.print(
        f"B-factor of ligand and its environment in {args.max_dist:.1f} Å in range [{min_bfactor:.2f}-{max_bfactor:.2f}]"
    )

    # attach to each ligand atom the average BFactor of the nearby protein atoms
    tag: int = oechem.OEGetTag("avg residue BFfactor")
    set_average_bfactor_of_nearby_protein_atoms(protein, ligand, tag, args.max_dist)

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

    main_frame = oedepict.OEImageFrame(
        image, args.width, args.height * 0.85, oedepict.OE2DPoint(0.0, 0.0)
    )
    legend_frame = oedepict.OEImageFrame(
        image,
        args.width,
        args.height * 0.15,
        oedepict.OE2DPoint(0.0, args.height * 0.85),
    )

    color_gradient = get_bfactor_color_gradient()

    opts = oedepict.OE2DMolDisplayOptions(
        main_frame.GetWidth(), main_frame.GetHeight(), oedepict.OEScale_AutoScale
    )
    opts.SetAtomColorStyle(oedepict.OEAtomColorStyle_WhiteMonochrome)
    opts.SetTitleLocation(oedepict.OETitleLocation_Hidden)

    depict_bfactor(main_frame, ligand, opts, color_gradient, tag, is_svg)
    depict_color_gradient(
        legend_frame, color_gradient, min_bfactor, max_bfactor, avg_bfactor
    )

    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(
    image: oedepict.OEImageBase,
    ligand: oechem.OEMolBase,
    opts: oedepict.OE2DMolDisplayOptions,
    color_gradient: oechem.OEColorGradientBase,
    tag: int,
    is_svg: bool,
) -> None:
    """Depict B-factor."""
    # prepare ligand for depiction

    oegrapheme.OEPrepareDepictionFrom3D(ligand)

    clear_coords, suppress_hydrogens = False, False
    prep_opts = oedepict.OEPrepareDepictionOptions(clear_coords, suppress_hydrogens)
    prep_opts.SetDepictOrientation(oedepict.OEDepictOrientation_Horizontal)
    oedepict.OEPrepareDepiction(ligand, prep_opts)

    arc_fxn = BFactorArcFxn(color_gradient, tag)
    for atom in ligand.GetAtoms():
        oegrapheme.OESetSurfaceArcFxn(ligand, atom, arc_fxn)
    opts.SetScale(oegrapheme.OEGetMoleculeSurfaceScale(ligand, opts))

    # render ligand and visualize BFactor

    disp = oedepict.OE2DMolDisplay(ligand, opts)

    if is_svg:
        font = oedepict.OEFont(
            oedepict.OEFontFamily_Default,
            oedepict.OEFontStyle_Default,
            14,
            oedepict.OEAlignment_Center,
            oechem.OEBlack,
        )
        for atom_disp in disp.GetAtomDisplays():
            atom = atom_disp.GetAtom()
            if not oechem.OEHasResidue(atom):
                continue
            res = oechem.OEAtomGetResidue(atom)
            hover_text = f"bfactor={res.GetBFactor():.2f}"
            oedepict.OEDrawSVGHoverText(disp, atom_disp, hover_text, font)

    color_bfactor = ColorLigandAtomByBFactor(color_gradient)
    oegrapheme.OEAddGlyph(disp, color_bfactor, oechem.OEIsTrueAtom())

    oegrapheme.OEDraw2DSurface(disp)

    oedepict.OERenderMolecule(image, disp)


def depict_color_gradient(
    image: oedepict.OEImageBase,
    color_gradient: oechem.OEColorGradientBase,
    min_bfactor: float,
    max_bfactor: float,
    avg_bfactor: float,
) -> None:
    """Depicts color gradient."""
    opts = oegrapheme.OEColorGradientDisplayOptions()
    opts.SetColorStopPrecision(1)
    opts.AddMarkedValue(avg_bfactor)
    opts.SetBoxRange(min_bfactor, max_bfactor)

    oegrapheme.OEDrawColorGradient(image, color_gradient, opts)


def get_average_bfactor(protein: oechem.OEMolBase, ligand: oechem.OEMolBase) -> float:
    """Calculate the average b-factor for the all ligand and protein atoms."""
    num_atoms, sum_bfactor = 0, 0.0
    for mol in [protein, ligand]:
        for atom in mol.GetAtoms():
            if not oechem.OEHasResidue(atom):
                continue
            res = oechem.OEAtomGetResidue(atom)
            sum_bfactor += res.GetBFactor()
            num_atoms += 1
    return sum_bfactor / num_atoms


class NotHydrogenOrWater(oechem.OEUnaryAtomPred):
    """Predicate used to identify heady atoms but ignore OH2."""

    def __call__(self, atom: oechem.OEAtomBase) -> bool:
        """Evaluate atom."""
        if atom.GetAtomicNum() == oechem.OEElemNo_H:
            return False
        if not oechem.OEHasResidue(atom):
            return False

        water_pred = oechem.OEIsWater()
        return not water_pred(atom)


def get_min_and_max_bfactor(
    protein: oechem.OEMolBase, ligand: oechem.OEMolBase, max_distance: float
) -> tuple[float, float]:
    """Calculate the range of the b-factor."""
    min_bfactor, max_bfactor = float("inf"), float("-inf")

    # ligand atoms

    for atom in ligand.GetAtoms(oechem.OEIsHeavy()):
        if not oechem.OEHasResidue(atom):
            continue
        res = oechem.OEAtomGetResidue(atom)
        min_bfactor = min(min_bfactor, res.GetBFactor())
        max_bfactor = max(max_bfactor, res.GetBFactor())

    # protein atoms close to ligand atoms
    consider_bfactor = NotHydrogenOrWater()

    nn = oechem.OENearestNbrs(protein, max_distance)
    for lig_atom in ligand.GetAtoms(oechem.OEIsHeavy()):
        for neigh in nn.GetNbrs(lig_atom):
            prot_atom = neigh.GetBgn()

            if consider_bfactor(prot_atom):
                res = oechem.OEAtomGetResidue(prot_atom)
                min_bfactor = min(min_bfactor, res.GetBFactor())
                max_bfactor = max(max_bfactor, res.GetBFactor())

    return min_bfactor, max_bfactor


def set_average_bfactor_of_nearby_protein_atoms(
    protein: oechem.OEMolBase, ligand: oechem.OEMolBase, tag: int, max_distance: float
) -> None:
    """Set average b-factor on protein atoms close to ligand."""
    consider_bfactor = NotHydrogenOrWater()
    nn = oechem.OENearestNbrs(protein, max_distance)
    for ligand_atom in ligand.GetAtoms(oechem.OEIsHeavy()):
        sum_bfactor = 0.0
        neighs = []
        for neigh in nn.GetNbrs(ligand_atom):
            pro_atom = neigh.GetBgn()
            if consider_bfactor(pro_atom):
                res = oechem.OEAtomGetResidue(pro_atom)
                sum_bfactor += res.GetBFactor()
                neighs.append(pro_atom)

        avg_bfactor = 0.0
        if len(neighs) > 0:
            avg_bfactor = sum_bfactor / len(neighs)
        ligand_atom.SetDoubleData(tag, avg_bfactor)


def get_bfactor_color_gradient() -> oechem.OELinearColorGradient:
    """Initialise color gradient used to visualize b-factor values."""
    color_gradient = oechem.OELinearColorGradient()
    color_gradient.AddStop(oechem.OEColorStop(0.0, oechem.OEDarkBlue))
    color_gradient.AddStop(oechem.OEColorStop(10.0, oechem.OELightBlue))
    color_gradient.AddStop(oechem.OEColorStop(25.0, oechem.OEYellowTint))
    color_gradient.AddStop(oechem.OEColorStop(50.0, oechem.OERed))
    color_gradient.AddStop(oechem.OEColorStop(100.0, oechem.OEDarkRose))
    return color_gradient


class BFactorArcFxn(oegrapheme.OESurfaceArcFxnBase):
    """Surface drawer around ligand."""

    def __init__(self, color_gradient: oechem.OEColorGradientBase, tag: int) -> None:
        """Initialize."""
        oegrapheme.OESurfaceArcFxnBase.__init__(self)
        self._color_gradient = color_gradient
        self._tag = tag

    def __call__(
        self, image: oedepict.OEImageBase, arc: oegrapheme.OESurfaceArc
    ) -> bool:
        """Draw arc."""
        atom_disp = arc.GetAtomDisplay()
        if atom_disp is None or not atom_disp.IsVisible():
            return False

        atom = atom_disp.GetAtom()
        if atom is None:
            return False

        avg_residue_bfactor = atom.GetDoubleData(self._tag)
        if avg_residue_bfactor == 0.0:
            return True
        color = self._color_gradient.GetColorAt(avg_residue_bfactor)

        pen = oedepict.OEPen(color, color, oedepict.OEFill_Off, 5.0)

        center = arc.GetCenter()
        bgn_angle, end_angle = arc.GetBgnAngle(), arc.GetEndAngle()
        radius = arc.GetRadius()

        oegrapheme.OEDrawDefaultSurfaceArc(
            image, center, bgn_angle, end_angle, radius, pen
        )

        return True

    def CreateCopy(self):  # noqa: ANN201, N802
        """Copy constructor."""
        return BFactorArcFxn(self._color_gradient, self._tag).__disown__()


class ColorLigandAtomByBFactor(oegrapheme.OEAtomGlyphBase):
    """Class used to color ligand atoms based on their b-factor."""

    def __init__(self, color_gradient: oechem.OEColorGradientBase) -> None:
        """Initialize."""
        oegrapheme.OEAtomGlyphBase.__init__(self)
        self._color_gradient = color_gradient

    def RenderGlyph(  # noqa: N802
        self, disp: oedepict.OE2DMolDisplay, atom: oechem.OEAtomBase
    ) -> bool:
        """Highlight atom."""
        atom_disp = disp.GetAtomDisplay(atom)
        if atom_disp is None or not atom_disp.IsVisible():
            return False

        if not oechem.OEHasResidue(atom):
            return False

        res = oechem.OEAtomGetResidue(atom)
        bfactor = res.GetBFactor()
        color = self._color_gradient.GetColorAt(bfactor)

        pen = oedepict.OEPen(color, color, oedepict.OEFill_On, 1.0)
        radius = disp.GetScale() / 3.0

        layer = disp.GetLayer(oedepict.OELayerPosition_Below)
        circle_style = oegrapheme.OECircleStyle_Default
        oegrapheme.OEDrawCircle(layer, circle_style, atom_disp.GetCoords(), radius, pen)
        return True

    def CreateCopy(self):  # noqa: ANN201, N802
        """Copy constructor."""
        return ColorLigandAtomByBFactor(self._color_gradient).__disown__()


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 get_average_bfactor function is used to calculate the average B-factor for the whole complex. When a protein-ligand complex is read from a PDB file, an OEResidue object is attached to each atom as generic data. This OEResidue object stores a number of fields specific to processing macromolecules including the B-factor.

def get_average_bfactor(protein: oechem.OEMolBase, ligand: oechem.OEMolBase) -> float:
    """Calculate the average b-factor for the all ligand and protein atoms."""
    num_atoms, sum_bfactor = 0, 0.0
    for mol in [protein, ligand]:
        for atom in mol.GetAtoms():
            if not oechem.OEHasResidue(atom):
                continue
            res = oechem.OEAtomGetResidue(atom)
            sum_bfactor += res.GetBFactor()
            num_atoms += 1
    return sum_bfactor / num_atoms

The get_min_and_max_bfactor function is utilized to calculate the minimum and maximum B-factors at the region of the active site. Therefore, only protein atoms that are close to the ligand are considered (as defined by the maxdistance parameter). These nearby protein atoms are identified by using the OENearestNbrs class. The user defined atom predicate, NotHydrogenOrWater, is also used to ignore any hydrogen atoms or water molecules.

def get_min_and_max_bfactor(
    protein: oechem.OEMolBase, ligand: oechem.OEMolBase, max_distance: float
) -> tuple[float, float]:
    """Calculate the range of the b-factor."""
    min_bfactor, max_bfactor = float("inf"), float("-inf")

    # ligand atoms

    for atom in ligand.GetAtoms(oechem.OEIsHeavy()):
        if not oechem.OEHasResidue(atom):
            continue
        res = oechem.OEAtomGetResidue(atom)
        min_bfactor = min(min_bfactor, res.GetBFactor())
        max_bfactor = max(max_bfactor, res.GetBFactor())

    # protein atoms close to ligand atoms
    consider_bfactor = NotHydrogenOrWater()

    nn = oechem.OENearestNbrs(protein, max_distance)
    for lig_atom in ligand.GetAtoms(oechem.OEIsHeavy()):
        for neigh in nn.GetNbrs(lig_atom):
            prot_atom = neigh.GetBgn()

            if consider_bfactor(prot_atom):
                res = oechem.OEAtomGetResidue(prot_atom)
                min_bfactor = min(min_bfactor, res.GetBFactor())
                max_bfactor = max(max_bfactor, res.GetBFactor())

    return min_bfactor, max_bfactor

The NotHydrogenOrWater atom predicate returns true if the given atom is not a hydrogen and not part of a water molecule.

class NotHydrogenOrWater(oechem.OEUnaryAtomPred):
    """Predicate used to identify heady atoms but ignore OH2."""

    def __call__(self, atom: oechem.OEAtomBase) -> bool:
        """Evaluate atom."""
        if atom.GetAtomicNum() == oechem.OEElemNo_H:
            return False
        if not oechem.OEHasResidue(atom):
            return False

        water_pred = oechem.OEIsWater()
        return not water_pred(atom)

Grapheme TK’s 2D molecule surface depiction style is used to visualize the B-factor of the protein around the ligand. The set_average_bfactor_of_nearby_protein_atoms function below loops over the ligand atoms and calculates the average B-factor of nearby protein atoms. This average B-factor will later be used to determine the color of the arc drawn by the BFactorArcFxn class.

def set_average_bfactor_of_nearby_protein_atoms(
    protein: oechem.OEMolBase, ligand: oechem.OEMolBase, tag: int, max_distance: float
) -> None:
    """Set average b-factor on protein atoms close to ligand."""
    consider_bfactor = NotHydrogenOrWater()
    nn = oechem.OENearestNbrs(protein, max_distance)
    for ligand_atom in ligand.GetAtoms(oechem.OEIsHeavy()):
        sum_bfactor = 0.0
        neighs = []
        for neigh in nn.GetNbrs(ligand_atom):
            pro_atom = neigh.GetBgn()
            if consider_bfactor(pro_atom):
                res = oechem.OEAtomGetResidue(pro_atom)
                sum_bfactor += res.GetBFactor()
                neighs.append(pro_atom)

        avg_bfactor = 0.0
        if len(neighs) > 0:
            avg_bfactor = sum_bfactor / len(neighs)
        ligand_atom.SetDoubleData(tag, avg_bfactor)

The linear color gradient used in the depiction is initialized by the get_bfactor_color_gradient function. The color gradient does not depend on the depicted complex, thereby allowing easy comparison of B-factor images.

def get_bfactor_color_gradient() -> oechem.OELinearColorGradient:
    """Initialise color gradient used to visualize b-factor values."""
    color_gradient = oechem.OELinearColorGradient()
    color_gradient.AddStop(oechem.OEColorStop(0.0, oechem.OEDarkBlue))
    color_gradient.AddStop(oechem.OEColorStop(10.0, oechem.OELightBlue))
    color_gradient.AddStop(oechem.OEColorStop(25.0, oechem.OEYellowTint))
    color_gradient.AddStop(oechem.OEColorStop(50.0, oechem.OERed))
    color_gradient.AddStop(oechem.OEColorStop(100.0, oechem.OEDarkRose))
    return color_gradient

The ColorLigandAtomByBFactor class is used to draw the colored circles below the ligand atoms indicating their B-factor. The color gradient used in this class is generated by the get_bfactor_color_gradient function.

class ColorLigandAtomByBFactor(oegrapheme.OEAtomGlyphBase):
    """Class used to color ligand atoms based on their b-factor."""

    def __init__(self, color_gradient: oechem.OEColorGradientBase) -> None:
        """Initialize."""
        oegrapheme.OEAtomGlyphBase.__init__(self)
        self._color_gradient = color_gradient

    def RenderGlyph(  # noqa: N802
        self, disp: oedepict.OE2DMolDisplay, atom: oechem.OEAtomBase
    ) -> bool:
        """Highlight atom."""
        atom_disp = disp.GetAtomDisplay(atom)
        if atom_disp is None or not atom_disp.IsVisible():
            return False

        if not oechem.OEHasResidue(atom):
            return False

        res = oechem.OEAtomGetResidue(atom)
        bfactor = res.GetBFactor()
        color = self._color_gradient.GetColorAt(bfactor)

        pen = oedepict.OEPen(color, color, oedepict.OEFill_On, 1.0)
        radius = disp.GetScale() / 3.0

        layer = disp.GetLayer(oedepict.OELayerPosition_Below)
        circle_style = oegrapheme.OECircleStyle_Default
        oegrapheme.OEDrawCircle(layer, circle_style, atom_disp.GetCoords(), radius, pen)
        return True

    def CreateCopy(self):  # noqa: ANN201, N802
        """Copy constructor."""
        return ColorLigandAtomByBFactor(self._color_gradient).__disown__()

The BFactorArcFxn class draws the arc segments around the depicted ligand. The color of each arc represents the average B-factor of the protein atoms around the corresponding ligand atom. These values are calculated by the set_average_bfactor_of_nearby_protein_atoms function and attached to the ligand atoms as generic data. The color gradient used in this class is generated by the get_bfactor_color_gradient function.

class BFactorArcFxn(oegrapheme.OESurfaceArcFxnBase):
    """Surface drawer around ligand."""

    def __init__(self, color_gradient: oechem.OEColorGradientBase, tag: int) -> None:
        """Initialize."""
        oegrapheme.OESurfaceArcFxnBase.__init__(self)
        self._color_gradient = color_gradient
        self._tag = tag

    def __call__(
        self, image: oedepict.OEImageBase, arc: oegrapheme.OESurfaceArc
    ) -> bool:
        """Draw arc."""
        atom_disp = arc.GetAtomDisplay()
        if atom_disp is None or not atom_disp.IsVisible():
            return False

        atom = atom_disp.GetAtom()
        if atom is None:
            return False

        avg_residue_bfactor = atom.GetDoubleData(self._tag)
        if avg_residue_bfactor == 0.0:
            return True
        color = self._color_gradient.GetColorAt(avg_residue_bfactor)

        pen = oedepict.OEPen(color, color, oedepict.OEFill_Off, 5.0)

        center = arc.GetCenter()
        bgn_angle, end_angle = arc.GetBgnAngle(), arc.GetEndAngle()
        radius = arc.GetRadius()

        oegrapheme.OEDrawDefaultSurfaceArc(
            image, center, bgn_angle, end_angle, radius, pen
        )

        return True

    def CreateCopy(self):  # noqa: ANN201, N802
        """Copy constructor."""
        return BFactorArcFxn(self._color_gradient, self._tag).__disown__()

depict_bfactor is the main function that performs the depiction of the ligand and draws the molecule surface around it to represent its environment.

1. First, the depiction coordinates of the ligand are generated by calling the OEPrepareDepictionFrom3D function. Then the OEPrepareDepictionOptions class and OEPrepareDepiction function of the OEDepict TK are utilized to orient the generated 2D coordinates horizontally.

2. The BFactorArcFxn arc drawing functor is added to each ligand (see also step 7)

3. Scaling of the depicted ligand is reduced by using the OEGetMoleculeSurfaceScale function to make sure that arcs representing the molecule surface of the ligand will not be drawn outside the image.

4. The molecule display for the ligand is constructed.

5. If the output image file format is svg, then the OEDrawSVGHoverText function can be used to associate each display atom with a text of its B-factor value. The generated svg image will be interactive i.e. the B-factor texts will be displayed when the mouse is hovered over the ligand atoms in the image.

6. The OEAddGlyph function is used to draw a glyph (in this case a circle colored by the B-factor) underneath each ligand atom as defined in the ColorLigandAtomByBFactor class.

7. The OEDraw2DSurface function is called that will invoke the arc drawing functor set on each ligand atom in step 2. resulting in the rendering of the molecule into the display.

8. Finally, the OERenderMolecule function renders the molecule display into the image).

def depict_bfactor(
    image: oedepict.OEImageBase,
    ligand: oechem.OEMolBase,
    opts: oedepict.OE2DMolDisplayOptions,
    color_gradient: oechem.OEColorGradientBase,
    tag: int,
    is_svg: bool,
) -> None:
    """Depict B-factor."""
    # prepare ligand for depiction

    oegrapheme.OEPrepareDepictionFrom3D(ligand)

    clear_coords, suppress_hydrogens = False, False
    prep_opts = oedepict.OEPrepareDepictionOptions(clear_coords, suppress_hydrogens)
    prep_opts.SetDepictOrientation(oedepict.OEDepictOrientation_Horizontal)
    oedepict.OEPrepareDepiction(ligand, prep_opts)

    arc_fxn = BFactorArcFxn(color_gradient, tag)
    for atom in ligand.GetAtoms():
        oegrapheme.OESetSurfaceArcFxn(ligand, atom, arc_fxn)
    opts.SetScale(oegrapheme.OEGetMoleculeSurfaceScale(ligand, opts))

    # render ligand and visualize BFactor

    disp = oedepict.OE2DMolDisplay(ligand, opts)

    if is_svg:
        font = oedepict.OEFont(
            oedepict.OEFontFamily_Default,
            oedepict.OEFontStyle_Default,
            14,
            oedepict.OEAlignment_Center,
            oechem.OEBlack,
        )
        for atom_disp in disp.GetAtomDisplays():
            atom = atom_disp.GetAtom()
            if not oechem.OEHasResidue(atom):
                continue
            res = oechem.OEAtomGetResidue(atom)
            hover_text = f"bfactor={res.GetBFactor():.2f}"
            oedepict.OEDrawSVGHoverText(disp, atom_disp, hover_text, font)

    color_bfactor = ColorLigandAtomByBFactor(color_gradient)
    oegrapheme.OEAddGlyph(disp, color_bfactor, oechem.OEIsTrueAtom())

    oegrapheme.OEDraw2DSurface(disp)

    oedepict.OERenderMolecule(image, disp)

Usage

See Download section to download the script.

> bfactor2img --help

Visualizing 1A1B.pdb complex shown in Figure 1

> bfactor2img --complex 1A1B.pdb --image image.svg

Discussion

The B-factor (temperature factor) is a measure of how much an atom vibrates around the position specified in the PDB file. Atoms at side-chain termini are expected to exhibit more freedom of movement than main-chain atoms. Similarly, ligand atoms exposed to solvent can have more freedom of movement (see 1A4K example below).

B-factors can indicate not only the mobility of the atoms but they can also reveal where there are errors in model building. Visualization helps to make an instant judgment on the quality of the model.

PDB: 1A4K	PDB: 1ATL

Related Python Cookbook Recipes

• Visualizing Protein-Ligand B-factor Map

• Visualizing Protein-Ligand B-factor Heat Map

• 🔄 Visualizing Protein-Ligand Contacts

• 🔄 Visualizing Protein-Ligand Interactions

• Visualizing Protein-Ligand Unpaired Interactions

• Visualizing Protein-Ligand Maps

See also in OEChem TK manual

Theory

• Generic Data chapter

• Biopolymers chapter

API

• OEAtomGetResidue function

• OEGetResidueIndex function

• OEHasResidue function

• OELinearColorGradient class

• OENearestNbrs class

• OEResidue class

See also in OEDepict TK manual

Theory

• Molecule Depiction chapter

• Generating Interactive SVG images chapter

API

• OE2DMolDisplay class

• OE2DMolDisplayOptions class

• OEDrawSVGHoverText function

• OEFont class

• OEImage class

• OEPrepareDepictionOptions class

• OEPrepareDepiction function

• OERenderMolecule function

See also in Grapheme^TM TK manual

Theory

• Annotating Atoms and Bonds chapter

• Drawing a Molecule Surface chapter

API

• OEAddGlyph function

• OEAtomGlyphBase base class

• OEDraw2DSurface function

• OEDrawCircle function

• OEDrawDefaultSurfaceArc function

• OEGetMoleculeSurfaceScale function

• OEPrepareDepictionFrom3D function

• OESetSurfaceArcFxn function

• OESurfaceArc class

• OESurfaceArcFxnBase base class