This task comes from working with researchers who have a desktop tool which displays a large number (a few hundred) structures on a scrolling page. The structures are the database hits from searching with a given SMARTS query, and the common SMARTS substructures are aligned in each depictions. This makes it easier for them to identify the differences between the structures.

They want to replace it with something newer, and they would like more control over how the structure is displayed. (For example, some people want the depiction tool to first follow corporate guidelines on depiction layout. Doing that is outside the scope of this task.)


Use the first 16 structures of the benzodiazepine SD file to make a 4x4 grid of depictions as a single image. The first structure is in the upper-left corner, the second is to its right, and so on. Each depiction should include the title field of the corresponding record, which in this case is the PubChem identifier.

Use "[#7]~1~[#6]~[#6]~[#7]~[#6]~[#6]~2~[#6]~[#6]~[#6]~[#6]~[#6]12" as the common SMARTS substructure. This is the fused ring of the benzodiazepine system but without bond type or atom aromaticity information. Use the first molecule as the reference depiction. All other depictions must have the depiction of their common substructure aligned to the reference.

The final image should be 400x400 pixels, which means each structure depiction should be 100x100 pixels. If the toolkit doesn't support multiple depictions in the same image then depict each to a different image. A desktop tool would need a way to write the image data to a GUI widget, but for now just write the result to a PNG or GIF file.

The point here is to show how to do the 2D structure alignment. The other details are less important. Use your own preferences for coloring scheme (eg, color on white background or on black background). You may assume the input structure has the correct MDL chemistry. (This is only relevant for toolkits which support multiple chemistry models.)


Aligned depictions with the OpenEye toolkits.

Output from OpenEye/Python

from openeye.oechem import *
from openeye.oedepict import *
def make_view():
    view = OEDepictView(CELL_WIDTH, CELL_HEIGHT)
    view.SetBackColor(0, 0, 0)
    view.SetForeColor(255, 255, 255)
    return view
# This assumes the input was in proper SD format already
def prepare_mol_for_depiction(mol):
    if mol.GetDimension() == 3:
core_subsearch = OESubSearch("[#7]~1~[#6]~[#6]~[#7]~[#6]~[#6]~2~[#6]~[#6]~[#6]~[#6]~[#6]12")
ifs = oemolistream("benzodiazepine.sdf.gz")
img = OE8BitImage(CELL_WIDTH*4, CELL_HEIGHT*4)
view = make_view()
at_first_structure = True
for record_num, mol in enumerate(ifs.GetOEGraphMols()):
    for match in core_subsearch.Match(mol):
        if at_first_structure:
            # save the first structure's depiction coordinates
            coords = [mol.GetCoords(atom) for atom in match.GetTargetAtoms()]
            # Align to the first structure's coordinates
            fixed = [0] * mol.GetMaxAtomIdx()
            for atom, coord in zip(match.GetTargetAtoms(), coords):
                mol.SetCoords(atom, coord)
                fixed[atom.GetIdx()] = 1
            OEDepictFixedCoordinates(mol, fixed)
        view.RenderImage(img, at_first_structure,
                         record_num//4 * CELL_WIDTH, record_num %4 * CELL_HEIGHT)
        at_first_structure = False
        # Only align the first match
    # Stop after 16 structures
    if record_num == 15:
# Save the final image
ofs = oeofstream("align_oe.gif")
OEWriteGIF(ofs, img)


Aligned depictions with the RDKit toolkit.

Output from RDKit/Python

from rdkit import Chem,Geometry
from rdkit.Chem import AllChem
# read in the first 16 molecules:
suppl = Chem.SDMolSupplier('benzodiazepine.sdf')
first16 = [suppl[x] for x in range(16)]
# do the substructure matching:
patt = Chem.MolFromSmarts("[#7]~1~[#6]~[#6]~[#7]~[#6]~[#6]~2~[#6]~[#6]~[#6]~[#6]~[#6]12")
matchVs= [x.GetSubstructMatch(patt) for x in first16]
# compute the reference coordinates:
coords = [first16[0].GetConformer().GetAtomPosition(x) for x in matchVs[0]]
coords2D = [Geometry.Point2D(pt.x,pt.y) for pt in coords]
# now generate coords for the other molecules using that reference:
for molIdx in range(1,16):
    mol = first16[molIdx]
    for i,coord in enumerate(coords2D):
        coordDict[matchVs[molIdx][i]] = coord
# generate 100x100 images for each molecule:
from rdkit.Chem import Draw
imgs = [Draw.MolToImage(x,size=(100,100)) for x in first16]
# and use PIL to combine them:
import Image
img ='RGBA',(400,400))
for i in range(16): img.paste(imgs[i],(100*(i%4),100*(i//4),100*(i%4+1),100*(i//4+1)))'combined.rdkit.png')
Community content is available under CC-BY-SA unless otherwise noted.