[gmx-users] Percentage of H-bonds

[Carla Jamous]

Hi Justin,

Please can you explain this particular comment in your perl script: "# There
should now be $nres lines left in the file
# The HB map for the last index is written first (top-down in .xpm file)
#   * Element 0 is index $nres, element 1 is $nres-1, etc."

Is the perl script reading the index file beginning from the bottom? If so,
why is that? Because in the manual, it says that in the matrix, the ordering
is identical to that in the index file.

Thank you.

Carla

On Mon, Oct 11, 2010 at 4:31 PM, Justin A. Lemkul <jalemkul at vt.edu> wrote:

>
> 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.
>
> -Justin
>
>
>
> #!/usr/bin/perl
>
> #
> # plot_hbmap.pl - plot the probability of finding a particular hydrogen
> bond
> # based on several input files:
> #   1. coordinate file (for atom naming) - MUST be a .pdb file with NO
> CHAIN IDENTIFIERS
> #   2. hbmap.xpm
> #   3. hbond.ndx (modified to contain only the atom numbers in the
> [hbonds...] section, nothing else)
> #
>
> use strict;
>
> unless(@ARGV) {
>    die "Usage: perl $0 -s structure.pdb -map hbmap.xpm -index hbond.ndx\n";
> }
>
> # define input hash
> my %args = @ARGV;
>
> # input variables
> my $map_in;
> my $struct;
> my $ndx;
>
> if (exists($args{"-map"})) {
>    $map_in = $args{"-map"};
> } else {
>    die "No -map specified!\n";
> }
>
> if (exists($args{"-s"})) {
>    $struct = $args{"-s"};
> } else {
>    die "No -s specified!\n";
> }
>
> if (exists($args{"-index"})) {
>    $ndx = $args{"-index"};
> } else {
>    die "No -index specified!\n";
> }
>
> # open the input
> open(MAP, "<$map_in") || die "Cannot open input map file!\n";
> my @map = <MAP>;
> close(MAP);
>
> open(STRUCT, "<$struct") || die "Cannot open input coordinate file!\n";
> my @coord = <STRUCT>;
> close(STRUCT);
>
> open(NDX, "<$ndx") || die "Cannot open input index file!\n";
> my @index = <NDX>;
> close(NDX);
>
> # determine number of HB indices and frames
> my $nres = 0;
> my $nframes = 0;
> for (my $i=0; $i<scalar(@map); $i++) {
>    if ($map[$i] =~ /static char/) {
>        my $res_line = $map[$i+1];
>        my @info = split(" ", $res_line);
>
>        $nframes = $info[0];
>        my @nframes = split('', $nframes);
>        shift(@nframes);    # get rid of the "
>        $nframes = join('', @nframes);
>
>        $nres = $info[1];
>    }
> }
>
> print "Processing the map file...\n";
> print "There are $nres HB indices.\n";
> print "There are $nframes frames.\n";
>
> # initialize hashes for later output writing
> # counter $a holds the HB index from hbond.ndx
> my %hbonds;
> for (my $a=0; $a<$nres; $a++) {
>    $hbonds{$a+1} = 0;
> }
>
> # donor/acceptor hashes for bookkeeping purposes
> my %donors;
> for (my $b=1; $b<=$nres; $b++) {
>    $donors{$b} = 0;
> }
>
> my %acceptors;
> for (my $c=1; $c<=$nres; $c++) {
>    $acceptors{$c} = 0;
> }
>
> # clean up the output - up to 18 lines of comments, etc.
> splice(@map, 0, 18);
>
> # remove any "x-axis" or "y-axis" lines
> for (my $n=0; $n<scalar(@map); $n++) {
>    if (($map[$n] =~ /x-axis/) || ($map[$n] =~ /y-axis/)) {
>            shift(@map);
>            $n--;
>    }
> }
>
> # There should now be $nres lines left in the file
> # The HB map for the last index is written first (top-down in .xpm file)
> #   * Element 0 is index $nres, element 1 is $nres-1, etc.
>
> for (my $i=$nres; $i>=1; $i--) {
>    # There will be $nframes+2 elements in @line (extra two are " at
> beginning
>    # and end of the line)
>    # Establish a conversion factor and split the input lines
>    my $j = $nres - $i;
>    my @line = split('', $map[$j]);
>
>    # for each index, write to hash
>    for (my $k=1; $k<=($nframes+1); $k++) {
>        if ($line[$k] =~ /o/) {
>            $hbonds{$i}++;
>        }
>    }
> }
>
> print "Processing the index file...\n";
>
> # Open up the index file and work with it
> for (my $n=0; $n<$nres; $n++) {
>    my @line = split(" ", $index[$n]);
>    $donors{$n+1} = $line[0];
>    $acceptors{$n+1} = $line[2];
> }
>
>
> # some arrays for donor and acceptor atom names
> my @donor_names;
> my @donor_resn;
> my @acceptor_names;
> my @acceptor_resn;
>
> # Open up the structure file and work with it
> print "Processing coordinate file...\n";
> foreach $_ (@coord) {
>    my @line = split(" ", $_);
>    my $natom = $line[1];
>    my $name = $line[2];
>    my $resn = $line[3];
>    my $resnum = $line[4];
>
>    if ($line[0] =~ /ATOM/) {
>        unless ($resn =~ /SOL/) {
>            for (my $z=1; $z<=$nres; $z++) {
>                if ($donors{$z} == $natom) {
>                    $donor_names[$z] = $name;
>                    $donor_resn[$z] = join('', $resn, $resnum);
>                } elsif ($acceptors{$z} == $natom) {
>                    $acceptor_names[$z] = $name;
>                    $acceptor_resn[$z] = join('', $resn, $resnum);
>                }
>            }
>        }
>    }
> }
>
> # open a single output file for writing
> open(OUT, ">>summary_HBmap.dat") || die "Cannot open output file!\n";
> printf(OUT "%10s\t%10s\t%10s\t%10s\%10s\n", "#    Donor", " ", "Acceptor",
> " ", "% Exist.");
>
> for (my $o=1; $o<=$nres; $o++) {
>    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));
> }
>
> close(OUT);
>
> exit;
>
>
>
>
> Carla Jamous wrote:
>
>> Hi everyone,
>>
>> 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?
>>
>> Thank you.
>>
>> Cheers,
>> Carla
>>
>>
> --
> ========================================
>
> Justin A. Lemkul
> Ph.D. Candidate
> ICTAS Doctoral Scholar
> MILES-IGERT Trainee
> Department of Biochemistry
> Virginia Tech
> Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
>
> ========================================
> --
> gmx-users mailing list    gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the www interface
> or send it to gmx-users-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20101026/e1ef52bf/attachment.html>