|
|
| (8 intermediate revisions by the same user not shown) |
| Line 5: |
Line 5: |
| | | | |
| | ==I. Introduction== | | ==I. Introduction== |
| | + | |
| | + | Yaping |
| | | | |
| | ==II. Structural Preparation== | | ==II. Structural Preparation== |
| Line 10: |
Line 12: |
| | | | |
| | ====Antechamber, Parmchk, tLeap==== | | ====Antechamber, Parmchk, tLeap==== |
| | + | |
| | + | Omar, Katie |
| | | | |
| | ==III. Simulation using pmemd== | | ==III. Simulation using pmemd== |
| Line 15: |
Line 19: |
| | ====PMEMD==== | | ====PMEMD==== |
| | | | |
| − | =
| + | Agatha, Beilei |
| | | | |
| | ==IV. Simulation Analysis== | | ==IV. Simulation Analysis== |
| Line 21: |
Line 25: |
| | ===Ptraj=== | | ===Ptraj=== |
| | | | |
| − | You should create another work directory for this step (004.ptraj, for example).
| + | Lauren, Haoyue |
| − | '''1.''' At first we want to concatenate the two 1ns trajectories together, stripping off the waters, and creating a .strip-file as an output. Below is the input file which will allow us to do so.
| |
| − | | |
| − | '''ptraj.1.in'''
| |
| − | trajin ../003.pmemd/10md.trj 1 1000 1
| |
| − | trajin ../003.pmemd/11md.trj 1 1000 1
| |
| − | trajout 4TKG.trj.strip nobox
| |
| − | strip :WAT
| |
| − | | |
| − | The two sets of numbers ''1 1000 1'' give the input information about which frames are used for the Ptraj. The first two numbers ''1'' and ''1000'' specify the starting and ending snapshots from the trajectory file. The ending number of the snapshot doesn't need to be accurate because if you actually don't have enough snapshots in your trajectory file, Ptraj will read up to the last one you have. The last number ''1'' specifies the frequency of the snapshot saved, in this case, we are saving every frame of the trajectory file. And the last line of the input file will take away all the water molecules.
| |
| − | | |
| − | As the input file is prepared we can launch the first analysis as follows:
| |
| − | ptraj ../002.tleap/4TKG.com.wat.leap.parm ptraj.1.in > ptraj.1.log
| |
| − | | |
| − | As the output you will obtain '''4TKG.trj.strip''' file which contains ''2000'' snapshots of the trajectory.
| |
| − | | |
| − | '''2.''' Later on we want to compare the output file just obtained to our reference file '''4TKG.com.gas.leap.rst7''', using the following input file.
| |
| − | | |
| − | '''ptraj.2.in'''
| |
| − | trajin 4TKG.trj.strip 1 2000 1
| |
| − | trajout 4TKG.com.trj.stripfit
| |
| − | reference ../002.tleap/4TKG.com.gas.leap.rst7
| |
| − | rms reference out 4TKG.rmsd.CA.txt :1-206@CA
| |
| − | | |
| − | Since we have just concatenated the two trajectories, we will have ''2000'' snapshots in '''4TKG.trj.strip'''. The third line in the input specifies the reference file, we have taken away all the water molecules during the first step, hence we are working here with the gas phase complex. The last line says we are calculating the rmsd for alpha carbon number ''1'' to ''206''.
| |
| − | | |
| − | And as we have this file filled out, we can run this step:
| |
| − | ptraj ../002.tleap/4TKG.com.gas.leap.parm ptraj.2.in > ptraj.2.log
| |
| − | | |
| − | The output file '''4TKG.rmsd.CA.txt''' wil contain two columns, the first one is the number of frame and the second one stands for rmsd value.
| |
| − | | |
| − | '''3.''' Afterwards we will generate a similar file for ligand, using the following input file.
| |
| − | | |
| − | '''ptraj.3.in'''
| |
| − | trajin 4TKG.com.trj.stripfit 1 2000 1
| |
| − | reference ../002.tleap/4TKG.com.gas.leap.rst7
| |
| − | rms reference out 4TKG.lig.rmsd.txt :207@C*,N*,O*,S* nofit
| |
| − | | |
| − | And then run this step:
| |
| − | ptraj ../002.tleap/4TKG.com.gas.leap.parm ptraj.3.in > ptraj.3.log
| |
| − | | |
| − | By doing this we will compare the trajectory file to '''4TKG.com.gas.leap.rst7''' as well, but working with the ligand instead of the receptor (we specified that by pointing that we want to calculate the rmsd for carbon, nitrogen, oxygen and sulfur in residue ''207''.
| |
| − | | |
| − | The last two steps are to obtain energetic information about the system. To do this we take a trajectory file of the gas phase complex '''4TKG.com.trj.stripfit''', and want to create two more trajectory files containing the information on only receptor and only ligand correspondingly.
| |
| − | | |
| − | '''4.''' At this step we consider receptor only. The input file is provided below:
| |
| − | | |
| − | '''ptraj.4.in'''
| |
| − | trajin 4TKG.com.trj.stripfit 1 2000 1
| |
| − | trajout 4TKG.rec.trj.stripfit
| |
| − | strip :207
| |
| − | | |
| − | And then run this step:
| |
| − | ptraj ../002.tleap/4TKG.com.gas.leap.parm ptraj.4.in > ptraj.4.log
| |
| − | | |
| − | '''5.''' And the same procedure for the ligand, with the following input file:
| |
| − | | |
| − | '''ptraj.5.in'''
| |
| − | trajin 4TKG.com.trj.stripfit 1 2000 1
| |
| − | trajout 4TKG.lig.trj.stripfit
| |
| − | strip :1-206
| |
| − | | |
| − | And then run this step:
| |
| − | ptraj ../002.tleap/4TKG.com.gas.leap.parm ptraj.5.in > ptraj.5.log
| |
| − | | |
| − | '''6.''' ''(optional)'' Visualization
| |
| − | | |
| − | As we've gone through all these steps, the analysis is done. If you want to visualize the trajectories, you first need to copy the trajectory files to Herbie like this, for example (being a level above '''004.PTRAJ''' directory):
| |
| − | scp 004.ptraj/ your_username@herbie.mathlab.sunysb.edu:~AMS536/amber_tutorial/004.ptraj
| |
| − | | |
| − | Now, launch VMD, then open one of the .prm7 files in 002.tleap. If you want to visualize the whole complex in gas state, you can open 4TKG.com.gas.leap.prm7 with AMBER7 Parm from 002.tleap and then 4TKG.com.trj.stripfit with AMBER coordinates from 004.ptraj. With these files, you can look at the real-time movement of the complex in the gas state. You can repeat this procedure to observe the real-time movement of the complex in the water state. Just open 4TKG.com.wat.leap.prm7 instead of 4TKG.com.gas.leap.prm7.
| |
| | | | |
| | ===MM-GBSA Energy Calculation=== | | ===MM-GBSA Energy Calculation=== |
| − | Molecular Mechanics-Generalized Born Surface Area (MM-GBSA) is a great method to calculate or estimate relative binding affinity of a ligand(s) to a receptor. The binding energy calculated from this method are also known as free energies of binding, where the more negative values indicate stronger binding. For this section, the topology files for the ligand, receptor and complex are needed.
| |
| − |
| |
| − | Create a new directory:
| |
| − | mkdir 005.MMGBSA
| |
| − | Create an input file name
| |
| − | vim gb.rescore.in
| |
| − | Enter the following into the input file:
| |
| − | Single point GB energy calc
| |
| − | &cntrl
| |
| − | ntf = 1, ntb = 0, ntc = 2,
| |
| − | idecomp= 0,
| |
| − | igb = 5, saltcon= 0.00,
| |
| − | gbsa = 2, surften= 1.0,
| |
| − | offset = 0.09, extdiel= 78.5,
| |
| − | cut = 99999.0, nsnb = 99999,
| |
| − | imin = 5, maxcyc = 1, ncyc = 0,
| |
| − | /
| |
| − | Create a tcsh/bash/csh script (run.sander.rescore.csh) with the following information:
| |
| − | #! /bin/tcsh
| |
| − | #PBS -l nodes=1:ppn=1
| |
| − | #PBS -l walltime=48:00:00
| |
| − | #PBS -o zzz.qsub.out
| |
| − | #PBS -e zzz.qsub.err
| |
| − | #PBS -V
| |
| − | #PBS -N mmgbsa
| |
| − |
| |
| − | set workdir = /nfs/user03/kbelfon/amber_tutorial/005.mmgbsa
| |
| − |
| |
| − | cd $workdir
| |
| − |
| |
| − | sander -O -i gb.rescore.in \
| |
| − | -o gb.rescore.out.com \
| |
| − | -p ../002.tleap/4TKG.com.gas.leap.prm7 \
| |
| − | -c ../002.tleap/4TKG.com.gas.leap.rst7 \
| |
| − | -y ../004.ptraj/4TKG.com.trj.stripfit \
| |
| − | -r restrt.com \
| |
| − | -ref ../002.tleap/4TKG.com.gas.leap.rst7 \
| |
| − | -x mdcrd.com \
| |
| − | -inf mdinfo.com
| |
| − |
| |
| − | sander -O -i gb.rescore.in \
| |
| − | -o gb.rescore.out.lig \
| |
| − | -p ../002.tleap/4TKG.lig.gas.leap.prm7 \
| |
| − | -c ../002.tleap/4TKG.lig.gas.leap.rst7 \
| |
| − | -y ../004.ptraj/4TKG.lig.trj.stripfit \
| |
| − | -r restrt.lig \
| |
| − | -ref ../002.tleap/4TKG.lig.gas.leap.rst7 \
| |
| − | -x mdcrd.lig \
| |
| − | -inf mdinfo.lig
| |
| − |
| |
| − | sander -O -i gb.rescore.in \
| |
| − | -o gb.rescore.out.test.rec \
| |
| − | -p ../002.tleap/4TKG.rec.gas.leap.prm7 \
| |
| − | -c ../002.tleap/4TKG.rec.gas.leap.rst7 \
| |
| − | -y ../004.ptraj/4TKG.rec.trj.stripfit \
| |
| − | -r restrt.rec \
| |
| − | -ref ../002.tleap/4TKG.rec.gas.leap.rst7 \
| |
| − | -x mdcrd.rec \
| |
| − | -inf mdinfo.rec
| |
| − |
| |
| − | exit
| |
| − |
| |
| − | Execute this script on the seawulf cluster or machine(s) of your preference
| |
| − | qsub run.sander.rescore.csh
| |
| − | Three output files will be generated once the job is completed:
| |
| − | '''gb.rescore.out.com''', '''gb.rescore.out.lig''', and '''gb.rescore.out.rec'''
| |
| − | These files represent the single point energy calculation results for the complex (.com), the ligand (.lig) and the receptor (.rec). The energy will be output by the program Sander for each frame specified in the input file. The final results for one frame in one of the three files should look as the following:
| |
| − |
| |
| − | FINAL RESULTS
| |
| − |
| |
| − |
| |
| − |
| |
| − | NSTEP ENERGY RMS GMAX NAME NUMBER
| |
| − | 1 5.9132E+03 2.0005E+01 1.2640E+02 C 159
| |
| − |
| |
| − | BOND = 661.8980 ANGLE = 1751.7992 DIHED = 2581.7692
| |
| − | VDWAALS = -1696.6585 EEL = -13958.9335 EGB = -3125.9524
| |
| − | 1-4 VDW = 747.0185 1-4 EEL = 7750.8118 RESTRAINT = 0.0000
| |
| − | ESURF = 11201.4791
| |
| − | minimization completed, ENE= 0.59132314E+04 RMS= 0.200047E+02
| |
| − | '''Extracting Data from MM-GBSA calculation and calculating Free energy of Binding'''
| |
| − | From the output files above:
| |
| − |
| |
| − | VDWAALS = ΔGvdw
| |
| − |
| |
| − | EELS = ΔGcoul
| |
| − |
| |
| − | EGB = ΔGpolar
| |
| − |
| |
| − | SASA = ESURF
| |
| − |
| |
| − | With this information ΔGnonpolar can be solved using equation(1):
| |
| − |
| |
| − | ΔGnonpolar = SASA*0.00542 + 0.92 (1)
| |
| − |
| |
| − | Once ΔGnonpolar is solved then ΔGmmgbsa can be determined by equation(2):
| |
| − |
| |
| − | ΔGmmgbsa = ΔGvdw + ΔGcoul + ΔGpolar + ΔGnonpolar (2)
| |
| − |
| |
| − | Solve equation 2 and 3 using the extracted information from all three output files. So therefore you should have ΔGmmgbsa for the complex, ligand and receptor
| |
| − |
| |
| − | Finally ΔΔGbind can be calculated using equation (3):
| |
| − |
| |
| − | ΔΔGbind = ΔGmmgbsa,complex – (ΔGmmgbsa,lig + ΔGmmgbsa,rec) (3)
| |
| − |
| |
| − | Plot your ΔΔGbind and examine the plot for changes in the ligand position and the ΔΔGbind. Also, you should calculate the mean and standard deviation for your ΔΔGbind.
| |
| − |
| |
| − | The following script (get.mmgbsa.sh) can be used to extract the energy from the three output files obtained above and calculate ΔΔGbind:
| |
| − |
| |
| − | #! /bin/bash
| |
| − | # by Haoquan
| |
| − | echo com lig rec > namelist
| |
| − | LIST=`cat namelist`
| |
| − | for i in $LIST ; do
| |
| − | grep VDWAALS gb.rescore.out.$i | awk '{print $3}' > $i.vdw
| |
| − | grep EGB gb.rescore.out.$i | awk '{print $9}' > $i.polar
| |
| − | grep EELS gb.rescore.out.$i | awk '{print $6}' > $i.coul
| |
| − | grep ESURF gb.rescore.out.$i | awk '{print $3 * 0.00542 + 0.92}' > $i.surf
| |
| − | paste -d " " $i.vdw $i.polar $i.surf $i.coul | awk '{print $1 + $2 + $3 + $4}' > data.$i
| |
| − | rm $i.*
| |
| − | done
| |
| − | paste -d " " data.com data.lig data.rec | awk '{print $1 - $2 - $3}' > data.all
| |
| − | for ((j=1; j<=`wc -l data.all | awk '{print $1}'`; j+=1)) do
| |
| − | echo $j , >> time
| |
| − | done
| |
| − | paste -d " " time data.all > MMGBSA_vs_time.dat
| |
| − | rm namelist time data.*
| |
| − |
| |
| − | Execute this script:
| |
| − | bash get.mmgbsa.sh
| |
| − |
| |
| − | A text file called MMGBSA_vs_time.dat with x and y values separated by a space and comma should be created. Use XMGRACE to plot this dat file using the following command in Linux:
| |
| | | | |
| − | xmgrace MMGBSA_vs_time.dat
| + | Monaf |
| | | | |
| | ==V. Frequently Encountered Problems== | | ==V. Frequently Encountered Problems== |
In this tutorial, we will learn how to run a molecular dynamics simulation of a protein-ligand complex. We will then post-process that simulation by calculating structural fluctuations (with RMSD) and free energies of binding (MM-GBSA).
