DISI - User contributions [en]

FlexPepDock

2013-05-17T17:54:26Z

Oliv Eidam: /* 1. Create peptide coordinates in Pymol */

=FlexPepDock: Wanna dock a peptide? Use FlexPepDock!=

[[File:Rec w arg327.png|thumb|right|alt=X-ray structure of Neurotensin receptor. The C-terminus of NT(8-13) was predicted to interact with Arg327, which is depicted with spheres. |X-ray structure of Neurotensin receptor. The C-terminus of NT(8-13) was predicted to interact with Arg327, which is depicted with spheres.]]

Interested in peptide docking? Use FlexPepDock, which is a peptide docking software implemented in Rosetta. The easiest way to do that is to use the online server: 
http://flexpepdock.furmanlab.cs.huji.ac.il

If you want to run FlexPepDock locally, follow the steps below. 
In my example I docked the hexapeptide NT(8-13) with the sequence RRPYIL into the neurotensin receptor. I did this prospectively, without having looked at the bound peptide conformation before doing this. All I knew was that the C-terminal carboxylate interacts with Arg327 ([http://www.jbc.org/content/275/1/328 Barroso et al, JBC, 2000]) and that the peptide adopts an extended conformation ([http://www.pnas.org/content/100/19/10706 Luca et al, PNAS, 2003]).

==1. Create peptide coordinates in Pymol==

[[File:Rec w input peptides.png|thumb|right|alt=Backbone representation of manually docked starting peptides (orange and green). |Backbone representation of manually docked starting peptides (orange and green).]]

You can generate peptide coordinates using Pymol (Build->Residue->Alanine). If you want an extended peptide, it is best to generate coordinates for poly-Ala first, and then mutate to your peptide using the mutation wizard (Wizard->Mutagenesis->Mutate to Arg->Apply). 
Save molecule.
Do NOT add a amine at the N-terminus and do NOT add a C-terminal carboxylate to the pymol coordinates!
'''Important note:''' FlexPepDock does not change your Omega angles, so make sure that a proline is in trans if that's what you want!!

==2. Dock peptide coordinates in Pymol to generate a starting model for FlexPepDock ==

You need a rough starting peptide model for FlexPepDock. The peptide starting model should be within 5 Angstroem RMSD to the native structure. You can manually dock your peptide created above in Pymol into to the receptor by switching to Editing Mode, and dragging/rotating the peptide by holding the Shift-Button and Middle and Left-Mouse button, respectively. Don't worry about clashes of the peptide with the protein: the most important thing is that the backbone is within 5 Anstroem of the native structure.

Save the docked peptide coordinates and '''add them at the below of the PDB coordinates''' of the apo receptor structure. It is important to add them below the receptor coordinates, separated by a TER statement. Also, there should be no END statements.

==3. Run prepack.sh==

Run prepack.sh (see code below):

./prepack.sh NTS1_rrpyil_input.pdb

prepack.sh generates pNTS1_rrpyil_input_0001.pdb as output: the protein is protonated and a N-terminal amine and a C-terminal carboxylate have been added to the peptide.

'''Code of prepack.sh:'''

#!/bin/csh

#$1 is the start.pdb

set arch = `uname -p`
if ( $arch == 'x86_64') then
~londonir/rosetta/rosetta_source/bin/FlexPepDocking.static.linuxgccrelease -s $1 -database /raid1/people/londonir/rosetta/rosetta_database -native $1 -flexpep_prepack -ex1 -ex2aro -unboundrot $1 > log.prepack
else
echo "I can only run on a x86_64 system..."
endif

==4. Run FlexPepDock on Cluster using submit_fpdock.sh==

[[File:Rec w output peptides.png|thumb|right|alt=Backbone representation of best scoring output peptides (orange and green). |Backbone representation of best scoring output peptides (orange and green).]]

Run FlexPepDock on Cluster using submit_fpdock.sh (see code below):

./submit_fpdock.sh NTS1_rrpyil_input_0001.pdb 200

submit_fpdock.sh takes two arguments: the protonated protein-peptide input model from prepack.sh (NTS1_rrpyil_input_0001.pdb) and the second argument (200) is the number of models you want to generate - 200 is a good number. submit_fpdock.sh calls a script called single_fpdock.sh (see code below), which executes the actual peptide docking. The peptide docking is pretty fast: it takes about 2 minutes per model.

'''Code of submit_fpdock.sh:'''

#!/bin/csh

#$1 is the start file (start.pdb)
#$2 is the number of jobs to send (each will attempt nstruct=1)

foreach i (`seq 1 $2`)
qsub -l arch=lx24-amd64 -q all.q -cwd -e ./error -o ./out -v arg1=$i,arg2=$1 ./single_fpdock.sh
end

'''Code of single_fpdock.sh:'''

#!/bin/csh

#arg1 is the run prefix number
#arg2 is the start.pdb (and native.pdb)

set arch = `uname -p`
if ( $arch == 'x86_64') then
~londonir/rosetta/rosetta_source/bin/FlexPepDocking.static.linuxgccrelease -s $arg2 -database /raid1/people/londonir/rosetta/rosetta_database -native $arg2 -pep_refine -ex1 -ex2aro -use_input_sc -nstruct 1 -unboundrot $arg2 -out:prefix $arg1'.' > log.$arg1
endif

==5. Analyze output models==

[[File:Pose1 input output vs xray.png|thumb|right|alt=Comparison of native peptide conformation (slate) to best scoring output model (orange). The backbone of the starting peptide is depicted in yellow. |Comparison of native peptide conformation (slate) to best scoring output model (orange). The backbone of the starting peptide is depicted in yellow.]]

FlexPepDock outputs a pdb, a score file (.sc) and a log file (log) for each model requested. The score files is a tab-delimited file with many scores and other interesting numbers like Interface buried surface area (I_bsa), number of hydrogen bonds (I_hb), score for the Interface (I_sc) and a score for the peptide (pep_sc). Your main interest is most likely the total score (total_score).

You can print the total scores of each model using this script:

#!/bin/csh
# prints total scores for each model

foreach file (*.sc)
echo $file":" `awk 'NR==3 {print $2}' $file`
end

In my prospective peptide docking, the best scoring output model starting from the orange pose had a total score of -324, which was better than the score of the green starting pose (-316).

==6. Calculate RMSD values compared to native structure==

FlexPepDock automatically outputs RMSD values of the models compared to the input structure (rmsALL, rmsBB etc.). But if you know the native structure, you are probably interested how the models compare to the native structure.

You can calculate RMSD values for each model compared to the native structure using the script rescore_batch.csh.
First, you have to generate a list of pdb IDs for which you want to generate RMSDs compared to the native structure:

ls -1 *.pdb > list_of_pdbs.txt

Then run rescore_batch.csh (see code below):
./rescore_batch.csh

In my prospective prediction, the orange output model had a backbone RMSD of 1.9 Angstroem (3.4 A over HA all atoms), which is a freaking good prediction. The the backbones align perfectly over the last four amino acids, and the two arginines at the N-terminus are actually (and surprisingly) not perfectly ordered in the crystal structure. 
As a side note: a few models had RMSD values as low as 1.2 A, but scored worse.

'''Code of rescore_batch.csh:'''

#!/bin/csh

#-s starting
#-native: reference for rmsd calc
# -l list file

~londonir/rosetta/rosetta_source/bin/FlexPepDocking.static.linuxgccrelease -l list_of_pdbs.txt -database /raid1/people/londonir/rosetta/rosetta_database -native 4GRV_AB_ATOM.pdb -out:prefix rescore_list'.' -flexpep_score_only -out:pdb false > log.rescore_list

#create list with pdbs
# ls -1 *.pdb > list_of_pdbs.txt

This script takes a few minutes to run and outputs a file called log.rescore_list. For each model, you can print the scores, backbone RMSDs and model name, sorted by backbone RMSD, by typing this command:

awk '{print $2, $41, $53}' rescore_list.score.sc | sort -rgk 2

=Acknowlegements=
Many thanks to Nir for compiling Rosetta on our server, the introduction into FlexPepDock and the fp* scripts.

Dock Sampling

2013-05-08T21:48:29Z

Oliv Eidam:

=Pose Sampling in DOCK=

Don't get right DOCK pose? Well, this happens a lot... Docking is a complicated process, but it can be broken down to two main processes: scoring and sampling. This page is about sampling (in DOCK).

To investigate why you don't get the right DOCK pose (in the case you know the answer from a crystal structure), you can consider the following things:

1. Check if the bioactive conformation is sampled in docking. This is described here: http://wiki.uoft.bkslab.org/index.php/Screen3d
 
2. Dock the bioactive conformation rigidly. This is described here: http://wiki.uoft.bkslab.org/index.php/Multimol2db.py
 
3. Run scoreopt on the bioactive conformation to get a score for the crystallographic pose of your ligand. This is described here: http://wiki.uoft.bkslab.org/index.php/Scoreopt
 
4. Check if a bioactive orientation (or an orientation close to it is actually sampled in the docking process. This is what is described on this page.

==Bioactive orientation sampled?==

[[File:Nz14 xtal vs dock.png|thumb|right|alt=Prediction by DOCK depicted in yellow, actual crystal structure depicted in green. |Prediction by DOCK depicted in yellow, actual crystal structure depicted in green. ]]

I want to make the following example: I predicted a fragment binding to the oxyanion hole, but when we determined the crystal structure, the fragment bound to the so-called distal site of AmpC (see picture on the right). In this case I was wondering whether a orientation similar to the crystallographic pose was actually sampled during the docking process.

'''It requires two steps to address this question:''' 
1. Use Niu Huang's ligand clustering tool (described here: http://wiki.uoft.bkslab.org/index.php/Dock_Ligand_Clustering) to output all orientations sampled during the docking process. 
2. Calculate RMSD values for all orientations compared to the crystal structure pose. You can do this easily and fast using DOCK 6.6 as described here: http://wiki.uoft.bkslab.org/index.php/Calculate_DOCK6_RMSD 

'''In this process you will re-dock your ligand and there are (at least) two caveats to this process:''' 
1. Niu Huang's ligand clustering tool is a branch of DOCK 3.6 and therefore (most likely) not the same DOCK binary which you have been using in your docking. 
2. If you re-dock your ligand starting from a new db file for your ligand (because your ligand was part of a big library and the original db-file is no more available which is often the case when you docked a big ZINC library), then you won't dock the same Omega conformations and the partial charges generated by Amsol (and therefore also the desolvation energies) will also be changed. 

That said: it is the best approximation at this point.

===1. Output all orientations sampled during docking===

To output all orientations sampled during docking, you can add the ~30 lines of code at the bottom of your original INDOCK file (see here how it works: http://wiki.uoft.bkslab.org/index.php/Dock_Ligand_Clustering), and re-dock the ligand into the original grids.
Make sure you set pose_clustering_method to 0 in your INDOCK for ligand clustering (to ensure all poses get written out).

This will output a file called test.eel1.clu, which is a eel1 file with all poses sampled during the docking process. If you want to know how many poses were sampled, type:

grep -c family test.eel1.clu

The answer was 2293 for the small fragment discussed here.

===2. Calculate RMSD values for all orientations compared to the crystal structure pose===

[[File:Nz14 xtal vs dock vs sampled.png|thumb|right|alt=Prediction by DOCK depicted in yellow, actual crystal structure depicted in green, orientation with lowest RMSD value (1.3A) depicted in magenta. |Prediction by DOCK depicted in yellow, actual crystal structure depicted in green, orientation with lowest RMSD value (1.3A) depicted in magenta. ]]

To calculate RMSD values for all orientations generated above (test.eel1.clu) using DOCK 6.6. you have to do execute two programs: 

'''1. Convert test.eel1.clu to a mol2 file without hydrogens''' (you need the same number of atoms to calculate RMSD values and your crystal structure likely does not have hydrogens). You can use the script eel12pdb_no_H.csh (see code below) to do this:

./eel12pdb_no_H.csh test.eel1.clu

This will output a file called test.eel1.clu_no_H.mol2.

'''Code of eel12pdb_no_H.csh:'''

#!/bin/sh -f

# awk script to convert eel1 file to pdb
# run:
# ./eel12pdb.csh file.eel1

awk '{
if ($1=="ATOM")
print substr($0,0,56)"1.00 20.00 "substr($0,14,1);
else if ($1=="TER")
print $0"\n""END";
else
print $0;
}' $1 > $1.pdb

# Remove hydrogens
awk '($3!~"H")' $1.pdb > $1_no_H.pdb

# convert to mol2 format
convert.py --i $1_no_H.pdb --o $1_no_H.mol2

'''2. Calculate RMSD values for each orientation using DOCK 6.6''' as described here: http://wiki.uoft.bkslab.org/index.php/Calculate_DOCK6_RMSD

calc_dock6_rmsd.csh test.eel1.clu_no_H.mol2 nz14_xtal_pose.mol2 | grep RMSDh | sort -rgk 3

The result of this analysis was that there about a dozen orientations sampled with a RMSD < 2 Angstroem. The closest orientation had a RMSD value of 1.3 Angstroem compared to the crystal structure and is depicted on the right in magenta. Its DOCK score was -30, which was much lower than than the best scoring orientation (DOCK score: -80).

===3. Optional: calculate electrostatics, vdw and desolvation for a DOCK orientation ===

One little problem of the ligand clustering DOCK binary is that it outputs only the total DOCK score for each orientation sampled. And for some strange reason it prefers not to output the desolvation and atom type columns in test.eel1.clu. But you may be interested in the individual terms (electrostatics, vdw and desolvation) contributing to the total DOCK score of an individual orientation!

To get the individual terms of a orientation output in test.eel1.clu, you can use the the script append_solv2eel1_file.csh (see code below).
Before running this script, you have to copy a individual orientation from test.eel1.clu into a new file. In the example below I copied the fifth highest ranking orientation. Copy the three REMARK lines, all ATOM lines, and the TER statement. 
Then run:

./append_solv2eel1_file.csh test.eel1.clu5

This creates a file called test.eel1.clu5_solv, with desolvation energies and atom types extracted from test.eel1, which is output for the highest ranking orientation. 
You can then run [[Scoreopt]]:

/raid5/people/mysinger/code/mud/trunk/doscoreopt.csh test.eel1.clu5_solv ../../../grids > test.eel1.clu5_solv.scoreopt.log

The output file test.eel1.clu5_solv.scores summarizes then the individual terms (electrostatics, vdw and desolvation) contributing to the total DOCK score.

Code of append_solv2eel1_file.csh:

#!/bin/bash

awk '{
if ($1=="ATOM")
print substr($0,65,10);
if ($1=="TER")
print $0;

}' test.eel1 > solvation.txt

awk '{
if (($1=="ATOM") || ($1=="TER"))
print $0;
}' $1 > $1_ATOM

awk '
BEGIN {
# load array with contents f1.txt
while ( getline < "solvation.txt" > 0 )
{
f1_counter++
f1[f1_counter] = $0
}
}
{
if ($1=="ATOM")
print $0,"",f1[NR];
else
print $0;
} ' $1_ATOM > $1_solv

=Acknowledgements=
Thanks to RGC for lots of docking advice and [[Multimol2db.py]]. Thanks to MM for [[Scoreopt]]. Thanks to Niu Huang and his group for [[Dock Ligand Clustering]]: it's awesome!

Dock Sampling

2013-05-08T21:37:55Z

Oliv Eidam:

Dock Sampling

2013-05-08T21:00:35Z

Oliv Eidam: /* 2. Calculate RMSD values for all orientations compared to the crystal structure pose */

2013-05-08T19:23:27Z

Tools for protein and ligand analysis

2013-05-08T18:00:09Z

Oliv Eidam:

*[http://wiki.uoft.bkslab.org/index.php/Plop Plop]: minimize ligand and active site
*[http://wiki.uoft.bkslab.org/index.php/Omega.parm omega.parm]: change settings for ligand conformation generation
*[http://wiki.uoft.bkslab.org/index.php/screen3d screen3d]: align ligands in 3D
*[http://wiki.uoft.bkslab.org/index.php/Phenix Phenix]: calculate electron densities and display them in Pymol
*[http://wiki.uoft.bkslab.org/index.php/filter.py: filter.py]: calculate atoms distances of DOCK poses
*[http://wiki.uoft.bkslab.org/index.php/FlexPepDock FlexPepDock]: peptide docking
*[http://wiki.uoft.bkslab.org/index.php/Dock_Sampling DOCK Sampling]: pose sampling in DOCK
*[http://wiki.uoft.bkslab.org/index.php/Calculate_DOCK6_RMSD DOCK6 RMSD]: calculate ligand RMSDs

FlexPepDock

2013-05-08T00:09:17Z

Oliv Eidam: /* Acknowlegements */

=FlexPepDock: Wanna dock a peptide? Use FlexPepDock!=

[[File:Rec w arg327.png|thumb|right|alt=X-ray structure of Neurotensin receptor. The C-terminus of NT(8-13) was predicted to interact with Arg327, which is depicted with spheres. |X-ray structure of Neurotensin receptor. The C-terminus of NT(8-13) was predicted to interact with Arg327, which is depicted with spheres.]]

Interested in peptide docking? Use FlexPepDock, which is a peptide docking software implemented in Rosetta. The easiest way to do that is to use the online server: 
http://flexpepdock.furmanlab.cs.huji.ac.il

If you want to run FlexPepDock locally, follow the steps below. 
In my example I docked the hexapeptide NT(8-13) with the sequence RRPYIL into the neurotensin receptor. I did this prospectively, without having looked at the bound peptide conformation before doing this. All I knew was that the C-terminal carboxylate interacts with Arg327 ([http://www.jbc.org/content/275/1/328 Barroso et al, JBC, 2000]) and that the peptide adopts an extended conformation ([http://www.pnas.org/content/100/19/10706 Luca et al, PNAS, 2003]).

==1. Create peptide coordinates in Pymol==

[[File:Rec w input peptides.png|thumb|right|alt=Backbone representation of manually docked starting peptides (orange and green). |Backbone representation of manually docked starting peptides (orange and green).]]

You can generate peptide coordinates using Pymol (Build->Residue->Alanine). If you want an extended peptide, it is best to generate coordinates for poly-Ala first, and then mutate to your peptide using the mutation wizard (Wizard->Mutagenesis->Mutate to Arg->Apply). 
Save molecule.
Do NOT add a amine at the N-terminus and do NOT add a C-terminal carboxylate to the pymol coordinates!
'''Important note:''' FlexPepDock does not change your Psi angles, so make sure that a proline is in trans if that's what you want!!

==2. Dock peptide coordinates in Pymol to generate a starting model for FlexPepDock ==

You need a rough starting peptide model for FlexPepDock. The peptide starting model should be within 5 Angstroem RMSD to the native structure. You can manually dock your peptide created above in Pymol into to the receptor by switching to Editing Mode, and dragging/rotating the peptide by holding the Shift-Button and Middle and Left-Mouse button, respectively. Don't worry about clashes of the peptide with the protein: the most important thing is that the backbone is within 5 Anstroem of the native structure.

Save the docked peptide coordinates and '''add them at the below of the PDB coordinates''' of the apo receptor structure. It is important to add them below the receptor coordinates, separated by a TER statement. Also, there should be no END statements.

==3. Run prepack.sh==

Run prepack.sh (see code below):

./prepack.sh NTS1_rrpyil_input.pdb

prepack.sh generates pNTS1_rrpyil_input_0001.pdb as output: the protein is protonated and a N-terminal amine and a C-terminal carboxylate have been added to the peptide.

'''Code of prepack.sh:'''

#!/bin/csh

#$1 is the start.pdb

set arch = `uname -p`
if ( $arch == 'x86_64') then
~londonir/rosetta/rosetta_source/bin/FlexPepDocking.static.linuxgccrelease -s $1 -database /raid1/people/londonir/rosetta/rosetta_database -native $1 -flexpep_prepack -ex1 -ex2aro -unboundrot $1 > log.prepack
else
echo "I can only run on a x86_64 system..."
endif

==4. Run FlexPepDock on Cluster using submit_fpdock.sh==

[[File:Rec w output peptides.png|thumb|right|alt=Backbone representation of best scoring output peptides (orange and green). |Backbone representation of best scoring output peptides (orange and green).]]

Run FlexPepDock on Cluster using submit_fpdock.sh (see code below):

./submit_fpdock.sh NTS1_rrpyil_input_0001.pdb 200

submit_fpdock.sh takes two arguments: the protonated protein-peptide input model from prepack.sh (NTS1_rrpyil_input_0001.pdb) and the second argument (200) is the number of models you want to generate - 200 is a good number. submit_fpdock.sh calls a script called single_fpdock.sh (see code below), which executes the actual peptide docking. The peptide docking is pretty fast: it takes about 2 minutes per model.

'''Code of submit_fpdock.sh:'''

#!/bin/csh

#$1 is the start file (start.pdb)
#$2 is the number of jobs to send (each will attempt nstruct=1)

foreach i (`seq 1 $2`)
qsub -l arch=lx24-amd64 -q all.q -cwd -e ./error -o ./out -v arg1=$i,arg2=$1 ./single_fpdock.sh
end

'''Code of single_fpdock.sh:'''

#!/bin/csh

#arg1 is the run prefix number
#arg2 is the start.pdb (and native.pdb)

set arch = `uname -p`
if ( $arch == 'x86_64') then
~londonir/rosetta/rosetta_source/bin/FlexPepDocking.static.linuxgccrelease -s $arg2 -database /raid1/people/londonir/rosetta/rosetta_database -native $arg2 -pep_refine -ex1 -ex2aro -use_input_sc -nstruct 1 -unboundrot $arg2 -out:prefix $arg1'.' > log.$arg1
endif

==5. Analyze output models==

[[File:Pose1 input output vs xray.png|thumb|right|alt=Comparison of native peptide conformation (slate) to best scoring output model (orange). The backbone of the starting peptide is depicted in yellow. |Comparison of native peptide conformation (slate) to best scoring output model (orange). The backbone of the starting peptide is depicted in yellow.]]

FlexPepDock outputs a pdb, a score file (.sc) and a log file (log) for each model requested. The score files is a tab-delimited file with many scores and other interesting numbers like Interface buried surface area (I_bsa), number of hydrogen bonds (I_hb), score for the Interface (I_sc) and a score for the peptide (pep_sc). Your main interest is most likely the total score (total_score).

You can print the total scores of each model using this script:

#!/bin/csh
# prints total scores for each model

foreach file (*.sc)
echo $file":" `awk 'NR==3 {print $2}' $file`
end

In my prospective peptide docking, the best scoring output model starting from the orange pose had a total score of -324, which was better than the score of the green starting pose (-316).

==6. Calculate RMSD values compared to native structure==

FlexPepDock automatically outputs RMSD values of the models compared to the input structure (rmsALL, rmsBB etc.). But if you know the native structure, you are probably interested how the models compare to the native structure.

You can calculate RMSD values for each model compared to the native structure using the script rescore_batch.csh.
First, you have to generate a list of pdb IDs for which you want to generate RMSDs compared to the native structure:

ls -1 *.pdb > list_of_pdbs.txt

Then run rescore_batch.csh (see code below):
./rescore_batch.csh

In my prospective prediction, the orange output model had a backbone RMSD of 1.9 Angstroem (3.4 A over HA all atoms), which is a freaking good prediction. The the backbones align perfectly over the last four amino acids, and the two arginines at the N-terminus are actually (and surprisingly) not perfectly ordered in the crystal structure. 
As a side note: a few models had RMSD values as low as 1.2 A, but scored worse.

'''Code of rescore_batch.csh:'''

#!/bin/csh

#-s starting
#-native: reference for rmsd calc
# -l list file

~londonir/rosetta/rosetta_source/bin/FlexPepDocking.static.linuxgccrelease -l list_of_pdbs.txt -database /raid1/people/londonir/rosetta/rosetta_database -native 4GRV_AB_ATOM.pdb -out:prefix rescore_list'.' -flexpep_score_only -out:pdb false > log.rescore_list

#create list with pdbs
# ls -1 *.pdb > list_of_pdbs.txt

This script takes a few minutes to run and outputs a file called log.rescore_list. For each model, you can print the scores, backbone RMSDs and model name, sorted by backbone RMSD, by typing this command:

awk '{print $2, $41, $53}' rescore_list.score.sc | sort -rgk 2

=Acknowlegements=
Many thanks to Nir for compiling Rosetta on our server, the introduction into FlexPepDock and the fp* scripts.

FlexPepDock

2013-05-07T23:46:35Z

Oliv Eidam:

=FlexPepDock: Wanna dock a peptide? Use FlexPepDock!=

[[File:Rec w arg327.png|thumb|right|alt=X-ray structure of Neurotensin receptor. The C-terminus of NT(8-13) was predicted to interact with Arg327, which is depicted with spheres. |X-ray structure of Neurotensin receptor. The C-terminus of NT(8-13) was predicted to interact with Arg327, which is depicted with spheres.]]

Interested in peptide docking? Use FlexPepDock, which is a peptide docking software implemented in Rosetta. The easiest way to do that is to use the online server: 
http://flexpepdock.furmanlab.cs.huji.ac.il

If you want to run FlexPepDock locally, follow the steps below. 
In my example I docked the hexapeptide NT(8-13) with the sequence RRPYIL into the neurotensin receptor. I did this prospectively, without having looked at the bound peptide conformation before doing this. All I knew was that the C-terminal carboxylate interacts with Arg327 ([http://www.jbc.org/content/275/1/328 Barroso et al, JBC, 2000]) and that the peptide adopts an extended conformation ([http://www.pnas.org/content/100/19/10706 Luca et al, PNAS, 2003]).

==1. Create peptide coordinates in Pymol==

[[File:Rec w input peptides.png|thumb|right|alt=Backbone representation of manually docked starting peptides (orange and green). |Backbone representation of manually docked starting peptides (orange and green).]]

You can generate peptide coordinates using Pymol (Build->Residue->Alanine). If you want an extended peptide, it is best to generate coordinates for poly-Ala first, and then mutate to your peptide using the mutation wizard (Wizard->Mutagenesis->Mutate to Arg->Apply). 
Save molecule.
Do NOT add a amine at the N-terminus and do NOT add a C-terminal carboxylate to the pymol coordinates!
'''Important note:''' FlexPepDock does not change your Psi angles, so make sure that a proline is in trans if that's what you want!!

==2. Dock peptide coordinates in Pymol to generate a starting model for FlexPepDock ==

You need a rough starting peptide model for FlexPepDock. The peptide starting model should be within 5 Angstroem RMSD to the native structure. You can manually dock your peptide created above in Pymol into to the receptor by switching to Editing Mode, and dragging/rotating the peptide by holding the Shift-Button and Middle and Left-Mouse button, respectively. Don't worry about clashes of the peptide with the protein: the most important thing is that the backbone is within 5 Anstroem of the native structure.

Save the docked peptide coordinates and '''add them at the below of the PDB coordinates''' of the apo receptor structure. It is important to add them below the receptor coordinates, separated by a TER statement. Also, there should be no END statements.

==3. Run prepack.sh==

Run prepack.sh (see code below):

./prepack.sh NTS1_rrpyil_input.pdb

prepack.sh generates pNTS1_rrpyil_input_0001.pdb as output: the protein is protonated and a N-terminal amine and a C-terminal carboxylate have been added to the peptide.

'''Code of prepack.sh:'''

#!/bin/csh

#$1 is the start.pdb

set arch = `uname -p`
if ( $arch == 'x86_64') then
~londonir/rosetta/rosetta_source/bin/FlexPepDocking.static.linuxgccrelease -s $1 -database /raid1/people/londonir/rosetta/rosetta_database -native $1 -flexpep_prepack -ex1 -ex2aro -unboundrot $1 > log.prepack
else
echo "I can only run on a x86_64 system..."
endif

==4. Run FlexPepDock on Cluster using submit_fpdock.sh==

[[File:Rec w output peptides.png|thumb|right|alt=Backbone representation of best scoring output peptides (orange and green). |Backbone representation of best scoring output peptides (orange and green).]]

Run FlexPepDock on Cluster using submit_fpdock.sh (see code below):

./submit_fpdock.sh NTS1_rrpyil_input_0001.pdb 200

submit_fpdock.sh takes two arguments: the protonated protein-peptide input model from prepack.sh (NTS1_rrpyil_input_0001.pdb) and the second argument (200) is the number of models you want to generate - 200 is a good number. submit_fpdock.sh calls a script called single_fpdock.sh (see code below), which executes the actual peptide docking. The peptide docking is pretty fast: it takes about 2 minutes per model.

'''Code of submit_fpdock.sh:'''

#!/bin/csh

#$1 is the start file (start.pdb)
#$2 is the number of jobs to send (each will attempt nstruct=1)

foreach i (`seq 1 $2`)
qsub -l arch=lx24-amd64 -q all.q -cwd -e ./error -o ./out -v arg1=$i,arg2=$1 ./single_fpdock.sh
end

'''Code of single_fpdock.sh:'''

#!/bin/csh

#arg1 is the run prefix number
#arg2 is the start.pdb (and native.pdb)

set arch = `uname -p`
if ( $arch == 'x86_64') then
~londonir/rosetta/rosetta_source/bin/FlexPepDocking.static.linuxgccrelease -s $arg2 -database /raid1/people/londonir/rosetta/rosetta_database -native $arg2 -pep_refine -ex1 -ex2aro -use_input_sc -nstruct 1 -unboundrot $arg2 -out:prefix $arg1'.' > log.$arg1
endif

==5. Analyze output models==

[[File:Pose1 input output vs xray.png|thumb|right|alt=Comparison of native peptide conformation (slate) to best scoring output model (orange). The backbone of the starting peptide is depicted in yellow. |Comparison of native peptide conformation (slate) to best scoring output model (orange). The backbone of the starting peptide is depicted in yellow.]]

FlexPepDock outputs a pdb, a score file (.sc) and a log file (log) for each model requested. The score files is a tab-delimited file with many scores and other interesting numbers like Interface buried surface area (I_bsa), number of hydrogen bonds (I_hb), score for the Interface (I_sc) and a score for the peptide (pep_sc). Your main interest is most likely the total score (total_score).

You can print the total scores of each model using this script:

#!/bin/csh
# prints total scores for each model

foreach file (*.sc)
echo $file":" `awk 'NR==3 {print $2}' $file`
end

In my prospective peptide docking, the best scoring output model starting from the orange pose had a total score of -324, which was better than the score of the green starting pose (-316).

==6. Calculate RMSD values compared to native structure==

FlexPepDock automatically outputs RMSD values of the models compared to the input structure (rmsALL, rmsBB etc.). But if you know the native structure, you are probably interested how the models compare to the native structure.

You can calculate RMSD values for each model compared to the native structure using the script rescore_batch.csh.
First, you have to generate a list of pdb IDs for which you want to generate RMSDs compared to the native structure:

ls -1 *.pdb > list_of_pdbs.txt

Then run rescore_batch.csh (see code below):
./rescore_batch.csh

In my prospective prediction, the orange output model had a backbone RMSD of 1.9 Angstroem (3.4 A over HA all atoms), which is a freaking good prediction. The the backbones align perfectly over the last four amino acids, and the two arginines at the N-terminus are actually (and surprisingly) not perfectly ordered in the crystal structure. 
As a side note: a few models had RMSD values as low as 1.2 A, but scored worse.

'''Code of rescore_batch.csh:'''

#!/bin/csh

#-s starting
#-native: reference for rmsd calc
# -l list file

~londonir/rosetta/rosetta_source/bin/FlexPepDocking.static.linuxgccrelease -l list_of_pdbs.txt -database /raid1/people/londonir/rosetta/rosetta_database -native 4GRV_AB_ATOM.pdb -out:prefix rescore_list'.' -flexpep_score_only -out:pdb false > log.rescore_list

#create list with pdbs
# ls -1 *.pdb > list_of_pdbs.txt

This script takes a few minutes to run and outputs a file called log.rescore_list. For each model, you can print the scores, backbone RMSDs and model name, sorted by backbone RMSD, by typing this command:

awk '{print $2, $41, $53}' rescore_list.score.sc | sort -rgk 2

=Acknowlegements=
Many thanks to Nir for the introduction into FlexPepDock.

FlexPepDock

2013-05-07T23:41:51Z

Oliv Eidam:

=FlexPepDock: Wanna dock a peptide? Use FlexPepDock!=

[[File:Rec w arg327.png|thumb|right|alt=X-ray structure of Neurotensin receptor. The C-terminus of NT(8-13) was predicted to interact with Arg327, which is depicted with spheres. |X-ray structure of Neurotensin receptor. The C-terminus of NT(8-13) was predicted to interact with Arg327, which is depicted with spheres.]]

Interested in peptide docking? Use FlexPepDock, which is a peptide docking software implemented in Rosetta. The easiest way to do that is to use the online server: 
http://flexpepdock.furmanlab.cs.huji.ac.il

If you want to run FlexPepDock locally, follow the steps below. 
In my example I docked the hexapeptide NT(8-13) with the sequence RRPYIL into the neurotensin receptor. I did this prospectively, without having looked at the bound peptide conformation before doing this. All I knew was that the C-terminal carboxylate interacts with Arg327 ([http://www.jbc.org/content/275/1/328 Barroso et al, JBC, 2000]) and that the peptide adopts an extended conformation ([http://www.pnas.org/content/100/19/10706 Luca et al, PNAS, 2003]).

==1. Create peptide coordinates in Pymol==

[[File:Rec w input peptides.png|thumb|right|alt=Backbone representation of manually docked starting peptides (orange and green). |Backbone representation of manually docked starting peptides (orange and green).]]

You can generate peptide coordinates using Pymol (Build->Residue->Alanine). If you want an extended peptide, it is best to generate coordinates for poly-Ala first, and then mutate to your peptide using the mutation wizard (Wizard->Mutagenesis->Mutate to Arg->Apply). 
Save molecule.
Do NOT add a amine at the N-terminus and do NOT add a C-terminal carboxylate to the pymol coordinates!
'''Important note:''' FlexPepDock does not change your Psi angles, so make sure that a proline is in trans if that's what you want!!

==2. Dock peptide coordinates in Pymol to generate a starting model for FlexPepDock ==

You need a rough starting peptide model for FlexPepDock. The peptide starting model should be within 5 Angstroem RMSD to the native structure. You can manually dock your peptide created above in Pymol into to the receptor by switching to Editing Mode, and dragging/rotating the peptide by holding the Shift-Button and Middle and Left-Mouse button, respectively. Don't worry about clashes of the peptide with the protein: the most important thing is that the backbone is within 5 Anstroem of the native structure.

Save the docked peptide coordinates and '''add them at the below of the PDB coordinates''' of the apo receptor structure. It is important to add them below the receptor coordinates, separated by a TER statement. Also, there should be no END statements.

==3. Run prepack.sh==

Run prepack.sh (see code below):

./prepack.sh NTS1_rrpyil_input.pdb

prepack.sh generates pNTS1_rrpyil_input_0001.pdb as output: the protein is protonated and a N-terminal amine and a C-terminal carboxylate have been added to the peptide.

'''Code of prepack.sh:'''

#!/bin/csh

#$1 is the start.pdb

set arch = `uname -p`
if ( $arch == 'x86_64') then
~londonir/rosetta/rosetta_source/bin/FlexPepDocking.static.linuxgccrelease -s $1 -database /raid1/people/londonir/rosetta/rosetta_database -native $1 -flexpep_prepack -ex1 -ex2aro -unboundrot $1 > log.prepack
else
echo "I can only run on a x86_64 system..."
endif

==4. Run FlexPepDock on Cluster using submit_fpdock.sh==

[[File:Rec w output peptides.png|thumb|right|alt=Backbone representation of best scoring output peptides (orange and green). |Backbone representation of best scoring output peptides (orange and green).]]

Run FlexPepDock on Cluster using submit_fpdock.sh (see code below):

./submit_fpdock.sh NTS1_rrpyil_input_0001.pdb 200

submit_fpdock.sh takes two arguments: the protonated protein-peptide input model from prepack.sh (NTS1_rrpyil_input_0001.pdb) and the second argument (200) is the number of models you want to generate - 200 is a good number. submit_fpdock.sh calls a script called single_fpdock.sh (see code below), which executes the actual peptide docking. The peptide docking is pretty fast: it takes about 2 minutes per model.

'''Code of submit_fpdock.sh:'''

#!/bin/csh

#$1 is the start file (start.pdb)
#$2 is the number of jobs to send (each will attempt nstruct=1)

foreach i (`seq 1 $2`)
qsub -l arch=lx24-amd64 -q all.q -cwd -e ./error -o ./out -v arg1=$i,arg2=$1 ./single_fpdock.sh
end

'''Code of single_fpdock.sh:'''

#!/bin/csh

#arg1 is the run prefix number
#arg2 is the start.pdb (and native.pdb)

set arch = `uname -p`
if ( $arch == 'x86_64') then
~londonir/rosetta/rosetta_source/bin/FlexPepDocking.static.linuxgccrelease -s $arg2 -database /raid1/people/londonir/rosetta/rosetta_database -native $arg2 -pep_refine -ex1 -ex2aro -use_input_sc -nstruct 1 -unboundrot $arg2 -out:prefix $arg1'.' > log.$arg1
endif

==5. Analyze output models==

[[File:Pose1 input output vs xray.png|thumb|right|alt=Comparison of native peptide conformation (slate) to best scoring output model (orange). The backbone of the starting peptide is depicted in yellow. |Comparison of native peptide conformation (slate) to best scoring output model (orange). The backbone of the starting peptide is depicted in yellow.]]

FlexPepDock outputs a pdb, a score file (.sc) and a log file (log) for each model requested. The score files is a tab-delimited file with many scores and other interesting numbers like Interface buried surface area (I_bsa), number of hydrogen bonds (I_hb), score for the Interface (I_sc) and a score for the peptide (pep_sc). Your main interest is most likely the total score (total_score).

You can print the total scores of each model using this script:

#!/bin/csh
# prints total scores for each model

foreach file (*.sc)
echo $file":" `awk 'NR==3 {print $2}' $file`
end

In my prospective peptide docking, the best scoring output model starting from the orange pose had a total score of -324, which was better than the score of the green starting pose (-316).

==6. Calculate RMSD values compared to native structure==

FlexPepDock automatically outputs RMSD values of the models compared to the input structure (rmsALL, rmsBB etc.). But if you know the native structure, you are probably interested how the models compare to the native structure.

You can calculate RMSD values for each model compared to the native structure using the script rescore_batch.csh.
First, you have to generate a list of pdb IDs for which you want to generate RMSDs compared to the native structure:

ls -1 *.pdb > list_of_pdbs.txt

Then run rescore_batch.csh (see code below):
./rescore_batch.csh

The orange output model had a backbone RMSD of 1.9 Angstroem (3.4 A over HA all atoms), which is a freaking good prediction. A few models had RMSD values as low as 1.2 A, but scored worse.

'''Code of rescore_batch.csh:'''

#!/bin/csh

#-s starting
#-native: reference for rmsd calc
# -l list file

~londonir/rosetta/rosetta_source/bin/FlexPepDocking.static.linuxgccrelease -l list_of_pdbs.txt -database /raid1/people/londonir/rosetta/rosetta_database -native 4GRV_AB_ATOM.pdb -out:prefix rescore_list'.' -flexpep_score_only -out:pdb false > log.rescore_list

#create list with pdbs
# ls -1 *.pdb > list_of_pdbs.txt

This script takes a few minutes to run and outputs a file called log.rescore_list. For each model, you can print the scores, backbone RMSDs and model name, sorted by backbone RMSD, by typing this command:

awk '{print $2, $41, $53}' rescore_list.score.sc | sort -rgk 2

=Acknowlegements=
Many thanks to Nir for the introduction into FlexPepDock.

File:Pose1 input output vs xray.png

2013-05-07T23:23:39Z

Oliv Eidam:

File:Rec w output peptides.png

2013-05-07T23:22:54Z

Oliv Eidam:

File:Rec w input peptides.png

2013-05-07T23:22:27Z

Oliv Eidam:

2013-05-07T18:16:58Z

Oliv Eidam:

FlexPepDock

2013-05-07T00:43:28Z

Oliv Eidam:

FlexPepDock

2013-05-07T00:25:54Z

Oliv Eidam:

FlexPepDock

2013-05-07T00:14:38Z

Oliv Eidam: Created page with "=FlexPepDock: Wanna dock a peptide? Use FlexPepDock!= Interested in peptide docking? Use FlexPepDock, which is a peptide docking software implemented in Rosetta. The easiest ..."

Calculate DOCK6 RMSD

2013-05-06T23:37:23Z

Oliv Eidam:

==Calculate ligand RMSDs using DOCK6==

Calculating ligand RMSDs is a tricky thing, but DOCK6 does this quite well.
All you need is two ligands in mol2 format (with the same number of atoms and same Sybil atom types) and a DOCK6 binary.

Then you run the following command:

~eidamo/work/scripts/calc_dock6_rmsd.csh xtal.mol2 dock.mol2

calc_dock6_rmsd.csh (see code below) outputs three lines:

HA_RMSDs: 4.2878 
HA_RMSDh: 2.7406 
HA_RMSDm: 0.9625

The most meaningful line is the one in the middle (HA_RMSDh), which outputs the RMSD using the Hungarian Algorithm (for more information: http://dock.compbio.ucsf.edu/DOCK_6/dock6_manual.htm)

===calc_dock6_rmsd.csh===

#!/bin/csh -f
if ($#argv != 2) then
echo "Usage: $0 mol.mol2 ref.mol2"
echo "calculates rmsd between 2 molecules using hungarian algorithm implemented in DOCK 6"
endif

setenv dock66_dir "/raid1/people/tbalius/zzz.programs/dock6_2012-10-09.stonybrook"

echo "ligand_atom_file $1" > temp.in
echo "limit_max_ligands no" >> temp.in
echo "skip_molecule no" >> temp.in
echo "read_mol_solvation no" >> temp.in
echo "calculate_rmsd yes" >> temp.in
echo "use_rmsd_reference_mol yes" >> temp.in
echo "rmsd_reference_filename $2" >> temp.in
echo "use_database_filter no" >> temp.in
echo "orient_ligand no" >> temp.in
echo "use_internal_energy no" >> temp.in
echo "flexible_ligand no" >> temp.in
echo "bump_filter no" >> temp.in
echo "score_molecules no" >> temp.in
echo "atom_model all" >> temp.in
echo "vdw_defn_file $dock66_dir/parameters/vdw_AMBER_parm99.defn" >> temp.in
echo "flex_defn_file $dock66_dir/parameters/flex.defn" >> temp.in
echo "flex_drive_file $dock66_dir/parameters/flex_drive.tbl" >> temp.in
echo "ligand_outfile_prefix output" >> temp.in
echo "write_orientations no" >> temp.in
echo "num_scored_conformers 1" >> temp.in
echo "rank_ligands no" >> temp.in

$dock66_dir/bin/dock6 -i temp.in > dock6.log

grep "HA_RMSD" output_scored.mol2

==Acknowledgments==
Thanks to Trent for pointing out this great new feature in DOCK 6.6.

Phenix

2013-05-06T23:34:17Z

Oliv Eidam:

==Calculate electron densities using Phenix==

Here are a couple of useful phenix commands to calculate electron densities. In docking, it is often useful to have a look at electron densities, especially for proteins and ligands determined at low resolution.
Also, you may want to display them in Pymol: here[link] is described how you generate a electron density map in CCP4 format map, which is required for Pymol.

===1. Download PDB coordinates and structure factors from PDB.ORG===

[[File:Pdb sf download.jpg|thumb|right|alt=Where to download files from PDB entry.|Where to download files from PDB entry.]] To calculate a electron density, you need to download PDB coordinates and structure factors (Miller indices and measured intensities) from the PDB. On pdb.org, there is a dropdown menu on the right top corner of each entry site where PDB and structure factors can be downloaded (see image on the right).

===2. Calculate electron densities from deposited structure factors===

You can calculate an electron density (CCP4 format for Pymol and mtz format for Coot) using a PDB and structure factors using Phenix. 
If you want a 2Fo-Fc map, simply type:

phenix.maps 4DJH.pdb 4djh-sf.cif

If you want a Fo-Fc difference density map from deposited structure factors, type:

phenix.maps 4DJH.pdb 4djh-sf.cif maps.map.map_type=Fo-F

===3. Calculate a simulated annealing (SA) omit map===

Sometimes it is interesting to see how much density of a ligand is left upon removal of the ligand (or protein residues or water molecules) from your model. 
You can calculate a simulated annealing (SA) omit map using Phenix like this:

# remove ligands and water molecules
grep -v HETATM 4DJH.pdb > 4DJH_no_HETATM.pdb

# Run Phenix refinement including simulated annealing
phenix.refine 4djh-sf.cif 4DJH_no_HETATM.pdb simulated_annealing=true

# Generate a CCP4 map for Pymol:
phenix.maps maps.map.map_type=Fo-Fc 4DJH_no_HETATM.pdb 4DJH_no_HETATM_refine_001.mtz

===4. Display electron densities in Pymol===

Pymol wants CCP4 maps as input. To generate those, we use Phenix as described above. 

[[File:Sa omit map2.png|thumb|upright=1.35|right|alt=SA omit map around kappa ligand in green, 2Fo-Fc map in cyan. The ligand and water molecules show up nicely, although they were omitted during refinement. |SA omit map around kappa ligand in green, 2Fo-Fc map in cyan. The ligand and water molecules show up nicely, although they were omitted during refinement.]]
Here is a Pymol script to show density around the ligand and it's surrounding residues:

load 4DJH.pdb, kappa_opiod_rec
load 4DJH_2mFo-DFc_map.ccp4
load 4DJH_no_HETATM_Fo-Fc_map.ccp4, SA_omit

select lig, resn JDC AND chain A
cmd.show("sticks" ,"lig")
select around_lig, lig expand 10

isomesh fo-fc_map, SA_omit, 3.0, around_lig
isomesh 2fo-fc_map, 4DJH_2mFo-DFc_map, 1.0, around_lig
cmd.color(5,"2fo-fc_map")
cmd.color(3,"fo-fc_map")

orient lig

set_view (\
-0.552325428, -0.811329126, -0.191545352,\
0.833266199, -0.543944836, -0.098800458,\
-0.024021570, -0.214180335, 0.976500750,\
0.001199273, -0.000215724, -52.703464508,\
4.263263702, -21.972831726, 57.798912048,\
48.832965851, 56.566112518, -20.000000000 )

set ray_shadow, off

==Acknowledgments==
Thanks to the Phenix and Pymol developers.

Screen3d

2013-05-06T21:48:02Z

Oliv Eidam:

Omega.parm

2013-05-06T21:44:37Z

Oliv Eidam: