# YASARA MACRO
# TOPIC:       5. Structure prediction
# TITLE:       Docking a covalently bound ligand to a receptor
# REQUIRES:    Structure
# AUTHOR:      Elmar Krieger
# LICENSE:     GPL
# DESCRIPTION: This macro runs VINA or AutoDock to predict the structure of a covalent 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,
# 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'

# 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 (alternatively you can provide
# the ligand as a *.yob file and fix certain dihedral angles only).
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). Ignored if you provide a scene file.
flexres='' 

# If receptor and ligand structures are provided separately, you need to select the
# two receptor and ligand hydrogen atoms that disappear when the covalent bond is formed.
# Either by adapting 'receptorsel' and 'ligandsel' below or by clicking 'Edit > Select':
receptorsel=''
ligandsel=''

# If receptor and ligand are provided covalently bound, you don't need to do anything as
# long as the ligand is the only covalently bound hetgroup. If there are more hetgroups
# or if the 'ligand' is a peptide itself, then you need to select the ligand.
# Either by adapting 'ligandsel' below or by clicking 'Edit > Select':
#ligandsel='Res PNM 800'

# Set to 1 if you also want to save PDB files of all clusters
pdbsaved=0

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

# 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

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

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

# Sanity checks
RequireVersion 19.1.27
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='start','receptor','ligand'

Clear
# Do we already have a start conformation scene of the complex, with covalently bound ligand?
scefilename='(MacroTarget)_start.sce'
scene = FileSize (scefilename)
if !scene
  # No, load PDB or YOb files of complex start conformation or separate receptor and ligand
  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),Model=1
        break
    if struct=='start' 
      if (struct)
        break
    elif !(struct)
      RaiseError 'The (struct) structure was not found at (MacroTarget)_(struct).* Please follow the instructions in the "Recipes" section exactly, especially when setting the macro target'
  if struct!='start'
    # Receptor and ligand loaded separately
    # Get the receptor/ligand hydrogens that disappear upon bond formation
    for struct in 'receptor','ligand'
      # Look for a permanently selected atom
      atomlist() = ListAtom selected Obj ((struct))
      if count atomlist
        (struct)sel='(join atomlist)'
      if (struct)sel==''
        RaiseError 'The (struct) hydrogen that disappears upon bond formation must be selected, either by providing "(struct)sel" in the macro, or via Edit > Select'
      atomlist() = ListAtom ((struct)sel)
      if count atomlist>1
        RaiseError 'Exactly one (struct) hydrogen that disappears upon bond formation must be selected, but (count atomlist) were found'
      element = ElementAtom ((struct)sel)
      if element!=1
        RaiseError 'Exactly one (struct) hydrogen that disappears upon bond formation must be selected, but atom ((struct)sel) is not a hydrogen'
      atomlist() = ListAtom all with bond to ((struct)sel)
      if count atomlist!=1
        RaiseError 'Exactly one heavy atom must be bound to the (struct) hydrogen that disappears upon bond formation must be selected, but (count atomlist) atoms were found'
      (struct)heavy = atomlist1
      # Pick a bond length of 1.5 A
      Distance ((struct)heavy),((struct)sel),Set=1.5
    # Did the user fix any receptor atoms and thus choose additional flexible side-chains?
    fixedatoms = CountAtom fixed Obj (receptor)
    if !fixedatoms
      # Fix receptor except ligand and the binding side-chain
      FixObj (receptor)
      FreeRes Atom (receptorheavy) Obj (receptor)
      FixAtom Backbone Obj (receptor)
    if flexres!=''
      FreeAtom Res (flexres) and Sidechain Obj (receptor) and not Res Ala Pro
    # Superpose the two atom pairs ligandhyd-ligandheavy and receptorheavy-receptorhyd
    SupOrderedAtom (ligandsel),(receptorheavy),(ligandheavy),(receptorsel)
    Cell Auto,Extension=10,Shape=Cuboid,Obj (ligand) or !fixed
    JoinObj (ligand),(receptor)
    AddBond (receptorheavy),(ligandheavy)
    DelAtom (receptorsel) or (ligandsel)
    start=receptor
  else
    # Starting structure is a complex
    receptor=start
    IdentifyLigand
    Cell Auto,Extension=10,Shape=Cuboid,Res Atom (ligandsel) or !fixed
  # Save the scene, for example to copy it for another project with slightly different ligand
  SaveSce (scefilename)
else
  # Yes, a scene is present. Here we cannot determine anymore what the ligand is, could be any flexible side-chain
  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
  if Objects>2
    ShowMessage 'Your scene contains (Objects) objects, while only the receptor and the docking cell are expected.'
    Wait ContinueButton
    HideMessage
  fixedatoms = CountAtom fixed
  if !fixedatoms
    # User provided a scene without fixed atoms, now we really need the ligand
    IdentifyLigand
# Warn if the cell is too large
x,y,z = Cell
if 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 object to make sure they can be identified later
NameObj (receptor),Receptor
# The segment name C*** is used to tag ligand Conformations, and cannot be used for the receptor
hit = ListAtom Segment C??? Obj Receptor
if hit
  MarkAtom (hit)
  segname = SegAtom (hit)
  RaiseError 'Segment names starting with "C" are unfortunately reserved to identify ligand 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'
# Docking is done without periodic boundaries
Boundary Wall
Longrange None  
# Perform the docking
Experiment Docking
  Method AutoDockLGA
  ReceptorObj (receptor)
  # The ligand is covalently bound, so there is no ligand object
  Runs (runs)
  ClusterRMSD (rmsdmin)
  # The number of cluster members to save individually, by default only the first member of each cluster is saved
  ClusterMembers 1
  # 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
  # Uncomment below to keep temporary files in the current working directory:
  # TmpFileID 3upo_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 'Covalent docking result analysis'
print '================================'
print
print '(runs) docking runs of the covalently bound ligand yielded the following results,'
print 'sorted by energy [more positive energies indicate stronger binding]'
print
print 'Run |Bind.energy[kcal/mol]|Dissoc. constant [pM]'
print '----+---------------------+---------------------'
clusters=0
for i=001 to runs
  result = DockingResult i,'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 after superposing on the receptor]' 
print
print 'Clu |Bind.energy[kcal/mol]|Dissoc. constant [pM]'
print '----+---------------------+---------------------'
for i=001 to clusters
  DelAll
  LoadYOb (MacroTarget)_(i)
  # Convert clusters to PDB format if requested
  if pdbsaved
    SavePDB 1,(MacroTarget)_(i)
  # Determine SegmentID of the ligand that seeds this cluster
  segnamelist(i) = SegAtom Segment C???
  # Ligands have SegmentID C001, C002 etc.
  result = DockingResult i,'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:'
LoadSce (MacroTarget)
for i=001 to clusters
  print '(MacroTarget)_(i).yob',Convert=No
StopLog
HideMessage
Console On
# 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 and binding energy, extracted from atomic BFactor and Property.
# 'num' is a sequential number, 'ligsel' selects the ligand (plus flexible side-chains)
def DockingResult num,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
  result='(000+num) |         (000000.0000+bindnrg) | (dissconst)'
  return result

# IDENTIFY THE LIGAND AND FIX ATOMS AS NEEDED
# ===========================================
def IdentifyLigand
  global ligandsel,flexres
  
  atomlist() = ListAtom selected
  if count atomlist
    ligandsel='(join atomlist)'
  if ligandsel==''
    # Identify ligand
    for connection in 'bond','arrow'
      atomlist() = ListAtom Hetgroup with (connection) to Protein
      if count atomlist
        break
    if !count atomlist
      RaiseError 'Could not find the ligand, i.e. a hetgroup bound to the protein. Please select the ligand with "ligandsel" in this macro or by clicking Edit > Select'
    if count atomlist>1
      reslist() = ListRes Atom (join atomlist),Format="RESNAME RESNUM"
      RaiseError 'Found more than one hetgroup bound to the protein: (reslist). Either delete unneeded hetgroups, or select the correct ligand with "ligandsel" in this macro or by clicking Edit > Select'
    ligandsel='(join atomlist)'
  # Also fix atoms as needed
  fixedatoms = CountAtom fixed
  if !fixedatoms
    FixAll
    FreeRes all with bond to Atom (ligandsel) or Atom (ligandsel)
    FixAtom Backbone and not (ligandsel)
  if flexres!=''
    FreeAtom Res (flexres) Sidechain and not Res Ala Pro