Visualizing Torsional Angle Distribution
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
Ingredients
|
Difficulty level
![../_images/chilly.png](../_images/chilly.png)
![../_images/chilly.png](../_images/chilly.png)
![../_images/chilly.png](../_images/chilly.png)
Download
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.
1class IsRotatableOrMacroCycleBond(oechem.OEUnaryBondPred):
2 """
3 Identifies rotatable bonds and single bonds in macro-cycles.
4 """
5 def __call__(self, bond):
6 """
7 :type mol: oechem.OEBondBase
8 :rtype: boolean
9 """
10 if bond.GetOrder() != 1:
11 return False
12 if bond.IsAromatic():
13 return False
14
15 isrotor = oechem.OEIsRotor()
16 if isrotor(bond):
17 return True
18
19 if oechem.OEBondGetSmallestRingSize(bond) >= 10:
20 return True
21
22 return False
1def get_dihedrals(mol, itag):
2 """
3 Iterates over rotatable bonds and identifies their dihedral
4 atoms. These atoms are added to the molecule in a group
5 using the given tag.
6
7 :type mol: oechem.OEMol
8 :type itag: int
9 :return: Number of dihedral angles identified
10 :rtype: int
11 """
12 nrdihedrals = 0
13 for bond in mol.GetBonds(IsRotatableOrMacroCycleBond()):
14 atomB = bond.GetBgn()
15 atomE = bond.GetEnd()
16
17 neighB = None
18 neighE = None
19
20 for atom in atomB.GetAtoms(oechem.OEIsHeavy()):
21 if atom != atomE:
22 neighB = atom
23 break
24 for atom in atomE.GetAtoms(oechem.OEIsHeavy()):
25 if atom != atomB:
26 neighE = atom
27 break
28
29 if neighB is None or neighE is None:
30 continue
31
32 atomorder = [neighB, atomB, atomE, neighE]
33 bondorder = [mol.GetBond(neighB, atomB), bond, mol.GetBond(neighE, atomE)]
34
35 if neighB.GetIdx() < neighE.GetIdx():
36 atomorder.reverse()
37 bondorder.reverse()
38
39 atoms = oechem.OEAtomVector(atomorder)
40 bonds = oechem.OEBondVector(bondorder)
41
42 nrdihedrals += 1
43 mol.NewGroup(itag, atoms, bonds)
44
45 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
1def set_dihedral_histograms(mol, itag, nrbins):
2 """
3 Iterates over the dihedral groups and bins the torsional
4 angles for each conformation. The histogram data is then
5 attached to the groups with the given tag.
6
7 :type mol: oechem.OEMol
8 :type itag: int
9 :type nrbins: int
10 """
11
12 angleinc = 360.0 / float(nrbins)
13
14 for group in mol.GetGroups(oechem.OEHasGroupType(itag)):
15 atoms = oechem.OEAtomVector()
16 for atom in group.GetAtoms():
17 atoms.append(atom)
18 histogram = [0] * nrbins
19
20 for conf in mol.GetConfs():
21 rad = oechem.OEGetTorsion(conf, atoms[0], atoms[1], atoms[2], atoms[3])
22 deg = math.degrees(rad)
23 deg = (deg + 360.0) % 360.0
24 binidx = int(math.floor((deg / angleinc)))
25 histogram[binidx] += 1
26
27 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.
1def depict_dihedrals(image, dimage, mol, refmol, opts, itag, nrbins, colorg):
2 """
3 Highlights the dihedral atoms of a torsion and the depicts the
4 corresponding dihedral angle histogram when hovering over
5 the center of the torsion on the molecule display.
6
7 :type image: oedepict.OEImageBase
8 :type dimage: oedepict.OEImageBase
9 :type mol: oechem.OEMol
10 :type refmol: oechem.OEMol
11 :type opts: oedepict.OE2DMolDisplayOptions
12 :type itag: int
13 :type nrbins: int
14 :type oechem.OEColorGradientBase
15 """
16
17 nrconfs = mol.NumConfs()
18 center = oedepict.OEGetCenter(dimage)
19 radius = min(dimage.GetWidth(), dimage.GetHeight()) * 0.40
20
21 draw_dihedral_circle(dimage, center, radius, nrbins, nrconfs)
22
23 suppressH = True
24 oegrapheme.OEPrepareDepictionFrom3D(mol, suppressH)
25 if refmol is not None:
26 oegrapheme.OEPrepareDepictionFrom3D(refmol, suppressH)
27
28 disp = oedepict.OE2DMolDisplay(mol, opts)
29
30 dihedrals = []
31 ref_dihedrals = []
32 centers = []
33 agroups = []
34 dgroups = []
35
36 nrdihedrals = 0
37 for group in mol.GetGroups(oechem.OEHasGroupType(itag)):
38
39 uniqueid = uuid.uuid4().hex
40 agroup = image.NewSVGGroup("torsion_area_" + uniqueid)
41 dgroup = image.NewSVGGroup("torsion_data_" + uniqueid)
42 oedepict.OEAddSVGHover(agroup, dgroup)
43
44 dihedrals.append(group)
45 if refmol is not None:
46 ref_dihedrals.append(get_reference_dihedral(group, refmol, itag))
47
48 centers.append(get_center(disp, group))
49 agroups.append(agroup)
50 dgroups.append(dgroup)
51 nrdihedrals += 1
52
53 for didx in range(0, nrdihedrals):
54
55 image.PushGroup(dgroups[didx])
56
57 dihedral = dihedrals[didx]
58 abset = oechem.OEAtomBondSet(dihedral.GetAtoms(), dihedral.GetBonds())
59 draw_highlight(image, disp, abset)
60 dihedral_histogram = dihedral.GetData(itag)
61 draw_dihedral_histogram(dimage, dihedral_histogram, center, radius, nrbins, nrconfs)
62
63 if refmol is not None and ref_dihedrals[didx] is not None:
64 ref_dihedral = ref_dihedrals[didx]
65 draw_reference_dihedral(dimage, ref_dihedral, itag, center, radius)
66
67 image.PopGroup(dgroups[didx])
68
69 clearbackground = True
70 oedepict.OERenderMolecule(image, disp, not clearbackground)
71
72 markpen = oedepict.OEPen(oechem.OEBlack, oechem.OEWhite, oedepict.OEFill_On, 1.0)
73 farpen = oedepict.OEPen(oechem.OEBlack, oechem.OERed, oedepict.OEFill_Off, 2.0)
74
75 angleinc = 360.0 / float(nrbins)
76
77 for didx in range(0, nrdihedrals):
78
79 image.PushGroup(agroups[didx])
80
81 dihedral = dihedrals[didx]
82 dihedral_histogram = dihedral.GetData(itag)
83 flexibility = determine_flexibility(dihedral_histogram)
84 color = colorg.GetColorAt(flexibility)
85 markpen.SetBackColor(color)
86
87 markradius = disp.GetScale() / 8.0
88 image.DrawCircle(centers[didx], markradius, markpen)
89
90 if refmol is not None and ref_dihedrals[didx] is not None:
91 ref_dihedral = ref_dihedrals[didx]
92 if get_closest_dihedral_angle(mol, dihedral, ref_dihedral, itag) > angleinc:
93 image.DrawCircle(centers[didx], markradius, farpen)
94
95 radius = disp.GetScale() / 4.0
96 image.DrawCircle(centers[didx], radius, oedepict.OESVGAreaPen)
97
98 image.PopGroup(agroups[didx])
Usage
Usage
dihedral2img.py
and acyclovi.sdf
multi-conformer molecule file
The following command will generate the image shown in Figure 1:
prompt > python3 dihedral2img.py -in acyclovi.sdf -out acyclovi.svg
Command Line Parameters
Simple parameter list
input/output options
-in : Input filename of a multi conformer molecule
-out : Output filename of the generated image
-ref : Input filename of reference molecule
visualization options
-flexibility : Visualize dihedral angle flexibility
-nrbins : Number of bins in the dihedral angle histogram
Discussion
OpenEye’s Omega TK can be used to generate diverse sets of low-energy conformations.
Usage
The following commands will generate a multi-conformer file
for molecule acyclovi
:
prompt > echo "c1nc2c(=O)[nH]c(nc2n1COCCO)N acyclovi" > acyclovi.ism
prompt > python3 omega.py acyclovi.ism acyclovi.sdf
Visualizing Torsion Flexibility
By using the -flexibility parameter, the flexibility of the torsions angles can be visualized using a color gradient. Torsions with high flexibility are colored red, while black color indicates restrained torsion angles.
The flexibility of the torsion is determined with the following rudimentary method:
1def determine_flexibility(histogram):
2 """
3 Returns the simple estimation of torsion flexibility by counting the
4 number of non-zero bins in the torsional histogram.
5
6 :type histogram: list(int)
7 """
8
9 nr_non_zero_bins = sum([1 for x in histogram if x > 0]) * 2
10 return nr_non_zero_bins
Usage
dihedral2img.py
and
penicillin.sdf
multi-conformer molecule file
The following command will generate the image shown in Figure 2:
prompt > python3 dihedral2img.py -in penicillin.sdf -out penicillin.svg -flexibility
The Figure 2 shows that the preferred angle of the amide bond in the molecule is around 180° coloring the corresponding bond circle black, while the other bonds have more flexibility. The utilized color gradient can revealed by hovering over the “Legend” label.
hover over any rotatable bond in the molecule (marked with a circle)
Figure 2. Example of visualizing torsion flexibility
Visualizing Torsional Angle Distribution with Reference
By using the -ref parameter, the torsion angle of a reference molecule (such as experimental conformation) can visualized in the generated image.
Usage
dihedral2img.py
and
0QI.sdf
multi-conformer molecule file with corresponding
0QI.pdb
reference
molecule file
The following command will generate the image shown in Figure 3:
prompt > python3 dihedral2img.py -in 0QI.sdf -ref 0QI.pdb -out 0QI.svg
When hovering over a rotatable bond in Figure 3 the corresponding torsional angle of the reference molecule is depicted in the middle of the radial histogram. Red circle is depicted around the bond marker when the reference angle is not close to any generated torsional angles. This allows to make instant judgment whether the generated conformations can reproduce an experimentally determined conformation.
hover over any rotatable bond in the molecule (marked with a circle)
Figure 3. Example of visualizing torsion flexibility with reference molecule
See also in OEChem TK manual
API
OEAtomBondSet class
OEBondGetSmallestRingSize function
OEExponentColorGradient class
OEGetTorsion function
OEGroupBase class
OEHasGroupType predicate
OEIsRotor predicate
OEUnaryBondPred class
See also in OEDepict TK manual
Theory
Molecule Depiction chapter
API
OE2DMolDisplay class
OE2DMolDisplayOptions class
OE2DPoint class
OEAddInteractiveIcon function
OEAddSVGHover function
OEDrawLegendLayout function
OEDrawSVGHoverText function
OEDrawTextToCenter function
OEFont class
OEImage class
OEImageFrame class
OELegendLayout class
OELegendLayoutOptions class
OEPen class
OERenderMolecule function
OESVGAreaPen object
See also in GraphemeTM TK manual
API
OEColorGradientLabel class
OEDrawColorGradient function
OEPrepareDepictionFrom3D function