This tutorial describes how to apply the 'flexible side chain' covalent docking as described in the paper:

   "Covalent docking using autodock: Two-point attractor and flexible side chain methods"
   http://www.ncbi.nlm.nih.gov/pubmed/26103917

If you use this method, please cite the paper.
Happy dockings!

==============
REQUIREMENTS
==============
- Python  >=2.6
- OpenBabel Python bindings
   (On Debian/Ubuntu derivatives, it can be installed with "apt-get install python-openbabel")
- Basic experience with docking with AutoDock 

==============
INSTALLATION 
==============
- Copy the directory 'adcovalent' in a location that is accessible from the user (i.e., '~/local/adcovalent').

==============
USAGE
==============

The preparation code is used to generate the geometry of the ligand covalently bound to the
residue. Once the ligand initial coordinate is generated, the PDBQT files of the receptor and
flexible ligand need to be generated using prepare_flexreceptor4.py script in AutoDockTools.

The directory '3upo_test' contains an example case of a ligand ready to processed.
The following examples assume a terminal opened in that location.


1. Generate the alignment
--------------------------
The input ligand structure must be modeled with a portion of the alkylated residue already present, 
i.e. for a ligand bound to a serine:

         LIGAND-OG-CB-CA-RESIDUE  [MOL1]

the ligand structure to process should be:
 
         LIGAND-O-C [MOL2]

where C and O are the two atoms that are going to be used for the alignment. More in general:

        0. Initial configuration

                    B
                   /
           [L]----A      I---J--[R]
                        /
                       X

        1. Translate [L] to overlap A->I

                    B
                   /
           [L]----A  ->  I---J--[R]
                        /
                       X

                          B
                         /
                 [L]----AI---J--[R]
                        /
                       X

        2. Rotate [L] around the normal of plane (B,A,J) to overlap B->J
              
                [L]
                  \
                   AI--BJ--[R]
                   / 
                  X 

The alignment could be skipped if the ligand is already attached to the receptor structure
(i.e., a re-docking experiment from a PDB is going to be perfoemd).

To perform the alignment, a receptor file and a residue name (Chain:ResNum) must
be specified, then a ligand file and an alignment criterion.

The ligand alignment can be defined in two ways. The first is by specifying 
the atom indices in the ligand file (i.e., first and second atoms):

    python  ~/local/adcovalent/prepareCovalent.py --ligand ligand.mol2 \
                  --ligindices 1,2\
                  --receptor 3upo_protein.pdb\
                  --residue B:SER222\
                  --outputfile ligcovalent.pdb

or, in a more general approach, by specifying a SMARTS pattern and the indices
of the alignment atoms in the pattern itself:

    python  ~/local/adcovalent/prepareCovalent.py --ligand ligand.mol2 \
                  --ligindices 4,3\
                  --ligsmart "C(=O)-O-C"\
                  --receptor 3upo_protein.pdb\
                  --residue B:SER222\
                  --outputfile ligcovalent.pdb

2. Generate PDBQT files
--------------------------
The standard PDBQT files for AutoDock need to be generated for both the receptor structure. The standard ADT scripts
are going to be used, assuming that $MGLROOT is the path where ADT has been installed.

*** NOTE: if a covalent ligand is already bound to the residue in the protein structure, it must be removed! ***

    $MGLROOT/bin/pythonsh $MGLROOT/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_receptor4.py -r 3upo_protein.pdb -A hydrogens
       [ this generates the file "3upo_protein.pdbqt" ]

and the covalent ligand that has been aligned:

    $MGLROOT/bin/pythonsh $MGLROOT/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_receptor4.py -r ligcovalent.pdb 
       [ this generates the file "ligcovalent.pdbqt" ]


3. Generate flexible/rigid PDBQT files
---------------------------------------
The PDBQT files are going to be used to generate the rigid and flexible components that are going to be used
for the docking. 
First, the receptor is processed to extract the rigid part by specifying which residue is going to be made flexible:

    $MGLROOT/bin/pythonsh $MGLROOT/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_flexreceptor4.py -r 3upo_protein.pdbqt -s 3upo_protein:B:SER222
       [ this generates the files "3upo_protein_rigid.pdbqt" and  "3upo_protein_flex.pdbqt" ]

The same thing is going to be done for the processed ligand:

    $MGLROOT/bin/pythonsh $MGLROOT/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_flexreceptor4.py -r ligcovalent.pdbqt -s ligcovalent:B:SER222
       [ this generates the files "ligcovalent_rigid.pdbqt" and  "ligcovalent_flex.pdbqt" ]

NOTE: if the ligand has not been prepared in step 1, make sure that all atoms in the ligand have the poper residue id
corresponding to the residue modified (i.e., Ser222A).


4. Generate GPF and DPF files
------------------------------
The parameter files for the actual calculation are going to be generated. The GPF for AutoGrid is generated with:

    $MGLROOT/bin/pythonsh $MGLROOT/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_gpf4.py -r 3upo_protein_rigid.pdbqt\
                -x ligcovalent_flex.pdbqt\
                -l ligcovalent_flex.pdbqt\
                -y -I 20\
                -o 3upo_priotein.gpf

The DPF for AutoDock is generated with:

    $MGLROOT/bin/pythonsh $MGLROOT/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_dpf4.py -r 3upo_protein_rigid.pdbqt\
                -x ligcovalent_flex.pdbqt\
                -l ligcovalent_flex.pdbqt\
                -o ligcovalent_3upo_protein.dpf\
                -p move='empty'

This command instructs the script to dock the covalent ligand as a flexible residue and ignore any 'true' ligand ("move='empty'"). To do this,
an empty file must be created, and it can be done with the following Unix command:
    
    touch empty

Also, in order to include the ligand atoms in the calculation of RMSD and clustering, the following parameters need to be added:

    rmsdatoms all

which will include flexible residue (and ligand) atoms in the RMSD calculation.
 Finally, the DPF file must be manually edited to set the appropriate energy model so that the docking score corresponds to the interaction between
the flexible residue (the ligand) and the rigid receptor. For this, the entry:
   
    unbound_model bound

must be replaced with:

    unbound_energy 0.0


5. Run AutoGrid and AutoDock
------------------------------
Both programs can be executed using the standard procedure:

    autogrid4 -p 3upo_priotein.gpf -l 3upo_priotein.glg
    autodock4 -p ligcovalent_3upo_protein.dpf -l  ligcovalent_3upo_protein.dlg



The output generated at each step of the process is available in the directory '3upo_test/output'.

