Thank you Justin,<br><br>you Perl script has been a great help!<br><br>Cheers,<br>Carla<br><br><div class="gmail_quote">On Mon, Oct 11, 2010 at 4:31 PM, Justin A. Lemkul <span dir="ltr"><<a href="mailto:jalemkul@vt.edu">jalemkul@vt.edu</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;"><br>
I wrote a Perl script to do a similar task (appended below). Perhaps it will be useful to you. I hope it works; I had to hack out some things that were specific to my needs and have done only limited testing.<br>
<br>
-Justin<br>
<br>
<br>
<br>
#!/usr/bin/perl<br>
<br>
#<br>
# <a href="http://plot_hbmap.pl" target="_blank">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 IDENTIFIERS<br>
# 2. hbmap.xpm<br>
# 3. hbond.ndx (modified to contain only the atom numbers in the [hbonds...] section, nothing else)<br>
#<br>
<br>
use strict;<br>
<br>
unless(@ARGV) {<br>
die "Usage: perl $0 -s structure.pdb -map hbmap.xpm -index hbond.ndx\n";<br>
}<br>
<br>
# define input hash<br>
my %args = @ARGV;<br>
<br>
# input variables<br>
my $map_in;<br>
my $struct;<br>
my $ndx;<br>
<br>
if (exists($args{"-map"})) {<br>
$map_in = $args{"-map"};<br>
} else {<br>
die "No -map specified!\n";<br>
}<br>
<br>
if (exists($args{"-s"})) {<br>
$struct = $args{"-s"};<br>
} else {<br>
die "No -s specified!\n";<br>
}<br>
<br>
if (exists($args{"-index"})) {<br>
$ndx = $args{"-index"};<br>
} else {<br>
die "No -index specified!\n";<br>
}<br>
<br>
# open the input<br>
open(MAP, "<$map_in") || die "Cannot open input map file!\n";<br>
my @map = <MAP>;<br>
close(MAP);<br>
<br>
open(STRUCT, "<$struct") || die "Cannot open input coordinate file!\n";<br>
my @coord = <STRUCT>;<br>
close(STRUCT);<br>
<br>
open(NDX, "<$ndx") || die "Cannot open input index file!\n";<br>
my @index = <NDX>;<br>
close(NDX);<br>
<br>
# 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);<br>
<br>
$nframes = $info[0];<br>
my @nframes = split('', $nframes);<br>
shift(@nframes); # get rid of the "<br>
$nframes = join('', @nframes);<br>
<br>
$nres = $info[1];<br>
}<br>
}<br>
<br>
print "Processing the map file...\n";<br>
print "There are $nres HB indices.\n";<br>
print "There are $nframes frames.\n";<br>
<br>
# 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>
}<br>
<br>
# donor/acceptor hashes for bookkeeping purposes<br>
my %donors;<br>
for (my $b=1; $b<=$nres; $b++) {<br>
$donors{$b} = 0;<br>
}<br>
<br>
my %acceptors;<br>
for (my $c=1; $c<=$nres; $c++) {<br>
$acceptors{$c} = 0;<br>
}<br>
<br>
# clean up the output - up to 18 lines of comments, etc.<br>
splice(@map, 0, 18);<br>
<br>
# 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>
}<br>
<br>
# 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.<br>
<br>
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]);<br>
<br>
# 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>
}<br>
<br>
print "Processing the index file...\n";<br>
<br>
# 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>
}<br>
<br>
<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;<br>
<br>
# 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];<br>
<br>
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>
}<br>
<br>
# 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", " ", "% Exist.");<br>
<br>
for (my $o=1; $o<=$nres; $o++) {<br>
printf(OUT "%10s\t%10s\t%10s\t%10s\t%10.3f\n", $donor_resn[$o], $donor_names[$o], $acceptor_resn[$o], $acceptor_names[$o], (($hbonds{$o}/$nframes)*100));<br>
}<br>
<br>
close(OUT);<br>
<br>
exit;<div><div></div><div class="h5"><br>
<br>
<br>
<br>
Carla Jamous wrote:<br>
<blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">
Hi everyone,<br>
<br>
I tried to analyze the H-bonds in my trajectory with g-hbond and I analysed the xpm and ndx file. But now I need to know the percentage of existence of each hbond during my trajectory. Is there a way to do it with a command line? Or is there a program (someone told me there are python programs for analysis of gromacs trajectories) to extract this information from the .xpm file?<br>
<br>
Thank you.<br>
<br>
Cheers,<br>
Carla<br>
<br>
</blockquote>
<br></div></div>
-- <br>
========================================<br>
<br>
Justin A. Lemkul<br>
Ph.D. Candidate<br>
ICTAS Doctoral Scholar<br>
MILES-IGERT Trainee<br>
Department of Biochemistry<br>
Virginia Tech<br>
Blacksburg, VA<br>
jalemkul[at]<a href="http://vt.edu" target="_blank">vt.edu</a> | (540) 231-9080<br>
<a href="http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin" target="_blank">http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin</a><br>
<br>
========================================<br><font color="#888888">
-- <br></font><div><div></div><div class="h5">
gmx-users mailing list <a href="mailto:gmx-users@gromacs.org" target="_blank">gmx-users@gromacs.org</a><br>
<a href="http://lists.gromacs.org/mailman/listinfo/gmx-users" target="_blank">http://lists.gromacs.org/mailman/listinfo/gmx-users</a><br>
Please search the archive at <a href="http://www.gromacs.org/Support/Mailing_Lists/Search" target="_blank">http://www.gromacs.org/Support/Mailing_Lists/Search</a> before posting!<br>
Please don't post (un)subscribe requests to the list. Use the www interface or send it to <a href="mailto:gmx-users-request@gromacs.org" target="_blank">gmx-users-request@gromacs.org</a>.<br>
Can't post? Read <a href="http://www.gromacs.org/Support/Mailing_Lists" target="_blank">http://www.gromacs.org/Support/Mailing_Lists</a><br>
</div></div></blockquote></div><br>