1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 """
24 Various structure related calculations.
25 """
26
27 import Numeric as N
28 from Biskit import molUtils as molU
29 import Biskit.tools as T
30
32 """
33 Collect a list with all potential hydrogen bonds in model.
34
35 @param model: PDBModel for which
36 @type model: PDBModel
37
38 @return: a list of potential hydrogen bonds containing a lists
39 with donor index, acceptor index, distance and angle.
40 @rtype: [ int, int, float, float ]
41 """
42 hbond_lst = []
43 donors = molU.hbonds['donors']
44 accept = molU.hbonds['acceptors']
45
46
47 d_ind = []
48 for res , aList in donors.items():
49 for a in aList:
50 if a in molU.hydrogenSynonyms.keys():
51 aList.append( molU.hydrogenSynonyms[a] )
52
53 d_ind += model.filterIndex( residue_name=res, name=aList )
54
55
56 a_ind = []
57 for res , aList in accept.items():
58 a_ind += model.filterIndex( residue_name=res, name=aList )
59
60
61 for d in d_ind:
62 d_xyz = model.xyz[d]
63 d_nr = model.atoms[d]['residue_number']
64 d_cid = model.atoms[d]['chain_id']
65 d_segi = model.atoms[d]['segment_id']
66
67 for a in a_ind:
68 a_xyz = model.xyz[a]
69 a_nr = model.atoms[a]['residue_number']
70 a_cid = model.atoms[a]['chain_id']
71 a_segi = model.atoms[a]['segment_id']
72
73 dist = N.sqrt( sum( (d_xyz - a_xyz)**2 ) )
74
75
76
77 if dist < 3.0 and not\
78 ( d_nr == a_nr and d_cid == a_cid and d_segi == a_segi ):
79
80
81 d_xyz_cov = xyzOfNearestCovalentNeighbour( d, model )
82 a_xyz_cov = xyzOfNearestCovalentNeighbour( a, model )
83 d_vec = d_xyz_cov - d_xyz
84 a_vec = a_xyz - a_xyz_cov
85
86 d_len = N.sqrt( sum( (d_vec)**2 ) )
87 a_len = N.sqrt( sum( (a_vec)**2 ) )
88
89 da_dot = N.dot( d_vec, a_vec)
90
91 angle = 180 - N.arccos( da_dot / (d_len * a_len) )*180/N.pi
92
93 if hbondCheck( angle, dist ):
94 hbond_lst += [[ d, a, dist, angle ]]
95
96 return hbond_lst
97
98
100 """
101 A h-bond is longer that 2.4A but shorter that 3.5A, measuring
102 from the nitrogen (donor) to the acceptor. The corresponding
103 length from the hydrogen (donor) to the acceptor (used here)
104 is 1.4 and 2.5A.
105
106 - the optimal length is about 2.8A (1.8A measured from the hydrogen).
107 - long h-bond (>3.6A) should be quite linear angle >150
108 - short h-bond (<3.6A) could be quite bent, but not more than 90-110
109
110 @param angle: angle of bond to check
111 @type angle: float
112 @param length: lenght of bond to check (D-H...A)
113 @type length: float
114
115 @return: if the test is passed the cutoff angle for a
116 hydrogen bond of the given length is returned, else None.
117 @rtype: float OR None
118 """
119
120
121 if length < 1.4 or length > 2.9:
122 return None
123
124
125
126 len_1, ang_1 = 1.5, 90
127 len_2, ang_2 = 2.8, 150
128 slope = (ang_2-ang_1)/(len_2-len_1)
129 intersect = ang_1 - slope*len_1
130
131 cutoff_angle = slope*length + intersect
132 if angle < cutoff_angle:
133 return None
134 else:
135 return cutoff_angle
136
137
139 """
140 Closest atom in the same residue as atom with index i
141
142 @param model: PDBModel
143 @type model: PDBModel
144 @param i: atom index
145 @type i: int
146
147 @return: coordinates of the nearest atom
148 @rtype: [float, float, float]
149 """
150 resModel = model.filter( residue_number=model.atoms[i]['residue_number'] )
151 dist = N.sqrt( N.sum( (resModel.xyz - model.xyz[i])**2 , 1) )
152
153
154 dist[ N.argmin(dist) ] = 100.
155
156 pos_shortest = N.nonzero( dist == min(dist) )[0]
157
158 return resModel.xyz[ pos_shortest ]
159
160
161 -def fasta( m , start=0, stop=None ):
162 """
163 Extract fasta sequence from model.
164
165 @param m: model
166 @type m: PDBModel
167 @param start: first residue
168 @type start: int
169 @param stop: last residue
170 @type stop: int
171
172 @return: fasta formated sequence and PDB code
173 @rtype: string
174 """
175 if not stop:
176 stop = m.lenResidues()
177
178 s = m.sequence()[start:stop]
179
180 n_chunks = len( s ) / 80
181
182 result = ">%s \n" % m.pdbCode
183
184 for i in range(0, n_chunks+1):
185
186 if i * 80 + 80 < len( s ):
187 chunk = s[i * 80 : i * 80 + 80]
188 else:
189 chunk = s[i * 80 :]
190
191 result += chunk
192 return result, m.pdbCode
193
194
195
196
197
198
200 """
201 Test class
202 """
203
204 - def run( self, local=0 ):
205 """
206 run function test
207
208 @param local: transfer local variables to global and perform
209 other tasks only when run locally
210 @type local: 1|0
211
212 @return: icmCad values
213 @rtype: [float]
214 """
215 from Biskit import PDBModel
216
217
218 m = PDBModel( T.testRoot() + '/lig/1A19.pdb' )
219 m = m.compress( m.maskProtein() )
220
221 hb = hbonds( m )
222
223 xyz = xyzOfNearestCovalentNeighbour( 40, m )
224
225 if local:
226 print '\nThe nearest covalently attached atom to the'
227 print ' atom with index 40 has the coordinates:'
228 print xyz
229
230 print 'Potential h-bonds in model:'
231 print '(donor index, acceptor index, distance and angle)'
232 for h in hb:
233 print h
234
235 globals().update( locals() )
236
237 return N.sum(N.ravel(hb[3:5])) + N.sum(xyz)
238
239
241 """
242 Precalculated result to check for consistent performance.
243
244 @return: icmCad values
245 @rtype: [float]
246 """
247 return 2025.8997840075292 + 152.687011719
248
249
250
251 if __name__ == '__main__':
252
253 test = Test()
254
255 assert abs( test.run( local=1 ) - test.expected_result() ) < 1e-8
256