Visualizing Torsional Angle Distribution new

Problem

You want to generate an interactive image (in svg file format) that visualizes the distribution of dihedral angles of rotatable bond in a multi-conformer molecule. See example in Figure 1.

hover over any rotatable bond in the molecule (marked with a circle)

Figure 1. Example of visualizing torsional information

../_images/dihedral2img-acyclovi.svg

Ingredients

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

Difficulty level

../_images/chilly6.png ../_images/chilly6.png ../_images/chilly6.png

Download

Download code

dihedral2img.py

See also the Usage subsection.

Solution

The first step of generating the image is to identify the rotatable bonds in the input molecule using the IsRotatableOrMacroCycleBond bond predicate. The get_dihedrals function iterates over rotatable bonds and identifies their dihedral atoms. These dihedral atoms are stored on the molecule in a OEGroupBase object for further process.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
class IsRotatableOrMacroCycleBond(oechem.OEUnaryBondPred):
    """
    Identifies rotatable bonds and single bonds in macro-cycles.
    """
    def __call__(self, bond):
        """
        :type mol: oechem.OEBondBase
        :rtype: boolean
        """
        if bond.GetOrder() != 1:
            return False
        if bond.IsAromatic():
            return False

        isrotor = oechem.OEIsRotor()
        if isrotor(bond):
            return True

        if oechem.OEBondGetSmallestRingSize(bond) >= 10:
            return True

        return False
 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
def get_dihedrals(mol, itag):
    """
    Iterates over rotatable bonds and identifies their dihedral
    atoms. These atoms are added to the molecule in a group
    using the given tag.

    :type mol: oechem.OEMol
    :type itag: int
    :return: Number of dihedral angles identified
    :rtype: int
    """
    nrdihedrals = 0
    for bond in mol.GetBonds(IsRotatableOrMacroCycleBond()):
        atomB = bond.GetBgn()
        atomE = bond.GetEnd()

        neighB = None
        neighE = None

        for atom in atomB.GetAtoms(oechem.OEIsHeavy()):
            if atom != atomE:
                neighB = atom
                break
        for atom in atomE.GetAtoms(oechem.OEIsHeavy()):
            if atom != atomB:
                neighE = atom
                break

        if neighB is None or neighE is None:
            continue

        atoms = oechem.OEAtomVector()
        atoms.append(neighB)
        atoms.append(atomB)
        atoms.append(atomE)
        atoms.append(neighE)

        bonds = oechem.OEBondVector()
        bonds.append(mol.GetBond(neighB, atomB))
        bonds.append(bond)
        bonds.append(mol.GetBond(neighE, atomE))

        nrdihedrals += 1
        mol.NewGroup(itag, atoms, bonds)

    return nrdihedrals

After the dihedral atoms are identified, the set_dihedral_histograms function is used to iterate over the conformation of the molecule and calculate torsion angles using the OEGetTorsion function . These angles are binned to

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
def set_dihedral_histograms(mol, itag, angleinc, nrbins):
    """
    Iterates over the dihedral groups and bins the torsional
    angles for each conformation.

    :type mol: oechem.OEMol
    :type itag: int
    :type angleinc: float
    :type nrbins: int
    """
    for group in mol.GetGroups(oechem.OEHasGroupType(itag)):
        atoms = oechem.OEAtomVector()
        for atom in group.GetAtoms():
            atoms.append(atom)

        histogram = oechem.OEUIntVector(nrbins)
        for conf in mol.GetConfs():
            rad = oechem.OEGetTorsion(conf, atoms[0], atoms[1], atoms[2], atoms[3])
            deg = math.degrees(rad)
            deg = (deg + 360.0) % 360.0
            binidx = int(math.floor((deg / angleinc)))
            histogram[binidx] += 1

        group.SetData(itag, histogram)

The last step is to highlight the dihedral atoms when hovered over and depict the corresponding dihedral angle histogram. In order to achieve the hover effect in the generated SVG image, SVG groups are utilized (OESVGGroup) in the depict_dihedrals function. For each dihedral two groups are created. These groups are associated by calling the OEAddSVGHover function: while hovering the mouse over objects drawn inside the torsion_area_<id> the objects drawn in the torsion_data_<id> will be displayed. The group id must be unique amongst all the ids in the SVG image. Everything that is rendered between the OEImageBase.PushGroup and the corresponding OEImageBase.PopGroup methods is considered “inside” the group.

It is important that the molecule is rendered into the image after the dihedral angles are highlighted (after the second loop). As a result the highlight will appear underneath the molecule rather than on top of.

In the last loop of the depict_dihedrals function transparent circles are drawn (using OESVGAreaPen) in the middle of the the dihedral angle representing the hover areas in the interactive SVG image.

 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
def depict_dihedrals(image, dimage, mol, opts, itag, angleinc, nrbins):
    """
    Highlights the dihedral atoms and the depicts the corresponding dihedral
    angle histogram in hover

    :type image: oedepict.OEImageBase
    :type dimage: oedepict.OEImageBase
    :type mol: oechem.OEMol
    :type opts: oedepict.OE2DMolDisplayOptions
    :type itag: int
    :type angleinc: float
    :type nrbins: int
    """

    nrconfs = mol.NumConfs()
    center = oedepict.OEGetCenter(dimage)
    radius = min(dimage.GetWidth(), dimage.GetHeight()) * 0.40
    draw_dihedral_circle(dimage, center, radius, angleinc, nrbins, nrconfs)

    oegrapheme.OEPrepareDepictionFrom3D(mol, True)
    disp = oedepict.OE2DMolDisplay(mol, opts)

    dihedrals = []
    centers = []
    agroups = []
    dgroups = []

    nrdihedrals = 0
    for group in mol.GetGroups(oechem.OEHasGroupType(itag)):

        uniqueid = uuid.uuid4().hex
        agroup = image.NewSVGGroup("torsion_area_" + uniqueid)
        dgroup = image.NewSVGGroup("torsion_data_" + uniqueid)
        oedepict.OEAddSVGHover(agroup, dgroup)

        dihedrals.append(group)
        centers.append(get_center(disp, group))
        agroups.append(agroup)
        dgroups.append(dgroup)
        nrdihedrals += 1

    for didx in range(0, nrdihedrals):

        image.PushGroup(dgroups[didx])

        dihedral = dihedrals[didx]
        abset = oechem.OEAtomBondSet(dihedral.GetAtoms(), dihedral.GetBonds())
        draw_highlight(image, disp, abset)
        draw_dihedral_histogram(dimage, dihedral.GetData(itag), center, radius, angleinc, nrbins, nrconfs)

        image.PopGroup(dgroups[didx])

    clearbackground = True
    oedepict.OERenderMolecule(image, disp, not clearbackground)

    markpen = oedepict.OEPen(oechem.OEBlack, oechem.OEWhite, oedepict.OEFill_On, 1.0)
    for didx in range(0, nrdihedrals):

        image.PushGroup(agroups[didx])

        markradius = disp.GetScale() / 10.0
        image.DrawCircle(centers[didx], markradius, markpen)
        radius = disp.GetScale() / 4.0
        image.DrawCircle(centers[didx], radius, oedepict.OESVGAreaPen)

        image.PopGroup(agroups[didx])

Usage

Download code

dihedral2img.py script and acyclovi.sdf multi-conformer molecule file

Usage

prompt > python3 dihedral2img.py -in acyclovi.sdf  -out acyclovi.svg

Discussion

OpenEye’s Omega TK toolkit can be used to generate diverse sets of low-energy conformations.

Download code

omega.py

Usage

prompt > echo "c1nc2c(=O)[nH]c(nc2n1COCCO)N acyclovi" > acyclovi.ism
prompt > omega.py acyclovi.ism acyclovi.sdf

See also in OEChem TK manual

API

See also in OEDepict TK manual

Theory

API

See also in GraphemeTM TK manual

API