# YASARA MACRO
# TOPIC:       5. Structure prediction
# TITLE:       Docking a ligand to a receptor locally
# REQUIRES:    Structure
# AUTHOR:      Elmar Krieger
# LICENSE:     GPL
# DESCRIPTION: This macro samples local ligand conformations in a user-supplied complex. An analysis log file is written at the end.

# Parameter section - adjust as needed
# ====================================

# 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'

# Number of docking runs
runs = 250

# 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 for rigid docking
rigid = 1

# Forcefield to use for docking
ForceField Amber99
# Docking is done without periodic boundaries
Boundary Wall
Longrange None  

# 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 runs>999
  RaiseError "This macro can handle at most 999 docking runs"

structlist='receptor','ligand'

# Load receptor and ligand
Clear
# Do we already have a complex scene?
scefilename='(MacroTarget)_complex.sce'
scene = FileSize (scefilename)
if !scene
  # No, load PDB or YOB files of receptor and ligand.
  # Receptor and ligand coordinates are assumed to match
  for struct in structlist
    (struct)=0
    for type in 'yob','pdb'
      filename='(MacroTarget)_(struct).(type)'
      exists = FileSize (filename)
      if exists 
        (struct) = Load(type) (filename),Center=No
        break
    if not (struct)
      RaiseError 'The (struct) was not found at (filename)'
  # Orient the scene, create a docking cell that includes the ligand. 
  NiceOriAll
  Zoom Steps=1
  RemoveObj (receptor)
  Cell Auto,Extension=7
  AddObj (receptor)
  # Energy-minimize the receptor-ligand interaction
  FixAtom Protein Backbone
  Experiment Minimization
  Experiment On
  Wait ExpEnd
  # Save for reuse later
  SaveSce (scefilename)
else
  # Yes, a scene is present
  LoadSce (scefilename)
  Boundary Wall
  # The receptor is the object containing the first atom
  receptor = ListObj Atom 1
  # The ligand is the object containing the last atom
  ligand = ListObj Atom (Atoms)
  if ligand==receptor
    RaiseError 'Receptor and ligand must be present as separate objects'
  
# Warn if the cell is too large
x,y,z = Cell
if x>47 or y>47 or z>47
  ShowMessage "A cell axis is longer than 47 A, which may reduce docking accuracy. Consider focusing on the active site."
  Wait ContinueButton
  HideMessage
  
# Rename receptor and ligand objects to make sure they can be identified later
NameObj (receptor),Receptor
NameObj (ligand),Ligand

if rigid
  FixAtom Obj Ligand

# Perform the docking with local sampling
Experiment Docking
  Method AutoDockLS
  ReceptorObj (receptor)
  LigandObj (ligand)
  Runs (runs)
  ClusterRMSD (rmsdmin)
  ResultFile (MacroTarget)001
  # If you want to keep all temporary files, uncomment below
  # TmpFileID YourChoice
Experiment On
Wait ExpEnd
  
# Save a scene with receptor and all ligand conformations
SaveSce (MacroTarget)

# Save a log file with an analysis
Console Off
RecordLog (MacroTarget)
print 'Docking result analysis'
print '======================='
print
print '(runs) docking runs of the ligand object (ligand) to the receptor object (receptor) yielded the following results,'
print 'sorted by binding energy: [More positive values indicates stronger binding, and negative values mean no binding]'
print
print 'Run |Bind.energy[kcal/mol]|Bind.constant [pM]| Contacting receptor residues'
print '----+---------------------+------------------+-----------------------------'
clusters=0
for i=001 to runs
  bindnrg = BFactorAtom Segment C(i)
  bindconst = PropAtom Segment C(i) 
  if bindnrg<=0
    bindconst='            None'
  else
    bindconst=0000000000000.00+bindconst
  reslist() = ListRes Obj (receptor) with distance<4.5 from Segment C(i),Format='RESNAME MOLNAME RESNUM'
  print '(i) |         (000000.0000+bindnrg) | (bindconst) | (join reslist)'
  # Count the clusters
  exists = FileSize (MacroTarget)(i).yob
  if exists
    clusters=clusters+1
print
print 'After clustering the (runs) runs, the following (clusters) distinct complex conformations were found:'
print '[They all differ by at least (rmsdmin) A heavy atom RMSD]' 
print
print 'Cluster|Bind.energy[kcal/mol]|Bind.constant [pM]| Contacting receptor residues'
print '-------+---------------------+------------------+-----------------------------'
for i=001 to clusters
  DelAll
  LoadYOb (MacroTarget)(i)
  # Identify the last molecule in the soup, which is the ligand
  ligand = ListMol Atom (atoms)
  bindnrg = BFactorMol (ligand)
  bindconst = PropMol (ligand) 
  if bindnrg<=0
    bindconst='            None'
  else
    bindconst=0000000000000.00+bindconst
  reslist() = ListRes not Mol (ligand) with distance<4.5 from Mol (ligand),Format='RESNAME MOLNAME RESNUM'
  print '(i)    |         (000000.0000+bindnrg) | (bindconst) | (join reslist)'
print 
print 'The results of the (runs) runs have been combined in a single scene at (MacroTarget).sce',Convert=No
print 'The (clusters) clusters have been saved as:'
for i=001 to clusters
  print '(MacroTarget)(i).yob',Convert=No
LoadSce (MacroTarget)
StopLog