Contando el número de enlaces de hidrógeno de múltiples archivos PDB

He estado tratando de encontrar una forma sistemática de contar la cantidad de enlaces de hidrógeno para múltiples archivos PDB.

DSSP me muestra el número total de enlaces de hidrógeno cuando hago dssp.exe 1A11.pdb, dándome el número total de 24(¿estoy en lo correcto, está bien?)

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

Pero cuando ejecuto este comando en R:

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

no parece tener el número de enlaces de hidrógeno. Muestra esto:

> 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

Entonces, mi pregunta: ¿el número total de enlaces de hidrógeno está oculto en algún lugar de esta salida o en algún otro lugar (o no existe)?

Si logro contar la cantidad de enlaces de hidrógeno en un solo PDB usando Bio3D y R, entonces puedo crear un ciclo e iterar sobre todos los archivos PDB restantes.

¡Gracias por todo!

Respuestas (2)

El número de enlaces de hidrógeno en realidad no se puede decir a partir del software, solo se puede inferir. Pero podría intentar usar STRIDE ( http://structure.usc.edu/stride/ ) que toma el nombre del archivo con un indicador -h para generar la cantidad de enlaces de hidrógeno. Luego podría escribir un pequeño script de shell que pasaría en cada archivo y almacenaría los datos que quisiera en el formato que quisiera.

¡Es cierto! Pensé que podría existir una forma directa de hacerlo :)

Esta es una solución inicial, pero se puede mejorar:

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)

Compilé una versión de Windows de Stride y la subí aquí: http://lcrserver.icmc.usp.br/~daniel/bin/stride.exe