# YASARA MACRO
# TOPIC:       5. Structure prediction
# TITLE:       Docking a ligand to a receptor
# REQUIRES:    Structure
# AUTHOR:      Elmar Krieger
# LICENSE:     GPL
# DESCRIPTION: This macro predicts the structure of a ligand-receptor complex. It can also continue a docking run that got interrupted. An analysis log file is written at the end.

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

# Number of docking runs (maximally 999, each run can take up to an hour)
runs = 25

# 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.
ForceField AMBER03

# 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 'Too many docking runs selected, (runs) would take forever'

structlist='receptor','ligand'

# Load receptor and ligand
Clear
# Do we already have a receptor scene?
scefilename='(MacroTarget)_receptor.sce'
scene = FileSize (scefilename)
if !scene
  # No, load PDB or YOB files 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)
  RemoveObj (ligand)
  Cell Auto,Extension=(ligandradius*2)
  AddObj (ligand)
else
  # Yes, a scene is present
  LoadSce (scefilename)
  # The receptor is the object containing the first atom
  receptor = ListObj Atom 1
  if Objects>2
    ShowMessage 'Your scene contains (Objects) objects, while only the receptor and the docking cell are expected.'
    Wait ContinueButton
  # Load the ligand
  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)'
  ligandradius = RadiusObj (ligand)  
# 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
# Docking is done without periodic boundaries
Boundary Wall
Longrange None  
# Rename receptor and ligand objects to make sure they can be identified later
NameObj (receptor),Receptor
NameObj (ligand),Ligand
# Move the ligand out of the cell, so that it does not block the view
cellpos() = PosObj SimCell
PosObj (ligand),(cellpos1),(cellpos2+y*0.5+ligandradius),(cellpos3)
# Perform rigid docking if requested
if rigid
  FixObj Ligand
# Perform the docking
Experiment Docking
  Method AutoDockLGA
  ReceptorObj (receptor)
  LigandObj (ligand)
  Runs (runs)
  ClusterRMSD (rmsdmin)
  # Result file number must be separated with '_', in case MacroTarget also ends with number
  ResultFile (MacroTarget)_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
  # Uncomment below to stop after the setup step (combine with TmpFileID above to use YASARA only for setup)
  # SetupOnly yes
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 'Global 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 energies indicate stronger binding, and negative energies mean no binding]'
print
print 'Run |Bind.energy[kcal/mol]|Inhib.constant Ki[pM]| Contacting receptor residues'
print '----+---------------------+---------------------+-----------------------------'
clusters=0
for i=001 to runs
  # Ligands have SegmentID C001, C002 etc.
  result = DockingResult i,'Obj (receptor)','Segment C(i)'
  print (result)
  # Count the clusters
  exists = FileSize (MacroTarget)_(i).yob
  if exists
    clusters=clusters+1
  if !(i%10)
    ShowMessage 'Analyzing docking results, (100*i/runs)% completed...'
    Wait 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 'Clu |Bind.energy[kcal/mol]|Inhib.constantKi[pM]| Contacting receptor residues'
print '----+---------------------+--------------------+-----------------------------'
for i=001 to clusters
  DelAll
  LoadYOb (MacroTarget)_(i)
  # Ligands have SegmentID C001, C002 etc.
  result = DockingResult i,'Segment !C???','Segment C???'
  print (result)
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
StopLog
HideMessage
Console On
LoadSce (MacroTarget)

# EXTRACT DOCKING RESULT
# ======================
# Returns a result string and binding energy, extracted from atomic BFactor and Property.
# 'num' is a sequential number, 'recsel' selects the receptor, 'ligsel' the ligand.
def DockingResult num,recsel,ligsel
  # Get the binding energy from the B-factor...
  bindnrg = BFactorAtom (ligsel)
  # ...and the binding constant from the atomic property.
  inhibconst = PropAtom (ligsel)
  if bindnrg<=0
    inhibconst='               None'
  else
    inhibconst= 00000000000000.0000+inhibconst
  # 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+num) |         (000000.0000+bindnrg) | (inhibconst) | (reslist)'
  return result