1
2
3
4
5
6
7
8
9 """
10 rdkit - A Cinfony module for accessing the RDKit from CPython
11
12 Global variables:
13 Chem and AllChem - the underlying RDKit Python bindings
14 informats - a dictionary of supported input formats
15 outformats - a dictionary of supported output formats
16 descs - a list of supported descriptors
17 fps - a list of supported fingerprint types
18 forcefields - a list of supported forcefields
19 """
20
21 import os
22
23 from rdkit import Chem
24 from rdkit.Chem import AllChem, Draw
25 from rdkit.Chem import Descriptors
26
27 _descDict = dict(Descriptors.descList)
28
29 import rdkit.DataStructs
30 import rdkit.Chem.MACCSkeys
31 import rdkit.Chem.AtomPairs.Pairs
32 import rdkit.Chem.AtomPairs.Torsions
33
34
35 try:
36 import Tkinter as tk
37 import Image as PIL
38 import ImageTk as PILtk
39 except:
40 PILtk = None
41
42
43 try:
44 import aggdraw
45 from rdkit.Chem.Draw import aggCanvas
46 except ImportError:
47 aggdraw = None
48
49 fps = ['rdkit', 'layered', 'maccs', 'atompairs', 'torsions', 'morgan']
50 """A list of supported fingerprint types"""
51 descs = _descDict.keys()
52 """A list of supported descriptors"""
53
54 _formats = {'smi': "SMILES",
55 'can': "Canonical SMILES",
56 'mol': "MDL MOL file",
57 'mol2': "Tripos MOL2 file",
58 'sdf': "MDL SDF file",
59 'inchi':"InChI",
60 'inchikey':"InChIKey"}
61 _notinformats = ['can', 'inchikey']
62 _notoutformats = ['mol2']
63 if not Chem.INCHI_AVAILABLE:
64 _notinformats += ['inchi']
65 _notoutformats += ['inchi', 'inchikey']
66
67 informats = dict([(_x, _formats[_x]) for _x in _formats if _x not in _notinformats])
68 """A dictionary of supported input formats"""
69 outformats = dict([(_x, _formats[_x]) for _x in _formats if _x not in _notoutformats])
70 """A dictionary of supported output formats"""
71
72 _forcefields = {'uff': AllChem.UFFOptimizeMolecule}
73 forcefields = _forcefields.keys()
74 """A list of supported forcefields"""
77 """Iterate over the molecules in a file.
78
79 Required parameters:
80 format - see the informats variable for a list of available
81 input formats
82 filename
83
84 You can access the first molecule in a file using the next() method
85 of the iterator:
86 mol = readfile("smi", "myfile.smi").next()
87
88 You can make a list of the molecules in a file using:
89 mols = list(readfile("smi", "myfile.smi"))
90
91 You can iterate over the molecules in a file as shown in the
92 following code snippet:
93 >>> atomtotal = 0
94 >>> for mol in readfile("sdf", "head.sdf"):
95 ... atomtotal += len(mol.atoms)
96 ...
97 >>> print atomtotal
98 43
99 """
100 if not os.path.isfile(filename):
101 raise IOError, "No such file: '%s'" % filename
102 format = format.lower()
103
104
105
106 if format=="sdf":
107 iterator = Chem.SDMolSupplier(filename)
108 def sdf_reader():
109 for mol in iterator:
110 yield Molecule(mol)
111 return sdf_reader()
112 elif format=="mol":
113 def mol_reader():
114 yield Molecule(Chem.MolFromMolFile(filename))
115 return mol_reader()
116 elif format=="mol2":
117 def mol_reader():
118 yield Molecule(Chem.MolFromMol2File(filename))
119 return mol_reader()
120 elif format=="smi":
121 iterator = Chem.SmilesMolSupplier(filename, delimiter=" \t",
122 titleLine=False)
123 def smi_reader():
124 for mol in iterator:
125 yield Molecule(mol)
126 return smi_reader()
127 elif format=='inchi' and Chem.INCHI_AVAILABLE:
128 def inchi_reader():
129 for line in open(filename, 'r'):
130 mol = Chem.inchi.MolFromInchi(line.strip())
131 yield Molecule(mol)
132 return inchi_reader()
133 else:
134 raise ValueError, "%s is not a recognised RDKit format" % format
135
137 """Read in a molecule from a string.
138
139 Required parameters:
140 format - see the informats variable for a list of available
141 input formats
142 string
143
144 Example:
145 >>> input = "C1=CC=CS1"
146 >>> mymol = readstring("smi", input)
147 >>> len(mymol.atoms)
148 5
149 """
150 format = format.lower()
151 if format=="mol":
152 mol = Chem.MolFromMolBlock(string)
153 elif format=="mol2":
154 mol = Chem.MolFromMol2Block(string)
155 elif format=="smi":
156 mol = Chem.MolFromSmiles(string)
157 elif format=='inchi' and Chem.INCHI_AVAILABLE:
158 mol = Chem.inchi.MolFromInchi(string)
159 else:
160 raise ValueError,"%s is not a recognised RDKit format" % format
161 if mol:
162 return Molecule(mol)
163 else:
164 raise IOError, "Failed to convert '%s' to format '%s'" % (
165 string, format)
166
168 """Represent a file to which *output* is to be sent.
169
170 Required parameters:
171 format - see the outformats variable for a list of available
172 output formats
173 filename
174
175 Optional parameters:
176 overwite -- if the output file already exists, should it
177 be overwritten? (default is False)
178
179 Methods:
180 write(molecule)
181 close()
182 """
183 - def __init__(self, format, filename, overwrite=False):
184 self.format = format
185 self.filename = filename
186 if not overwrite and os.path.isfile(self.filename):
187 raise IOError, "%s already exists. Use 'overwrite=True' to overwrite it." % self.filename
188 if format=="sdf":
189 self._writer = Chem.SDWriter(self.filename)
190 elif format=="smi":
191 self._writer = Chem.SmilesWriter(self.filename, isomericSmiles=True)
192 elif format in ('inchi', 'inchikey') and Chem.INCHI_AVAILABLE:
193 self._writer= open(filename, 'w')
194 else:
195 raise ValueError,"%s is not a recognised RDKit format" % format
196 self.total = 0
197
198 - def write(self, molecule):
199 """Write a molecule to the output file.
200
201 Required parameters:
202 molecule
203 """
204 if not self.filename:
205 raise IOError, "Outputfile instance is closed."
206 if self.format in ('inchi', 'inchikey'):
207 self._writer.write(molecule.write(self.format) +'\n')
208 else:
209 self._writer.write(molecule.Mol)
210 self.total += 1
211
213 """Close the Outputfile to further writing."""
214 self.filename = None
215 self._writer.flush()
216 del self._writer
217
219 """Represent an rdkit Molecule.
220
221 Required parameter:
222 Mol -- an RDKit Mol or any type of cinfony Molecule
223
224 Attributes:
225 atoms, data, formula, molwt, title
226
227 Methods:
228 addh(), calcfp(), calcdesc(), draw(), localopt(), make3D(), removeh(),
229 write()
230
231 The underlying RDKit Mol can be accessed using the attribute:
232 Mol
233 """
234 _cinfony = True
235
237 if hasattr(Mol, "_cinfony"):
238 a, b = Mol._exchange
239 if a == 0:
240 molecule = readstring("smi", b)
241 else:
242 molecule = readstring("mol", b)
243 Mol = molecule.Mol
244
245 self.Mol = Mol
246
247 @property
248 - def atoms(self): return [Atom(rdkatom) for rdkatom in self.Mol.GetAtoms()]
249 @property
251 @property
252 - def molwt(self): return Descriptors.MolWt(self.Mol)
253 @property
256
257 if "_Name" in self.data:
258 return self.data["_Name"]
259 else:
260 return ""
261 - def _settitle(self, val): self.Mol.SetProp("_Name", val)
262 title = property(_gettitle, _settitle)
263 @property
265 if self.Mol.GetNumConformers() == 0:
266 return (0, self.write("smi"))
267 else:
268 return (1, self.write("mol"))
269
271 """Add hydrogens."""
272 self.Mol = Chem.AddHs(self.Mol)
273
275 """Remove hydrogens."""
276 self.Mol = Chem.RemoveHs(self.Mol)
277
278 - def write(self, format="smi", filename=None, overwrite=False):
279 """Write the molecule to a file or return a string.
280
281 Optional parameters:
282 format -- see the informats variable for a list of available
283 output formats (default is "smi")
284 filename -- default is None
285 overwite -- if the output file already exists, should it
286 be overwritten? (default is False)
287
288 If a filename is specified, the result is written to a file.
289 Otherwise, a string is returned containing the result.
290
291 To write multiple molecules to the same file you should use
292 the Outputfile class.
293 """
294 format = format.lower()
295 if filename:
296 if not overwrite and os.path.isfile(filename):
297 raise IOError, "%s already exists. Use 'overwrite=True' to overwrite it." % filename
298 if format=="smi":
299 result = Chem.MolToSmiles(self.Mol, isomericSmiles=True, canonical=False)
300 elif format=="can":
301 result = Chem.MolToSmiles(self.Mol, isomericSmiles=True, canonical=True)
302 elif format=="mol":
303 result = Chem.MolToMolBlock(self.Mol)
304 elif format in ('inchi', 'inchikey') and Chem.INCHI_AVAILABLE:
305 result = Chem.inchi.MolToInchi(self.Mol)
306 if format == 'inchikey':
307 result = Chem.inchi.InchiToInchiKey(result)
308 else:
309 raise ValueError,"%s is not a recognised RDKit format" % format
310 if filename:
311 print >> open(filename, "w"), result
312 else:
313 return result
314
316 """Iterate over the Atoms of the Molecule.
317
318 This allows constructions such as the following:
319 for atom in mymol:
320 print atom
321 """
322 return iter(self.atoms)
323
326
328 """Calculate descriptor values.
329
330 Optional parameter:
331 descnames -- a list of names of descriptors
332
333 If descnames is not specified, all available descriptors are
334 calculated. See the descs variable for a list of available
335 descriptors.
336 """
337 if not descnames:
338 descnames = descs
339 ans = {}
340 for descname in descnames:
341 try:
342 desc = _descDict[descname]
343 except KeyError:
344 raise ValueError, "%s is not a recognised RDKit descriptor type" % descname
345 ans[descname] = desc(self.Mol)
346 return ans
347
348 - def calcfp(self, fptype="rdkit", opt=None):
349 """Calculate a molecular fingerprint.
350
351 Optional parameters:
352 fptype -- the fingerprint type (default is "rdkit"). See the
353 fps variable for a list of of available fingerprint
354 types.
355 opt -- a dictionary of options for fingerprints. Currently only used
356 for radius and bitInfo in Morgan fingerprints.
357 """
358 if opt == None:
359 opt = {}
360 fptype = fptype.lower()
361 if fptype=="rdkit":
362 fp = Fingerprint(Chem.RDKFingerprint(self.Mol))
363 elif fptype=="layered":
364 fp = Fingerprint(Chem.LayeredFingerprint(self.Mol))
365 elif fptype=="maccs":
366 fp = Fingerprint(Chem.MACCSkeys.GenMACCSKeys(self.Mol))
367 elif fptype=="atompairs":
368
369 fp = Chem.AtomPairs.Pairs.GetAtomPairFingerprintAsIntVect(self.Mol)
370 elif fptype=="torsions":
371
372 fp = Chem.AtomPairs.Torsions.GetTopologicalTorsionFingerprintAsIntVect(self.Mol)
373 elif fptype == "morgan":
374 info = opt.get('bitInfo', None)
375 radius = opt.get('radius', 4)
376 fp = Fingerprint(Chem.rdMolDescriptors.GetMorganFingerprintAsBitVect(self.Mol,radius,bitInfo=info))
377 else:
378 raise ValueError, "%s is not a recognised RDKit Fingerprint type" % fptype
379 return fp
380
381 - def draw(self, show=True, filename=None, update=False, usecoords=False):
382 """Create a 2D depiction of the molecule.
383
384 Optional parameters:
385 show -- display on screen (default is True)
386 filename -- write to file (default is None)
387 update -- update the coordinates of the atoms to those
388 determined by the structure diagram generator
389 (default is False)
390 usecoords -- don't calculate 2D coordinates, just use
391 the current coordinates (default is False)
392
393 Aggdraw or Cairo is used for 2D depiction. Tkinter and
394 Python Imaging Library are required for image display.
395 """
396 if not usecoords and update:
397 AllChem.Compute2DCoords(self.Mol)
398 usecoords = True
399 mol = Chem.Mol(self.Mol.ToBinary())
400 if not usecoords:
401 AllChem.Compute2DCoords(mol)
402
403 if filename:
404 Draw.MolToFile(mol, filename)
405 if show:
406 if not tk:
407 errormessage = ("Tkinter or Python Imaging "
408 "Library not found, but is required for image "
409 "display. See installation instructions for "
410 "more information.")
411 raise ImportError(errormessage)
412 img = Draw.MolToImage(mol)
413 root = tk.Tk()
414 root.title((hasattr(self, "title") and self.title)
415 or self.__str__().rstrip())
416 frame = tk.Frame(root, colormap="new", visual='truecolor').pack()
417 imagedata = PILtk.PhotoImage(img)
418 label = tk.Label(frame, image=imagedata).pack()
419 quitbutton = tk.Button(root, text="Close", command=root.destroy).pack(fill=tk.X)
420 root.mainloop()
421
422 - def localopt(self, forcefield = "uff", steps = 500):
423 """Locally optimize the coordinates.
424
425 Optional parameters:
426 forcefield -- default is "uff". See the forcefields variable
427 for a list of available forcefields.
428 steps -- default is 500
429
430 If the molecule does not have any coordinates, make3D() is
431 called before the optimization.
432 """
433 forcefield = forcefield.lower()
434 if self.Mol.GetNumConformers() == 0:
435 self.make3D(forcefield)
436 _forcefields[forcefield](self.Mol, maxIters = steps)
437
438 - def make3D(self, forcefield = "uff", steps = 50):
439 """Generate 3D coordinates.
440
441 Optional parameters:
442 forcefield -- default is "uff". See the forcefields variable
443 for a list of available forcefields.
444 steps -- default is 50
445
446 Once coordinates are generated, a quick
447 local optimization is carried out with 50 steps and the
448 UFF forcefield. Call localopt() if you want
449 to improve the coordinates further.
450 """
451 forcefield = forcefield.lower()
452 success = AllChem.EmbedMolecule(self.Mol)
453 if success == -1:
454 success = AllChem.EmbedMolecule(self.Mol,
455 useRandomCoords = True)
456 if success == -1:
457 raise Error, "Embedding failed!"
458 self.localopt(forcefield, steps)
459
461 """Represent an rdkit Atom.
462
463 Required parameters:
464 Atom -- an RDKit Atom
465
466 Attributes:
467 atomicnum, coords, formalcharge
468
469 The original RDKit Atom can be accessed using the attribute:
470 Atom
471 """
472
475 @property
477 @property
479 owningmol = self.Atom.GetOwningMol()
480 if owningmol.GetNumConformers() == 0:
481 raise AttributeError, "Atom has no coordinates (0D structure)"
482 idx = self.Atom.GetIdx()
483 atomcoords = owningmol.GetConformer().GetAtomPosition(idx)
484 return (atomcoords[0], atomcoords[1], atomcoords[2])
485 @property
487
489 if hasattr(self, "coords"):
490 return "Atom: %d (%.2f %.2f %.2f)" % (self.atomicnum, self.coords[0],
491 self.coords[1], self.coords[2])
492 else:
493 return "Atom: %d (no coords)" % (self.atomicnum)
494
496 """A Smarts Pattern Matcher
497
498 Required parameters:
499 smartspattern
500
501 Methods:
502 findall(molecule)
503
504 Example:
505 >>> mol = readstring("smi","CCN(CC)CC") # triethylamine
506 >>> smarts = Smarts("[#6][#6]") # Matches an ethyl group
507 >>> print smarts.findall(mol)
508 [(0, 1), (3, 4), (5, 6)]
509
510 The numbers returned are the indices (starting from 0) of the atoms
511 that match the SMARTS pattern. In this case, there are three matches
512 for each of the three ethyl groups in the molecule.
513 """
515 """Initialise with a SMARTS pattern."""
516 self.rdksmarts = Chem.MolFromSmarts(smartspattern)
517 if not self.rdksmarts:
518 raise IOError, "Invalid SMARTS pattern."
519
521 """Find all matches of the SMARTS pattern to a particular molecule.
522
523 Required parameters:
524 molecule
525 """
526 return molecule.Mol.GetSubstructMatches(self.rdksmarts)
527
529 """Store molecule data in a dictionary-type object
530
531 Required parameters:
532 Mol -- an RDKit Mol
533
534 Methods and accessor methods are like those of a dictionary except
535 that the data is retrieved on-the-fly from the underlying Mol.
536
537 Example:
538 >>> mol = readfile("sdf", 'head.sdf').next()
539 >>> data = mol.data
540 >>> print data
541 {'Comment': 'CORINA 2.61 0041 25.10.2001', 'NSC': '1'}
542 >>> print len(data), data.keys(), data.has_key("NSC")
543 2 ['Comment', 'NSC'] True
544 >>> print data['Comment']
545 CORINA 2.61 0041 25.10.2001
546 >>> data['Comment'] = 'This is a new comment'
547 >>> for k,v in data.iteritems():
548 ... print k, "-->", v
549 Comment --> This is a new comment
550 NSC --> 1
551 >>> del data['NSC']
552 >>> print len(data), data.keys(), data.has_key("NSC")
553 1 ['Comment'] False
554 """
558 if not key in self:
559 raise KeyError, "'%s'" % key
561 return self._mol.GetPropNames()
563 return [self._mol.GetProp(x) for x in self.keys()]
567 return iter(self.keys())
569 return iter(self.items())
571 return len(self.keys())
573 return self._mol.HasProp(key)
575 self._testforkey(key)
576 self._mol.ClearProp(key)
578 for key in self:
579 del self[key]
582 - def update(self, dictionary):
583 for k, v in dictionary.iteritems():
584 self[k] = v
586 self._testforkey(key)
587 return self._mol.GetProp(key)
589 self._mol.SetProp(key, str(value))
592
594 """A Molecular Fingerprint.
595
596 Required parameters:
597 fingerprint -- a vector calculated by one of the fingerprint methods
598
599 Attributes:
600 fp -- the underlying fingerprint object
601 bits -- a list of bits set in the Fingerprint
602
603 Methods:
604 The "|" operator can be used to calculate the Tanimoto coeff. For example,
605 given two Fingerprints 'a', and 'b', the Tanimoto coefficient is given by:
606 tanimoto = a | b
607 """
609 self.fp = fingerprint
611 return rdkit.DataStructs.FingerprintSimilarity(self.fp, other.fp)
613 if attr == "bits":
614
615 return list(self.fp.GetOnBits())
616 else:
617 raise AttributeError, "Fingerprint has no attribute %s" % attr
619 return ", ".join([str(x) for x in _compressbits(self.fp)])
620
622 """Compress binary vector into vector of long ints.
623
624 This function is used by the Fingerprint class.
625
626 >>> _compressbits([0, 1, 0, 0, 0, 1], 2)
627 [2, 0, 2]
628 """
629 ans = []
630 for start in range(0, len(bitvector), wordsize):
631 compressed = 0
632 for i in range(wordsize):
633 if i + start < len(bitvector) and bitvector[i + start]:
634 compressed += 2**i
635 ans.append(compressed)
636
637 return ans
638
639
640 if __name__=="__main__":
641 import doctest
642 doctest.testmod()
643