<div>Dear Justin</div>
<div> </div>
<div>thank you very much.</div>
<div> </div>
<div>what I used as HB.pl exactly is as follows:</div>
<div> </div>
<div>is there problem in that?</div>
<div> </div>
<div>#!/usr/bin/perl</div>
<div>#<br># <a href="http://plot_hbmap.pl">plot_hbmap.pl</a> - plot the probability of finding a particular hydrogen bond<br># based on several input files:<br># 1. coordinate file (for atom naming) - MUST be a .pdb file with NO CHAIN<br>
IDENTIFIERS<br># 2. hbmap.xpm<br># 3. hbond.ndx (modified to contain only the atom numbers in the [hbonds...] <br>section, nothing else)<br>#</div>
<div>use strict;</div>
<div>unless(@ARGV) {<br> die "Usage: perl $0 -s structure.pdb -map hbmap.xpm -index hbond.ndx\n";<br>}</div>
<div># define input hash<br>my %args = @ARGV;</div>
<div># input variables<br>my $map_in;<br>my $struct;<br>my $ndx;</div>
<div>if (exists($args{"-map"})) {<br> $map_in = $args{"-map"};<br>} else {<br> die "No -map specified!\n";<br>}</div>
<div>if (exists($args{"-s"})) {<br> $struct = $args{"-s"};<br>} else {<br> die "No -s specified!\n";<br>}</div>
<div>if (exists($args{"-index"})) {<br> $ndx = $args{"-index"};<br>} else {<br> die "No -index specified!\n";<br>}</div>
<div># open the input<br>open(MAP, "<$map_in") || die "Cannot open input map file!\n";<br>my @map = <MAP>;<br>close(MAP);</div>
<div>open(STRUCT, "<$struct") || die "Cannot open input coordinate file!\n";<br>my @coord = <STRUCT>;<br>close(STRUCT);</div>
<div>open(NDX, "<$ndx") || die "Cannot open input index file!\n";<br>my @index = <NDX>;<br>close(NDX);</div>
<div># determine number of HB indices and frames<br>my $nres = 0;<br>my $nframes = 0;<br>for (my $i=0; $i<scalar(@map); $i++) {<br> if ($map[$i] =~ /static char/) {<br> my $res_line = $map[$i+1];<br> my @info = split(" ", $res_line);</div>
<div> $nframes = $info[0];<br> my @nframes = split('', $nframes);<br> shift(@nframes); # get rid of the "<br> $nframes = join('', @nframes);</div>
<div> $nres = $info[1];<br> }<br>}</div>
<div>print "Processing the map file...\n";<br>print "There are $nres HB indices.\n";<br>print "There are $nframes frames.\n";</div>
<div># initialize hashes for later output writing<br># counter $a holds the HB index from hbond.ndx<br>my %hbonds;<br>for (my $a=0; $a<$nres; $a++) {<br> $hbonds{$a+1} = 0;<br>}</div>
<div># donor/acceptor hashes for bookkeeping purposes<br>my %donors;<br>for (my $b=1; $b<=$nres; $b++) {<br> $donors{$b} = 0;<br>}</div>
<div>my %acceptors;<br>for (my $c=1; $c<=$nres; $c++) {<br> $acceptors{$c} = 0;<br>}</div>
<div># clean up the output - up to 18 lines of comments, etc.<br>splice(@map, 0, 18);</div>
<div># remove any "x-axis" or "y-axis" lines<br>for (my $n=0; $n<scalar(@map); $n++) {<br> if (($map[$n] =~ /x-axis/) || ($map[$n] =~ /y-axis/)) {<br> shift(@map);<br> $n--;<br>
}<br>}</div>
<div># There should now be $nres lines left in the file<br># The HB map for the last index is written first (top-down in .xpm file)<br># * Element 0 is index $nres, element 1 is $nres-1, etc.</div>
<div>for (my $i=$nres; $i>=1; $i--) {<br> # There will be $nframes+2 elements in @line (extra two are " at beginning<br> # and end of the line)<br> # Establish a conversion factor and split the input lines<br>
my $j = $nres - $i;<br> my @line = split('', $map[$j]);</div>
<div> # for each index, write to hash<br> for (my $k=1; $k<=($nframes+1); $k++) {<br> if ($line[$k] =~ /o/) {<br> $hbonds{$i}++;<br> }<br> }<br>}</div>
<div>print "Processing the index file...\n";</div>
<div># Open up the index file and work with it<br>for (my $n=0; $n<$nres; $n++) {<br> my @line = split(" ", $index[$n]);<br> $donors{$n+1} = $line[0];<br> $acceptors{$n+1} = $line[2];<br>}</div>
<div><br># some arrays for donor and acceptor atom names<br>my @donor_names;<br>my @donor_resn;<br>my @acceptor_names;<br>my @acceptor_resn;</div>
<div># Open up the structure file and work with it<br>print "Processing coordinate file...\n";<br>foreach $_ (@coord) {<br> my @line = split(" ", $_);<br> my $natom = $line[1];<br> my $name = $line[2];<br>
my $resn = $line[3];<br> my $resnum = $line[4];</div>
<div> if ($line[0] =~ /ATOM/) {<br> unless ($resn =~ /SOL/) {<br> for (my $z=1; $z<=$nres; $z++) {<br> if ($donors{$z} == $natom) {<br> $donor_names[$z] = $name;<br>
$donor_resn[$z] = join('', $resn, $resnum);<br> } elsif ($acceptors{$z} == $natom) {<br> $acceptor_names[$z] = $name;<br> $acceptor_resn[$z] = join('', $resn, $resnum);<br>
}<br> }<br> }<br> }<br>}</div>
<div># open a single output file for writing<br>open(OUT, ">>summary_HBmap.dat") || die "Cannot open output file!\n";<br>printf(OUT "%10s\t%10s\t%10s\t%10s\%10s\n", "# Donor", " ", "Acceptor", " ",<br>
"% Exist.");</div>
<div>for (my $o=1; $o<=$nres; $o++) {<br> printf(OUT "%10s\t%10s\t%10s\t%10s\t%10.3f\n", $donor_resn[$o],<br>$donor_names[$o], $acceptor_resn[$o], $acceptor_names[$o],<br>(($hbonds{$o}/$nframes)*100));<br>
}</div>
<div>close(OUT);</div>
<div>exit;<br><br clear="all"><br>-- <br></div><pre style="FONT-FAMILY: arial,helvetica,sans-serif">Leila Karami<br>Ph.D. student of Physical Chemistry<br>K.N. Toosi University of Technology<br>Theoretical Physical Chemistry Group</pre>
<br>