# YASARA MACRO # TOPIC: 5. Structure prediction # TITLE: Showing the docking results interactively # REQUIRES: Structure # AUTHOR: Elmar Krieger and Kornel Ozvoldik # 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 RequireVersion 26.4.1 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???? Complex???? ligands = count ligandlist recconfs = CountMol Segment R??? complexes = CountObj Complex???? 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) and not Segment R??? # VINA supports docking groups of ligands groupligands = CountMol Segment C001 Obj (ligandlist1) receptorruns = receptorruns/groupligands # Collect all ligand energies for i=1 to ligands for j=001 to receptorruns if j<1000 ligseg='C(j)' else ligseg=j bindnrg = BFactorAtom Obj Ligand(i) Segment (ligseg) # Create a list of info strings for each run: -binding energy, ligand number, and segment ID runinfolist((i-1)*receptorruns+j)='(-bindnrg) (i) (ligseg) R(j)' 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 if !receptorruns RaiseError 'No ligands were found in your scene "(scefilename)", this file does not seem to have been created by one of the docking macros' for i=1 to receptorruns ligseg='C(000+i)' bindnrg = BFactorAtom Obj RecepFlex Segment (ligseg) # Create a list of info strings for each run: -binding energy, ligand number, and segment ID runinfolist(i)='(-bindnrg) (i) (ligseg) R(j)' runs=receptorruns # Sort the info strings using the negated 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,ligseg,recseg=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) Complex(cycle)' elif ligands # Normal docking recobj='Receptor' flxobj='RecFlex' ligobj='Ligand Complex' 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 using the current style SwitchObj !SimCell,off SwitchObj (flxobj) (ligobj),On if !recconfs SwitchObj (recobj),On # Set the current style for all atoms, hetgroup style is kept anyway bbstyle,scstyle=Style Style (bbstyle),(scstyle) # Hide what we don't need HideAtom not Segment (ligseg) (recseg) and Obj (flxobj) (ligobj) HideSecStrRes not Segment (ligseg) (recseg) and 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 (ligseg) 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 !(ligseg) 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 (ligseg) 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 Pi: intflags=intflags^8 ShowPlayer 'Cluster' Go (mainloop) # Ionic interactions Ionic: intflags=intflags^16 ShowPlayer 'Cluster' Go (mainloop)