cinfony.jchem

1 #-*. coding: utf-8 -*- 2 ## Copyright (c) 2012, Adrià Cereto-Massagué, Noel O'Boyle 3 ## All rights reserved. 4 ## 5 ## This file is part of Cinfony. 6 ## The contents are covered by the terms of the BSD license 7 ## which is included in the file LICENSE_BSD.txt. 8 9 """ 10 jchem - A Cinfony module for accessing ChemAxon's JChem from CPython and Jython 11 12 Global variables: 13 chemaxon - the underlying JChem Java library 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 import sys 21 import os 22 from glob import glob 23 24 if sys.platform[:4] == "java": 25 classpath = [] 26 if 'JCHEMDIR' in os.environ: 27 assert os.path.isdir(os.path.join(os.environ['JCHEMDIR'], 'lib')) 28 for jar in glob(os.path.join(os.path.join(os.environ['JCHEMDIR'],'lib'), '*.jar')): 29 classpath.append(jar) 30 31 if sys.platform[:4] == "java" or sys.platform[:3] == "cli": 32 import sys 33 sys.path = classpath + sys.path 34 import java, javax 35 import chemaxon 36 from chemaxon.util import MolHandler 37 #Exceptions are handled differently in jpype and jython. We need to wrap them: 38 MolExportException = chemaxon.marvin.io.MolExportException 39 MolFormatException = chemaxon.formats.MolFormatException 40 else: 41 from jpype import * 42 43 if not isJVMStarted(): 44 _jvm = os.environ['JPYPE_JVM'] 45 if _jvm[0] == '"': # Remove trailing quotes 46 _jvm = _jvm[1:-1] 47 _cp = os.pathsep.join(os.environ.get('CLASSPATH', '').split(os.pathsep)) 48 startJVM(_jvm, "-Djava.class.path=" + _cp) 49 50 chemaxon = JPackage("chemaxon") 51 MolHandler = chemaxon.util.MolHandler 52 try: 53 _testmol = MolHandler() 54 except TypeError: 55 raise ImportError, "jchem.jar file cannot be found." 56 57 # Exception wrappers for JPype 58 MolExportException = JavaException 59 MolFormatException = JavaException 60 61 _descset = set(['HAcc', 'HDon', 'Heavy', 'LogD', 'LogP', 'Mass', 'TPSA']) 62 _descset.update(dir(chemaxon.descriptors.scalars)) 63 descs = [cls for cls in _descset if hasattr(getattr(chemaxon.descriptors.scalars, cls),'generate') and cls != 'LogD'] + ['RotatableBondsCount'] 64 """A list of supported descriptors""" 65 fps = ['ecfp'] 66 """A list of supported fingerprint types""" 67 forcefields = ["mmff94"] 68 """A list of supported forcefields""" 69 70 informats = { 71 'smi': "SMILES" 72 ,'cxsmi': "ChemAxon exntended SMILES" 73 ,'mol': "MDL MOL" 74 ,'sdf': "MDL SDF" 75 ,'inchi': "InChI" 76 ,'cml': "Chemical Markup Language" 77 , 'mrv':'Marvin Documents' 78 , 'skc':'ISIS/Draw sketch file' 79 , 'cdx':'ChemDraw sketch file' 80 , 'cdxml':'ChemDraw sketch file' 81 , "name":"Common name" 82 , "peptide":"Aminoacid sequence" 83 , "sybyl":"Tripos SYBYL" 84 , "pdb":"PDB" 85 , "xyz":"XYZ" 86 , 'cube':'Gaussian cube' 87 , 'gout':'Gaussian output format' 88 } 89 """A dictionary of supported input formats""" 90 91 outformats = { 92 'smi': "SMILES" 93 ,'cxsmi': "ChemAxon exntended SMILES" 94 ,'mol': "MDL MOL" 95 ,'sdf': "MDL SDF" 96 ,'inchi': "InChI" 97 ,'inchikey': "InChIKey" 98 ,'cml': "CML" 99 , 'mrv':'Marvin Documents' 100 , 'skc':'ISIS/Draw sketch file' 101 , 'cdx':'ChemDraw sketch file' 102 , 'cdxml':'ChemDraw sketch file' 103 , "name":"Common name" 104 , "peptide":"Aminoacid sequence" 105 , "sybyl":"Tripos SYBYL" 106 , "pdb":"PDB" 107 , "xyz":"XYZ" 108 , 'cube':'Gaussian cube' 109 , 'gjf':'Gaussian input format' 110 } 111 """A dictionary of supported output formats"""

112 113 -def readfile(format, filename):

114 """Iterate over the molecules in a file. 115 116 Required parameters: 117 format - Ignored, but needed for compatibility with other cinfony 118 modules and also good for readability 119 filename 120 121 You can access the first molecule in a file using the next() method 122 of the iterator: 123 mol = readfile("smi", "myfile.smi").next() 124 125 You can make a list of the molecules in a file using: 126 mols = list(readfile("smi", "myfile.smi")) 127 128 You can iterate over the molecules in a file as shown in the 129 following code snippet: 130 >>> atomtotal = 0 131 >>> for mol in readfile("sdf", "head.sdf"): 132 ... atomtotal += len(mol.atoms) 133 ... 134 >>> print atomtotal 135 43 136 """ 137 if not os.path.isfile(filename): 138 raise IOError, "No such file: '%s'" % filename 139 if not format in outformats: 140 raise ValueError("%s is not a recognised JChem format" % format) 141 try: 142 mi = chemaxon.formats.MolImporter(filename) 143 mol = mi.read() 144 while mol: 145 mol.aromatize() 146 yield Molecule(mol) 147 mol = mi.read() 148 except chemaxon.formats.MolFormatException: 149 raise ValueError("%s is not a recognised JChem format" % format)

150

151 -def readstring(format, string):

152 """Read in a molecule from a string. 153 154 Required parameters: 155 format - Ignored, but needed for compatibility with other cinfony 156 modules and also good for readability 157 string 158 159 Example: 160 >>> input = "C1=CC=CS1" 161 >>> mymol = readstring("smi", input) 162 >>> len(mymol.atoms) 163 5 164 """ 165 format = format.lower() 166 if format not in informats: 167 raise ValueError("%s is not a recognised JChem format" % format) 168 try: 169 mh = MolHandler(string) 170 return Molecule(mh.molecule) 171 except MolFormatException, ex: 172 if sys.platform[:4] != "java": 173 #Jpype exception 174 ex = ex.message() 175 raise IOError, ex 176 else: 177 raise IOError("Problem reading the supplied string")

178

179 -class Outputfile(object):

180 """Represent a file to which *output* is to be sent. 181 182 Required parameters: 183 format - see the outformats variable for a list of available 184 output formats 185 filename 186 187 Optional parameters: 188 overwite -- if the output file already exists, should it 189 be overwritten? (default is False) 190 191 Methods: 192 write(molecule) 193 close() 194 """

195 - def __init__(self, format, filename, overwrite=False):

196 if ':' in format: 197 format, options = format.split(':') 198 if options: 199 options = ':' + options 200 else: 201 options = '' 202 self.format = format.lower() 203 self.filename = filename 204 if not overwrite and os.path.isfile(self.filename): 205 raise IOError, "%s already exists. Use 'overwrite=True' to overwrite it." % self.filename 206 if format in ("smi", 'cxsmi'): 207 if not options: 208 options = ':a-H' 209 out = chemaxon.formats.MolExporter.exportToFormat(self.Molecule,format +'les:a-H') 210 try: 211 self._writer = chemaxon.formats.MolExporter(filename, format + options) 212 except MolExportException, e: 213 raise ValueError(e) 214 self.total = 0 # The total number of molecules written to the file

215

216 - def write(self, molecule):

217 """Write a molecule to the output file. 218 219 Required parameters: 220 molecule 221 """ 222 if not self.filename: 223 raise IOError, "Outputfile instance is closed." 224 self._writer.write(molecule.Molecule) 225 self.total += 1

226

227 - def close(self):

228 """Close the Outputfile to further writing.""" 229 self.filename = None 230 self._writer.close()

231

232 -class Molecule(object):

233 """Represent a JChem Molecule. 234 235 Required parameters: 236 Molecule -- a JChem Molecule or any type of cinfony Molecule 237 238 Attributes: 239 atoms, data, exactmass, formula, molwt, title 240 241 Methods: 242 addh(), calcfp(), calcdesc(), draw(), removeh(), write() 243 244 The underlying JChem Molecule can be accessed using the attribute: 245 Molecule 246 The associated JChem MolHandler can be accessed using the attribute: 247 MolHandler 248 """ 249 _cinfony = True 250

251 - def __init__(self, Molecule):

252 253 if hasattr(Molecule, "_cinfony"): 254 a, b = Molecule._exchange 255 if a == 0: 256 mol = readstring("smi", b) 257 else: 258 mol = readstring("sdf", b) 259 Molecule = mol.Molecule 260 261 self.Molecule = Molecule 262 self.MolHandler = chemaxon.util.MolHandler(self.Molecule) 263 self.MolHandler.aromatize()

264 265 @property

266 - def atoms(self): return [Atom(atom) for atom in self.Molecule.atomArray]

267 @property

268 - def data(self): return MoleculeData(self)

269 @property

270 - def formula(self): return self.MolHandler.calcMolFormula()

271 @property

272 - def exactmass(self):

273 return self.MolHandler.calcMolWeightInDouble()

274 @property

275 - def molwt(self):

276 return self.MolHandler.calcMolWeight()

277 - def _gettitle(self): return self.Molecule.getName()

278 - def _settitle(self, val): self.Molecule.setName(val)

279 title = property(_gettitle, _settitle) 280 @property

281 - def _exchange(self):

282 if self.Molecule.dim > 1: 283 return (1, self.write("mol")) 284 else: 285 return (0, self.write("smi"))

286

287 - def __iter__(self):

288 """Iterate over the Atoms of the Molecule. 289 290 This allows constructions such as the following: 291 for atom in mymol: 292 print atom 293 """ 294 return iter(self.atoms)

295

296 - def __str__(self):

297 return self.write()

298

299 - def addh(self):

300 """Add hydrogens.""" 301 self.MolHandler.addHydrogens()

302

303 - def removeh(self):

304 """Remove hydrogens.""" 305 self.MolHandler.removeHydrogens()

306

307 - def write(self, format="smi", filename=None, overwrite=False):

308 """Write the molecule to a file or return a string. 309 310 Optional parameters: 311 format -- see the informats variable for a list of available 312 output formats (default is "smi") 313 filename -- default is None 314 overwite -- if the output file already exists, should it 315 be overwritten? (default is False) 316 317 If a filename is specified, the result is written to a file. 318 Otherwise, a string is returned containing the result. 319 320 To write multiple molecules to the same file you should use 321 the Outputfile class. 322 """ 323 if ':' in format: 324 format, options = format.split(':') 325 if options: 326 options = ':' + options 327 else: 328 options = '' 329 format = format.lower() 330 if format not in outformats: 331 raise ValueError("%s is not a recognised format" % format) 332 333 if filename is not None and not overwrite and os.path.isfile(filename): 334 raise IOError, "%s already exists. Use 'overwrite=True' to overwrite it." % filename 335 336 if format in ("smi", 'cxsmi'): 337 if not options: 338 options = ':a-H' 339 out = chemaxon.formats.MolExporter.exportToFormat(self.Molecule,format +'les' + options) 340 elif format == 'inchikey': 341 out = chemaxon.formats.MolExporter.exportToFormat(self.Molecule,'inchikey').replace('InChIKey=', '') 342 else: 343 out = chemaxon.formats.MolExporter.exportToFormat(self.Molecule,format + options) 344 if format == 'inchi': 345 out = out.split('AuxInfo=')[0] 346 if filename: 347 output = open(filename, "w") 348 print >> output, out 349 output.close() 350 return 351 else: 352 return out

353 354

355 - def calcfp(self, fp="ecfp"):

356 """Calculate a molecular fingerprint. 357 358 Optional parameters: 359 fptype -- the fingerprint type (default is "daylight"). See the 360 fps variable for a list of of available fingerprint 361 types. 362 """ 363 fp = fp.lower() 364 if fp in fps: 365 if fp == 'ecfp': 366 fp = chemaxon.descriptors.ECFP(ECFPConfiguration) 367 fp.generate(self.Molecule) 368 else: 369 raise ValueError, "%s is not a recognised fingerprint type" % fp 370 return Fingerprint(fp)

371

372 - def calcdesc(self, descnames=[]):

373 """Calculate descriptor values. 374 375 Optional parameter: 376 descnames -- a list of names of descriptors 377 378 If descnames is not specified, all available descriptors are 379 calculated. See the descs variable for a list of available 380 descriptors. 381 """ 382 if not descnames: 383 descnames = descs 384 ans = {} 385 for descname in descnames: 386 if descname not in descs: 387 raise ValueError, "%s is not a recognised descriptor type" % descname 388 if descname == 'RotatableBondsCount': 389 ta = chemaxon.calculations.TopologyAnalyser() 390 ta.setMolecule(self.Molecule) 391 ans[descname] = ta.rotatableBondCount() 392 else: 393 desc = getattr(chemaxon.descriptors.scalars, descname)('') 394 desc.generate(self.Molecule) 395 ans[descname] = desc.toFloatArray()[0] 396 return ans

397

398 - def make3D(self):

399 """Generate 3D coordinates. 400 401 Hydrogens are added, and a low energy conformer is found 402 using the MMFF94 forcefield. 403 """ 404 self.addh() 405 cp = chemaxon.marvin.calculations.ConformerPlugin() 406 cp.setMolecule(self.Molecule) 407 cp.setLowestEnergyConformerCalculation(True) 408 cp.setMMFF94Optimization(True) 409 success = cp.run() 410 optmol = cp.getMMFF94OptimizedStrucutre() 411 self.Molecule = optmol 412 self.MolHandler = chemaxon.util.MolHandler(self.Molecule) 413 self.MolHandler.aromatize()

414

415 - def draw(self, show=True, filename=None, update=False, 416 usecoords=False):

417 """Create a 2D depiction of the molecule. 418 """ 419 if not usecoords: 420 molecule = self.Molecule.clone() 421 molecule.setDim(0) 422 else: 423 molecule = self.Molecule 424 if update: 425 myMolecule = readstring("mol", Molecule(molecule).write("mol")) 426 self.Molecule = myMolecule.Molecule 427 self.MolHandler = myMolecule.MolHandler 428 bytearray = chemaxon.formats.MolExporter.exportToBinFormat(molecule, 'png') 429 if filename: 430 of = java.io.FileOutputStream(filename) 431 of.write(bytearray) 432 of.close() 433 if show: 434 source = java.io.ByteArrayInputStream(bytearray) 435 reader = javax.imageio.ImageIO.getImageReadersByFormatName('png').next() 436 iis = javax.imageio.ImageIO.createImageInputStream(source) 437 reader.setInput(iis, True) 438 param = reader.getDefaultReadParam() 439 image = reader.read(0, param) 440 frame = javax.swing.JFrame() 441 imageIcon = javax.swing.ImageIcon(image) 442 label = javax.swing.JLabel() 443 label.setIcon(imageIcon) 444 frame.getContentPane().add(label, java.awt.BorderLayout.CENTER) 445 frame.pack() 446 frame.setVisible(True) 447 frame.show()

448

449 450 -class Fingerprint(object):

451 """A Molecular Fingerprint. 452 453 Required parameters: 454 fingerprint -- a vector calculated by one of the fingerprint methods 455 456 Attributes: 457 fp -- the underlying fingerprint object 458 bits -- a list of bits set in the Fingerprint 459 460 Methods: 461 The "|" operator can be used to calculate the Tanimoto coeff. For example, 462 given two Fingerprints 'a', and 'b', the Tanimoto coefficient is given by: 463 tanimoto = a | b 464 """

465 - def __init__(self, fingerprint):

466 self.fp = fingerprint

467 - def __or__(self, other):

468 return 1 - self.fp.getTanimoto(other.fp)

469 - def __getattr__(self, attr):

470 if attr == "bits": 471 # Create a bits attribute on-the-fly 472 bs = self.fp.toBitSet() 473 bits = [-1] 474 while True: 475 setbit = bs.nextSetBit(bits[-1] + 1) 476 if setbit == -1: 477 break 478 bits.append(setbit) 479 return bits[1:] # Leave out the initial '-1' 480 else: 481 raise AttributeError, "Fingerprint has no attribute %s" % attr

482 - def __str__(self):

483 return ", ".join([str(x) for x in self.fp.toIntArray()])

484

485 -class Atom(object):

486 """Represent an Atom. 487 488 Required parameters: 489 Atom -- a JChem Atom 490 491 Attributes: 492 atomicnum, coords, formalcharge 493 494 The original JChem Atom can be accessed using the attribute: 495 Atom 496 """ 497

498 - def __init__(self, Atom):

499 self.Atom = Atom

500 501 @property

502 - def atomicnum(self): return self.Atom.getAtno()

503 @property

504 - def coords(self):

505 return (self.Atom.x, self.Atom.y, self.Atom.z)

506 @property

507 - def formalcharge(self):

508 return self.Atom.charge

509

510 - def __str__(self):

511 c = self.coords 512 return "Atom: %d (%.2f %.2f %.2f)" % (self.atomicnum, c[0], c[1], c[2])

513

514 -class Smarts(object):

515 """A Smarts Pattern Matcher 516 517 Required parameters: 518 smartspattern 519 520 Methods: 521 findall() 522 523 Example: 524 >>> mol = readstring("smi","CCN(CC)CC") # triethylamine 525 >>> smarts = Smarts("[#6][#6]") # Matches an ethyl group 526 >>> print smarts.findall(mol) 527 [(1, 2), (4, 5), (6, 7)] 528 """

529 - def __init__(self, smartspattern):

530 """Initialise with a SMARTS pattern.""" 531 self.search = chemaxon.sss.search.MolSearch() 532 smarts = MolHandler(smartspattern) 533 smarts.setQueryMode(True) 534 smarts.aromatize() 535 self.search.setQuery(smarts.molecule)

536

537 - def findall(self, molecule):

538 """Find all matches of the SMARTS pattern to a particular molecule. 539 540 Required parameters: 541 molecule 542 """ 543 self.search.setTarget(molecule.Molecule) 544 match = self.search.findAll() 545 result = [] 546 for i in xrange(len(match)): 547 result.append(tuple([n+1 for n in match[i]])) 548 return result

549

550 -class MoleculeData(object):

551 """Store molecule data in a dictionary-type object 552 553 Required parameters: 554 Molecule -- a JChem Molecule 555 556 Methods and accessor methods are like those of a dictionary except 557 that the data is retrieved on-the-fly from the underlying Molecule. 558 559 Example: 560 >>> mol = readfile("sdf", 'head.sdf').next() 561 >>> data = mol.data 562 >>> print data 563 {'Comment': 'CORINA 2.61 0041 25.10.2001', 'NSC': '1'} 564 >>> print len(data), data.keys(), data.has_key("NSC") 565 2 ['Comment', 'NSC'] True 566 >>> print data['Comment'] 567 CORINA 2.61 0041 25.10.2001 568 >>> data['Comment'] = 'This is a new comment' 569 >>> for k,v in data.iteritems(): 570 ... print k, "-->", v 571 Comment --> This is a new comment 572 NSC --> 1 573 >>> del data['NSC'] 574 >>> print len(data), data.keys(), data.has_key("NSC") 575 1 ['Comment'] False 576 """

577 - def __init__(self, Molecule):

578 self._data = Molecule.Molecule.properties()

579 - def _testforkey(self, key):

580 if not key in self: 581 raise KeyError, "'%s'" % key

582 - def keys(self):

583 return list(self._data.keys)

584 - def values(self):

585 return [self[k] for k in self._data.keys]

586 - def items(self):

587 return [(k, self[k]) for k in self._data.keys]

588 - def __iter__(self):

589 return iter(self.keys())

590 - def iteritems(self):

591 return iter(self.items())

592 - def __len__(self):

593 return len(self._data.keys)

594 - def __contains__(self, key):

595 return key in self.keys()

596 - def __delitem__(self, key):

597 self._testforkey(key) 598 self._data.setString(key, None)

599 - def clear(self):

600 for key in self: 601 del self[key]

602 - def has_key(self, key):

603 return key in self

604 - def update(self, dictionary):

605 for k, v in dictionary.iteritems(): 606 self[k] = v

607 - def __getitem__(self, key):

608 self._testforkey(key) 609 return self._data.get(key).propValue

610 - def __setitem__(self, key, value):

611 self._data.setString(key, str(value))

612 - def __repr__(self):

613 return dict(self.iteritems()).__repr__()

614 615 ECFPConfiguration = """<?xml version="1.0" encoding="UTF-8"?> 616 <ECFPConfiguration Version="0.1"> 617 618 <Parameters Length="1024" Diameter="4" Counts="no"/> 619 620 <IdentifierConfiguration> 621  622 <Property Name="AtomicNumber" Value="1"/> 623 <Property Name="HeavyNeighborCount" Value="1"/> 624 <Property Name="HCount" Value="1"/> 625 <Property Name="FormalCharge" Value="1"/> 626 <Property Name="IsRingAtom" Value="1"/> 627 628  629 <Property Name="ConnectionCount" Value="0"/> 630 <Property Name="Valence" Value="0"/> 631 <Property Name="Mass" Value="0"/> 632 <Property Name="MassNumber" Value="0"/> 633 <Property Name="HasAromaticBond" Value="0"/> 634 <Property Name="IsTerminalAtom" Value="0"/> 635 <Property Name="IsStereoAtom" Value="0"/> 636 </IdentifierConfiguration> 637 638 <StandardizerConfiguration Version="0.1"> 639 <Actions> 640 <Action ID="aromatize" Act="aromatize"/> 641 <RemoveExplicitH ID="RemoveExplicitH" Groups="target"/> 642 </Actions> 643 </StandardizerConfiguration> 644 645 <ScreeningConfiguration> 646 <ParametrizedMetrics> 647 <ParametrizedMetric Name="Tanimoto" ActiveFamily="Generic" Metric="Tanimoto" Threshold="0.5"/> 648 <ParametrizedMetric Name="Euclidean" ActiveFamily="Generic" Metric="Euclidean" Threshold="10"/> 649 </ParametrizedMetrics> 650 </ScreeningConfiguration> 651 652 </ECFPConfiguration> 653 """ 654 655 if __name__=="__main__": #pragma: no cover 656 mol = readstring("smi", "CC(=O)Cl") 657 mol.title = u"Adrià" 658 mol.draw() 659 660 for mol in readfile("sdf", "head.sdf"): 661 pass 662

Source Code for Module cinfony.jchem