Line 128: Line 128:
prop setparam E_GIF headerproperty E_NAME bgcolor white atomcolor element bonds 10
prop setparam E_GIF headerproperty E_NAME bgcolor white atomcolor element bonds 10
set fh [molfile open benzodiazepine.sdf.gz r hydrogens add]
set fh [molfile open benzodiazepine.sdf.gz r hydrogens add]
set dh [dataset create [molfile scan $fh {record <= 16} enslist]]
set dh [dataset create [molfile read $fh {} 16]]
set ss [ens create {[#7]~1~[#6]~[#6]~[#7]~[#6]~[#6]~2~[#6]~[#6]~[#6]~[#6]~[#6]12} smarts]
set ss [ens create {[#7]~1~[#6]~[#6]~[#7]~[#6]~[#6]~2~[#6]~[#6]~[#6]~[#6]~[#6]12} smarts]
foreach eh [dataset ens $dh] {
foreach eh [dataset ens $dh] {

Revision as of 14:53, November 12, 2010

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')


Aligned depictions with the Cactvs toolkit.

Output from Cactvs/Tcl

prop setparam D_GIF ncols 4 nrows 4 boxwidth 100 boxheight 100 filename aligned4x4.png format png padding 0
prop setparam E_GIF headerproperty E_NAME bgcolor white atomcolor element bonds 10
set fh [molfile open benzodiazepine.sdf.gz r hydrogens add]
set dh [dataset create [molfile read $fh {} 16]]
set ss [ens create {[#7]~1~[#6]~[#6]~[#7]~[#6]~[#6]~2~[#6]~[#6]~[#6]~[#6]~[#6]12} smarts]
foreach eh [dataset ens $dh] {
    match ss -align besteffort $ss $eh
dataset get $dh D_GIF
Community content is available under CC-BY-SA unless otherwise noted.