Count hydrogen bond number of all resides on receptor at each time range.

The ligand and receptor always form continuous hydrogen bonds. To investigate the occupancy of hydrogen bonds, I need to find data of following purpose.

1. Find the residues on the receptor that could form hydrogen bonds with ligand

2. Check the number of hydrogen bonds between each residue on receptor the whole ligand at each time range.

I used VMD and BASH script to do this.

Assuming in VMD the receptor residue number is 1 to 123. The residue number of ligand is 124 to 134. And the simulation has 40ns. Each frame is 10ps. The standard for hydrogen bonds are distance< 3.5A and angle> 130degree.

Build the following Tcl/Tk script for VMD.

###########hbond.tcl################

package require hbonds
mol new complex.prmtop type parm7 first 0 last -1 step 1 filebonds 1 autobonds 1 waitfor all
mol addfile complex.dcd type dcd first 0 last -1 step 1 filebonds 1 autobonds 1 waitfor all
set sel1 [atomselect top "residue 0 to 123"]
set sel2 [atomselect top "residue 124 to 124"]
hbonds -polar yes -sel1 $sel1 -sel2 $sel2 -dist 3.5 -ang 50 -plot no -writefile yes -outfile complex_hbond.dat -type all -detailout complex_hbond_details.txt
hbonds -polar yes -sel1 $sel1 -sel2 $sel2 -dist 3.5 -ang 50 -plot no -writefile yes -outfile complex_10_hbond.dat -type all -detailout complex_10_hbond_details.txt -frames 0:999
hbonds -polar yes -sel1 $sel1 -sel2 $sel2 -dist 3.5 -ang 50 -plot no -writefile yes -outfile complex_20_hbond.dat -type all -detailout complex_20_hbond_details.txt -frames 1000:1999
hbonds -polar yes -sel1 $sel1 -sel2 $sel2 -dist 3.5 -ang 50 -plot no -writefile yes -outfile complex_30_hbond.dat -type all -detailout complex_30_hbond_details.txt -frames 2000:2999
hbonds -polar yes -sel1 $sel1 -sel2 $sel2 -dist 3.5 -ang 50 -plot no -writefile yes -outfile complex_40_hbond.dat -type all -detailout complex_40_hbond_details.txt -frames 3000:3999
quit

###########################################

Run vmd -e hbond.tcl

The complex_hbond_details.txt contains all the residues formed hydrogen bond in the 40ns. I will pick the residue name from it and extract the occupancy in each 10ns files with a BASH script. And the complex_??_hbond_details.txt also contains the occupancy of each residue. I could use occupancy * 1000 to get the hydrogen bond number in the time range.

##############count_hbond.sh##########

((n=`cat $1 | wc -l`)) #total line of complex_hbond_details.txt

a=`cat $1 | sed -n 3,${n}p | sed s/\-/\ /g | awk '{print $1}'` #residue in the donor
b=`cat $1 | sed -n 3,${n}p | sed s/\-/\ /g | awk '{print $3}'` #residue in the acceptor
c=(`echo $a $b | tr ' ' '12' | sort | uniq | xargs echo`) #all residue and remove the duplicated terms
d=(`echo $a $b | tr ' ' '12' | cut -c 4- | sort -n | uniq | xargs echo`) # all the residue number, rand from small to large and remove deplicated

rm -f *ns.dat
#get the residue name on the receptor
for ((i=0;i<${#d[*]};i++));do
  if (( ${d[$i]} >=124 )); then
    break
  else
    e[$i]=`echo ${c[*]} | tr ' ' '\\12' | grep ${d[$i]}`
  fi
done
e=(`echo ${e[*]} | tr ' ' '\12' | sort | uniq`)
echo ${e[*]} | tr ' ' '12' > 0ns.dat

###get the number of hydrogen bonds of each residue on receptor in each 10ns
for i in `ls *_??_*txt`;do
  for ((j=0;j<${#e[*]};j++));do
    unset f
    ((sum=0))
    f=(`grep ${e[$j]} $i | awk '{print $3}' | rev | cut -c 2- | rev`)
      for ((p=0;p<${#f[*]};p++));do
        sum=`echo "${f[$p]}+$sum" | bc`
      done
    nn=`echo $i | cut -f2 -d"_"`
    sum=`echo "$sum * 10" | bc | cut -f1 -d"."`
    echo $sum >> ${nn}ns.dat
  done
done

i=`ls -rtx1 *ns.dat`
paste $i > total.dat
rm -r *ns.dat

unset c d e

##################################################

Run ./count_hbond.sh complex_hbond_details.txt

The total .dat is the output file. The format is below

#########total.dat format############################

residue1 hbond#(0-10ns) hbond#(10-20ns) hbond#(20-30ns) hbond#(30-40ns)

residue2 hbond#(0-10ns) hbond#(10-20ns) hbond#(20-30ns) hbond#(30-40ns)

########################################

Then you can use the data to draw a 2D map to indicate the hydrogen bond change during the time and explain the motif.

Advertisements

4 thoughts on “Count hydrogen bond number of all resides on receptor at each time range.

  1. Dear Chris,
    I have calculated all the residues of my target protein formed hydrogen bond with ligand in the 5ns simulation using cpptraj command. Out put file is pasted below;
    #Acceptor DonorH Donor Frames Frac AvgDist AvgAng
    SER_119@OG THA_533@H14 THA_533@N27 965 0.0386 2.9132 151.2571
    ASP_69@OD2 THA_533@H13 THA_533@N27 510 0.0204 2.8582 154.4103
    THA_533@N11 TYR_439@HH TYR_439@OH 310 0.0124 2.8903 151.9957
    THA_533@N27 TYR_118@HH TYR_118@OH 66 0.0026 2.9068 158.8886
    TYR_118@OH THA_533@H14 THA_533@N27 32 0.0013 2.9096 147.5846
    ASN_82@OD1 THA_533@H13 THA_533@N27 26 0.0010 2.9001 149.3425
    TYR_118@OH THA_533@H13 THA_533@N27 25 0.0010 2.9179 146.9786
    SER_119@OG THA_533@H13 THA_533@N27 15 0.0006 2.8861 145.6209
    TYR_331@OH THA_533@H13 THA_533@N27 6 0.0002 2.9221 153.1109
    THA_533@N11 TRP_81@HE1 TRP_81@NE1 2 0.0001 2.9397 136.8408
    ASP_69@OD1 THA_533@H13 THA_533@N27 1 0.0000 2.9049 149.3546

    I want to draw 2D picture of these interaction, like you have described in the following link
    https://sunxiaoquan.wordpress.com/2015/03/15/find-residues-of-hydrogen-bonds-per-picosecond-at-different-time-period/

    Can I use this data for plotting this 2D map or I have to follow your instructions?

    Like

  2. Dear chris,

    Can you please share your complex_hbond_details.txt file with me so that i can understand the difference between vmd and cpptraj hydrogen bonding occupancy. If both results are same then i can make 2D plot after running cpptraj instead of vmd command???

    Like

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s