# YASARA MACRO # TOPIC: 5. Structure prediction # TITLE: Docking a ligand to a receptor ensemble with flexible side-chains # REQUIRES: Structure # AUTHOR: Elmar Krieger # LICENSE: GPL # DESCRIPTION: This macro predicts the structure of a ligand-receptor complex. Receptor flexibility is considered by creating a receptor ensemble with alternative high-scoring solutions of the side-chain rotamer network. An analysis log file is written at the end. The macro can also continue a docking run that got interrupted, especially if the scene has not been moved or rotated manually during the docking. # Parameter section - adjust as needed, but NOTE that some changes only take effect # if you start an entirely new docking job, not if you continue an existing one. # ================================================================================= # You can either set the target structure by clicking on Options > Macro > Set target, # or by uncommenting the line below and specifying it directly. #MacroTarget='/home/myname/projects/docking/1sdf' # Size of the receptor ensemble with alternative side-chain rotamers receptors = 5 # Number of docking runs per receptor ensemble member # (maximally 999, each run can take up to an hour). # The total number of docking runs is receptors*receptorruns. receptorruns = 20 # Docking results usually cluster around certain hot spot conformations, # and the lowest energy complex in each cluster is saved. Two complexes belong to # different clusters if the ligand RMSD is larger than this minimum [A]: rmsdmin = 5.0 # Set to 1 to keep the ligand completely rigid rigid = 0 # Force field used for charge assignment. The actual scoring function it provided by AutoDock. chargefof='AMBER03' # The docking temperature. This currently affects the receptor flexibility only temp='298K' # 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" if receptors>999 or receptorruns>999 or receptors*receptorruns>9999 RaiseError 'Too many docking runs selected, (receptors)/(receptorruns) would take forever' structlist='receptor','ligand' Clear # Do we have a user-supplied ensemble? ensfilename='(MacroTarget)_receptor_ensemble.sce' scene = FileSize (ensfilename) if !scene # Do we already have an automatically generated receptor ensemble scene? ensfilename='(MacroTarget)_receptor_ensemble(receptors).sce' scene = FileSize (ensfilename) if !scene # Not yet. Do we have a receptor scene? scefilename='(MacroTarget)_receptor.sce' scene = FileSize (scefilename) if !scene # No, load PDB or YOB file of receptor and ligand for struct in structlist (struct)=0 for type in 'yob','pdb' filename='(MacroTarget)_(struct).(type)' exists = FileSize (filename) if exists (struct) = Load(type) (filename) break if not (struct) RaiseError 'The (struct) was not found at (filename)' # Orient the receptor, create a docking cell that includes the entire receptor # and also has enough empty space around to accomodate the ligand (which has to # be removed temporarily, since 'Cell Auto' encloses the entire soup). ligandradius = RadiusObj (ligand) NiceOriObj (receptor) # Ligand was only needed to determine the required cell size, will be loaded later DelObj (ligand) Cell Auto,Extension=(ligandradius*2) else # Load receptor scene LoadSce (scefilename) if Objects>2 # We can't continue, because the ensemble can only be created for a single object RaiseError 'Your scene contains (Objects) objects, while only the receptor and the docking cell are expected.' # Warn if the cell is too large x,y,z = Cell if x>48 or y>48 or z>48 ShowMessage "A cell axis is longer than 48 A, which may reduce docking accuracy. Consider focusing on the active site." Wait ContinueButton HideMessage # To create a receptor ensemble, the entire receptor must be inside the cell, # so we remember the current cell and add it again later x,y,z = Cell posorilist() = PosOriObj SimCell DelObj SimCell # Start an energy minimization but stop immediately. # This cleans the structure in case it was supplied without hydrogens. ForceField (chargefof) ShowMessage 'Checking if receptor is ready for docking...' Experiment Minimization Experiment On Experiment Off # Create a receptor ensemble with flexibility corresponding to the specified temperature ShowMessage 'Creating receptor ensemble with (receptors) members, each with alternative high-scoring side-chain conformations...' Style Ribbon,Stick ForceField YASARA2 Temp (temp) # Don't change side-chain rotamers of residues with pseudo-bonds to metals FixRes all with arrow to metal OptimizeAll Method=SCWALL,Structures=(receptors-1) FreeAll DelObj SimCell RenumberObj all,1 # Add original cell again Cell (x),(y),(z) PosOriObj SimCell,(posorilist) # Optimize the hydrogen bonding network of the ensemble members RemoveObj !SimCell for i=1 to receptors AddObj (i) ShowMessage 'Optimizing hydrogen bonding network of receptor ensemble member (i)/(receptors)...' OptHydAll Method=YASARA NameObj (i),Receptor(i) RemoveObj (i) AddObj !SimCell SaveSce (ensfilename) HideMessage # Load receptor ensemble scene LoadSce (ensfilename) # Docking is done without periodic boundaries Boundary Wall Longrange None ForceField (chargefof) # Count receptors, in case the user provided the ensemble receptors = CountObj !SimCell runs=receptors*receptorruns # Loop over the receptors in the ensemble RemoveObj !SimCell for i=1 to receptors AddObj (i) # Realign the simulation cell with the axes in case the user manually rotated it. alpha,beta,gamma = OriObj SimCell RotateAll Y=(-beta) RotateAll Z=(-gamma) RotateAll X=(-alpha) # Show progress (ShowMessage is used by the docking experiment) LabelAll 'Receptor ensemble member (i)/(receptors)',Height=0.3,Color=Yellow,Y=3.0 # Load the ligand for this receptor ensemble member ligand=0 for type in 'yob','pdb' filename='(MacroTarget)_ligand.(type)' exists = FileSize (filename) if exists ligand = Load(type) (filename) break if not ligand RaiseError 'The ligand was not found at (filename)' # Number the ligands sequentially NameObj (ligand),Ligand(i) # Fix ligand's internal degrees of freedom if requested if rigid FixObj Ligand(i) # Perform the docking Experiment Docking Method AutoDockLGA ReceptorObj (i) LigandObj Ligand(i) Runs (receptorruns) ClusterRMSD (rmsdmin) # Result file number must be separated with '_', in case MacroTarget also ends with number ResultFile (MacroTarget)_(i)_001 # Uncomment below to set the number of energy evaluations (AutoDock ga_num_evals): # DockPar ga_num_evals 25000000 # Uncomment below to add 1000 to the random number seed when more than 999 runs are needed # DockPar seed 1000 # Uncomment the two lines below to provide your own atom parameters (VdW radii etc.) # GridPar parameter_file /Path/To/Custom/AD4_parameters.dat # DockPar parameter_file /Path/To/Custom/AD4_parameters.dat # Uncomment below to keep temporary files in the current working directory: # TmpFileID 1adb_tmp Experiment On Wait ExpEnd UnlabelAll RemoveObj (i) Ligand(i) # Save a scene with all receptor and all ligand conformations AddObj all SaveSce (MacroTarget) # Collect the initial results Console off for i=1 to receptors for j=001 to receptorruns # To pass the docking results back to this macro, the docking experiment # stores the binding energy as the B-factor... run=(i-1)*receptorruns+j resultlist(run),bindnrglist(run) = DockingResult i,j,'Obj (i)','Obj Ligand(i) Segment C(j)' if !(run%10) ShowMessage 'Analyzing docking results, (100*run/runs)% completed...' Wait 1 # Sort the results by binding energy ShowMessage 'Sorting docking results...' Wait 1 SortResults 'resultlist','bindnrglist','None' # Save a log file with an initial analysis RecordLog (MacroTarget) print 'Ensemble docking result analysis' print '================================' print print 'The ligand was docked (receptorruns) times ["Run"] against each of the (receptors) receptors ["Rec"] in the ensemble, yielding the following (runs) results, sorted by binding energy:' print '[More positive energies indicate stronger binding, and negative energies mean no binding]' print print 'Rec | Run |Bind.energy[kcal/mol]|Binding constant[pM]| Contacting receptor residues' print '----+-----+---------------------+--------------------+-----------------------------' for i=1 to runs print (resultlist(i)) print # Delete previous resultlist resultlist()=0 # Now merge the clusters from the receptor ensemble members clusters=0 for i=1 to receptors for j=001 to receptorruns filename='(MacroTarget)_(i)_(j).yob' exists = FileSize (filename) if exists ShowMessage 'Merging clusters from receptor ensemble member (i) run (j), (clusters) found so far...' Wait 1 cluobj = LoadYOb (filename) # Combined clusters will be written at the end DelFile (filename) result,bindnrg = DockingResult i,j,'Obj (cluobj) Segment !C???','Obj (cluobj) Segment C???' added=1 if clusters # Get ligand RMSDs from all other clusters rmsdlist() = RMSDAtom Element !H Segment C??? Obj (cluobj),Element !H Segment C??? Obj (join cluobjlist) for k=1 to clusters if rmsdlist(k)bindnrglist(k) # Current one has better binding energy, delete old one resultlist(k)=result bindnrglist(k)=bindnrg SwapObj (cluobjlist(k)),(cluobj) DelObj (cluobj) break if added # Add a new cluster clusters=clusters+1 resultlist(clusters)=result bindnrglist(clusters)=bindnrg cluobjlist(clusters)=cluobj ShowMessage 'Sorting cluster results...' Wait 1 SortResults 'resultlist','bindnrglist','cluobjlist' # Save the unique clusters in order, mainly needed for the dock_play.mcr CenterObj (join cluobjlist) for i=1 to clusters SaveYOb (cluobjlist(i)),(MacroTarget)_(000+i).yob ShowMessage 'Saving sorted cluster (i)/(clusters)...' Wait 1 DelObj (join cluobjlist) # Print cluster log print 'After clustering the (runs) runs, the following (clusters) distinct complex conformations were found:' print '[They all differ by at least (rmsdmin) A heavy ligand atom RMSD, and have been saved as YOb files with terminal number "Num"]' print print 'Num | Rec | Run |Bind.energy[kcal/mol]|Binding constant[pM]| Contacting receptor residues' print '----+-----+-----+---------------------+--------------------+-----------------------------' for i=1 to clusters print '(000+i) | (resultlist(i))' Style Ribbon StopLog Console On HideMessage # EXTRACT DOCKING RESULT # ====================== # Returns a result string and binding energy, extracted from atomic BFactor and Property. # 'rec' is the receptor number in the ensemble, 'run' is the docking run for this receptor, # 'recsel' selects the receptor, 'ligsel' selects the ligand. def DockingResult rec,run,recsel,ligsel # Get the binding energy from the B-factor... bindnrg = BFactorAtom (ligsel) # ...and the binding constant from the atomic property. bindconst = PropAtom (ligsel) if bindnrg<=0 bindconst=' None' else bindconst= 0000000000000.0000+bindconst # Get receptor residues that contact the ligand reslist() = ListRes (recsel) with distance<4.0 from (ligsel),Format='RESNAME MOLNAME RESNUM' if !count reslist reslist='No residues in contact' else reslist=join reslist result='(000+rec) | (000+run) | (000000.0000+bindnrg) | (bindconst) | (reslist)' return result,bindnrg # SORT KEYLIST BY VALUES, HIGHEST FIRST # ===================================== # Lists are passed by reference, see Yanaconda doc section 'Calls to user-defined functions..' def SortResults keylist,vallist,objlist caller (keylist),(vallist),(objlist) elements=count (keylist) objs=count (objlist) do sorted=1 for i=1 to elements-1 if (vallist)(i)<(vallist)(i+1) swap=(vallist)(i) (vallist)(i)=(vallist)(i+1) (vallist)(i+1)=swap swap=(keylist)(i) (keylist)(i)=(keylist)(i+1) (keylist)(i+1)=swap if objs==elements # A list of objects is also provided swap=(objlist)(i) (objlist)(i)=(objlist)(i+1) (objlist)(i+1)=swap sorted=0 while not sorted