HEATMAP for intramolecular hydrogen bond analysis of a 10mer copolymer.

I am working on some simulations of certain copolymer. This copolymer could form intensive intramolecular hydrogen bonds. To observe how each residue interacts with other residues to form intramolecular hydrogen bonds, I choose HEATMAP to present the results.

hbond_intramolecular_occupancy

So the above figureĀ can show the average hydrogen bonds formed among different residue pairs on one copolymer during a period.

 

I did this in three step.

  1. Calculate intramolecular hydrogen bonds occupancy with VMD.
  2. Find pairs from the output files.
  3. Plot with GNUPLOT (version 5.0 patchlevel 1).

 

So easy. Mama will never worry about my study anymore ~~~~~~

———————-

The VMD script is below. My system only contains water and polymer, so I use “not water” to simplify my selection.

package require hbonds
mol new copolymer.prmtop type parm7 first 0 last -1 step 1 filebonds 1 autobonds 1 waitfor all
mol addfile trajectory.dcd type dcd first 1 last -1 step 1 filebonds 1 autobonds 1 waitfor all
set sel1 [atomselect top "not water"]
set sel2 [atomselect top "not water"]
hbonds -polar yes -sel1 $sel1 -sel2 $sel2 -dist 3.5 -ang 50 -plot no -writefile yes -type pair -outfile vmd_hbond_chain.dat -detailout vmd_hbond_chain_detail.dat
quit

 

———————-

The BASH to find pairs is below. This is 10mer copolymer. The names of residue in AMBER are MMH1 MMM2 MMM3 … MMM9 MME10. I chose perl to speed up the floating computing. It might save several minutes.

for i in {1..10};do echo \,${i};done | tr '\n' ' ' > hbond_intramolecular_occupancy.dat
echo >> hbond_intramolecular_occupancy.dat
a=()
for i in {1..10};do
 j=MMM${i}
 if (( $i == 1)); then j=MMH${i};fi
 if (( $i == 10)); then j=MME${i};fi
 for m in {1..10};do
  n=MMM${m}
  if (( $m == 1)); then n=MMH${m};fi
  if (( $m == 10)); then n=MME${m};fi
  a+=(`grep -E "${j}-.*${n}-" vmd_hbond_chain_detail.dat | awk '{print $3}' | cut -f1 -d"%" | perl -lne '$x += $_; END { print $x/100; }'`)
 done
 echo ${i} ${a[*]} | tr ' ' ',' >> hbond_intramolecular_occupancy.dat
 unset a
done
unset a

 

—————–

The GNUPLOT script is below. I got help from http://gnuplot.sourceforge.net/demo/heatmaps.html. And I did some slight modifications to fit my system.

reset
set grid front
set xrange [*:]
set terminal pngcairo
set encoding iso_8859_1
set xl "Acceptor Residue"
set yl "Donor Residue"
# Color runs from white to green
set palette rgbformula -7,2,-7
set cbrange [0:2]
set cblabel "Number of Hydrogen Bonds per Frame"
set cbtics

set output "hbond_intramolecular_occupancy.png"
set title "Intramolecular Hydrogen Bonds"
set datafile separator ","
plot 'hbond_intramolecular_occupancy.dat' matrix rowheaders columnheaders using 1:2:3 with image noti
set datafile separator
set output
Advertisements

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