1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
|
"""
Hbond extension
Find all hydrogen bonds as determined by the cutoffs specified.
Uses PDB2PQR to determine donors and acceptors, and displays
all available bonds to file.
The original bonding parameters were an angle of 20.0,
distance of 3.30, and using the old method for calculating
distance.
The original parameters for WHAT-IF output was an angle of
90.0, distance of 3.30, and using the old method for
calculating distance.
Authors: Todd Dolinsky, Michael J Bradley, Julie Mitchell, and Kyle Monson
"""
__date__ = "17 February 2006"
__author__ = "Todd Dolinsky, Michael J Bradley, Julie Mitchell, and Kyle Monson"
# NOTE: This extension edited and updated on 05 August 2008 by Michael J Bradley
# and again on 17 October by Kyle Monson.
# NOTE: The current defaults for hbonds used below were utilized in the following study:
#Bradley MJ, Chivers PT, Baker NA. Molecular dynamics simulation of the
#Escherichia coli NikR protein: Equilibrium conformational fluctuations reveal
#inter-domain allosteric communication pathways. Journal of Molecular Biology, 378,
#1155-1173, 2008. http://dx.doi.org/10.1016/j.jmb.2008.03.010
from src.utilities import distance, getAngle
from src.routines import Cells
from math import cos
import extensions
ANGLE_CUTOFF = 30.0 # A - D - H(D) angle
DIST_CUTOFF = 3.4 # D to A distance
def addExtensionOptions(extensionGroup):
"""
Add options to set output type, angle cutoff, distance cutoff, and distance calculating method.
"""
extensionGroup.add_option('--whatif', dest='whatif', action='store_true', default=False,
help='Change hbond output to WHAT-IF format.')
extensionGroup.add_option('--angle_cutoff', dest='angle_cutoff', type="float", action='store', default=ANGLE_CUTOFF,
help='Angle cutoff to use when creating hbond data (default %s)' % ANGLE_CUTOFF)
extensionGroup.add_option('--distance_cutoff', dest='distance_cutoff', type="float", action='store', default=DIST_CUTOFF,
help='Distance cutoff to use when creating hbond data (default %s)' % DIST_CUTOFF)
extensionGroup.add_option('--old_distance_method', dest='old_distance_method', action='store_true', default=False,
help='Use distance from donor hydrogen to acceptor to calculate distance used with --distance_cutoff.')
def usage():
return 'Print a list of hydrogen bonds to {output-path}.hbond'
#TODO: replace this with a ''.format() call.
def _residueString(residue, name):
return '%4d %-4s (%4d ) %s %-4s' % \
(residue.resSeq, residue.name, residue.resSeq, residue.chainID, name)
def create_hbond_output(routines, outfile, whatif=False,
angleCutoff=ANGLE_CUTOFF,
distanceCutoff=DIST_CUTOFF,
oldDistanceMethod=False):
routines.write("Printing hydrogen bond list...\n")
output = extensions.extOutputHelper(routines, outfile)
cellsize = int(distanceCutoff + 1.0 + 1.0)
protein = routines.protein
routines.setDonorsAndAcceptors()
routines.cells = Cells(cellsize)
routines.cells.assignCells(protein)
for donor in protein.getAtoms():
# Grab the list of donors
if not donor.hdonor:
continue
donorhs = []
for bond in donor.bonds:
if bond.isHydrogen():
donorhs.append(bond)
if donorhs == []:
continue
# For each donor, grab all acceptors
closeatoms = routines.cells.getNearCells(donor)
for acc in closeatoms:
if not acc.hacceptor:
continue
if donor.residue == acc.residue:
continue
#TODO: do we need to do this for plain hbond stuff?
if whatif and (donor.residue.chainID == acc.residue.chainID):
continue
# Do new style distance check
if not oldDistanceMethod:
dist = distance(donor.getCoords(), acc.getCoords())
if dist > distanceCutoff:
continue
for donorh in donorhs:
# Do old style distance check
if oldDistanceMethod:
dist = distance(donorh.getCoords(), acc.getCoords())
if dist > distanceCutoff:
continue
# Do angle check
angle = getAngle(acc.getCoords(), donor.getCoords(), donorh.getCoords())
if angle > angleCutoff:
continue
if whatif:
if (donor.tempFactor > 60.0):
continue
if (acc.tempFactor > 60.0):
continue
thisBstring='B' if donor.isBackbone() else 'S'
thatBstring='B' if acc.isBackbone() else 'S'
score= (1.7/dist) * cos(angle * 3.142 / 180.0)
output.write(_residueString(donor.residue, donor.name))
output.write('-> ')
output.write(_residueString(acc.residue, acc.name))
output.write('Sym= 1 Val= %6.3lf DA=%6.2f DHA=%6.2f (%s-%s)\n' %
(score, dist, angle, thisBstring, thatBstring))
# outfile.write("%4d %-4s (%4d ) %s %-4s-> %4d %-4s (%4d ) %s %-4sSym= 1 Val= %6.3lf DA=%6.2f DHA=%6.2f (%s-%s)\n" % \
# (donor.residue.resSeq,donor.residue.name,donor.residue.resSeq, donor.residue.chainID,donor.name,acc.residue.resSeq,acc.residue.name,acc.residue.resSeq, acc.residue.chainID,acc.name, score, dist, angle, thisBstring, thatBstring))
else:
s = "Donor: %s %s\tAcceptor: %s %s\tdist: %.2f\tAngle: %.2f\n" % \
(donor.residue, donor.name, acc.residue, acc.name, dist, angle)
output.write(s)
routines.write("\n")
def run_extension(routines, outroot, options):
"""
Print a list of hydrogen bonds.
Parameters
routines: A link to the routines object
outroot: The root of the output name
options: options object
"""
outname = outroot + ".hbond"
with open(outname, "w") as outfile:
create_hbond_output(routines, outfile, whatif=options.whatif,
angleCutoff=options.angle_cutoff,
distanceCutoff=options.distance_cutoff,
oldDistanceMethod=options.old_distance_method)
|
