[gmx-users] Percentage of H-bonds
Carla Jamous
carlajamous at gmail.com
Tue Oct 26 14:15:03 CEST 2010
Please I have a problem:
I concatenated 3 trajectories and got one long trajectory on which I ran the
perl script it gives me lets say 60% for hbond 1.
But when I look in each of the 3 trajectories alone, this hbond doesn't even
exist in 2 of them.
Do you have an idea where the problem might be because I checked if it's
well concatenated and it is. And I can't figure out why I have this error.
Thanks,
Carla
On Tue, Oct 26, 2010 at 2:08 PM, Justin A. Lemkul <jalemkul at vt.edu> wrote:
>
>
> Carla Jamous wrote:
>
>> 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.
>>
>>
> The order is the same, and the indices are mapped exactly as you would
> expect them. The .xpm file is written, however, from top down. So the
> final index in hbond.ndx is the top line in the matrix. As such, even
> though they are in the same "order," per se, the two files (.xpm and .ndx)
> have to be read in opposite directions.
>
> -Justin
>
> Thank you.
>>
>> Carla
>>
>>
>> On Mon, Oct 11, 2010 at 4:31 PM, Justin A. Lemkul <jalemkul at vt.edu<mailto:
>> 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 <http://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 <http://vt.edu> | (540) 231-9080
>>
>> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
>>
>> ========================================
>> -- gmx-users mailing list gmx-users at gromacs.org
>> <mailto: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
>> <mailto:gmx-users-request at gromacs.org>.
>>
>> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>>
>>
> --
> ========================================
>
> 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/52882024/attachment.html>
More information about the gromacs.org_gmx-users
mailing list