Find residues of hydrogen bonds per picosecond at different time period.

I was required to plot this figure. The hydrogen bonds formed between residues on receptor and ligand at different time period could show the kinetic change of the binding. My method is to use VMD to get the occupancy of hydrogen bonds and find the high occupancy receptor residues of the whole simulation. Then find the hydrogen bonds of these receptor residues at each 10ns time period. Dividing the the number of hydrogen bonds with the time to get the average hydrogen bonds per picosecond. And then use GNUPLOT to plot the figure.

I have pubilshed the blog how to find the hydrogen bonds of each time period and output the table of hydrogen bonds at https://sunxiaoquan.wordpress.com/2015/03/09/count-hydrogen-bond-number-of-all-resides-on-receptor-at-each-time-range/. The final output file is total.dat.

The format of total file is

residue1 hbond#(time period1)  hbond#(time period2) hbond#(time period3) ……

residue2 ……

 

Now I need to change the data format for gnuplot as below

Time residue1 residue2…

time1 hbond#  hbond#……

time2 hbond#  hbond#……

 

I use the following script to do this and generate a gnuplot script. Run ./plot.sh total.dat. The gnuplot script is hbond.gnu.

a=(`cat $1 | awk '{print $1}'`)
for ((i=0;i<${#a[*]};i++));do
        unset b
        ((j=$i+1))
        b=(`sed -n ${j}p $1 | cut -f2-`)
        for ((p=0;p<${#b[*]};p++));do
                ((c=`echo "${b[$p]} * 10 " | bc | awk '{printf "%f", $0}' | cut -f1 -d"."`))
                if (( $c >= 3 ));then
                        d[$i]=${a[$i]}
                        break
                fi
        done
done

d=(`echo ${d[*]}`)

rm -f big_occupancy.dat
for ((i=0;i<${#d[*]};i++));do
        grep ${d[$i]} $1 | tr '\t ' '12' > ${i}_temp.dat
done

i=`ls *_temp.dat | sort -n`
paste $i > big_occupancy_temp.dat
((p=0))
echo Time > time_gap1_temp.dat
((n=`cat big_occupancy_temp.dat | wc -l`))
for i in {0..200..10};do
        ((j=$i+10))
        echo ${i}\-${j}ns >> time_gap1_temp.dat
        ((i=$i+10))
done
sed -n 1,${n}p time_gap1_temp.dat > time_gap_temp.dat
paste time_gap_temp.dat big_occupancy_temp.dat > big_occupancy.dat
rm *_temp.dat

cat > hbond.gnu << \EOF
reset
set terminal pngcairo
set boxwidth 0.9 absolute
set style fill solid 1.00 border lt -1
set key outside right top vertical Right noreverse noenhanced autotitle nobox
set style histogram clustered gap 5 title textcolor lt -1
set datafile missing '-'
set style data histograms
set xtics border in scale 0,0 nomirror rotate by -45  autojustify
set xtics ()
set yl "Average Hydrogen Bonds Number (/ps)"
set title "Hydrogen bonds between ligand and protein"
EOF
name=`ls *dcd | cut -f1 -d"."`
((n=${#d[*]}+1))

cat >> hbond.gnu << EOF
set output "${name}.png"
p for [i=2:${n}] 'big_occupancy.dat' u i:xticlabels(1) ti col
set output
EOF

The figure looks like below example.

image

Advertisements

2 thoughts on “Find residues of hydrogen bonds per picosecond at different time period.

  1. […] The following script is still BASH. Because the running time is only several seconds, I didn’t use Python. Previously, I introduced a way to find residues on receptor to demonstrate the averaged hydrogen bonds at each time range at https://sunxiaoquan.wordpress.com/2015/03/09/count-hydrogen-bond-number-of-all-resides-on-receptor-at-each-time-range/ and https://sunxiaoquan.wordpress.com/2015/03/15/find-residues-of-hydrogen-bonds-per-picosecond-at-diffe&#8230;. […]

    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