Package AutoDockTools :: Module MoleculePreparation
[hide private]
[frames] | no frames]

Source Code for Module AutoDockTools.MoleculePreparation

   1  #!/usr/bin/env python2.3 
   2  # 
   3  # 
   4  # $Header: /opt/cvs/python/packages/share1.5/AutoDockTools/MoleculePreparation.py,v 1.49.2.3 2007/10/10 18:21:07 rhuey Exp $  
   5  # 
   6  # 
   7  # 
   8  import os, string 
   9  import Numeric, math  
  10   
  11  from MolKit import Read 
  12  from MolKit.torTree import TorTree 
  13  from MolKit.pdbWriter import PdbqWriter, PdbqsWriter, PdbqtWriter 
  14  import MolKit.molecule 
  15  import MolKit.protein 
  16  from MolKit.protein import ResidueSet 
  17   
  18  from MolKit.molecule import AtomSet, Bond, BondSet 
  19  from MolKit.chargeCalculator import GasteigerChargeCalculator  
  20  from MolKit.chargeCalculator import KollmanChargeCalculator 
  21  from MolKit.hydrogenBuilder import HydrogenBuilder 
  22  #import the Kollman charges dictionary to get keys it recognizes... 
  23  from MolKit.data.qkollua import q 
  24  from MolKit.bondSelector import RotatableBondSelector, AmideBondSelector, LeafBondSelector 
  25  from MolKit.bondSelector import GuanidiniumBondSelector 
  26   
  27   
  28   
  29  from AutoDockTools.atomTypeTools import NonpolarHydrogenMerger 
  30  from AutoDockTools.atomTypeTools import LonepairMerger 
  31  from AutoDockTools.atomTypeTools import SolvationParameterizer 
  32  from AutoDockTools.atomTypeTools import AromaticCarbonManager 
  33  from AutoDockTools.atomTypeTools import AutoDock4_AtomTyper 
  34   
  35  from AutoDockTools.AutoDockBondClassifier import AutoDockBondClassifier 
  36   
  37   
  38   
39 -class AutoDockMoleculePreparation:
40 """ 41 Base class for facades for preparing MolKit molecules to be used in an 42 AutoDock experiment. 43 44 This preparation involves: 45 1. optional cleanup 46 2. making atoms conform to AD atom format: 47 adding charges 48 merging nonpolarHydrogens (can be overridden) 49 merging lonepairsHydrogens (can be overridden) 50 adding solvation parameters 51 [optional: 52 adding Hydrogens] 53 3. writing an pdbq[s,t] outputfile 54 55 56 Collaborators: 57 in MolKit: 58 for adding charges: 59 KollmanChargeCalculator 60 [optional: GasteigerChargeCalculator] 61 for adding hydrogens: 62 HydrogenBuilder 63 in AutoDockTools: 64 atomTypeManager classes for conforming to AD atomtypes: 65 LonepairMerger 66 NonpolarHydrogenMerger 67 68 69 Preparation of a molecule 70 buildBonds if necessary 71 adds 'autodock_element field' to all Atoms 72 detects whether molecule isPeptide or not (for charge type to use) 73 74 75 Methods to be shared by derived classes 76 calcChargeError: calculation of totalCharge error 77 detectIsPeptide: returns boolean isPeptide 78 findNearest: calculation of atom nearest another (non-bonded) atom 79 """ 80 81 82 std_types = ['CYS', 'ILE', 'SER', 'VAL', 'GLN', 'LYS', 'ASN', 'PRO', 'THR', 'PHE', 'ALA', 'HIS', 'GLY', 'ASP', 'LEU', 'ARG', 'TRP', 'GLU', 'TYR', 'MET', 'HID', 'HSP', 'HIE', 'HIP', 'CYX', 'CSS'] 83 #what about nucleic acids??? 84
85 - def __init__(self, molecule, mode='', repairs='', 86 charges_to_add=None, cleanup='', 87 outputfilename=None, debug=False, 88 version=3, delete_single_nonstd_residues=False):
89 #debug = True 90 self.debug = debug 91 self.molecule = mol = molecule 92 self.mode = mode 93 self.lenNPHS = 0 94 #fix atom elements if necessary: 95 self.autodock_element_dict = {'Cl':'c', 'Fe': 'f', 'Br':'b'} 96 self.reverse_autodock_elements = {'c':'Cl', 'f':'Fe','b':'Br'} 97 self.version = version 98 self.delete_single_nonstd_residues = delete_single_nonstd_residues 99 if version==3: 100 revlist = ['c','f','b'] 101 reversedAts = AtomSet(filter(lambda x: x.element in revlist, mol.allAtoms)) 102 for at in reversedAts: 103 if debug: print 'restoring ', at.element, ' to ', 104 at.autodock_element = at.element 105 at.element = self.reverse_autodock_elements[at.element] 106 if debug: print at.element, ' for processing...' 107 # build bonds if necessary #ALWAYS DO THIS!!! 108 if not len(mol.allAtoms.bonds[0]): #what if not enough bonds 109 mol.buildBondsByDistance() 110 111 #PROCESS molecule 112 self.cleanup_type_list = string.split(cleanup, '_') 113 #if mode=='automatic': #DO this ANYWAY??? 114 #remove waters + chains composed of nonstd residues 115 self.cleanUpResidues(mol, self.cleanup_type_list) 116 117 # REPAIR: possibly add bonds to non-bonded atoms + add hydrogens 118 self.repair_type_list = string.split(repairs, '_') 119 #if mode=='automatic': #SHOULD THIS BE DONE ANYWAY??? 120 self.repairMol(mol, self.repair_type_list) 121 122 # CHARGES: ALWAYS add charges last after atoms are set but 123 # before merging nphs and lps 124 self.chargeType = charges_to_add 125 self.chargeError = 'ERROR' 126 if debug: print "charges_to_add=", charges_to_add 127 #if mode=='automatic': #SHOULD THIS BE DONE ANYWAY??? 128 self.addCharges(mol, charges_to_add) 129 if self.chargeType!='ERROR': 130 self.chargeError = self.calcChargeError() 131 #need a two step cleanup to do merges after adding bonds/hydrogens 132 #and AFTER adding charges 133 # CLEANUP: possibly merge nonpolarhydrogens and/or lonepairs 134 self.cleanUpAtomTypes(mol, self.cleanup_type_list) 135 #set autodock_element for all atoms 136 mol.allAtoms.autodock_element = mol.allAtoms.element 137 self.outputfilename = outputfilename 138 if debug: print "end of base class init"
139 140
141 - def repairMol(self, mol, repair_list):
142 #@@ what about metals? 143 # possibly check for non-bonded atoms and add bonds to them 144 # possibly add hydrogens 145 if 'bonds' in repair_list: 146 non_bonded_atoms = mol.allAtoms.get(lambda x: len(x.bonds)==0) 147 if non_bonded_atoms: 148 if self.debug: print "!!adding bonds to ", non_bonded_atoms.name, "!!" 149 for a in non_bonded_atoms: 150 nearestAt = self.findNearest(a, mol.allAtoms) 151 Bond(a, nearestAt, bondOrder=1, origin='AddedBond', check=0, 152 addBond=1) 153 if self.debug: print 'added bond:', a.bonds[0] 154 repair_list.remove('bonds') 155 # possibly add hydrogens 156 #MODE SWITCH 1: adding hydrogens ||IN USE|| 157 self.newHs = 0 158 hs = mol.allAtoms.get(lambda x: x.element=='H') 159 if 'hydrogens' in repair_list: 160 if self.debug: print "addinghydrogens!" 161 self.newHs = self.addHydrogens(mol) 162 if self.debug: print "added ", self.newHs, " hydrogens" 163 repair_list.remove('hydrogens') 164 elif not hs and 'checkhydrogens' in self.repair_type_list: 165 if self.debug: print "addinghydrogens from check hydrogens" 166 self.newHs = self.addHydrogens(mol) 167 if self.debug: print "added ", self.newHs, " hydrogens" 168 repair_list.remove('checkhydrogens')
169
170 - def addHydrogens(self, mol):
171 beforeLen = len(mol.allAtoms) 172 #NB: this adds all hydrogens (?!?) 173 HB = HydrogenBuilder() 174 HB.addHydrogens(mol) 175 afterLen = len(mol.allAtoms) 176 newHs = afterLen - beforeLen 177 #NB this adds the HYDROGENS as "ATOM" not HETATM 178 return newHs
179 180
181 - def addCharges(self, mol, charges_to_add):
182 #detect whether isPeptide first 183 debug = self.debug 184 if debug: 185 print "in addCharges:mol.name=", mol.name, " and charges_to_add=", charges_to_add 186 isPeptide = mol.isPeptide = self.detectIsPeptide() 187 if charges_to_add=='gasteiger': 188 chargeCalculator = self.chargeCalculator = GasteigerChargeCalculator() 189 if isPeptide: 190 print "adding gasteiger charges to peptide" 191 chargeCalculator = self.chargeCalculator = GasteigerChargeCalculator() 192 elif charges_to_add=='Kollman': 193 if not mol.isPeptide: 194 print "adding Kollman charges to non-peptide" 195 chargeCalculator = self.chargeCalculator = KollmanChargeCalculator() 196 #KEEP charges from file, if there are any 197 #MODE SWITCH 2: adding charges???? ||NOT IN USE|| 198 if debug: print "self.chargeType=", self.chargeType 199 if self.chargeType is None: 200 len_zero_charges = len(filter(lambda x: x.charge==0, 201 mol.chains.residues.atoms)) 202 if len_zero_charges==len(mol.allAtoms): 203 if debug: print "WARNING: all atoms had zero charges!" 204 chargeCalculator.addCharges(mol.allAtoms) 205 self.chargeType = mol.allAtoms[0].chargeSet 206 else: 207 len_charged = len(filter(lambda x:x.chargeSet!=None, 208 mol.chains.residues.atoms)) 209 if len_charged!=len(mol.chains.residues.atoms): 210 #charges could come from file or ? 211 if debug: print "WARNING: some atoms missing charges!" 212 self.chargeType = "ERROR" 213 if debug: print "no change to charges!" 214 self.chargeType = mol.allAtoms[0].chargeSet 215 else: 216 chargeCalculator.addCharges(mol.allAtoms) 217 if debug: print 'added ' + charges_to_add + 'charges to ', mol.name
218 219
220 - def cleanUpAtomTypes(self, mol, cleanup_list):
221 debug = self.debug 222 merge_nonpolar_hydrogens = 'nphs' in cleanup_list 223 if merge_nonpolar_hydrogens: 224 NPHM = self.NPHM = NonpolarHydrogenMerger() 225 lenNPHS = self.lenNPHS = NPHM.mergeNPHS(mol.allAtoms) 226 if debug: print 'merged ', lenNPHS, ' nonpolar hydrogens' 227 cleanup_list.remove('nphs') 228 merge_lonepairs = 'lps' in cleanup_list 229 if merge_lonepairs: 230 LPM = self.LPM = LonepairMerger() 231 #lps = filter(lambda x: x.element=='Xx' and \ 232 #(x.name[0]=='L' or x.name[1]=='L'), mol.allAtoms) 233 lenLPS = self.lenLPS = LPM.mergeLPS(mol.allAtoms) 234 if debug: print 'merged ', lenLPS, ' lonepairs' 235 cleanup_list.remove('lps') 236 if self.delete_single_nonstd_residues: 237 #delete non-standard residues from within a single chain 238 non_std = [] 239 for res in mol.chains.residues: 240 if res.type not in self.std_types: 241 non_std.append(res) 242 #delete non_std 243 names = "" 244 #remove all bonds to other atoms: 245 for r in non_std: 246 for a in r.atoms: 247 for b in a.bonds: 248 at2 = b.atom1 249 if at2==a: at2 = b.atom2 250 at2.bonds.remove(b) 251 a.bonds.remove(b) 252 names = names + r.parent.id +r.name + '_' 253 print "\'Deleting non-standard residues:" + names + " from " + mol.name 254 for res in non_std: 255 res.parent.remove(res) 256 del res 257 #fix allAtoms short cut 258 mol.allAtoms = mol.chains.residues.atoms 259 delete_alternateB = 'deleteAltB' in cleanup_list 260 if delete_alternateB: 261 if self.debug: print "deleting ", len(mol.allAtoms), " alternateB positions" 262 altB_ats = mol.allAtoms.get('*@B') 263 if not len(altB_ats): 264 return 265 for a in altB_ats: 266 res = a.parent 267 for b in a.bonds: 268 at2 = b.atom1 269 if at2==a: at2 = b.atom2 270 at2.bonds.remove(b) 271 a.bonds.remove(b) 272 res.remove(a) 273 del a 274 #remove any remaining @A 275 altA_ats = mol.allAtoms.get('*@A') 276 for a in altA_ats: 277 a.name = a.name[:-2] 278 mol.allAtoms = mol.chains.residues.atoms 279 if self.debug: print "resetting numbers to 1-", len(mol.allAtoms)+1 280 mol.allAtoms.number = range(1, len(mol.allAtoms)+1)
281 282 283
284 - def cleanUpResidues(self, mol, cleanup_list):
285 self.lenHOHS = 0 286 remove_waters = 'waters' in cleanup_list 287 if remove_waters: 288 hohs = mol.allAtoms.get(lambda x: x.parent.type=='HOH') 289 #this line added at Alex Perryman's suggestion: 290 hohs2 = mol.allAtoms.get(lambda x: x.parent.type=='WAT') 291 hohs = hohs + hohs2 292 if hohs: 293 #remove(hohs) 294 lenHOHS = self.lenHOHS = len(hohs) 295 for h in hohs: 296 for b in h.bonds: 297 c = b.atom1 298 if c==h: 299 c = b.atom2 300 c.bonds.remove(b) 301 h.bonds = BondSet() 302 res = h.parent 303 h.parent.remove(h) 304 if len(h.parent.children)==0: 305 res = h.parent 306 chain = res.parent 307 if self.debug: print "removing residue: ", res.name 308 chain.remove(res) 309 if len(chain.children)==0: 310 mol = chain.parent 311 if self.debug: print "removing chain", chain.id 312 mol.remove(chain) 313 del chain 314 del res 315 del h 316 #fix allAtoms short cut 317 mol.allAtoms = mol.chains.residues.atoms 318 cleanup_list.remove('waters') 319 320 remove_nonstdres = 'nonstdres' in cleanup_list 321 self.lenChainsRemoved = 0 322 if remove_nonstdres: 323 if len(mol.chains)>1: 324 chains_to_delete = [] 325 for c in mol.chains: 326 num_res = len(c.residues) 327 non_std = [] 328 for res in c.residues: 329 if res.type not in self.std_types: 330 non_std.append(res) 331 if len(c.residues)==len(non_std): 332 print "\'"+ c.name, "\' apparently composed of not std residues. Deleting " 333 chains_to_delete.append(c) 334 if len(chains_to_delete): 335 self.lenChainsRemoved = len(chains_to_delete) 336 for c in chains_to_delete: 337 #remove(c) #NB this would remove HOH and ligand from 1HSG 338 c.parent.remove(c) 339 del c 340 #fix allAtoms short cut 341 mol.allAtoms = mol.chains.residues.atoms 342 #if self.delete_single_nonstd_residues: 343 # #delete non-standard residues from within a single chain 344 # non_std = [] 345 # for res in mol.chains.residues: 346 # if res.type not in self.std_types: 347 # non_std.append(res) 348 # #delete non_std 349 # names = "" 350 # #remove all bonds to other atoms: 351 # for r in non_std: 352 # for a in r.atoms: 353 # for b in a.bonds: 354 # at2 = b.atom1 355 # if at2==a: at2 = b.atom2 356 # at2.bonds.remove(b) 357 # a.bonds.remove(b) 358 # names = names + r.parent.id +r.name + '_' 359 # print "\'Deleting non-standard residues:" + names + " from " + mol.name 360 # for res in non_std: 361 # res.parent.remove(res) 362 # del res 363 # #fix allAtoms short cut 364 # mol.allAtoms = mol.chains.residues.atoms 365 cleanup_list.remove('nonstdres')
366 367
368 - def calcChargeError(self):
369 s = Numeric.add.reduce(self.molecule.allAtoms.charge) 370 return min(math.ceil(s)-s, s-math.floor(s))
371 372
373 - def detectIsPeptide(self):
374 isPeptide = True 375 d = {} 376 for r in self.molecule.chains.residues: 377 d[r.type] = 0 378 for t in d.keys(): 379 if t not in q.keys(): 380 isPeptide = False 381 break 382 return isPeptide
383 384
385 - def findNearest(self, atom, bonded_atoms):
386 lenK = len(bonded_atoms) 387 lenC = 1 388 nonbonded_atoms = AtomSet([atom]) 389 c = Numeric.array(nonbonded_atoms.coords, 'f') 390 k = Numeric.array(bonded_atoms.coords, 'f') 391 bigK = Numeric.resize(k, (lenK, lenK, 3)) 392 c.shape = (1, 1, 3) 393 bigM = bigK[:1] 394 d = bigM - c 395 dSQ = d*d 396 dSQMAT = Numeric.sum(dSQ, 2) 397 mm = dSQMAT[0][0] 398 mindex = 0 399 for i in range(lenK): 400 if dSQMAT[0][i]<mm: 401 mm = dSQMAT[0][i] 402 mindex = i 403 #print "found closest atom %d at dist %f" %( mindex, mm) 404 return bonded_atoms[mindex]
405 406 407
408 -class ReceptorPreparation(AutoDockMoleculePreparation):
409 """ 410 Facade for preparing a MolKit molecule to be used as the receptor in an 411 AutoDock experiment derived from AutoDockMoleculePreparation class. 412 413 414 Receptor preparation involves: 415 1. adding solvation parameters 416 2. writing a pdbqs outputfile 417 418 Receptor preparation Collaborators extend base class collaborators: 419 in AutoDockTools: 420 in atomTypeManager classes for conforming to AD atomtypes: 421 SolvationParameterizer 422 in this file: 423 ReceptorWriter 424 425 """ 426
427 - def __init__(self, molecule, mode='automatic', repairs='checkhydrogens', 428 charges_to_add='Kollman', 429 cleanup='nphs_lps_waters_nonstdres', 430 outputfilename=None, debug=False, 431 preserved={}, 432 delete_single_nonstd_residues=False):
433 434 AutoDockMoleculePreparation.__init__(self, molecule, mode, repairs, 435 charges_to_add, cleanup, 436 debug=debug, 437 delete_single_nonstd_residues=False) 438 439 if len(preserved): 440 for atom, chargeList in preserved.items(): 441 atom._charges[chargeList[0]] = chargeList[1] 442 atom.chargeSet = chargeList[0] 443 if debug: print 'preserved charge on ', atom.name, ' charge=', atom.charge 444 445 # add solvation parameters here: 446 SP = self.SP = SolvationParameterizer() 447 unknownatoms = SP.addParameters(molecule.chains.residues.atoms) 448 if len(unknownatoms) and debug: 449 print len(unknownatoms), " unknown atoms: set solvation parameters to 0. 0." 450 451 self.writer = ReceptorWriter() 452 #MODE SWITCH 5: write outputfile ||IN USE|| 453 if mode=='automatic': 454 #write outputfile now without waiting ... 455 #fix this: what records should be written? 456 self.write(outputfilename) 457 self.outputfilename = outputfilename
458 459
460 - def summarize(self):
461 mol = self.molecule 462 msg = "setup "+ mol.name+ ":\n" 463 if self.newHs!=0: 464 msg = msg + " added %d new hydrogen(s)\n" %self.newHs 465 if self.chargeType!=None: 466 if self.chargeType in ['pdbq', 'mol2', 'pdbqt', 'pdbqs']: 467 msg = msg + " kept charges from " + self.chargeType + " file\n" 468 else: 469 msg = msg + " added %s charges\n" %self.chargeType 470 if self.lenNPHS: 471 msg = msg + " merged %d non-polar hydrogens\n" %self.lenNPHS 472 if self.lenLPS: 473 msg = msg + " merged %d lone pairs\n" %self.lenLPS 474 totalCh = Numeric.add.reduce(mol.allAtoms.charge) 475 if hasattr(self, 'chargeError') and abs(self.chargeError)>0.000005: 476 msg = msg + " total charge error = %6.4f\n"%self.chargeError 477 return msg
478 479
480 - def write(self, outputfilename=None):
481 #write to outputfile 482 mol = self.molecule 483 if not outputfilename and not self.outputfilename: 484 outputfilename = self.outputfilename = mol.name + ".pdbqs" 485 if self.debug: print 'using std outputfilename ', outputfilename 486 #should this be done at Molecule level or at Atom level? 487 #Atom level would allow multiple molecules in Receptor, at no price 488 self.writer.write(mol, outputfilename) 489 if self.debug: print "wrote ", mol.name, " to ", outputfilename
490 491 492
493 -class AD4ReceptorPreparation(AutoDockMoleculePreparation):
494 """ 495 Facade for preparing a MolKit molecule to be used as the receptor in an 496 AutoDock4 experiment derived from AutoDockMoleculePreparation class. 497 498 Receptor preparation involves: 499 1. adding autodock_element types 500 2. writing a pdbqt outputfile 501 502 Receptor preparation Collaborators extend base class collaborators: 503 in AutoDockTools: 504 in atomTypeManager classes for conforming to AD atomtypes: 505 AutoDock4_AtomTyper 506 in this file: 507 AD4ReceptorWriter 508 509 """ 510
511 - def __init__(self, molecule, mode='automatic', repairs='checkhydrogens', 512 charges_to_add='gasteiger', 513 cleanup='nphs_lps_waters_nonstdres', 514 outputfilename=None, debug=False, 515 version=4, preserved={}, 516 delete_single_nonstd_residues=False):
517 518 519 AutoDockMoleculePreparation.__init__(self, molecule, 520 mode=mode, repairs=repairs, 521 charges_to_add=charges_to_add, cleanup=cleanup, 522 outputfilename=outputfilename, debug=debug, 523 version=version, delete_single_nonstd_residues=delete_single_nonstd_residues) 524 525 #NEED TO TYPE ATOMS!! 526 try: 527 delattr(mol.allAtoms, 'autodock_element') 528 except: 529 pass 530 atomTyper = AutoDock4_AtomTyper() 531 atomTyper.setAutoDockElements(molecule, reassign=True) #4/15/05 catch here? 532 if len(preserved): 533 for atom, chargeList in preserved.items(): 534 atom._charges[chargeList[0]] = chargeList[1] 535 atom.chargeSet = chargeList[0] 536 if debug: print 'preserved charge on ', atom.name, ' charge=', atom.charge 537 538 self.writer = AD4ReceptorWriter() 539 #MODE SWITCH 5: write outputfile ||IN USE|| 540 if mode=='automatic': 541 #write outputfile now without waiting ... 542 #fix this: what records should be written? 543 self.write(outputfilename) 544 self.outputfilename = outputfilename
545 546
547 - def write(self, outputfilename=None):
548 #write to outputfile 549 mol = self.molecule 550 if not outputfilename and not self.outputfilename: 551 outputfilename = self.outputfilename = mol.name + ".pdbqt" 552 if self.debug: print 'using std outputfilename ', outputfilename 553 #should this be done at Molecule level or at Atom level? 554 #Atom level would allow multiple molecules in Receptor, at no price 555 self.writer.write(mol, outputfilename) 556 if self.debug: print "wrote ", mol.name, " to ", outputfilename
557 558 559 560
561 -class ReceptorWriter:
562 """ 563 ReceptorWriter class writes a receptor molecule which has a ReceptorPreparationObject to a file 564 """ 565
566 - def __init__(self, write_CONECT=False):
567 self.writer = PdbqsWriter() 568 #each atom must have charge and AtVol and AtSolPar 569 self.conditions = [lambda x: hasattr(x, 'AtVol'), \ 570 lambda x: hasattr(x, 'AtSolPar'), \ 571 lambda x: x.chargeSet is not None] 572 self.write_CONECT = write_CONECT
573 574
575 - def write(self, receptor, outputfile):
576 #should this be done at Molecule level or at Atom level? 577 receptor_atoms = receptor.chains.residues.atoms 578 len_receptor_atoms = len(receptor_atoms) 579 580 ##check each condition for receptor_atoms 581 for cond in self.conditions: 582 assert len(filter(cond, receptor_atoms))==len_receptor_atoms 583 584 outptr = open(outputfile, 'w') 585 records = ['ATOM'] 586 #if self.write_CONECT: 587 # records.append('CONECT') 588 #for at in receptor_atoms: 589 # self.writer.write_atom(outptr, at) 590 for ch in receptor.chains: 591 for at in ch.residues.atoms: 592 self.writer.write_atom(outptr, at) 593 rec = self.writer.defineTERRecord(at) 594 outptr.write(rec) 595 #put in the TER and END lines??? 596 #optional CONECT records??? 597 if self.write_CONECT: 598 self.writeCONECTRecords(receptor, outptr) 599 outptr.close()
600 601
602 - def writeCONECTRecords(self, fptr, mol):
603 atms = mol.allAtoms 604 for atm in atms: 605 rec = 'CONECT%5i'%atm.number 606 for b in atm.bonds: 607 a2 = b.atom1 608 if a2==atm: a2 = b.atom2 609 if not a2 in atms: continue #don't write intermolecular bonds 610 rec = rec + '%5i'%a2.number 611 rec = rec + '\n' 612 fptr.write(rec)
613 614 615 616
617 -class AD4ReceptorWriter(ReceptorWriter):
618 """ 619 AD4ReceptorWriter class writes a receptor molecule which has a ReceptorPreparationObject to a pdbqt file 620 """ 621
622 - def __init__(self):
623 ReceptorWriter.__init__(self) 624 self.writer = PdbqtWriter() 625 self.conditions = [lambda x: hasattr(x, 'autodock_element'), \ 626 lambda x: x.chargeSet is not None]
627 628 629
630 -class LigandPreparation(AutoDockMoleculePreparation):
631 632 """ 633 634 Facade for preparing a MolKit molecule to be used as the ligand in an 635 AutoDock experiment. 636 637 LigandPreparation involves: 638 1. managing definition of flexibility pattern in ligand 639 using metaphor of a Tree 640 +setting ROOT 641 +setting pattern of rotatable bonds 642 +setting TORSDOF, number of possible rotatable bonds- 643 number of bonds which only rotate hydrogens 644 optional: 645 -disallowing amide torsions 646 -disallowing peptidebackbone torsions 647 -disallowing guanidinium torsions 648 +/- toggling activity of any rotatable bond 649 2. writing an outputfile with keywords recognized by AutoDock 650 3. writing types list and number of active torsions to a dictionary file 651 652 LigandPreparation Collaborators: 653 in AutoDockTools: 654 for rotability pattern management: 655 RotatableBondManager 656 [which uses an AutoDockTools/AutoDockBondClassifier 657 for detecting possible flexibility patterns in bonds] 658 for managing cyclic aromatic carbons: 659 AromaticCarbonManager 660 (including changing cutoff from 7.5) 661 662 """ 663
664 - def __init__(self, molecule, mode='automatic', repairs='hydrogens_bonds', 665 charges_to_add='gasteiger', cleanup='nphs_lps', 666 allowed_bonds='backbone', #backbone on, amide + guanidinium off 667 root='auto', outputfilename=None, dict=None, debug=False, 668 check_for_fragments=False, bonds_to_inactivate=[], 669 inactivate_all_torsions=False, 670 version=3, limit_torsions=False, 671 delete_single_nonstd_residues=False):
672 # why aren't repairs, cleanup and allowed_bonds lists?? 673 # to run tests: use allowed_bonds = 'guanidinium' 674 if debug: print "LPO: allowed_bonds=", allowed_bonds 675 676 AutoDockMoleculePreparation.__init__(self, molecule, mode, repairs, 677 charges_to_add, cleanup, 678 debug=debug, version=version, 679 delete_single_nonstd_residues=False) 680 681 self.version = version 682 self.charge_error_tolerance = 0.000005 683 #FOR FLEXIBILITY MODEL: setup RotatableBondManager 684 #process bonds 685 #MODE SWITCH 4: disallow allowed_bonds ||NOT IN USE|| 686 self.RBM = RotatableBondManager(molecule, 687 allowed_bonds, root, debug=False, 688 check_for_fragments=check_for_fragments, 689 bonds_to_inactivate=bonds_to_inactivate) 690 if inactivate_all_torsions is True: 691 #in this case, make all active torsions inactive 692 self.RBM.set_all_torsions(False) 693 elif limit_torsions is not False: 694 self.RBM.limit_torsions(limit_torsions) 695 696 #detect aromatic cycle carbons and rename them for AutoDock3 697 rename = self.version==3 698 ACM = self.ACM = AromaticCarbonManager(rename=rename) 699 self.aromCs = self.ACM.setAromaticCarbons(molecule) 700 701 #optional output summary filename to write types and number of torsions 702 self.dict = dict 703 704 self.outputfilename = outputfilename 705 #if mode is 'automatic': write outputfile now 706 #without waiting to set torsions etc interactively 707 #MODE SWITCH 5: write outputfile ||IN USE|| 708 self.writer = LigandWriter() 709 if mode=='automatic': 710 self.write(outputfilename) 711 self.outputfilename = outputfilename
712 713
714 - def summarize(self):
715 mol = self.molecule 716 msg = "setup "+ mol.name+ ":\n" 717 if self.newHs!=0: 718 msg = msg + " added %d new hydrogen(s)\n" %self.newHs 719 if self.chargeType!=None: 720 if self.chargeType in ['pdbq', 'mol2']: 721 msg = msg + " kept charges from " + self.chargeType + " file\n" 722 else: 723 msg = msg + " added %s charges\n" %self.chargeType 724 if self.lenNPHS: 725 msg = msg + " merged %d non-polar hydrogens\n" %self.lenNPHS 726 if self.lenLPS: 727 msg = msg + " merged %d lone pairs\n" %self.lenLPS 728 if len(self.aromCs): 729 msg = msg + " found %d aromatic carbons\n" %len(self.aromCs) 730 totalCh = Numeric.add.reduce(mol.allAtoms.charge) 731 msg = msg + " detected %d rotatable bonds\n" %len(mol.possible_tors_bnds) 732 msg = msg + " set TORSDOF to %d\n"%mol.TORSDOF 733 #if hasattr(self, 'chargeError') and abs(self.chargeError)>0.000005: 734 if hasattr(self, 'chargeError') and abs(self.chargeError)>self.charge_error_tolerance: 735 msg = msg + " total charge error = %6.4f\n"%self.chargeError 736 return msg
737
738 - def set_autodock_element_types_for_writing(self, allAtoms):
739 sublist = self.autodock_element_dict.keys() 740 problem_ats = AtomSet(filter(lambda x: x.element in sublist, allAtoms)) 741 for at in problem_ats: 742 #there is a substitute 'name' for this atom 743 if self.debug: print "changed ", at.name, 744 newval = self.autodock_element_dict[at.element] 745 at.name = newval 746 #have to change element to get correct written output 747 #because the stupid pdbWriter writes the element not the name 748 at.element = newval 749 at.autodock_element = at.element 750 if self.debug: print "to ", at.name 751 return problem_ats
752 753
754 - def restore_autodock_element_types_after_writing(self, problem_ats):
755 for at in problem_ats: 756 #there is a substitute 'name' for this atom 757 if self.debug: print "changed back", at.name,"'s element " 758 newval = self.reverse_autodock_elements[at.element] 759 #have to change element to get correct written output 760 #because the stupid pdbWriter writes the element not the name 761 at.element = newval 762 if self.debug: print "to ", at.element
763 764
765 - def write(self, outputfilename=None):
766 #write ligand to outputfile 767 mol = self.molecule 768 if not hasattr(mol, 'ROOT'): 769 print 'must specify ROOT before writing' 770 return 'ERROR' 771 if hasattr(mol, 'torTree') and not hasattr(mol, 'torscount'): 772 mol.torscount = len(mol.torTree.torsionMap) 773 if not outputfilename: 774 stem = os.path.splitext(os.path.basename(mol.parser.filename))[0] 775 if self.version==4: 776 outputfilename = stem + ".pdbqt" 777 else: 778 outputfilename = stem + ".pdbq" 779 if self.debug: print 'using std outputfilename ', outputfilename 780 #if self.debug: print 'using std outputfilename ', outputfilename 781 #writer = LigandWriter() 782 #print "calling writer.write with ", mol.name, ' ', outputfilename 783 #change the names of Cl, Fe and Br atoms 784 problem_ats = self.set_autodock_element_types_for_writing(mol.allAtoms) 785 #if self.version==3: 786 #sublist = self.autodock_element_dict.keys() 787 #problem_ats = AtomSet(filter(lambda x: x.element in sublist, mol.allAtoms)) 788 #for at in problem_ats: 789 ##there is a substitute 'name' for this atom 790 #if self.debug: print "changed ", at.name, 791 #newval = self.autodock_element_dict[at.element] 792 #at.name = newval 793 ##have to change element to get correct written output 794 ##because the stupid pdbWriter writes the element not the name 795 #at.element = newval 796 #at.autodock_element = at.element 797 #if self.debug: print "to ", at.name 798 self.writer.write(mol, outputfilename) 799 self.outputfilename = outputfilename 800 #print 'back from writing' 801 if self.dict is not None: 802 if not os.path.exists(self.dict): 803 fptr = open(self.dict, 'a') 804 ostr = "summary = d = {}\n" 805 fptr.write(ostr) 806 else: 807 fptr = open(self.dict, 'a') 808 type_dict = {} 809 for a in self.molecule.allAtoms: 810 type_dict[a.autodock_element] = 1 811 atom_types = type_dict.keys() 812 atom_types.sort() 813 ostr = "d['" + self.molecule.name +"'] = {" 814 ostr = ostr + "'atom_types': [" 815 for t in atom_types[:-1]: 816 ostr = ostr + "'%s', "%t 817 ostr = ostr + "'%s' "%atom_types[-1] 818 ostr = ostr + "],\n\t\t\t'rbonds':" + str(self.molecule.torscount) 819 ostr = ostr + ",\n\t\t\t'zero_charge' : [" 820 zc_ats = self.molecule.allAtoms.get(lambda x: x.charge==0) 821 for at in zc_ats: 822 ostr = ostr + "'%s, '"%at.name 823 ostr = ostr + "],\n\t\t\t}\n" 824 fptr.write(ostr) 825 fptr.close() 826 self.restore_autodock_element_types_after_writing(problem_ats)
827 #if self.version==3: 828 #for at in problem_ats: 829 ##there is a substitute 'name' for this atom 830 #if self.debug: print "changed back", at.name,"'s element " 831 #newval = self.reverse_autodock_elements[at.element] 832 ##have to change element to get correct written output 833 ##because the stupid pdbWriter writes the element not the name 834 #at.element = newval 835 #if self.debug: print "to ", at.element 836 837
838 - def autoroot(self):
839 self.RBM.autoroot()
840 841
842 - def setroot(self, index):
843 self.RBM.setroot(index)
844 845
846 - def limit_torsions(self, val, type):
847 self.RBM.limit_torsions(val, type)
848 849
850 - def set_all_torsions(self, flag):
851 mol = self.molecule 852 if self.debug: print "LPO:set_all_torsions with flag=", flag 853 if self.debug: print "activeBonds=", len(filter(lambda x: x.activeTors, mol.allAtoms.bonds[0])) 854 self.RBM.set_all_torsions(flag) 855 if self.debug: print "2: activeBonds=", len(filter(lambda x: x.activeTors, mol.allAtoms.bonds[0]))
856 857
858 - def set_amide_torsions(self, flag):
859 if not len(self.molecule.amidebnds): 860 return 861 self.RBM.set_amide_torsions(flag)
862 863
864 - def set_ppbb_torsions(self, flag):
865 if not len(self.molecule.ppbbbnds): 866 return 867 self.RBM.set_peptidebackbone_torsions(flag)
868 869
870 - def set_guanidinium_torsions(self, flag):
871 if not len(self.molecule.guanidiniumbnds): 872 return 873 self.RBM.set_guanidinium_torsions(flag)
874 875
876 - def toggle_torsion(self,ind1, ind2):
877 self.RBM.toggle_torsion(ind1, ind2)
878 879
880 - def changePlanarityCriteria(self, cutoff):
881 oldAromCs = self.aromCs 882 self.aromCs = self.ACM.setAromaticCarbons(self.molecule, cutoff) 883 if self.debug: print "now there are ", len(self.aromCs), " aromaticCs"
884 885
886 - def set_carbon_names(self, atoms, type):
887 self.ACM.set_carbon_names(atoms, type) 888 #have to reset this field here 889 self.aromCs = AtomSet(filter(lambda x: (x.autodock_element=='A' and x.element=='C'), 890 self.molecule.allAtoms))
891 892
893 - def setAromaticCarbons(self):
894 self.aromCs = self.ACM.setAromaticCarbons(self.molecule) 895 return self.aromCs
896 897 #LigandPreparation 898
899 -class AD4LigandPreparation(LigandPreparation):
900
901 - def __init__(self, molecule, mode='automatic', repairs='checkhydrogens', 902 charges_to_add='gasteiger', cleanup='nphs_lps_waters_nonstdres', 903 allowed_bonds='backbone', #backbone on, amide + guanidinium off 904 root='auto', outputfilename=None, dict=None, debug=False, 905 check_for_fragments=False, bonds_to_inactivate=[], 906 inactivate_all_torsions=False, 907 version=4, typeAtoms=True, limit_torsion_number=False, 908 limit_torsions_type=None, limit_torsions=False, 909 delete_single_nonstd_residues=False):
910 911 912 913 #FIX THIS: what if molecule already has autodock_element set??? 914 LigandPreparation.__init__(self, molecule, mode='interactive', 915 repairs=repairs, charges_to_add=charges_to_add, 916 cleanup=cleanup, allowed_bonds=allowed_bonds, 917 root=root, outputfilename=