|
| 1 | + |
| 2 | +''' |
| 3 | +Calculate and display the dimensions of a protein. |
| 4 | +
|
| 5 | +
|
| 6 | +This is a first version, please use at your own risk! |
| 7 | +
|
| 8 | +REQUIREMENTS |
| 9 | +
|
| 10 | +numpy (http://numpy.scipy.org) that should be built into the newers versions of Pymol |
| 11 | +
|
| 12 | +(c) Pablo Guardado Calvo |
| 13 | +
|
| 14 | +Based on "inertia_tensor.py" (c) 2010 by Mateusz Maciejewski |
| 15 | +
|
| 16 | +License: MIT |
| 17 | +
|
| 18 | +''' |
| 19 | + |
| 20 | +from __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 | + |
| 61 | +from pymol import cmd, cgo |
| 62 | +from pymol.cgo import * |
| 63 | +import numpy |
| 64 | +from random import randint |
| 65 | + |
| 66 | +def 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 | + |
| 120 | +def 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 | + |
| 142 | +def 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 | + |
| 162 | +def 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 | + |
| 186 | +def 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 | + |
| 199 | +def 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 | + |
| 274 | +def 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 | + |
| 348 | +def draw_BB(selection): |
| 349 | + draw_AABB(selection) |
| 350 | + draw_IABB(selection) |
| 351 | + |
| 352 | +def draw_Protein_Dimensions(selection): |
| 353 | + draw_IABB(selection) |
| 354 | + cmd.hide("cgo") |
| 355 | + |
| 356 | +cmd.extend ("draw_Protein_Dimensions", draw_Protein_Dimensions) |
| 357 | +cmd.extend ("draw_BB", draw_BB) |
| 358 | + |
| 359 | + |
| 360 | + |
| 361 | + |
| 362 | + |
| 363 | + |
| 364 | + |
0 commit comments