Biskit.Dock.Complex

46 """ 47 Manage receptor and ligand structure, store additional information, 48 calculate some properties. 49 """ 50

51 - def __init__(self, rec_model=None,lig_model=None, ligMatrix=None,info={} ):

52 """ 53 @param rec_model: model of original receptor conformation 54 @type rec_model: PDBModel OR PCRModel 55 @param lig_model: model of original ligand conformation 56 @type lig_model: PDBModel OR PCRModel 57 @param ligMatrix: Numeric array 4 by 4, ligand transformation matrix 58 @type ligMatrix: matrix 59 @param info: optional dictionary with additional infos 60 e.g. {'eshape':-123.3, 'rms':12.2 } 61 @type info: dict 62 """ 63 self.rec_model = rec_model # PCRModel object for receptor 64 self.lig_model = lig_model # " for ligand 65 self.lig_transformed = None # " with transformed coordinates 66 self.pw_dist = None # cached pw atom distances rec x lig 67 68 self.info = { 'date':t.dateSortString() } ## default info record 69 self.info.update( info ) 70 71 self.ligandMatrix = ligMatrix 72 if self.ligandMatrix == None: 73 self.ligandMatrix = N.array([ [1, 0, 0, 0], 74 [0, 1, 0, 0], 75 [0, 0, 1, 0], 76 [0, 0, 0, 1],], N.Float32) 77 78 ## compressed by slim 79 self.contacts = None 80 81 ## version as of creation of this object 82 self.initVersion = self.version()

83 84

85 - def version( self ):

86 """ 87 Version of Dock.Complex 88 89 @return: version of class 90 @rtype: str 91 """ 92 return 'Complex $Revision: 2.17 $'

93 94

95 - def __getstate__( self ):

96 """ 97 Called before pickling. 98 """ 99 self.slim() 100 return self.__dict__

101 102

103 - def __getitem__( self, key ):

104 """ 105 @return: C.info[ key ] 106 @rtype: dict 107 """ 108 return self.info[ key ]

109 110

111 - def __setitem__( self, key, v ):

112 self.info[ key ] = v

113 114

115 - def __delitem__( self, key ):

116 del self.info[ key ]

117 118

119 - def __contains__( self, key ):

120 return key in self.info

121 122

123 - def __setstate__(self, state ):

124 """ 125 called for unpickling the object. 126 """ 127 self.__dict__ = state 128 self.ligandMatrix = N.array( self.ligandMatrix,N.Float32 ) 129 ## backwards compability 130 self.__defaults()

131 132

133 - def __defaults(self ):

134 """ 135 backwards compatibility to earlier pickled models 136 """ 137 self.pw_dist = getattr( self, 'pw_dist', None )

138 139

140 - def keys( self ):

141 """ 142 Keys of info dictionary. 143 144 @return: list of keys 145 @rtype: [str] 146 """ 147 return self.info.keys()

148 149

150 - def has_key( self, k ):

151 """ 152 Check if info dictionary has key. 153 154 @param k: key to test 155 @type k: str 156 157 @return: dict has key 158 @rtype: 1|0 159 """ 160 return self.info.has_key( k )

161 162

163 - def values( self, keys=[], default=None ):

164 """ 165 Use:: 166 values( keys=None, default=None ) -> list of info dict values 167 168 @param keys: give only values of these keys 169 @type keys: [str] 170 @param default: only used with keys!=[], default value for missing keys 171 @type default: any 172 173 @return: all values (default) or values of requested keys (ordered) 174 @rtype: [any] 175 """ 176 if not keys: 177 return self.info.values() 178 return [ self.get( k, default ) for k in keys ]

179 180

181 - def get( self, key, default=None):

182 """ 183 Use:: 184 C.get(k[,default]) -> C.info[k] if C.info.has_key(k), else 185 default defaults to None. 186 187 @param key: key to call 188 @type key: str 189 @param default: default value if dictionary doesn't have key 190 (default: None) 191 @type default: any 192 193 @return: dictionary value of key OR else default 194 @rtype: any OR None 195 """ 196 return self.info.get( key, default )

197 198

199 - def rec(self):

200 """ 201 Return PCRModel object. 202 203 @return: receptor model 204 @rtype: PCRModel 205 """ 206 return self.rec_model

207 208

209 - def lig(self, force=0, cache=1 ):

210 """ 211 Return ligand structure transformed. Use cached object, if 212 available. 213 214 @param force: force new transformation (default: 0) 215 @type force: 1|0 216 @param cache: cache transformed model (default: 1) 217 @type cache: 1|0 218 219 @return: ligand model 220 @rtype: PCRModel 221 """ 222 try: 223 lig = self.lig_transformed 224 except: 225 lig = None 226 227 if lig == None or force: 228 ## get transformation matrix and apply it 229 lig = self.lig_model.transform( *self.ligMatrix() ) 230 231 if cache: 232 self.lig_transformed = lig 233 234 return lig

235 236

237 - def model(self):

238 """ 239 Model with both receptor and ligand. 240 241 @return: single PDBModel with first rec then lig 242 @rtype: PCRModel 243 """ 244 return self.rec().concat( self.lig() )

245 246

247 - def ligMatrix(self):

248 """ 249 Return transformation matrix to transform original Ligand 250 PCRModel into docked conformation. 251 252 @return: tuple with rotation and transformation matrix 253 @rtype: array, vector 254 """ 255 a = self.ligandMatrix 256 ## return tuple of rotation matrix and translation vector 257 return (a[0:3,0:3], a[0:3, 3])

258 259

260 - def setLigMatrix( self, m ):

261 """ 262 Set a ligand translation/rotation matrix. 263 264 @param m: translation array 4x4 with rotation 265 @type m: array 266 """ 267 self.ligandMatrix = m 268 self.lig_transformed = None

269 270

271 - def getInfo(self):

272 """ 273 Return contents of the info dictionary. 274 275 @return: info dictionary 276 @rtype: dict 277 """ 278 return self.info

279 280

281 - def writeLigand(self, fname, ter=0 ):

282 """ 283 Write transformed ligand to pdb file. Remove TER record for 284 xplor usage (unless removeTer!=0). 285 286 @param fname: output filename 287 @type fname: str 288 @param ter: option for how to treat the TER statement:: 289 - 0, don't write any TER statements (default) 290 - 1, restore original TER statements 291 - 2, put TER between all detected chains 292 @type ter: 0|1|2 293 """ 294 pdb = self.lig() 295 pdb.writePdb( t.absfile( fname ), ter )

296 297

298 - def writeReceptor(self, fname, ter=0 ):

299 """ 300 Write receptor to pdb without TER statement (unless removeTer!=0). 301 302 @param fname: output file name 303 @type fname: str 304 @param ter: option for how to treat the TER statement:: 305 - 0, don't write any TER statements (default) 306 - 1, restore original TER statements 307 - 2, put TER between all detected chains 308 @type ter: 0|1|2 309 """ 310 self.rec().writePdb( t.absfile( fname), ter )

311 312

313 - def writeComplex(self, fname, ter=0):

314 """ 315 Write single pdb containing both receptor and ligand 316 separated by END but not TER (unless ter!=0). 317 318 @param fname: filename of pdb 319 @type fname: str 320 @param ter: option for how to treat the TER statement:: 321 - 0, don't write any TER statements (default) 322 - 1, restore original TER statements 323 - 2, put TER between all detected chains 324 """ 325 self.model().writePdb( t.absfile( fname ), ter )

326 327

328 - def slim(self):

329 """ 330 Remove coordinates and atoms of ligand and receptor from memory, 331 if they can be restored from file, compress contact matrix. 332 @note: CALLED BEFORE PICKLING 333 """ 334 self.lig_transformed = None 335 self.pw_dist = None 336 337 ## self.ligandMatrix = self.ligandMatrix.tolist() 338 339 if 'matrix' in self.info: 340 del self.info['matrix'] 341 342 ## compress contact matrix array 343 if self.contacts != None and \ 344 len(N.shape( self.contacts['result'] ) )==2: 345 m = self.contacts['result'] 346 self.contacts['shape'] = N.shape( m ) 347 348 self.contacts['result'] = N.nonzero( N.ravel( m ) ).astype(N.Int32)

349 350

351 - def contactsOverlap(self, ref, cutoff=None):

352 """ 353 Fraction of overlapping B{residue-residue} contacts between this and 354 reference complex. 355 356 @param ref: reference complex 357 @type ref: Complex 358 @param cutoff: maximal atom-atom distance, None .. previous setting 359 @type cutoff: float 360 361 @return: fraction of contacts shared between this and ref 362 (normalized to number of all contacts) 363 @rtype: float 364 """ 365 equal = N.logical_and(self.resContacts( cutoff=cutoff ), 366 ref.resContacts( cutoff=cutoff ) ) 367 total = N.logical_or( self.resContacts(cutoff), 368 ref.resContacts(cutoff) ) 369 370 return N.sum(N.sum( equal )) * 1.0 / N.sum(N.sum( total ))

371 372

373 - def contactsShared(self, reference, cutoff=None):

374 """ 375 Number of equal B{residue-residue} contacts in this and 376 reference complex. 377 378 @param reference: reference complex 379 @type reference: Complex 380 @param cutoff: cutoff for atom-atom contact to be counted 381 @type cutoff: float 382 @return: the number or residue-residue contacts that are common to 383 both this and reference:: 384 abs( N.sum( N.sum( contactMatrix_a - contactMatrix_b ))) 385 @rtype: int 386 """ 387 equality = N.logical_and(self.resContacts( cutoff=cutoff ), 388 reference.resContacts( cutoff=cutoff ) ) 389 return abs(N.sum(N.sum( equality )))

390 391

392 - def contactsDiff(self, ref, cutoff=None):

393 """ 394 Number of different B{residue-residue} contacts in this and 395 reference complex. 396 397 @param ref: to compare this one with 398 @type ref: Complex 399 @param cutoff: maximal atom-atom distance, None .. previous setting 400 @type cutoff: float 401 402 @return: number of contacts different in this and refererence complex. 403 @rtype: int 404 """ 405 both = N.logical_or( self.resContacts(cutoff), ref.resContacts(cutoff)) 406 return N.sum(N.sum(both)) - self.contactsShared( ref, cutoff )

407 408

409 - def fractionNativeContacts(self, ref, cutoff=None ):

410 """ 411 Fraction of native B{residue-residue} contacts. 412 413 @param ref: native complex 414 @type ref: Complex 415 @param cutoff: maximal atom-atom distance, None .. previous setting 416 @type cutoff: float 417 418 @return: fraction of native contacts 419 @rtype: float 420 """ 421 cont = self.resContacts( cutoff, refComplex=ref ) 422 ref_cont = ref.resContacts( cutoff ) 423 424 result = N.sum(N.sum( ref_cont * cont ))*1.0 425 return result / N.sum( N.sum( ref_cont ))

426 427

428 - def fractionNativeSurface(self, cont, contRef ):

429 """ 430 fraction of atoms/residues that are involved in B{any} contacts 431 in both complexes. 432 433 @param cont: contact matrix 434 @type cont: matrix 435 @param contRef: reference contact matrix 436 @type contRef: matrix 437 438 @return: (fractRec, fractLig), fraction of atoms/residues that 439 are involved in any contacts in both complexes 440 @rtype: (float, float) 441 442 """ 443 lig, ligRef = N.clip( N.sum(cont),0,1), N.clip( N.sum(contRef), 0,1) 444 rec = N.clip( N.sum(cont, 1),0,1) 445 recRef = N.clip( N.sum(contRef, 1), 0,1) 446 447 fLig = N.sum( N.logical_and( lig, ligRef )) *1./ N.sum( ligRef ) 448 fRec = N.sum( N.logical_and( rec, recRef )) *1./ N.sum( recRef ) 449 450 return (fRec, fLig)

451 452

453 - def rmsLig( self, ref ):

454 """ 455 Rms of ligand from reference ligand after superimposing the two 456 receptors. 457 458 @param ref: reference complex, must have identical atoms 459 @type ref: Complex 460 461 @return: ligand rmsd 462 @rtype: float 463 464 @note: not implemented 465 """ 466 pass

467 468

469 - def rmsInterface( self, ref, cutoff=4.5, fit=1 ):

470 """ 471 Rmsd between this and reference interface. The interface is 472 defined as any residue that has an atom which is within the 473 distance given by |cutoff| from its partner. 474 475 @param ref: reference complex 476 @type ref: Complex 477 @param cutoff: atom distance cutoff for interface residue definition 478 (default: 4.5) 479 @type cutoff: float 480 @param fit: least-squares fit before calculating the rms (default: 1) 481 @type fit: 1|0 482 483 @return: interface rmad 484 @rtype: float 485 """ 486 ## casting 487 this = self 488 if not ref.rec_model.equals( self.rec_model )[1] \ 489 or not ref.lig_model.equals( self.lig_model )[1]: 490 491 m_rec, m_rec_ref, m_lig, m_lig_ref = self.equalAtoms( ref ) 492 this = self.compress( m_rec, m_lig ) 493 ref = ref.compress( m_rec_ref, m_lig_ref ) 494 495 ## determine interface 496 contacts = ref.resContacts( cutoff ) 497 498 if_rec = ref.rec_model.res2atomMask( N.sum( contacts, 1 ) ) 499 if_lig = ref.lig_model.res2atomMask( N.sum( contacts, 0 ) ) 500 501 mask_interface = N.concatenate( (if_rec, if_lig) ) 502 mask_heavy = N.concatenate( (ref.rec().maskHeavy(), 503 ref.lig_model.maskHeavy()) ) 504 mask_interface = mask_interface * mask_heavy 505 506 ## rms 507 ref_model = ref.model() 508 this_model= this.model() 509 510 return ref_model.rms( this_model, mask_interface, fit=fit)

511 512

513 - def contactResPairs(self, cm=None ):

514 """ 515 Get list of residue names paired up in contacts. 516 517 @param cm: pre-calculated contact matrix (default: None) 518 @type cm: matrix 519 520 @return: list of tuples [('N','G'), ('P','C')..] 521 @rtype: [(str.str)..] 522 """ 523 if cm == None: 524 cm = self.resContacts() 525 526 seq_lig = self.lig().sequence() 527 seq_rec = self.rec().sequence() 528 529 result = [] 530 531 for i in range( 0, len(seq_rec) ): 532 recRes = seq_rec[i] 533 534 for j in range( 0, len(seq_lig)): 535 ligRes = seq_lig[j] 536 537 if cm[i,j]: 538 result += [(recRes, ligRes)] 539 540 return result

541 542

543 - def contactResDistribution( self, cm=None ):

544 """ 545 Count occurrence of residues in protein-protein interface. 546 547 @param cm: pre-calculated contact matrix (default: None) 548 @type cm: matrix 549 550 @return: dict {'A':3, 'C':1, .. } (20 standard amino acids) 551 @rtype: dict 552 """ 553 if cm == None: 554 cm = self.resContacts() 555 556 ## get mask for residues involved in contacts 557 maskLig = N.sum( cm ) 558 maskRec = N.sum( N.transpose( cm )) 559 560 ## get sequence of contact residues only 561 seqLig = N.compress( maskLig, self.lig().sequence() ) 562 seqRec = N.compress( maskRec, self.rec().sequence() ) 563 seq = ''.join( seqLig ) + ''.join(seqRec) ## convert back to string 564 565 ## count occurrence of letters 566 result = {} 567 for aa in molUtils.allAA(): 568 result[aa] = seq.count( aa ) 569 570 return result

571 572

573 - def contactTypeDistribution( self, cm = None ):

574 """ 575 Count occurrence of residue types aromatic (a), charged(c), 576 polar(p), else/unpolar(u). 577 578 @param cm: pre-calculated contact matrix (default: None) 579 @type cm: matrix 580 581 @return: dict {'a':4, 'c':5, 'p':12, 'u':4} 582 @rtype: dict 583 """ 584 resDistr = self.contactResDistribution( cm ) 585 586 result = {'a':0, 'c':0, 'p':0, 'u':0} 587 588 for r in resDistr.keys(): 589 result[ molUtils.resType( r )[0] ] += resDistr[ r ] 590 591 ## normalize 592 s = N.sum( result.values() ) 593 for r in result.keys(): 594 result[r] = result[r]*1.0 / s 595 596 return result

597 598

599 - def resPairCounts(self, cm=None ):

600 """ 601 Count how often (all possible) res-res combinations occur 602 in the given contacts. 603 604 @param cm: pre-calculated contact matrix (default: None) 605 @type cm: matrix 606 607 @return: dict {'AA':3,'AC':0, .. 'WW':2, 'WY':0 } 608 @rtype: dict 609 """ 610 resPairs = self.contactResPairs( cm ) 611 allAA = molUtils.allAA() 612 613 allPairs = {} 614 615 ## create dictionary entries for each possible pair of residues 616 for i in range(0, len(allAA)): 617 618 for j in range(i, len( allAA ) ): 619 key = t.sortString( allAA[j]+allAA[i] ) 620 allPairs[key] = 0 621 622 ## count occurrence of each pair 623 for pair in resPairs: 624 625 key = t.sortString( pair[0] + pair[1] ) 626 allPairs[ key ] += 1 627 628 return allPairs

629 630

631 - def loadResContacts( self ):

632 """ 633 Uncompress residue contact matrix if necessary. 634 635 @return: dict with contact matrix and parameters OR None 636 @rtype: dict OR None 637 """ 638 ## Backwards compatibility 639 if self.contacts != None and type( self.contacts ) == str: 640 self.contacts = t.Load( self.contacts ) 641 EHandler.warning("loading old-style pickled contacts.") 642 return self.contacts 643 644 ## New, uncompression from list of indices into raveled array 645 if self.contacts != None and \ 646 len( N.shape( self.contacts['result'])) == 1: 647 648 try: 649 lenRec, lenLig = self.contacts['shape'] 650 except: 651 EHandler.warning("uncompressing contacts without shape") 652 lenRec = self.rec().lenResidues() 653 lenLig = self.lig().lenResidues() 654 655 m = N.zeros( lenRec * lenLig ) 656 N.put( m, self.contacts['result'], 1 ) 657 658 self.contacts['result'] = N.reshape( m, (lenRec, lenLig) ) 659 660 return self.contacts

661 662

663 - def resContacts(self, cutoff=4.5, maskRec=None, maskLig=None, 664 refComplex=None, force=0, cache=1, cache_pw=0 ):

665 """ 666 Matrix of all residue - residue contacts between receptor and 667 ligand. Result is cached. 668 669 @param cutoff: float/int, cutoff in \AA for atom-atom contact to be 670 counted ( default 4.5; if None, last one used or 4.5) 671 @type cutoff: float 672 @param maskRec: atom mask, receptor atoms to consider 673 @type maskRec: [1|0] 674 @param maskLig: atom mask, ligand atoms to consider 675 @type maskLig: [1|0] 676 @param force: re-calculate even if cached matrix is available 677 (default: 0) 678 @type force: 0|1 679 @param cache: cache result (default: 1) 680 @type cache: 0|1 681 @param cache_pw: cache pairwise atom distances (default:0) 682 @type cache_pw: 0|1 683 @param refComplex: if <> None, 'reshape' contact matrix, so that 684 residue positions match residue positions in 685 matrices from castComplex. Useful if calling 686 complex is a crystallized reference structure 687 with slight sequence variations from the docked 688 receptor and ligand. 689 @type refComplex: Complex 690 691 @return: residue contact matrix, 692 2-D array(residues_receptor x residues_ligand) of 0 or 1 693 @rtype: array 694 """ 695 if cutoff == None: 696 if self.contacts != None: 697 cutoff = self.contacts['cutoff'] 698 else: 699 cutoff = 4.5 700 701 if not force and self.loadResContacts() != None \ 702 and self.contacts['cutoff'] == cutoff \ 703 and self.contacts['maskRec'] == maskRec \ 704 and self.contacts['maskLig'] == maskLig: 705 706 result = self.contacts['result'] 707 708 else: 709 result = self.__resContacts( cutoff, maskRec, maskLig, cache_pw ) 710 711 if cache: 712 self.contacts = {} 713 self.contacts['cutoff'] = cutoff 714 self.contacts['maskRec'] = maskRec 715 self.contacts['maskLig'] = maskLig 716 self.contacts['result'] = result 717 718 # delete/insert rows or columns to match sequence of reference complex 719 if refComplex != None: 720 result = self.__alignMatrixDimension(result, 721 self.rec_model.sequence(), 722 refComplex.rec_model.sequence() ) 723 result = self.__alignMatrixDimension(result, 724 self.lig_model.sequence(), 725 refComplex.lig_model.sequence(), 726 1) 727 728 return result

729 730

731 - def __alignMatrixDimension(self, cm, thisSeq, castSeq, axis=0):

732 """ 733 Correct one dimension of contactMatrix by inserting and deleting 734 columns, so that it can be later compared to contact matrices based 735 on slightly different sequences. 736 737 @param cm: contact matrix, 2D matrix of residue contacts 738 recceptor x ligand sequence 739 @type cm: array 740 @param thisSeq: AA sequence of this dimension of the contactMatrix 741 @type thisSeq: string 742 @param castSeq: AA sequence of this dimension in the other contact 743 @type castSeq: string 744 @param axis: which dimension to adapt (0=receptor, 1=ligand) 745 @type axis: 1|0 746 747 @return: contact matrix with residue contacts compatible to refSeq. 748 @rtype: 2D array 749 """ 750 # compare the two sequences 751 seqdiff = SequenceMatcher(None, thisSeq, castSeq) 752 seqDiff = seqdiff.get_opcodes() 753 ## print seqDiff 754 755 # decide which dimension to work on 756 if not axis: 757 cm = N.transpose( cm ) 758 759 seqCount = 0 # keep track of sequence length changes 760 i=0 761 762 for list in seqDiff: 763 764 # remove the column corresponding to the deletion in the 765 # docked sequence 766 if str( seqDiff[i][0] ) == 'delete': 767 768 # separate matrix into before and after deletion 769 matrixSeg1 = cm[ :, : seqDiff[i][1] + seqCount ] 770 matrixSeg2 = cm[ :, seqDiff[i][2] + seqCount : ] 771 # concatenate part 772 cm = N.concatenate( ( matrixSeg1, matrixSeg2 ), 1) 773 seqCount = seqCount + seqDiff[i][1] - seqDiff[i][2] 774 775 # inserts zeros in the column where there is a insertion in the 776 # docked sequence 777 if str( seqDiff[i][0] ) == 'insert': 778 779 # create a matrix to be inserted 780 insertZeros= seqDiff[i][4] - seqDiff[i][3] 781 insertColumns = N.array( [ [0] * insertZeros ] * N.size(cm,0) ) 782 # separate matrix into before and after insertion 783 matrixSeg1 = cm[ :, : seqDiff[i][1] + seqCount ] 784 matrixSeg2 = cm[ :, seqDiff[i][2] + seqCount : ] 785 # concatenate parts with the zero matrix 786 cm = N.concatenate( (matrixSeg1,insertColumns,matrixSeg2), 1) 787 seqCount = seqCount + seqDiff[i][4] - seqDiff[i][3] 788 789 i=i+1 790 791 if not axis: 792 return N.transpose( cm ) 793 return cm

794 795

796 - def __unmaskedMatrix( self, contacts, rec_mask, lig_mask ):

797 """ 798 Map contacts between selected rec and lig atoms back to all atoms 799 matrix. 800 801 @param contacts: contact matrix, array sum_rec_mask x sum_lig_mask 802 @type contacts: array 803 @param rec_mask: atom mask 804 @type rec_mask: [1|0] 805 @param lig_mask: atom mask 806 @type lig_mask: [1|0] 807 808 @return: atom contact matrix, array N_atoms_rec x N_atoms_lig 809 @rtype: array 810 """ 811 l_rec = len( self.rec_model ) 812 l_lig = len( self.lig_model ) 813 814 ## map contacts back to all atoms matrix 815 r = N.zeros( l_rec * l_lig ) 816 rMask = N.ravel( N.outerproduct( rec_mask, lig_mask ) ) 817 818 ## (Optimization: nonzero is time consuming step) 819 N.put( r, N.nonzero( rMask ), N.ravel( contacts ) ) 820 821 return N.resize( r, (l_rec, l_lig))

822 823

824 - def __atomContacts(self, cutoff, rec_mask, lig_mask, cache):

825 """ 826 Intermolecular distances below cutoff after applying the two masks. 827 828 @param cutoff: cutoff for B{atom-atom} contact in \AA 829 @type cutoff: float 830 @param rec_mask: atom mask 831 @type rec_mask: [1|0] 832 @param lig_mask: atom mask 833 @type lig_mask: [1|0] 834 @param cache: cache pairwise atom distance matrix 835 @type cache: 1|0 836 837 @return: atom contact matrix, array sum_rec_mask x sum_lig_mask 838 @rtype: array 839 """ 840 ## get atom coordinats as array 3 x all_atoms 841 rec_xyz = self.rec().getXyz() 842 lig_xyz = self.lig().getXyz() 843 844 ## get pair-wise distances -> atoms_rec x atoms_lig 845 dist = getattr( self, 'pw_dist', None ) 846 if dist is None or \ 847 N.shape( dist ) != ( N.sum(rec_mask), N.sum(lig_mask) ): 848 dist = self.__pairwiseDistances(N.compress( rec_mask, rec_xyz, 0), 849 N.compress( lig_mask, lig_xyz, 0) ) 850 if cache: 851 self.pw_dist = dist 852 853 ## reduce to 1 (distance < cutoff) or 0 -> n_atoms_rec x n_atoms_lig 854 return N.less( dist, cutoff )

855 856

857 - def atomContacts( self, cutoff=4.5, rec_mask=None, lig_mask=None, cache=0, 858 map_back=1 ):

859 """ 860 Find all inter-molecular B{atom-atom} contacts between rec and lig 861 862 @param cutoff: cutoff for atom - atom contact in \AA 863 @type cutoff: float 864 @param rec_mask: atom mask (default: all heavy) 865 @type rec_mask: [1|0] 866 @param lig_mask: atom mask (default: all heavy) 867 @type lig_mask: [1|0] 868 @param cache: cache pairwise atom distance matrix (default: 0) 869 @type cache: 1|0 870 @param map_back: map masked matrix back to matrix for all atoms 871 @type map_back: 1|0 872 873 @return: atom contact matrix, Numpy array N(atoms_lig) x N(atoms_rec) 874 @rtype: array 875 """ 876 if lig_mask == None: 877 lig_mask = self.lig().maskHeavy() 878 879 if rec_mask == None: 880 rec_mask = self.rec().maskHeavy() 881 882 contacts = self.__atomContacts( cutoff, rec_mask, lig_mask, cache ) 883 884 if not map_back: 885 ## contact matrix after masking rec and lig 886 return contacts 887 888 return self.__unmaskedMatrix( contacts, rec_mask, lig_mask )

889 890

891 - def __resContacts(self, cutoff, maskRec=None, maskLig=None, cache=0 ):

892 """ 893 Find all inter-molecule B{residue-residue} contacts between receptor 894 and ligand. A contact between A and B is set if any heavy atom of 895 A is within |cutoff| A of B. 896 897 @param cutoff: distance cutoff in \AA 898 @type cutoff: float 899 @param maskRec: atom mask (default: all heavy) 900 @type maskRec: [1|0] 901 @param maskLig: atom mask (default: all heavy) 902 @type maskLig: [1|0] 903 @param cache: cache pairwise atom distance matrix to pw_dist 904 (default:0) 905 @type cache: 1|0 906 907 @return: residue contact matrix, 2-D Numpy 908 N.array(residues_receptor x residues_ligand) where 909 1-contact, 0-no contact 910 @rtype: array 911 """ 912 ## get contact matrix atoms_rec x atoms_lig 913 c = self.atomContacts( cutoff, maskRec, maskLig, cache ) 914 915 ## convert atoms x atoms to residues x residues matrix 916 return self.__atom2residueMatrix( c )

917 918

919 - def __atom2residueMatrix( self, m ):

920 """ 921 Reduce binary matrix of n x k atoms to binary matrix of i x j residues. 922 923 @param m: atom contact matrix, 924 array n x k with 1(contact) or 0(no contact) 925 @type m: array 926 927 @return: residue contact matrix, 928 2-D numpy array(residues_receptor x residues_ligand) 929 @rtype: array 930 """ 931 recInd = N.concatenate((self.rec().resIndex(), 932 [ self.rec().lenAtoms()] )) 933 ligInd = N.concatenate((self.lig_model.resIndex(), 934 [ self.lig_model.lenAtoms() ] )) 935 936 residueMatrix = N.zeros(( len(recInd)-1, len(ligInd)-1 ), N.Int) 937 938 for r in range( len(recInd)-1 ): 939 940 for l in range( len(ligInd)-1 ): 941 942 res2res = m[ int(recInd[r]):int(recInd[r+1]), 943 int(ligInd[l]):int(ligInd[l+1]) ] 944 945 if N.any( res2res ): 946 residueMatrix[r, l] = 1 947 948 return residueMatrix

949 950

951 - def equalAtoms( self, ref ):

952 """ 953 Compare two complexes' atom content of receptor and ligand. 954 Apply to SORTED models without HETATOMS. Coordinates are not checked. 955 956 @param ref: reference complex 957 @type ref: Complex 958 959 @return: (mask_rec, mask_lig, mask_rec_ref, mask_lig_ref), 960 4 atom masks for all equal atoms 961 @rtype: [1|0], [1|0], [1|0], [1|0] 962 963 @note: in some rare cases m1.equalAtoms( m2 ) gives a different 964 result than m2.equalAtoms( m1 ). This is due to the used 965 SequenceMatcher class. 966 """ 967 m_rec, m_rec_ref = self.rec_model.equalAtoms( ref.rec_model ) 968 m_lig, m_lig_ref = self.lig_model.equalAtoms( ref.lig_model ) 969 return m_rec, m_lig, m_rec_ref, m_lig_ref

970 971

972 - def compareAtoms( self, ref ):

973 """ 974 Compare atom content and oder of rec and lig of 2 complexes. 975 Get list of all atom positions of rec and lig in both complexes that 976 have a matching atom (by residue and atom name) in the other complex. 977 Indices for this complex are returned in the right oder to match the 978 atom order of ref. I.e., in contrast to equalAtoms, castAtoms can 979 equalize two complexes where atoms are ordered differently within the 980 residues. 981 982 @param ref: reference complex 983 @type ref: Complex 984 985 @return: (ind_rec, ind_lig, ind_rec_ref, ind_lig_ref), 986 4 lists of indices 987 @rtype: [int], [int], [int], [int] 988 """ 989 i_rec, i_rec_ref = self.rec_model.compareAtoms( ref.rec_model ) 990 i_lig, i_lig_ref = self.lig_model.compareAtoms( ref.lig_model ) 991 992 return i_rec, i_lig, i_rec_ref, i_lig_ref

993 994

995 - def take( self, rec_pos, lig_pos ):

996 """ 997 Get copy of this complex with given atoms of rec and lig. 998 999 @param rec_pos: receptor indices to take 1000 @type rec_pos: [int] 1001 @param lig_pos: ligand indices to take 1002 @type lig_pos: [int] 1003 1004 @return: new complex 1005 @rtype: Complex 1006 """ 1007 r = self.__class__() 1008 r.lig_model = self.lig_model.take( lig_pos ) 1009 r.rec_model = self.rec_model.take( rec_pos ) 1010 r.info = deepcopy( self.info ) 1011 1012 if self.pw_dist: 1013 r.pw_dist = N.take( self.pw_dist, rec_pos, 1 ) 1014 r.pw_dist = N.take( r.pw_dist, lig_pos ) 1015 1016 r.ligandMatrix = copy( self.ligandMatrix ) 1017 1018 ## todo: take cached contacts as well 1019 1020 return r

1021 1022

1023 - def compress( self, rec_mask, lig_mask ):

1024 """ 1025 Compress complex using a rec and lig mask. 1026 1027 @param rec_mask: atom mask 1028 @type rec_mask: [1|0] 1029 @param lig_mask: atom mask 1030 @type lig_mask: [1|0] 1031 1032 @return: compressed complex 1033 @rtype: Complex 1034 """ 1035 return self.take( N.nonzero( rec_mask ), N.nonzero( lig_mask ) )

1036 1037

1038 - def contPairScore(self, cutoff=6.0):

1039 """ 1040 Score interaction surface residue pairs. 1041 Info on Scoring matrix see L{Biskit.molUtils} 1042 1043 @param cutoff: CB-CB distance cutoff for defining a contact 1044 (default: 6.0) 1045 @type cutoff: float 1046 1047 @return: score 1048 @rtype: float 1049 """ 1050 score = 0 1051 cm = self.resContacts(cutoff,self.rec().maskCB(), self.lig().maskCB(), 1052 cache=0 ) 1053 1054 pairFreq = self.resPairCounts(cm) 1055 1056 for pair in pairFreq: 1057 score += pairFreq[pair]*molUtils.pairScore[pair] 1058 1059 return score

1060 1061

1062 - def interfaceArea( self, profiles=0, log=StdLog(), verbose=1 ):

1063 """ 1064 Calculate the difference between the surface area of the 1065 complex vs. its free components in square angstrom. 1066 1067 @param profiles: option to return the lig, rec and com profiles 1068 rather than the value (both absolute and relative 1069 values are returned) 1070 @type profiles: 1|0 (default: 0) 1071 @param log: log file [STDOUT] 1072 @type log: Biskit.LogFile 1073 @param verbose: give progress report [1] 1074 @type verbose: bool | int 1075 1076 @return: AS area, MS area OR 1077 a dictionary of lig, rec and com profiles 1078 @rtype: (float, float) or dict 1079 """ 1080 rcom = self.rec() 1081 lcom = self.lig() 1082 ccom = self.model() 1083 result = {} 1084 1085 def getSurface( model, key ): 1086 if verbose: 1087 log.write("Calculating SurfaceRacer data for %s..."%key) 1088 d = PDBDope( model ) 1089 d.addSurfaceRacer( probe=1.4 ) 1090 if verbose: 1091 log.writeln('Done.') 1092 result['%s_AS'%key] = model.profile('AS') 1093 result['%s_MS'%key] = model.profile('MS') 1094 result['%s_relAS'%key] = model.profile('relAS') 1095 result['%s_relMS'%key] = model.profile('relMS')

1096 1097 getSurface( rcom, 'rec' ) 1098 getSurface( lcom, 'lig' ) 1099 getSurface( ccom, 'com' ) 1100 1101 if not profiles: 1102 return sum(result['rec_AS']) + sum(result['lig_AS']) - \ 1103 sum(result['com_AS']),\ 1104 sum(result['rec_MS']) + sum(result['lig_MS']) - \ 1105 sum(result['com_MS']) 1106 else: 1107 return result

Source Code for Module Biskit.Dock.Complex