Prepare a reference structure

The INPHARMA program uses the PDB atom numbering for the calculations of the theoretical INPHARMA NOEs. It is therefore necessary to properly prepare the PDB files prior to the calculations. These steps models missing residues or atoms, and re-maps the atom numbering such that it is consistent between difference PDB files.

Process the structure file

As a reference, we will download a PDB file (PDB id 3DND) containing ligand (with id LL1).

pdb="3DNE"
ligand="LL1"

# Fetch the PDB file
wget -Oxray.pdb http://www.pdb.org/pdb/files/${pdb}.pdb

To fetch the ligand you can use the following PyMOL script. Note: This script will also align (superimpose) the x-ray structure onto a reference PDB file (all structures your prepare for INPHARMA should be aligned on the binding site).

# Extract the ligand 
cat<<EOF>align.pml
  # Load the xray PDB file
  load xray.pdb

  # A reference to align to
  load reference.amber.pdb, reference

  align xray, reference
  save xray_fit.pdb, xray

  select ligand, resname LL1
  save ligand.pdb, ligand
EOF

# Run the PyMOL script
pymol align.pml

# Extract ATOM records only
grep "^ATOM" xray_fit.pdb > xray_ATOMs.pdb

At this point your should have the following files in your directory:

xray_ATOMs.pdb

raw PDB file containing the receptor structure ONLY. Should be aligned to a reference PDB file.

ligand.pdb

raw PDB file of the ligand ONLY. Should be in the binding pocket of the receptor.

Modeller

The next step is to ensure that our receptor structure is complete (no missing residues/atoms, and the correct sequence). For this we will use two Modeller-script (Python).

The first script aligns the sequence of the structure (xray_ATOMs.pdb) and a sequence file. The sequence file should contain the protein sequence in Modeller format (named target_sequence.ali).

Requirement: Modeller should be installed and working.

>P1;target
sequence:target:::::::0.00: 0.00
HLDQFERIKTLGTGSFGRVMLVKHKETGNHYAMKILDKQKVVKLKQIEHTLNEKRILQAV
NFPFLVKLEFSFKDNSNLYMVMEYVPGGEMFSHLRRIGRFSEPHARFYAAQIVLTFEYLH
SLDLIYRDLKPENLLIDQQGYIQVTDFGFAKRVKGRTWTLCGTPEYLAPEIILSKGYNKA
VDWWALGVLIYEMAAGYPPFFADQPIQIYEKIVSGKVRFPSHFSSDLKDLLRNLLQVDLT
KRFGNLKNGVNDIKNHKWFATTDWIAIYQRKVEAPFIPKFKGPGDTSNFDDYEEEEIR*

from modeller import *

env = environ()
aln = alignment(env)
mdl = model(env, file='xray_ATOMs.pdb', model_segment=('FIRST:A','LAST:A'))
aln.append_model(mdl, align_codes='XRAY', atom_files='xray_ATOMs.pdb')
aln.append(file='target_sequence.ali', align_codes='target')
aln.align2d()

aln.write(file='alignment.ali', alignment_format='PIR')
aln.write(file='alignment.pap', alignment_format='PAP')

Run the sequence alignment with the following command:

python 01_align.py

Next script makes a new receptor PDB file based on the sequence alignment generated above:

from modeller import *
from modeller.automodel import *

log.verbose()
env = environ()

# directories for input atom files
env.io.atom_files_directory = ['.']

a = automodel(env, alnfile='alignment.ali',
              knowns='XRAY', sequence='target',
              assess_methods=(assess.DOPE, assess.GA341))

a.md_level = None 
a.rand_method = None
a.max_var_iteration = 1

a.starting_model = 1
a.ending_model = 1
a.make()

Run the script with the following commands, and re-name the output PDB file:

python 02_build_model.py
mv target.B99990001.pdb receptor_model.pdb

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

receptor_model.pdb

a PDB file of the receptor built with Modeller without missing residues, and with the correct sequence.

Amber

Finally, we will use Amber to add hydrogen atoms to our model structures. Using Amber will also allow for energy minimization, and molecular dynamics simulations, which can come in handy.

The first step is to generate parameters for the ligand. To do this, the ligand must be correctly protonated. There are several options to this. Babel is one of them:

babel -h -ipdb ligand.pdb -opdb ligand_H.pdb

Alternatively, you can download the ligand (already protonated) with this command:

wget -Oligand_H.pdb http://ligand-expo.rcsb.org/reports/${ligand:0:1}/${ligand}/${ligand}_${model}.pdb

Note: Protonating ligands might not be trivial, and there is no guarantee that these programs yields reasonable protonation.

Once hydrogens are (correctly) added to the ligand, we will use Antechamber to build ligand parameters:

antechamber -i ligand_H.pdb -fi pdb -o prepi -fo prepi -c bcc -s 2 -nc 0 
parmchk -i prepi -f prepi -o frcmod

Amber paramter/topolgy, and coordinate files can be generated with the following Leap script:

cat<<EOF> build
  source leaprc.ff10 
  source leaprc.gaff
  loadAmberPrep prepi
  loadamberparams frcmod

  set default PBradii mbondi2

  receptor = loadpdb receptor_model.pdb
  ligand = loadpdb ligand.pdb

  saveamberparm receptor receptor.prmtop receptor.inpcrd
  saveamberparm ligand ligand.prmtop ligand.inpcrd

  complex = combine { receptor ligand }
  saveamberparm complex complex.prmtop complex.inpcrd
EOF 

tleap -f build

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

receptor.prmtop receptor.inpcrd

Amber topology and coordinate files of the receptor.

complex.prmtop complex.inpcrd

Amber topology and coordinate files of the receptor-ligand complex.

ligand.prmtop ligand.inpcrd

Amber topology and coordinate files of the ligand.

Convert back to PDB files

# Convert to PDB-files
ambpdb -p complex.prmtop < complex.inpcrd > complex.amber.pdb
ambpdb -p receptor.prmtop < receptor.inpcrd > receptor.amber.pdb

# Convert to SYBYL mol2-files
top2mol2 -p receptor.prmtop -c receptor.inpcrd -o receptor.sybyl.mol2
top2mol2 -p ligand.prmtop -c ligand.inpcrd -o ligand.sybyl.mol2

# Fetch the ligand, and save to a seperate PDB file
grep ${ligand} complex.amber.pdb > ligand.amber.pdb

Energy minimization

To get rid of unreasonable hydrogen-hydrogen distances, we will perform a energy minimization using Amber:

CPUS=4
ligand="LL1"
ligand_resi="299"

cat<<EOF>min.in
short minimization
 &cntrl
  imin=1, maxcyc=1000, ncyc=700,
  cut=15., rgbmax=12., ntb=0,
  igb=2, saltcon=0.2,
  ntr=1, restraintmask=':${ligand_resi}',
  restraint_wt=0.5,
  ntpr=100
/
EOF

$MPI_HOME/bin/mpirun -np $CPUS $AMBERHOME/bin/sander.MPI -O -i min.in -o min.out \
               -p complex.prmtop -c complex.inpcrd -r min.rst -ref complex.inpcrd

# Convert to PDB files
ambpdb -p complex.prmtop < min.rst > min.amb.pdb
grep ${ligand} min.amb.pdb > ${ligand}.amb.pdb

Convert complexed structure to a receptor structure:

cat<<EOF>stripLigand.ptraj
  trajin min.rst
  strip :${ligand_resi}
  trajout receptor.min.amb.rst restart
EOF

ptraj complex.prmtop stripLigand.ptraj
mv receptor.min.amb.rst.1 receptor.min.amb.rst

ambpdb -p receptor.prmtop < receptor.min.amb.rst > receptor.min.amb.pdb

# Convert to SYBYL mol2-files
top2mol2 -p receptor.prmtop -c receptor.min.amb.rst -o receptor.sybyl.mol2

Summary

This description should provide the steps needed to build the reference PDB files needed for INPHARMA experiments. The INPHARMA program will depend on the atom numbering of the files complex.amber.pdb, receptor.amber.pdb, and ligand.amber.pdb.