Visualizing Protein-Ligand B-factor Heat Map

Problem

You want to visualize the B-factor of an protein-ligand complex in a heat-map style. See example in Figure 1. In the sequence chain view (on the left) each residue is colored by the average B-factors of its atoms. When hovering over a residue in the chain sequence view, the corresponding residue will appear in the residue display window (on the right). The atoms of these displayed residues are colored according to their B-factors.

The utilized color gradient is depicted at the bottom right corner of the image. Residues and atoms with high B-factor (i.e. high flexibility) are colored red, while blue color indicates low B-factors.

hover the mouse over the sequence chain to display residue B-factor

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

../_images/bfactorheatmap2img-pdb1a1b.svg

Ingredients

  • OEChem TK - cheminformatics toolkit (including OEBio TK)

  • OEDepict TK - molecule depiction toolkit

  • Grapheme TK - molecule and property visualization toolkit

Difficulty level

../_images/chilly.png ../_images/chilly.png ../_images/chilly.png

Download

Download code

bfactorheatmap2img.py

See also Usage subsection.

Solution

The depict_bfactor_heatmap function shows how to generate the color coded sequence view with hover mode: The depict_bfactor_heatmap function takes the group of residues that are split into equal sizes by the split_residues function. These residues are then visualized in a sequence view:

  1. First the parameters are calculated that determine the layout of the sequence view and the font that is used to render information (lines 14-33).

  2. Then iterating over the pre-chunked residue groups the label for each sequence line is rendered (lines 34-44).

  3. The inside loop iterates over the residues of each group and rendering a rectangle for each residue when the color of the rectangle indicates the average B-factors of the residue atoms (lines 34-44).

    For each residue, two SVG groups are generated to accomplish the hover effect. The rectangle representing the residue in the sequence view is drawn into a group called area group (See also OEAddSVGHover function). When hovering the mouse over objects drawn inside the area group the objects drawn in the second group, called target group, will pop up. In this example, both the residue label box and the residue molecule (rendered by depict_residue) are drawn inside the target group (lines 60-76). Everything that is rendered between the OEImageBase.PushGroup and the corresponding OEImageBase.PopGroup methods is considered “inside” a group.

 1def depict_bfactor_heatmap(image, residue_frame, protein, residue_groups,
 2                           max_residues_per_line, colorg):
 3    """
 4    Depicts the residues and the b-factor heat-map
 5
 6    :type image: oedepict.OEImageBase
 7    :type residue_frame: oedepict.OEImageBase
 8    :type protein: oechem.OEMolBase
 9    :type residue_groups: list[list[oechem.OEHierResidue]]
10    :type max_residues_per_line: int
11    :type colorg: oechem.OELineraColorGradient
12    """
13
14    nr_residue_groups = len(residue_groups)
15
16    linegap = image.GetHeight() * 0.85 / (nr_residue_groups + 1)
17    linegap = min(linegap, 50.0)
18    bgn_line = oedepict.OE2DPoint(80, linegap / 2.0)
19    end_line = oedepict.OE2DPoint(image.GetWidth() - 20, linegap / 2.0)
20    line_offset = oedepict.OE2DPoint(0, linegap)
21
22    linelength = get_distange(bgn_line, end_line)
23
24    res_box_width = linelength / max_residues_per_line
25    res_box_height = min(linegap / 1.5, res_box_width * 4.0)
26
27    res_hgroup_font = oedepict.OEFont(oedepict.OEFontFamily_Default, oedepict.OEFontStyle_Default,
28                                      10, oedepict.OEAlignment_Left, oechem.OEBlack)
29    if nr_residue_groups > 20:
30        res_hgroup_font.SetSize(8)
31
32    res_font = oedepict.OEFont(oedepict.OEFontFamily_Default, oedepict.OEFontStyle_Bold,
33                               9, oedepict.OEAlignment_Center, oechem.OEBlack)
34
35    for gidx, res_group in enumerate(residue_groups):
36
37        image.DrawLine(bgn_line, end_line, oedepict.OELightGreyPen)
38
39        label = get_residue_group_label(res_group)
40        image.DrawText(bgn_line + oedepict.OE2DPoint(-70, res_hgroup_font.GetSize() * 0.33), label, res_hgroup_font)
41
42        box_offset = oedepict.OE2DPoint(res_box_width, 0)
43        tr_box = oedepict.OE2DPoint(bgn_line) - oedepict.OE2DPoint(0, res_box_height / 2.0)
44        bl_box = oedepict.OE2DPoint(bgn_line) + oedepict.OE2DPoint(res_box_width, res_box_height / 2.0)
45
46        for res in res_group:
47
48            bgroup, lgroup = add_residue_svg_groups(image)
49
50            avg_value = get_average_bfactor(res)
51            color = colorg.GetColorAt(avg_value)
52
53            # area group
54
55            image.PushGroup(bgroup)
56            pen = oedepict.OEPen(color, color, oedepict.OEFill_On, 1.0, oedepict.OEStipple_NoLine)
57            image.DrawRectangle(tr_box, bl_box, pen)
58            image.PopGroup(bgroup)
59
60            # target group
61
62            image.PushGroup(lgroup)
63
64            label_box_width, label_box_height = 60, 30
65            label_offset = (tr_box + bl_box) / 2.0
66            label_offset -= oedepict.OE2DPoint(label_box_width / 2.0, label_box_height * 1.2 + res_box_height / 2.0)
67            lframe = oedepict.OEImageFrame(image, label_box_width, label_box_height, label_offset)
68
69            draw_residue_info_label_box(lframe, color)
70
71            res_label = get_residue_label(res.GetOEResidue())
72            center = oedepict.OEGetCenter(lframe)
73            lframe.DrawText(center - oedepict.OE2DPoint(0, 3), res_label, res_font)
74            lframe.DrawText(center + oedepict.OE2DPoint(0, 6), "avg: {:.2f}".format(avg_value), res_font)
75
76            depict_residue(residue_frame, protein, res, res_label, colorg)
77
78            image.PopGroup(lgroup)
79
80            tr_box += box_offset
81            bl_box += box_offset
82
83        bgn_line += line_offset
84        end_line += line_offset
85
86    draw_x_axis(image, bgn_line, end_line, max_residues_per_line)

The depict_residue function shows how to depict residues with color coded atom B-factors:

  1. First residue is extracted from the protein, 2D coordinates are generated and and molecule display is created with monochrome style (lines 12-26).

  2. Then the R-group atoms of the residue are identified and set to be invisible (lines 28-35).

  3. Any bond connected to a R-group atom is then marked with a zigzag line to (lines 37-51).

  4. Finally, circle is drawn underneath each atom of the molecule display with color that indicates the corresponding B-factor (lines 53-65).

 1def depict_residue(image, protein, res, res_label, colorg):
 2    """
 3    Depicts the given residue with b-factor.
 4
 5    :type image: oedepict.OEImageBase
 6    :type protein: oechem.OEMolBase
 7    :type res: oechem.OEHierResidue
 8    :type res_label: string
 9    :type colorg: oechem.OELineraColorGradient
10    """
11
12    resmol = oechem.OEGraphMol()
13    respred = oechem.OEAtomIsInResidue(res.GetOEResidue())
14    adjustHCount, RGroup = False, True
15    oechem.OESubsetMol(resmol, protein, respred, adjustHCount, RGroup)
16    resmol.SetTitle(res_label)
17
18    popts = oedepict.OEPrepareDepictionOptions()
19    popts.SetDepictOrientation(oedepict.OEDepictOrientation_Vertical)
20    oedepict.OEPrepareDepiction(resmol, popts)
21
22    dopts = oedepict.OE2DMolDisplayOptions(image.GetWidth(), image.GetHeight(), oedepict.OEScale_AutoScale)
23    dopts.SetScale(oedepict.OEGetMoleculeScale(resmol, dopts) * 0.9)
24    dopts.SetAtomColorStyle(oedepict.OEAtomColorStyle_WhiteMonochrome)
25    disp = oedepict.OE2DMolDisplay(resmol, dopts)
26
27    rgroups = oechem.OEAtomBondSet()
28    for adisp in disp.GetAtomDisplays():
29        if not adisp.IsVisible():
30            continue
31        atom = adisp.GetAtom()
32        if atom.GetAtomicNum() == 0 and atom.GetMapIdx() > 0:
33            adisp.SetVisible(False)
34            rgroups.AddAtom(atom)
35
36    diagonal = True
37    zigzag_glyph = oegrapheme.OEBondGlyphZigZag(oedepict.OEBlackPen, 0.25,
38                                                oedepict.OELayerPosition_Above, diagonal)
39
40    tpen = oedepict.OEPen(oechem.OEWhite, oechem.OEColor(0, 0, 0, 1), oedepict.OEFill_On, 1.0)
41    for bdisp in disp.GetBondDisplays():
42        bond = bdisp.GetBond()
43        atomB = bond.GetBgn()
44        atomE = bond.GetEnd()
45
46        if rgroups.HasAtom(atomB):
47            bdisp.SetBgnPen(tpen)
48            oegrapheme.OEAddGlyph(disp, zigzag_glyph, oechem.OEHasBondIdx(bond.GetIdx()))
49        if rgroups.HasAtom(atomE):
50            bdisp.SetEndPen(tpen)
51            oegrapheme.OEAddGlyph(disp, zigzag_glyph, oechem.OEHasBondIdx(bond.GetIdx()))
52
53    for adisp in disp.GetAtomDisplays():
54        if not adisp.IsVisible():
55            continue
56        atom = adisp.GetAtom()
57        if not oechem.OEHasResidue(atom):
58            continue
59        res = oechem.OEAtomGetResidue(atom)
60        bfactor = res.GetBFactor()
61        color = colorg.GetColorAt(bfactor)
62        pen = oedepict.OEPen(color, color, oedepict.OEFill_On, 1.0)
63        radius_scale = 1.2
64        circle_glyph = oegrapheme.OEAtomGlyphCircle(pen, oegrapheme.OECircleStyle_Default, radius_scale)
65        oegrapheme.OEAddGlyph(disp, circle_glyph, oechem.OEHasAtomIdx(atom.GetIdx()))
66
67    oedepict.OERenderMolecule(image, disp)

The split_residues function iterates over the chains of the OEHierView and splits the residues into equal sizes.

 1def split_residues(hierview, nr_residues):
 2    """
 3    Splits the residues of the chains of the protein-ligand complex
 4    into equal sizes.
 5
 6    :type hierview: oechem.OEHierView
 7    :type nr_residues: int
 8    """
 9    residue_groups = []
10    for chain in hierview.GetChains():
11        for frag in chain.GetFragments():
12            residues = [res for res in frag.GetResidues()]
13            residue_slices = [residues[i:i+nr_residues] for i in range(0, len(residues), nr_residues)]
14            for r in residue_slices:
15                residue_groups.append(r)
16    return residue_groups

Usage

Usage

bfactorheatmap2img.py

The following commands will generate the image shown in Figure 1:

prompt > wget https://files.rcsb.org/download/1a1b.pdb
prompt > python3 bfactorheatmap2img.py -complex 1a1b.pdb -out 1a1b.svg

Command Line Parameters

Simple parameter list
    input/output options
      -complex : Input filename of the protein-ligand complex
      -out : Output filename of the generated image

    residue options
      -ignore-isolated : Removes isolated atoms prior to depiction
      -ignore-water : Removes water molecule prior to depiction
      -max-residues : Controls the maximum number of residues depicted in each
                      line

Discussion

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

See also in OEChem TK manual

Theory

API

See also in OEDepict TK manual

Theory

API

See also in GraphemeTM TK manual

API