0
$\begingroup$

I've been trying to figure out a systematic way to count the number of hydrogen bonds for multiple PDB files.

DSSP shows me the total number of hydrogen bonds when I make dssp.exe 1A11.pdb, giving me the total number of 24 (am I correct, ok?)

24 96.0   TOTAL NUMBER OF HYDROGEN BONDS OF TYPE O(I)-->H-N(J)  , SAME NUMBER PER 100 RESIDUES

But when I run this command in R:

> require(biod3d)
> x=dssp(read.pdb("1A11"), full=T)

it doesn't look to have the number of hydrogen bonds. It shows this:

> x$hbonds
     BP1 BP2 NH-O.1   E1 O-HN.1   E2 NH-O.2   E3 O-HN.2   E4 ChainBP1 ChainBP2 Chain1 Chain2 Chain3 Chain4
1-A   NA  NA     NA  0.0      5 -2.1     NA  0.0      4 -0.3     <NA>     <NA>   <NA>      A   <NA>      A
2-A   NA  NA      3 -0.2      6 -1.5      4 -0.2      7 -0.2     <NA>     <NA>      A      A      A      A
3-A   NA  NA      4 -0.2      7 -1.8      5 -0.2      2 -0.2     <NA>     <NA>      A      A      A      A
4-A   NA  NA      1 -0.3      8 -1.5      6 -0.2      7 -0.6     <NA>     <NA>      A      A      A      A
5-A   NA  NA      1 -2.1      9 -1.9      6 -0.2     10 -0.2     <NA>     <NA>      A      A      A      A
6-A   NA  NA      2 -1.5     10 -1.3      1 -0.3      5 -0.2     <NA>     <NA>      A      A      A      A
7-A   NA  NA      3 -1.8     11 -1.3      4 -0.6      5 -0.2     <NA>     <NA>      A      A      A      A
8-A   NA  NA      4 -1.5     12 -1.2      3 -0.2      6 -0.2     <NA>     <NA>      A      A      A      A
9-A   NA  NA      5 -1.9     13 -1.1      4 -0.2      7 -0.2     <NA>     <NA>      A      A      A      A
10-A  NA  NA      6 -1.3     14 -2.3      5 -0.2     13 -0.4     <NA>     <NA>      A      A      A      A
11-A  NA  NA      7 -1.3     15 -1.8     12 -0.2     16 -0.3     <NA>     <NA>      A      A      A      A
12-A  NA  NA      8 -1.2     16 -1.8     13 -0.2     11 -0.2     <NA>     <NA>      A      A      A      A
13-A  NA  NA      9 -1.1     17 -2.4     10 -0.4     18 -0.5     <NA>     <NA>      A      A      A      A
14-A  NA  NA     10 -2.3     18 -1.6     15 -0.2     13 -0.2     <NA>     <NA>      A      A      A      A
15-A  NA  NA     11 -1.8     19 -1.8     17 -0.2     20 -0.5     <NA>     <NA>      A      A      A      A
16-A  NA  NA     12 -1.8     20 -0.8     17 -0.3     14 -0.2     <NA>     <NA>      A      A      A      A
17-A  NA  NA     13 -2.4     21 -1.9     19 -0.2     16 -0.3     <NA>     <NA>      A      A      A      A
18-A  NA  NA     14 -1.6     22 -2.0     13 -0.5     21 -1.0     <NA>     <NA>      A      A      A      A
19-A  NA  NA     15 -1.8     23 -2.0     20 -0.3     24 -0.3     <NA>     <NA>      A      A      A      A
20-A  NA  NA     16 -0.8     24 -1.5     15 -0.5     19 -0.3     <NA>     <NA>      A      A      A      A
21-A  NA  NA     17 -1.9     25 -1.0     18 -1.0     19 -0.2     <NA>     <NA>      A      A      A      A
22-A  NA  NA     18 -2.0     25 -0.3     23 -0.2     20 -0.2     <NA>     <NA>      A      A      A      A
23-A  NA  NA     19 -2.0     22 -0.2     18 -0.3     21 -0.2     <NA>     <NA>      A      A      A      A
24-A  NA  NA     20 -1.5     23 -0.2     19 -0.3     22 -0.2     <NA>     <NA>      A      A      A      A
25-A  NA  NA     21 -1.0     23 -0.2     22 -0.3     24 -0.2     <NA>     <NA>      A      A      A      A

So, my question: Is the total number of hydrogen bonds hidden somewhere in this output or somewhere else (or doesn't exist)?

If I manage to count the number of hydrogen bonds in a single PDB using Bio3D and R, then I can create a loop and iterate over all remainder PDB files.

Thanks for all!

$\endgroup$
1
$\begingroup$

The number of hydrogen bonds cannot actually be said from software only inferred. But you could try using STRIDE (http://structure.usc.edu/stride/) which takes the file name with a -h flag to output the number of hydrogen bonds. You could then write a small shell script which would pass in each file and store the data you wanted in whatever format you liked.

$\endgroup$
  • $\begingroup$ It's true! I thought it could exist a direct way to do that :) $\endgroup$ – Daniel Bonetti Nov 22 '14 at 19:46
0
$\begingroup$

This an initial solution, but can be improved:

out = system2("f:/stride/stride.exe", "f:/stride/1a11.pdb -h", stdout=T)
nh  = strsplit(out[grepl("HBT", out)], "  ")[[1]][2]
print(nh) # Output: 26 (hydrogen bonds)

I compiled an Windows version of Stride and uploaded here: http://lcrserver.icmc.usp.br/~daniel/bin/stride.exe

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.