· 7 years ago · Jan 04, 2019, 08:44 PM
1
2'''
3Calculate and display the dimensions of a protein.
4
5
6This is a first version, please use at your own risk!
7
8REQUIREMENTS
9
10numpy (http://numpy.scipy.org) that should be built into the newers versions of Pymol
11
12(c) Pablo Guardado Calvo
13
14Based on "inertia_tensor.py" (c) 2010 by Mateusz Maciejewski
15
16License: MIT
17
18'''
19
20from __future__ import print_function
21
22__author__ = 'Pablo Guardado Calvo'
23__version__ = '0.2'
24__email__ = 'pablo.guardado (at) gmail.com'
25__date__ = '13/08/2015'
26__modification_date__ = '05/02/2018'
27__modification_reason__ = 'Error in the code produced sometimes inverted structures'
28
29###########################################################################################################################################################
30# USAGE
31#
32# The idea behing this script is to calculate an aproximate minimal bounding box to extract the cell dimensions of a protein. To calculate the minimal bounding
33# is not trivial and usually the Axis Aligned Bounding Box (AABB) does not show up the real dimensions of the protein. This script calculates the inertia tensor
34# of the object, extract the eigenvalues and use them to rotate the molecule (using as rotation matrix the transpose of the eigenvalues matrix). The result is that
35# the molecule is oriented with the inertia axis aligned with the cartesian axis. A new Bounding Box is calculated that is called Inertia Axis Aligned Bounding Box
36#(IABB), whose volume is always lower than AABB volume, and in many cases will correspond with the lowest volume. Of course, maybe it exists another Bounding Box
37# with a lower volume (the minimal Bounding Box).
38#
39# As always with these type of things, you have to use at your own risk. I did not try all the possible combinations, but if you find a bug, do
40# not hesitate to contact me (pablo.guardado (at) gmail.com) or try to modify the code for yourself to correct it.
41#
42# To load the script just type:
43#
44# run path-to-the-script/Draw_Protein_Dimensions.py
45#
46# or if you want something more permanent add the previous line to your .pymolrc file
47#
48# The script works just typing:
49#
50# draw_Protein_Dimensions selection
51#
52# This will draw the cell dimensions of your selection based on a IABB. It also generates the IABB box and the inertia axis, you just need to do "show cgo" to display them.
53#
54# You could also try:
55#
56# draw_BB selection
57#
58# This will draw the AABB and IABB boxes with their cell dimensions and show in the command line their volumes, you can compare both of them.
59############################################################################################################################################################
60
61from pymol import cmd, cgo
62from pymol.cgo import *
63import numpy
64from random import randint
65
66def matriz_inercia(selection):
67 '''
68 DESCRIPTION
69
70 The method calculates the mass center, the inertia tensor and the eigenvalues and eigenvectors
71 for a given selection. Mostly taken from inertia_tensor.py
72 '''
73
74 model = cmd.get_model(selection)
75 totmass = 0.0
76 x,y,z = 0,0,0
77 for a in model.atom:
78 m = a.get_mass()
79 x += a.coord[0]*m
80 y += a.coord[1]*m
81 z += a.coord[2]*m
82 totmass += m
83 global cM
84 cM = numpy.array([x/totmass, y/totmass, z/totmass])
85
86
87 I = []
88 for index in range(9):
89 I.append(0)
90
91 for a in model.atom:
92 temp_x, temp_y, temp_z = a.coord[0], a.coord[1], a.coord[2]
93 temp_x -= x
94 temp_y -= y
95 temp_z -= z
96
97 I[0] += a.get_mass() * (temp_y**2 + temp_z**2)
98 I[1] -= a.get_mass() * temp_x * temp_y
99 I[2] -= a.get_mass() * temp_x * temp_z
100 I[3] -= a.get_mass() * temp_x * temp_y
101 I[4] += a.get_mass() * (temp_x**2 + temp_z**2)
102 I[5] -= a.get_mass() * temp_y * temp_z
103 I[6] -= a.get_mass() * temp_x * temp_z
104 I[7] -= a.get_mass() * temp_y * temp_z
105 I[8] += a.get_mass() * (temp_x**2 + temp_y**2)
106
107 global tensor
108 tensor = numpy.array([(I[0:3]), (I[3:6]), (I[6:9])])
109
110 global autoval, autovect, ord_autoval, ord_autovect
111 autoval, autovect = numpy.linalg.eig(tensor)
112 auto_ord = numpy.argsort(autoval)
113 ord_autoval = autoval[auto_ord]
114 ord_autovect_complete = autovect[:, auto_ord].T
115 ord_autovect = numpy.around(ord_autovect_complete, 3)
116
117 return ord_autoval
118
119
120def draw_inertia_axis(selection):
121 '''
122 DESCRIPTION
123
124 This method draw the inertia axis calculated with the method matriz_inercia.
125 '''
126
127 matriz_inercia(selection)
128 axis1 = ord_autovect[0]
129 x1, y1, z1 = cM[0], cM[1], cM[2]
130 x2, y2, z2 = cM[0]+50*axis1[0], cM[1]+50*axis1[1], cM[2]+50*axis1[2]
131 eje1 = [cgo.CYLINDER, x1, y1, z1, x2, y2, z2, 0.6, 1, 0, 0, 1, 0, 0, 0.0]
132 cmd.load_cgo(eje1, 'Inertia_Axis1')
133 axis2 = ord_autovect[1]
134 x3, y3, z3 = cM[0]+40*axis2[0], cM[1]+40*axis2[1], cM[2]+40*axis2[2]
135 eje1 = [cgo.CYLINDER, x1, y1, z1, x3, y3, z3, 0.6, 1, 0.5, 0, 1, 0.5, 0, 0.0]
136 cmd.load_cgo(eje1, 'Inertia_Axis2')
137 axis4 = ord_autovect[2]
138 x4, y4, z4 = cM[0]+30*axis4[0], cM[1]+30*axis4[1], cM[2]+30*axis4[2]
139 eje1 = [cgo.CYLINDER, x1, y1, z1, x4, y4, z4, 0.6, 1, 1, 0, 1, 1, 0, 0.0]
140 cmd.load_cgo(eje1, 'Inertia_Axis3')
141
142def translacion_cM(selection):
143 '''
144 DESCRIPTION
145
146 Translate the center of mass of the molecule to the origin.
147 '''
148 model = cmd.get_model(selection)
149 totmass = 0.0
150 x,y,z = 0,0,0
151 for a in model.atom:
152 m = a.get_mass()
153 x += a.coord[0]*m
154 y += a.coord[1]*m
155 z += a.coord[2]*m
156 totmass += m
157 cM = numpy.array([x/totmass, y/totmass, z/totmass])
158 trans_array = ([1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, -cM[0], -cM[1], -cM[2], 1])
159 model_trans = cmd.transform_selection(selection, trans_array)
160
161
162def rotacion_orig(selection):
163 '''
164 DESCRIPTION
165
166 Find the proper rotation matrix, i.e. the transpose of the matrix formed by the eigenvectors of the inertia tensor
167 '''
168
169 translacion_cM(selection)
170 matriz_inercia(selection)
171 global transf, transf_array, ord_autovect_array, transf_array_print
172 ord_autovect_array = numpy.array([[ord_autovect[0][0], ord_autovect[0][1], ord_autovect[0][2]],
173 [ord_autovect[1][0], ord_autovect[1][1], ord_autovect[1][2]],
174 [ord_autovect[2][0], ord_autovect[2][1], ord_autovect[2][2]]])
175 if numpy.linalg.det(ord_autovect_array) < 0.:
176 ord_autovect_array = numpy.array([[ord_autovect[2][0], ord_autovect[2][1], ord_autovect[2][2]],
177 [ord_autovect[1][0], ord_autovect[1][1], ord_autovect[1][2]],
178 [ord_autovect[0][0], ord_autovect[0][1], ord_autovect[0][2]]])
179 transf = numpy.transpose(ord_autovect_array)
180 transf_array = numpy.array([transf[0][0], transf[0][1], transf[0][2], 0,
181 transf[1][0], transf[1][1], transf[1][2], 0,
182 transf[2][0], transf[2][1], transf[2][2], 0,
183 0, 0, 0, 1])
184
185
186def transformar(selection):
187 '''
188 DESCRIPTION
189
190 Rotate the molecule and draw the inertia axis.
191 '''
192
193 rotacion_orig(selection)
194 model_rot = cmd.transform_selection(selection, transf_array, homogenous=0, transpose=1);
195 draw_inertia_axis(selection)
196
197
198
199def draw_AABB(selection):
200 """
201 DESCRIPTION
202 For a given selection, draw the Axes Aligned bounding box around it without padding. Code taken and modified from DrawBoundingBox.py.
203
204 """
205
206 AA_original = selection + "_original"
207 model_orig = cmd.create(AA_original, selection)
208
209 ([min_X, min_Y, min_Z],[max_X, max_Y, max_Z]) = cmd.get_extent(AA_original)
210
211 print("The Axis Aligned Bounding Box (AABB) dimensions are (%.2f, %.2f, %.2f)" % (max_X-min_X, max_Y-min_Y, max_Z-min_Z))
212 print("The Axis Aligned Bounding Box (AABB) volume is %.2f A3" % ((max_X-min_X)*(max_Y-min_Y)*(max_Z-min_Z)))
213
214
215 min_X = min_X
216 min_Y = min_Y
217 min_Z = min_Z
218 max_X = max_X
219 max_Y = max_Y
220 max_Z = max_Z
221
222 boundingBox = [
223 LINEWIDTH, float(2),
224
225 BEGIN, LINES,
226 COLOR, float(1), float(1), float(0),
227
228 VERTEX, min_X, min_Y, min_Z,
229 VERTEX, min_X, min_Y, max_Z,
230 VERTEX, min_X, max_Y, min_Z,
231 VERTEX, min_X, max_Y, max_Z,
232 VERTEX, max_X, min_Y, min_Z,
233 VERTEX, max_X, min_Y, max_Z,
234 VERTEX, max_X, max_Y, min_Z,
235 VERTEX, max_X, max_Y, max_Z,
236 VERTEX, min_X, min_Y, min_Z,
237 VERTEX, max_X, min_Y, min_Z,
238 VERTEX, min_X, max_Y, min_Z,
239 VERTEX, max_X, max_Y, min_Z,
240 VERTEX, min_X, max_Y, max_Z,
241 VERTEX, max_X, max_Y, max_Z,
242 VERTEX, min_X, min_Y, max_Z,
243 VERTEX, max_X, min_Y, max_Z,
244 VERTEX, min_X, min_Y, min_Z,
245 VERTEX, min_X, max_Y, min_Z,
246 VERTEX, max_X, min_Y, min_Z,
247 VERTEX, max_X, max_Y, min_Z,
248 VERTEX, min_X, min_Y, max_Z,
249 VERTEX, min_X, max_Y, max_Z,
250 VERTEX, max_X, min_Y, max_Z,
251 VERTEX, max_X, max_Y, max_Z,
252
253 END
254 ]
255
256 p0 = '_0' + str(randint(0, 100))
257 p1 = '_1' + str(randint(0, 100))
258 p2 = '_2' + str(randint(0, 100))
259 p3 = '_3' + str(randint(0, 100))
260 cmd.pseudoatom (pos=[min_X, min_Y, min_Z], object=p0)
261 cmd.pseudoatom (pos=[min_X, min_Y, max_Z], object=p1)
262 cmd.pseudoatom (pos=[min_X, max_Y, min_Z], object=p2)
263 cmd.pseudoatom (pos=[max_X, min_Y, min_Z], object=p3)
264 cmd.distance(None, p0, p3)
265 cmd.distance(None, p0, p2)
266 cmd.distance(None, p0, p1)
267 cmd.hide("nonbonded")
268
269 boxName = "box_AABB_" + str(randint(0, 100))
270 cmd.load_cgo(boundingBox,boxName)
271 return boxName
272
273
274def draw_IABB(selection):
275 """
276 DESCRIPTION
277 For a given selection, draw the Inertia Axes Aligned bounding box around it without padding. Code taken and modified from DrawBoundingBox.py.
278
279 """
280
281 transformar(selection)
282
283 ([minX, minY, minZ],[maxX, maxY, maxZ]) = cmd.get_extent(selection)
284
285 print("The Inertia Axis Aligned Bounding Box (IABB) dimensions are (%.2f, %.2f, %.2f)" % (maxX-minX, maxY-minY, maxZ-minZ))
286 print("The Inertia Axis Aligned Bounding Box (IABB) volume is %.2f A3" % ((maxX-minX)*(maxY-minY)*(maxZ-minZ)))
287
288
289 minX = minX
290 minY = minY
291 minZ = minZ
292 maxX = maxX
293 maxY = maxY
294 maxZ = maxZ
295
296 boundingBox = [
297 LINEWIDTH, float(2),
298
299 BEGIN, LINES,
300 COLOR, float(1), float(0), float(0),
301
302 VERTEX, minX, minY, minZ,
303 VERTEX, minX, minY, maxZ,
304 VERTEX, minX, maxY, minZ,
305 VERTEX, minX, maxY, maxZ,
306 VERTEX, maxX, minY, minZ,
307 VERTEX, maxX, minY, maxZ,
308 VERTEX, maxX, maxY, minZ,
309 VERTEX, maxX, maxY, maxZ,
310 VERTEX, minX, minY, minZ,
311 VERTEX, maxX, minY, minZ,
312 VERTEX, minX, maxY, minZ,
313 VERTEX, maxX, maxY, minZ,
314 VERTEX, minX, maxY, maxZ,
315 VERTEX, maxX, maxY, maxZ,
316 VERTEX, minX, minY, maxZ,
317 VERTEX, maxX, minY, maxZ,
318 VERTEX, minX, minY, minZ,
319 VERTEX, minX, maxY, minZ,
320 VERTEX, maxX, minY, minZ,
321 VERTEX, maxX, maxY, minZ,
322 VERTEX, minX, minY, maxZ,
323 VERTEX, minX, maxY, maxZ,
324 VERTEX, maxX, minY, maxZ,
325 VERTEX, maxX, maxY, maxZ,
326
327 END
328 ]
329
330 p4 = '_4' + str(randint(0, 100))
331 p5 = '_5' + str(randint(0, 100))
332 p6 = '_6' + str(randint(0, 100))
333 p7 = '_7' + str(randint(0, 100))
334 cmd.pseudoatom (pos=[minX, minY, minZ], object=p4)
335 cmd.pseudoatom (pos=[minX, minY, maxZ], object=p5)
336 cmd.pseudoatom (pos=[minX, maxY, minZ], object=p6)
337 cmd.pseudoatom (pos=[maxX, minY, minZ], object=p7)
338 cmd.distance(None, p4, p7)
339 cmd.distance(None, p4, p6)
340 cmd.distance(None, p4, p5)
341 cmd.hide("nonbonded")
342
343 boxName = "box_IABB_" + str(randint(0, 100))
344 cmd.load_cgo(boundingBox,boxName)
345 return boxName
346
347
348def draw_BB(selection):
349 draw_AABB(selection)
350 draw_IABB(selection)
351
352def draw_Protein_Dimensions(selection):
353 draw_IABB(selection)
354 cmd.hide("cgo")
355
356cmd.extend ("draw_Protein_Dimensions", draw_Protein_Dimensions)
357cmd.extend ("draw_BB", draw_BB)