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

Ingredients

  • OEChem TK - cheminformatics toolkit (including OEBio TK)
  • OEDepict TK - molecule depiction toolkit
  • Grapheme TK - molecule and property visualization toolkit

Difficulty level

../_images/chilly7.png ../_images/chilly7.png ../_images/chilly7.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.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
def depict_bfactor_heatmap(image, residue_frame, protein, residue_groups,
                           max_residues_per_line, colorg):
    """
    Depicts the residues and the b-factor heat-map

    :type image: oedepict.OEImageBase
    :type residue_frame: oedepict.OEImageBase
    :type protein: oechem.OEMolBase
    :type residue_groups: list[list[oechem.OEHierResidue]]
    :type max_residues_per_line: int
    :type colorg: oechem.OELineraColorGradient
    """

    nr_residue_groups = len(residue_groups)

    linegap = image.GetHeight() * 0.85 / (nr_residue_groups + 1)
    linegap = min(linegap, 50.0)
    bgn_line = oedepict.OE2DPoint(80, linegap / 2.0)
    end_line = oedepict.OE2DPoint(image.GetWidth() - 20, linegap / 2.0)
    line_offset = oedepict.OE2DPoint(0, linegap)

    linelength = get_distange(bgn_line, end_line)

    res_box_width = linelength / max_residues_per_line
    res_box_height = min(linegap / 1.5, res_box_width * 4.0)

    res_hgroup_font = oedepict.OEFont(oedepict.OEFontFamily_Default, oedepict.OEFontStyle_Default,
                                      10, oedepict.OEAlignment_Left, oechem.OEBlack)
    if nr_residue_groups > 20:
        res_hgroup_font.SetSize(8)

    res_font = oedepict.OEFont(oedepict.OEFontFamily_Default, oedepict.OEFontStyle_Bold,
                               9, oedepict.OEAlignment_Center, oechem.OEBlack)

    for gidx, res_group in enumerate(residue_groups):

        image.DrawLine(bgn_line, end_line, oedepict.OELightGreyPen)

        label = get_residue_group_label(res_group)
        image.DrawText(bgn_line + oedepict.OE2DPoint(-70, res_hgroup_font.GetSize() * 0.33), label, res_hgroup_font)

        box_offset = oedepict.OE2DPoint(res_box_width, 0)
        tr_box = oedepict.OE2DPoint(bgn_line) - oedepict.OE2DPoint(0, res_box_height / 2.0)
        bl_box = oedepict.OE2DPoint(bgn_line) + oedepict.OE2DPoint(res_box_width, res_box_height / 2.0)

        for res in res_group:

            bgroup, lgroup = add_residue_svg_groups(image)

            avg_value = get_average_bfactor(res)
            color = colorg.GetColorAt(avg_value)

            # area group

            image.PushGroup(bgroup)
            pen = oedepict.OEPen(color, color, oedepict.OEFill_On, 1.0, oedepict.OEStipple_NoLine)
            image.DrawRectangle(tr_box, bl_box, pen)
            image.PopGroup(bgroup)

            # target group

            image.PushGroup(lgroup)

            label_box_width, label_box_height = 60, 30
            label_offset = (tr_box + bl_box) / 2.0
            label_offset -= oedepict.OE2DPoint(label_box_width / 2.0, label_box_height * 1.2 + res_box_height / 2.0)
            lframe = oedepict.OEImageFrame(image, label_box_width, label_box_height, label_offset)

            draw_residue_info_label_box(lframe, color)

            res_label = get_residue_label(res.GetOEResidue())
            center = oedepict.OEGetCenter(lframe)
            lframe.DrawText(center - oedepict.OE2DPoint(0, 3), res_label, res_font)
            lframe.DrawText(center + oedepict.OE2DPoint(0, 6), "avg: {:.2f}".format(avg_value), res_font)

            depict_residue(residue_frame, protein, res, res_label, colorg)

            image.PopGroup(lgroup)

            tr_box += box_offset
            bl_box += box_offset

        bgn_line += line_offset
        end_line += line_offset

    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).
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
def depict_residue(image, protein, res, res_label, colorg):
    """
    Depicts the given residue with b-factor.

    :type image: oedepict.OEImageBase
    :type protein: oechem.OEMolBase
    :type res: oechem.OEHierResidue
    :type res_label: string
    :type colorg: oechem.OELineraColorGradient
    """

    resmol = oechem.OEGraphMol()
    respred = oechem.OEAtomIsInResidue(res.GetOEResidue())
    adjustHCount, RGroup = False, True
    oechem.OESubsetMol(resmol, protein, respred, adjustHCount, RGroup)
    resmol.SetTitle(res_label)

    popts = oedepict.OEPrepareDepictionOptions()
    popts.SetDepictOrientation(oedepict.OEDepictOrientation_Vertical)
    oedepict.OEPrepareDepiction(resmol, popts)

    dopts = oedepict.OE2DMolDisplayOptions(image.GetWidth(), image.GetHeight(), oedepict.OEScale_AutoScale)
    dopts.SetScale(oedepict.OEGetMoleculeScale(resmol, dopts) * 0.9)
    dopts.SetAtomColorStyle(oedepict.OEAtomColorStyle_WhiteMonochrome)
    disp = oedepict.OE2DMolDisplay(resmol, dopts)

    rgroups = oechem.OEAtomBondSet()
    for adisp in disp.GetAtomDisplays():
        if not adisp.IsVisible():
            continue
        atom = adisp.GetAtom()
        if atom.GetAtomicNum() == 0 and atom.GetMapIdx() > 0:
            adisp.SetVisible(False)
            rgroups.AddAtom(atom)

    diagonal = True
    zigzag_glyph = oegrapheme.OEBondGlyphZigZag(oedepict.OEBlackPen, 0.25,
                                                oedepict.OELayerPosition_Above, diagonal)

    tpen = oedepict.OEPen(oechem.OEWhite, oechem.OEColor(0, 0, 0, 1), oedepict.OEFill_On, 1.0)
    for bdisp in disp.GetBondDisplays():
        bond = bdisp.GetBond()
        atomB = bond.GetBgn()
        atomE = bond.GetEnd()

        if rgroups.HasAtom(atomB):
            bdisp.SetBgnPen(tpen)
            oegrapheme.OEAddGlyph(disp, zigzag_glyph, oechem.OEHasBondIdx(bond.GetIdx()))
        if rgroups.HasAtom(atomE):
            bdisp.SetEndPen(tpen)
            oegrapheme.OEAddGlyph(disp, zigzag_glyph, oechem.OEHasBondIdx(bond.GetIdx()))

    for adisp in disp.GetAtomDisplays():
        if not adisp.IsVisible():
            continue
        atom = adisp.GetAtom()
        if not oechem.OEHasResidue(atom):
            continue
        res = oechem.OEAtomGetResidue(atom)
        bfactor = res.GetBFactor()
        color = colorg.GetColorAt(bfactor)
        pen = oedepict.OEPen(color, color, oedepict.OEFill_On, 1.0)
        radius_scale = 1.2
        circle_glyph = oegrapheme.OEAtomGlyphCircle(pen, oegrapheme.OECircleStyle_Default, radius_scale)
        oegrapheme.OEAddGlyph(disp, circle_glyph, oechem.OEHasAtomIdx(atom.GetIdx()))

    oedepict.OERenderMolecule(image, disp)

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
def split_residues(hierview, nr_residues):
    """
    Splits the residues of the chains of the protein-ligand complex
    into equal sizes.

    :type hierview: oechem.OEHierView
    :type nr_residues: int
    """
    residue_groups = []
    for chain in hierview.GetChains():
        for frag in chain.GetFragments():
            residues = [res for res in frag.GetResidues()]
            residue_slices = [residues[i:i+nr_residues] for i in range(0, len(residues), nr_residues)]
            for r in residue_slices:
                residue_groups.append(r)
    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