Performing a hydrogen-bond analysis
An analysis of the hydrogen-bonding network can provide useful structural information in many situations. For example, it can be used to classify minima (and hence colour disconnectivity graphs) or analyse which contacts are important in protein/ligand binding and identify possible mutations to make. With AMBER, this is done using the ptraj tool. You can find a full listing of ptraj functionality in the AMBER tools manual AMBER Tools here.
Running a hydrogen-bond analysis is broken into three stages, preparing the ptraj input, running the analysis and analysing the output. This tutorial will cover the first two stages, as analysing the output is highly system and purpose specific!
Preparing the ptraj input
A ptraj input file for hydrogen bond analysis is broken into three parts:
- reading in the structure to be analysed
- defining donor and acceptor atoms
- running the hbond analysis
We are going to write a script which will then be passed to 'ptraj' to run the analysis. In this example, we will call it hbond.in.
Reading in the structure to be analysed
ptraj is able to read in both PDB and AMBER .rst format files. Here we will use PDB files. To read a PDB file into ptraj, we simply use the trajin command:
#-- READING IN A SINGLE FILE trajin A_Netherlands_602_2009_DG222_0_SIA26GAL14NAG_OME_resp_md1.pdb
When doing a hydrogen-bond analysis, you need to decide whether you are interested in the bonding within a single structure (for classifying minima for example), or average properties of many minima (looking at hydrogen-bond occupation). You can read as many PDB files in as you'd like, just by adding additional trajin lines. For example:
#-- READING IN 10 FILES trajin A_Netherlands_602_2009_DG222_0_SIA26GAL14NAG_OME_resp_md1.pdb trajin A_Netherlands_602_2009_DG222_0_SIA26GAL14NAG_OME_resp_md2.pdb trajin A_Netherlands_602_2009_DG222_0_SIA26GAL14NAG_OME_resp_md3.pdb trajin A_Netherlands_602_2009_DG222_0_SIA26GAL14NAG_OME_resp_md4.pdb trajin A_Netherlands_602_2009_DG222_0_SIA26GAL14NAG_OME_resp_md5.pdb trajin A_Netherlands_602_2009_DG222_0_SIA26GAL14NAG_OME_resp_md6.pdb trajin A_Netherlands_602_2009_DG222_0_SIA26GAL14NAG_OME_resp_md7.pdb trajin A_Netherlands_602_2009_DG222_0_SIA26GAL14NAG_OME_resp_md8.pdb trajin A_Netherlands_602_2009_DG222_0_SIA26GAL14NAG_OME_resp_md9.pdb trajin A_Netherlands_602_2009_DG222_0_SIA26GAL14NAG_OME_resp_md10.pdb
Defining donor and acceptor atoms
You need to tell ptraj which atoms to consider as donors and acceptors. This is done using the mask syntax that is common in AMBER. If you're working on a protein, comprised entirely of standard amino acids, with no custom ligand involved, you simply need to define the standard donor-acceptor pairs for the sidechains and backbone:
#-- DONORS from standard amino acid sidechains donor mask :GLN@OE1 donor mask :GLN@NE2 donor mask :ASN@OD1 donor mask :ASN@ND2 donor mask :TYR@OH donor mask :ASP@OD1 donor mask :ASP@OD2 donor mask :GLU@OE1 donor mask :GLU@OE2 donor mask :SER@OG donor mask :THR@OG1 donor mask :HIS@ND1 donor mask :HIE@ND1 donor mask :HID@NE2 #-- ACCEPTORS from standard amino acid sidechains acceptor mask :ASN@ND2 :ASN@HD21 acceptor mask :ASN@ND2 :ASN@HD22 acceptor mask :TYR@OH :TYR@HH acceptor mask :GLN@NE2 :GLN@HE21 acceptor mask :GLN@NE2 :GLN@HE22 acceptor mask :TRP@NE1 :TRP@HE1 acceptor mask :LYS@NZ :LYS@HZ1 acceptor mask :LYS@NZ :LYS@HZ2 acceptor mask :LYS@NZ :LYS@HZ3 acceptor mask :SER@OG :SER@HG acceptor mask :THR@OG1 :THR@HG1 acceptor mask :ARG@NH2 :ARG@HH21 acceptor mask :ARG@NH2 :ARG@HH22 acceptor mask :ARG@NH1 :ARG@HH11 acceptor mask :ARG@NH1 :ARG@HH12 acceptor mask :ARG@NE :ARG@HE acceptor mask :HIS@NE2 :HIS@HE2 acceptor mask :HIE@NE2 :HIE@HE2 acceptor mask :HID@ND1 :HID@HD1 acceptor mask :HIP@ND1,NE2 :HIP@HE2,HD1 #-- Backbone donors and acceptors # N-H for prolines do not exist so need to exclude :PRO@N from the acceptor mask donor mask '@O' acceptor mask '!:PRO & @N' '@H' #Terminal residues have different atom names - OXT is only present in termini so can use always donor mask @OXT #This assumes there is only one protein segment i.e. if there was a terminal BEFORE the ligand #you would need to add more residue numbers e.g. 'acceptor mask :1,20@N :1,20@H1' etc acceptor mask :1@N :1@H1 acceptor mask :1@N :1@H2 acceptor mask :1@N :1@H3
DONOR ACCEPTORH ACCEPTOR atom# :res@atom atom# :res@atom atom# :res@atom %occupied distance angle | 3969 :255@O1A | 1481 :98@HG1 1480 :98@OG1 | 100.00 2.629 ( 0.10) 16.33 ( 4.50)
From the AMBER Tools manual (AMBER 11 version)
hbond [<dataset name>] [out <filename>] <mask> [angle <cut>] [dist <cut>] [avgout <avgfilename>] Search for hydrogen bond donor and acceptor atoms in the region specified by <mask> (currently following the simplistic criterion that “hydrogen bonds are FON”, i.e., hydrogens bonded to F, O, and N atoms are considered) and calculate the number of hydrogen bonds formed for each frame, writing the results to the file specified by “out”. Hydrogen bonds are considered to have the form: Acceptor ... Hydrogen-Donor and are determined via the distance between the heavy atoms using a cutoff of 3.0 Å (or the value specified by “dist”) and the angle between the acceptor, hydrogen, and donor atoms using a cutoff of 135° (or the value specified by “angle”)