# YASARA MACRO # TOPIC: 5. Structure prediction # TITLE: Showing the docking results interactively # REQUIRES: Structure # AUTHOR: Elmar Krieger # LICENSE: GPL # DESCRIPTION: This macro provides an interactive player to cycle through the docking poses and clusters # You can either set the target structure by clicking on Options > Macro > Set target, # by providing it as command line argument (see docs at Essentials > The command line), # or by uncommenting the line below and specifying it directly. #MacroTarget '/home/myname/projects/docking/1sdf' # Set playback waittime (1 is maximum playback speed) waittime=12 # Flags of interactions to show by default (combined with |): # 1=H-Bonds, 2=Hydrophobic, 4=Pi-Pi, 8=Cation-Pi, 16=Ionic intflags=1|2|16 # Normally no change required below this point # ============================================ # Sanity checks if MacroTarget=='' RaiseError "This macro requires a target. Either edit the macro file or click Options > Macro > Set target to choose a target structure" Console off intnamelist='Hydrophobic','PiPi','CationPi','Ionic' # PART 1 - DOCKING POSE ANALYSIS DockingPoses: # Load the result scene Clear scefilename='(MacroTarget).sce' scene = FileSize (scefilename) if not scene # No scene present, presumable dock_runensemble was used with scesaved=0, proceed to cluster playback. receptors=999 Go DockingClusters LoadSce (scefilename) # Determine if this was a single receptor&ligand, a receptor ensemble, or a virtual screening receptorlist() = ListObj Receptor??? receptors = count receptorlist ligandlist() = ListObj Ligand???? ligands = count ligandlist if !receptors RaiseError 'No receptor object with name Receptor... was found, please make sure to run the dock_run macro first' if ligands receptorruns = CountMol Obj (ligandlist1) # Collect all ligand energies for i=1 to ligands for j=001 to receptorruns if j<1000 seg='C(j)' else seg=j bindnrg = BFactorAtom Obj Ligand(i) Segment (seg) # Create a list of info strings for each run: -binding energy, ligand number, and segment ID runinfolist((i-1)*receptorruns+j)='(-bindnrg) (i) (seg)' if receptors>1 # If there is more than one receptor (dock_runensemble), then there is just one ligand ligands=1 runs=receptors*receptorruns*ligands else # Covalent docking doesn't have a ligand object if receptors>1 RaiseError 'Only one receptor expected for covalent docking' receptorruns = CountMol Segment C??? Obj RecepFlex for i=1 to receptorruns seg='C(000+i)' bindnrg = BFactorAtom Obj RecepFlex Segment (seg) # Create a list of info strings for each run: -binding energy, ligand number, and segment ID runinfolist(i)='(-bindnrg) (i) (seg)' runs=receptorruns # Sort the info strings using the reversed binding energy (best binders first) runinfolist()=sort runinfolist if runs!=count runinfolist RaiseError 'Number of runs differ [(runs) and (count runinfolist)], please contact support' # Create the user interface ShowPlayer 'Pose' run=0 step=1 delay=waittime mainloop='PoseMainLoop' # Loop through the docking poses PoseMainLoop: while 1 if run<0 or run>=runs Go RewindToStart # Determine docking cycle (receptor ensemble member or virtual screening ligand) # and segment names for current run _,cycle,seg=split runinfolist(run+1) if receptors>1 # Receptor ensemble recobj='Receptor(cycle)' flxobj='RecFlex(cycle)' ligobj='Ligand(cycle)' elif ligands>1 # Virtual screening recobj='Receptor' flxobj='RecFlex(cycle)' ligobj='Ligand(cycle)' elif ligands # Normal docking recobj='Receptor' flxobj='RecFlex' ligobj='Ligand' else # Covalent docking recobj='Receptor' flxobj='RecepFlex' ligobj=flxobj flexobjects = CountObj (flxobj) if !flexobjects # No flexible receptor part, try RecepFlex instead of RecFlex (needed for screening) _,_,_,chrlist() = crack flxobj flxobj='Recep'+fuse chrlist # Show the corresponding receptor and ligand SwitchObj !SimCell,off SwitchObj (recobj) (flxobj) (ligobj),On HideObj (flxobj) (ligobj) HideSecStrObj (flxobj) (ligobj) ShowAtom Segment (seg) Obj (flxobj) (ligobj) ShowSecStrRes Segment (seg) Obj (flxobj) (ligobj) # Always hide the flexible side-chains in the receptor fixedatoms = CountAtom fixed Obj (recobj) if fixedatoms # Get a list of partly free residues in recobj reslist() = ListRes !fixed Obj (recobj) # Show their sidechains ShowAtom SCwithBB Res (join reslist) Obj (recobj) # And hide whatever is not fixed (to show RecepFlex atoms instead) HideAtom !fixed Obj (recobj) # For covalent docking, there may be free backbone atoms, don't show secondary structure there HideSecStrRes !fixed Backbone Obj (recobj) # Change the marked atoms in ligand and flexible sidechains to the currently visible ones atomlist() = MarkAtom for i=1 to 4 atom = ListAtom (atomlist(i)),Format='ATOMNAME Res RESNAME RESNUM Mol MOLNAME Obj OBJNUM visible' atomlist(i) = ListAtom (atom) MarkAtom (atomlist),Zoom=1 # Update the player UpdatePlayer run,runs,'Segment (seg) Obj (ligobj)' # Proceed to next pose run=run+step if step # Wait a bit Wait (delay) else # Pause mode, wait until next click Wait forever # PART 2 - DOCKING CLUSTER ANALYSIS DockingClusters: Clear # Determine cluster format for format in 'yob','pdb' exists = FileSize (MacroTarget)_001.(format) if exists break # Count the clusters for runs=000 to 9999 exists = FileSize (MacroTarget)_(runs+1).(format) if !exists break if runs<1 RaiseError "No docking clusters have been found. This macro must be applied to a docking project created with one of the 'dock_run*' macros" # Create the user interface ShowPlayer 'Cluster' run=0 step=1 delay=waittime posorilist=0,0,50,0,0,0 mainloop='ClusterMainLoop' # Loop through the docking clusters ClusterMainLoop: while 1 DelObj (join intnamelist) if run<0 or run>=runs Go RewindToStart obj = Load(format) (MacroTarget)_(001+run) # Remove line below to display ligand names for screening runs made with YASARA versions < 15.10 UnlabelAll CenterAtom CA Obj (obj) if Objects==1 Pos 0,0,50 Style Ribbon,Stick BallRes Segment C??? else TransferObj 2,1,Local=Keep CopyStyleObj 1,2 SwapObj 1,2 DelObj 2 FreeAll # Show the various interactions if intflags&1 ShowHBoAtom Segment C???,Extend=Yes else HideHBoAll for i=1 to count intnamelist if intflags&(1<=100000000 dissconst='Gigantic' elif dissconst>=100000 dissconst=0+dissconst else dissconst=0.000+dissconst # Display run, binding energy and dissociation constant FillRect X=170,Y=44,Width=150,Height=70 PosText X=170,Y=46,Justify=Left Print '(1+run)/(0+runs)\n(0.00+bindnrg) kcal/mol\n(dissconst) pM' # Display first 20 characters of ligand name (compound) FillRect X=100,Y=115,Width=220,Height=18 compound = CompoundMol (ligandsel) if compound=='' compound='unnamed' compound()=crack compound PosText X=100,Y=115,Justify=Left len=count compound # Maximum reached with uppercase ALCOHOL DEHYDROGENASE if len>15 len=15 for i=1 to len Print '(compound(i))' if runs>1 # Show the '>' at the color gradient FillRect X=(gpos1-20),Y=(gpos2),Width=15,Height=260 PosText X=(gpos1-20),Y=(gpos2+240*run/(runs-1)) Print '>' # HANDLE THE PLAYER BUTTONS # ========================= # Press the Play button Play: step=1 delay=waittime Go (mainloop) # Press the Pause button Pause: step=0 delay=waittime Go (mainloop) # Press the step forward button StepForward: run=run+1 Go (mainloop) # Press the step backward button StepBackward: run=run-1 Go (mainloop) # Press the fast forward button FastForward: step=1 delay=1 Go (mainloop) # Press the rewind to start button RewindToStart: run=0 if step<0 Go Play Go (mainloop) # Click the color gradient on the right ClickGradient: run=(runs-1)*GoPar/10 Go (mainloop) # Save a pose in PDB format PDB: format='pdb' Go SavePose # Save a pose in YOb format YOb: format='yob' SavePose: # Ask for PDB/YOb filename to save _,filename = ShowWin Type=FileSelection,Title='Choose filename to save current complex in (format) format', MultipleSelections=No,Filename='*.(format)' # Duplicate the receptor ligand objects objlist() = DuplicateObj (recobj) (ligobj) # Keep only the current ligand pose DelAtom Segment !(seg) Obj (objlist2) # Transfer atom coordinates of flexible residues to receptor reslist() = ListRes Obj (flxobj),Format='RESNAME RESNUM Mol MOLNAME' for res in reslist # Get atom names in PDB3 format, so that equivalent hydrogens are numbered atomlist() = ListAtom Res (res) Segment (seg) Obj (flxobj) namelist() = ListAtom (join atomlist),Format='ATOMNAMEPDB3' poslist() = PosAtom (join atomlist) for i=1 to count atomlist PosAtom PDB3 (namelist(i)) Res (res) Obj (objlist1),(poslist(i*3-2)),(poslist(i*3-1)),(poslist(i*3)) # Recreate aliphatic hydrogens of flexible residues in receptor, docking doesn't yield coordinates for them DelAtom Element H Res (res) Obj (objlist1) with bond to Element C AddHydRes (res) Obj (objlist1),Update=No # Join to single object and add aliphatic hydrogens again JoinObj (objlist2),(objlist1) FreeObj (objlist1) # Save and delete Save(format) (objlist1),(filename) DelObj (objlist1) Go (mainloop) # Toggle the interactions during cluster visualization # Hydrogen bonds HBonds: intflags=intflags^1 ShowPlayer 'Cluster' Go (mainloop) # Hydrophobic interactions Hydrophobic: intflags=intflags^2 ShowPlayer 'Cluster' Go (mainloop) # Pi/Pi interactions PiPi: intflags=intflags^4 ShowPlayer 'Cluster' Go (mainloop) # Cation/Pi interactions CationPi: intflags=intflags^8 ShowPlayer 'Cluster' Go (mainloop) # Ionic interactions Ionic: intflags=intflags^16 ShowPlayer 'Cluster' Go (mainloop)