Protein-ligand docking

The calculation of theoretical INPHARMA NOEs is (very) sensitive to hydrogen-hydrogen distances. It is therefore necessary to perform energy minimization after protein-ligand docking. It is moreover important that the atom numbering is consistent throughout the docking modes, and that all protein-ligand complexes are properly aligned on the binding-site.

These steps are illustrated below using docking with Glide and minimization with Amber.

Glide docking

With Glide we will use the prepwizard to generate the maestro files (.mae) needed for the docking. We will dock to the receptor structure generated in step Prepare a reference structure.

Prepare the receptor structures:

# Prepare receptor for Glide docking
reduce -Trim receptor.min.amb.pdb > receptor.pdb
prepwizard -WAIT receptor.pdb receptor.mae

# Prepare ligand for Glide docking
reduce -Trim ligand.amber.pdb > ligand.pdb
prepwizard -WAIT ligand.pdb ligand.mae

cat<<EOF>glide_grid.in
 USECOMPMAE YES
 INNERBOX 10, 10, 10
 ACTXRANGE 30.000000
 ACTYRANGE 30.000000
 ACTZRANGE 30.000000
 GRID_CENTER 6.989,  10.156,   2.855
 OUTERBOX 30.000000, 30.000000, 30.000000
 ENTRYTITLE receptor
 RECEP_FILE receptor.mae
 GRIDFILE glide-grid.zip
EOF


cat<<EOF>glide_dock.in
 USECOMPMAE YES
 RINGCONFCUT 2.500000
 POSES_PER_LIG 50
 GRIDFILE glide-grid.zip
 LIGANDFILE ligand.mae
EOF

# Make the grid
glide -WAIT glide_grid.in

# Run docking
glide -WAIT glide_dock.in

# Convert to PDB files
pdbconvert  -imae glide_dock_pv.maegz -opdb glide_dock.pdb

This step should provide you with following file(s):

glide_dock-1.pdb, glide_dock-2.pdb, ...

A set of docking modes in PDB format

Minimization of docking modes

After docking, minimization of each docking mode should be carried out. Below is a bash script which can be used:

#!/bin/bash

# Number of CPUs to use for minimization
CPUS=4

# Directory for reference files
reference_dir="../01_build"

# Residue name of ligand
ligand_name="LL1"

# Directory where minimized files should be dumped
output_dir="pdbs_min"

# Define workdir as an anchor
base=`pwd`

# Make the output directory
mkdir $output_dir

# Fetch files needed for minimization
cp ${reference_dir}/frcmod .
cp ${reference_dir}/prepi .
ln -s ${reference_dir}/receptor.min.amb.pdb .

# 500 steps minimization should be enough
cat<<EOF>min.in
short minimization
 &cntrl
  imin=1, maxcyc=500, ncyc=350,
  cut=15., rgbmax=12., ntb=0,
  igb=2, saltcon=0.2,
  ntpr=10,
/
EOF

# Make the Leap-input file
cat<<EOF>build
source leaprc.ff10 
source leaprc.gaff
loadAmberPrep prepi
loadamberparams frcmod
           
set default PBradii mbondi2
                  
receptor = loadpdb receptor.pdb
ligand = loadpdb ligand.pdb
                           
complex = combine { receptor ligand }
saveamberparm complex complex.prmtop complex.inpcrd
EOF


# The first file is the receptor - rename it
if [ -s glide_dock-1.pdb ]; then
        mv glide_dock-1.pdb glide_dock-receptor.pdb
fi

# Rename the docking modes
if [ -s "glide_dock-2.pdb" ]; then
        i=1
        for i in $(seq 1 200) 
        do
                let i=${i}+1
                file=glide_dock-${i}.pdb
                if [ ! -s $file ]; then
                        break
                else
                        let j=${i}-1
                        n=$(printf %03d $j)
                        mv glide_dock-${i}.pdb glide_dock-${n}.pdb
                fi
        done
fi

# Find all docked PDB files
pdbfiles=`ls glide_dock-???.pdb`

# Loop the files - minimize each one of them
i=0
for file in ${pdbfiles}
  do

  cd ${base}

  # Make a temp directory for each docking mode
  let i=${i}+1
  n=$(printf %03d $i)
  tmpdir=_tmp_${n}

  mkdir $tmpdir
  cd $tmpdir
  workdir=`pwd`
  echo "Current workdir is "${workdir}

  ## Copy required files to the dir
  cp ../$file ligand.glide.pdb
  ln -s ../receptor.min.amb.pdb receptor.pdb
  ln -s ../prepi .
  ln -s ../frcmod .
  ln -s ../build . 
  ln -s ../min.in .

  reduce -Trim ligand.glide.pdb > ligand.noH.pdb
  grep $ligand_name ligand.noH.pdb > ligand.pdb

  outfile=${base}/${output_dir}/${file}

cat <<eof >launch_me.sh
  # Generate Amber files
  \$AMBERHOME/bin/tleap -f build
  
  # Run Minimization
  \$MPI_HOME/bin/mpirun -np $CPUS \$AMBERHOME/bin/pmemd.MPI -O -i min.in -o min.out \ 
                  -p complex.prmtop -c complex.inpcrd -r min.rst
  # Convert to PDB format
  \$AMBERHOME/bin/ambpdb -p complex.prmtop < min.rst > min.amb.pdb

  # Copy to output directory
  cp min.amb.pdb ${outfile}
eof

  chmod +x launch_me.sh
  if [ `hostname` == "submaster" ]; then
      bsub -n ${CPUS} -R "span[hosts=1]" ./launch_me.sh
  else
      echo "Running local version - This will take some time..."
      ./launch_me.sh
  fi
  
  echo ${outfile}
  sleep 1;
done

# Finally delete all temp directories
rm -rf _tmp_???

This step should provide you with a directory:

pdbs_min

Containing a set of minimized docking modes

Check atom numbering

If everything went fine, you should now have a set of minimized docking modes with the same atom numbering as your reference protein-ligand complex.

grep LL1 ../01_build/complex.amber.pdb | awk '{print $2"\t"$3"\t"$4}' > a
grep LL1 pdbs_min/glide_dock-001.pdb | awk '{print $2"\t"$3"\t"$4}' > b

diff a b