# YASARA MACRO
# TOPIC:       5. Structure prediction
# TITLE:       Docking a ligand to a receptor locally
# REQUIRES:    Structure
# AUTHOR:      Elmar Krieger and Kornel Ozvoldik
# LICENSE:     GPL
# DESCRIPTION: This macro runs VINA, AutoDock or Boltz to sample local ligand conformations in a user-supplied complex, possibly on the GPU if selected. 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,
# 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'

# Docking method, either VINALS (LS=Local Search), VINASO (SO=Score Only), AutoDockLS, or BoltzLS.
# Note that VINALS yields only a single result. If you want an ensemble of ligand conformations,
# choose method 'AutoDockLS' or BoltzLS below, but note that these do huge sampling steps and may
# even turn your ligand around.
# VINASO is not recommended, since a minimization should be done before calculating the energy.
# If left empty, a window will be shown to choose (GUI) or VINA will be used (console mode)
method=''

# Number of docking runs. VINA local search always yields a single result only.
runs=5

# 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 (alternatively you can provide
# the ligand as a *.yob file and fix certain dihedral angles only). Not supported by Boltz.
rigid=0

# A selection of receptor residues whose side-chains should be kept flexible, e.g.
# flexres = 'Lys 91 Leu 100', or (if the receptor is not a monomer)
# flexres = 'Res Lys 91 Mol A or Res Leu 100 Mol B'.
# Alternatively you can provide the receptor as a *.sce or *.yob file with fixed
# atoms, which gives better control (e.g. you can keep only part of a side-chain
# or even a terminal backbone flexible). Not supported by Boltz.
flexres='' 

# YASARA commands follow below, no '=' is used
# --------------------------------------------

# Uncomment below to set a certain random seed
# RandomSeed 1000

# Normally no change required below this point
# ============================================

# Do we have a target?
if MacroTarget==''
  RaiseError "This macro requires a target. Either edit the macro file or click Options > Macro > Set target to choose a target structure"

# Do we have a method?
if method==''
  if ConsoleMode
    # No, use VINA as default in console mode
    method='VINALS'
  else
    # No, ask user interactively
    methodnum,rmsdmin,runs,rigid =
      ShowWin Type=Custom,Title="Select options for local docking",Width=600,Height=400,
              Widget=Text,        X= 20,Y= 48,Text="Docking method (use GPU at Options > Processors > Set compute GPU)",
              Widget=RadioButtons,Options=4,Default=2,
                                  X= 30,Y= 80,Text="_A_utoDock",
                                  X= 30,Y=120,Text="_V_INA (CPU only)",
                                  X= 30,Y=160,Text="_V_INA (score only)",
                                  X= 30,Y=200,Text="_B_oltz",
              Widget=Text,        X=250,Y=130,Text="Cluster poses with minimum",
              Widget=NumberInput, X=250,Y=150,Text="R_M_SD",Default=(rmsdmin),Min=0,Max=20,
              Widget=Text,        X= 30,Y=256,Text="Docking",
              Widget=NumberInput, X= 30,Y=276,Text="_r_uns",Default=(runs),Min=1,Max=999,
              Widget=CheckBox,    X=250,Y=300,Text="_K_eep ligand rigid (not with Boltz)",Default=(rigid),
              Widget=Button,      X=548,Y=348,Text="_O_ K"
    methodlist='AutoDockLS','VINALS','VINASO','BoltzLS'
    method=methodlist(methodnum)

# Sanity checks
if runs>999
  RaiseError "This macro can handle at most 999 docking runs"
if method=='Boltz'
  if rigid
    RaiseError "The docking method 'Boltz' does not support rigid ligands. Set rigid=0 at the beginning of the macro"
  if flexres!=''
    RaiseError "The docking method 'Boltz' does not distinguish between fixed and flexible side-chains. Set flexres='' at the beginning of the macro"

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','sdf','mol2'
      filename='(MacroTarget)_(struct).(type)'
      exists = FileSize (filename)
      if exists and (type=='yob' or type=='pdb' or struct=='ligand')
        if type=='yob'
          (struct) = LoadYOb (filename)
        else
          (struct) = Load(type) (filename),Center=No
        # Make sure coordinate systems match
        TransferObj ((struct)),(receptor)
        break
    if not (struct)
      RaiseError 'The (struct) was not found at (MacroTarget)_(struct).* Please follow the instructions in the "Recipes" section exactly, especially when setting the macro target'
  # Orient the scene, create a docking cell that includes the ligand. 
  NiceOriAll
  ZoomAll Steps=1
  Cell Auto,Extension=7,Shape=Cuboid,Obj (ligand)
  ForceField NOVA,SetPar=Yes
  Boundary Wall
  # Remember which atoms are fixed, and additionally fix the backbone
  fixedlist() = ListAtom fixed
  FixAtom Protein Backbone Obj (receptor)
  # Energy-minimize the receptor-ligand interaction
  Experiment Minimization
  Experiment On
  Wait ExpEnd
  Free
  if count fixedlist
    FixAtom (join fixedlist)
  # Save for reuse later
  SaveSce (scefilename)
else
  # Yes, a scene is present
  LoadSce (scefilename)
  # Verify that the cell is present
  simcell = CountObj SimCell
  if !simcell
    RaiseError 'If you provide a scene, it must contain a simulation cell, but none was found in (scefilename)'
  # 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'
# Force field used for charge assignment (AutoDock only, not VINA) and to determine which ligand bonds can rotate.
# Do not change to an in vacuo force field (NOVA), all other force fields should work too (not tested).
ForceField AMBER03
# Docking is done without periodic boundaries
Boundary Wall
Longrange None
# Warn if the cell is too large
x,y,z = Cell
if method!='BoltzLS' and (x>96 or y>96 or z>96)
  ShowMessage "A cell axis is longer than 96 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
# Keep selected side-chains flexible
if flexres!=''
  FixObj Receptor
  FreeAtom Res (flexres) and Sidechain Obj Receptor
  FixRes Cys Atom SG with bond to Atom SG or Ala Pro and Obj Receptor
# The segment name C*** and R*** are used to tag ligand and Receptor Conformations, and cannot be used for the receptor
hit = ListAtom Segment C??? R??? Obj Receptor
if hit
  MarkAtom (hit)
  segname = SegAtom (hit)
  RaiseError 'Segment names starting with "C" or "R" are unfortunately reserved to identify ligand and receptor conformations, please click "Edit > Name > Segment", double-click object (receptor) on the right, and click OK to set an empty segment name'  
# The molecule name 'l' is reserved for the ligand
hit = ListAtom Mol l
if hit
  RaiseError 'The receptor must not contain a molecule named "l", please click Edit > Name > Molecule and try again'
# Perform rigid docking if requested
if rigid
  FixAtom Obj Ligand
# Perform the docking with local sampling
Experiment Docking
  Method (method)
  ReceptorObj (receptor)
  LigandObj (ligand)
  Runs (runs)
  ClusterRMSD (rmsdmin)
  # Filename to save cluster members, number must be separated with '_', in case MacroTarget also ends with number.
  # If PDB files are needed too, replace .yob with .pdb to get both
  ResultFile (MacroTarget)_001.yob
  # Uncomment below to set the number of energy evaluations (AutoDock ga_num_evals):
  # DockPar ga_num_evals 25000000
  # 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
  # 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)

# Clear surface environment to calculate contact surface
RemoveEnvAll
SurfPar Molecular=Gaussian
# Must be done during simulation so that the contact surface between different objects can be calculated
RemoveObj RecepFlex
Sim Pause
Sim On
# Save a log file with an analysis
Console Off
RecordLog (MacroTarget)
print 'Local docking result analysis'
print '============================='
print
print '(runs) (method) 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 '"Con.Surf" is the molecular contact surface.'
print
print 'Run |Bind.energy[kcal/mol]|Dissoc. constant [pM]| Con.Surf[A^2] | Contacting receptor residues'
print '----+---------------------+---------------------+---------------+-----------------------------'
clusters=0
for i=001 to runs
  # Ligands have SegmentID C001, C002 etc.
  recsel='Obj (receptor)'
  if method=='BoltzLS'
    recsel='Segment R(i)'
  result = DockingResult i,recsel,'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
Sim Off    
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]|Dissoc. constant [pM]| Con.Surf[A^2] | 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 
rndseed = RandomSeed
print 'The current random seed is (rndseed).'
fof = ForceField
print 'Point charges and dihedral barriers were obtained from the (fof) force field.'
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)
# Exit YASARA if this macro was provided as command line argument in console mode and not included from another macro
if runWithMacro and ConsoleMode and !IndentationLevel
  Exit

# EXTRACT DOCKING RESULT
# ======================
# Returns a result string. '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 dissociation constant from the atomic property.
  dissconst = PropAtom (ligsel)
  if bindnrg<=0
    dissconst='               None'
  else
    dissconst= 00000000000000.0000+dissconst
  # Calculate contact surface in both directions to average result
  consurf1=ConSurfAtom (recsel),(ligsel)
  consurf2=ConSurfAtom (ligsel),(recsel)
  # Get receptor residues that contact the ligand
  reslist() = ListRes (recsel) with distance<4.0 from (ligsel),Format='MOLNAME RESNAME RESNUM'
  if !count reslist
    reslist='No residues in contact'
  else
    reslist=join reslist
  result='(000+num) |         (000000.0000+bindnrg) | (dissconst) |     (000000.00+mean consurf) | (reslist)'
  return result