1# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- 2# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 3# 4# MDAnalysis --- https://www.mdanalysis.org 5# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors 6# (see the file AUTHORS for the full list of names) 7# 8# Released under the GNU Public Licence, v2 or any higher version 9# 10# Please cite your use of MDAnalysis in published work: 11# 12# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, 13# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. 14# MDAnalysis: A Python package for the rapid analysis of molecular dynamics 15# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th 16# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy. 17# 18# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. 19# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. 20# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 21# 22# Python implementation of FORTRAN HELANAL 23# [Bansal et al, J Biomol Struct Dyn. 17 (2000), 811] 24# Copyright (c) 2009 Benjamin Hall <benjamin.hall@bioch.ox.ac.uk> 25# Published under the GNU General Public License v2 or higher 26# 27# Copyright (c) 2011 Oliver Beckstein <orbeckst@gmail.com> 28# Integrated into MDAnalysis and NumPy-fied 29 30 31""" 32HELANAL --- analysis of protein helices 33======================================= 34 35:Author: Benjamin Hall <benjamin.a.hall@ucl.ac.uk>, Oliver Beckstein, Xavier Deupi 36:Year: 2009, 2011, 2013 37:License: GNU General Public License v2 (or higher) 38 39The :mod:`MDAnalysis.analysis.helanal` module is a Python implementation of the 40HELANAL_ algorithm [Bansal2000]_ in `helanal.f`_, which is also available 41through the `HELANAL webserver`_. 42 43Please cite the paper [Bansal2000]_ (and possibly [Kumar1996]_ and 44[Kumar1998]_) in published work when using 45:mod:`~MDAnalysis.analysis.helanal.helanal_trajectory` or 46:mod:`~MDAnalysis.analysis.helanal.helanal_main`. 47 48HELANAL_ quantifies the geometry of helices in proteins on the basis of their Cα 49atoms alone. It can extract the helices from the structure files and then 50characterises the overall geometry of each helix as being linear, curved or 51kinked, in terms of its local structural features, viz. local helical twist and 52rise, virtual torsion angle, local helix origins and bending angles between 53successive local helix axes. Even helices with large radius of curvature are 54unambiguously identified as being linear or curved. The program can also be 55used to differentiate a kinked helix and other motifs, such as helix-loop-helix 56or a helix-turn-helix (with a single residue linker) with the help of local 57bending angles. In addition to these, the program can also be used to 58characterise the helix start and end as well as other types of secondary 59structures. 60 61.. _HELANAL: http://www.ccrnp.ncifcrf.gov/users/kumarsan/HELANAL/helanal.html 62.. _Helanal webserver: http://nucleix.mbu.iisc.ernet.in/helanal/helanal.shtml 63.. `helanal.f`: http://www.webcitation.org/5y1RpVJtF 64.. _`helanal.f`: http://www.ccrnp.ncifcrf.gov/users/kumarsan/HELANAL/helanal.f 65 66Background 67---------- 68 69From the HELANAL_ home page: 70 71HELANAL_ can be used to characterize the geometry of helices with a minimum 9 72residues. The geometry of an alpha helix is characterized by computing local 73helix axes and local helix origins for four contiguous C-Alpha atoms, using the 74procedure of Sugeta and Miyazawa [Sugeta1967]_ and sliding this window over 75the length of the helix in steps of one C-alpha atom. 76 77The angles between successive local helix axes can identify *local bends* or 78*kinks* as well as occurrence of *smooth curvature* in the helix. A matrix, whose 79elements *M(I, J)* are the *bending angles* between local helix axes *I* and *J*, 80is obtained to get an idea about the overall geometry of the helix. 81 82*Unit twist* and *unit height* of the alpha helix are also computed to analyze the 83uniformity of the helix. The *local helix origins* trace out the path described 84by the helix in three dimensional space. The local helix origins are reoriented 85in *X-Y* plane and the reoriented points are used to fit a circle as well as a 86line, by least squares method. Based on the relative goodness of line and 87circle fit to local helix origins, the helix is *classified as being linear or 88curved*. A helix is classified as being *kinked*, if at least one local bending 89angle in the middle of the helix is greater than 20 degrees. 90 91 92References 93---------- 94 95.. [Sugeta1967] Sugeta, H. and Miyazawa, T. 1967. General method for 96 calculating helical parameters of polymer chains from bond lengths, bond 97 angles and internal rotation angles. *Biopolymers* 5 673 - 679 98 99.. [Kumar1996] Kumar, S. and Bansal, M. 1996. Structural and sequence 100 characteristics of long alpha-helices in globular proteins. *Biophysical 101 Journal* 71(3):1574-1586. 102 103.. [Kumar1998] Kumar, S. and Bansal, M. 1998. Geometrical and sequence 104 characteristics of alpha helices in globular proteins. *Biophysical Journal* 105 75(4):1935-1944. 106 107.. [Bansal2000] Bansal M, Kumar S, Velavan R. 2000. 108 HELANAL - A program to characterise helix geometry in proteins. 109 *J Biomol Struct Dyn.* 17(5):811-819. 110 111Functions 112--------- 113 114.. autofunction:: helanal_trajectory 115.. autofunction:: helanal_main 116 117""" 118from __future__ import print_function, division, absolute_import 119from six.moves import range, zip 120 121import os 122 123import numpy as np 124 125import MDAnalysis 126from MDAnalysis import FinishTimeException 127from MDAnalysis.lib.log import ProgressMeter 128from MDAnalysis.lib import mdamath 129 130import logging 131logger = logging.getLogger("MDAnalysis.analysis.helanal") 132 133def center(coordinates): 134 """Return the geometric center (centroid) of the coordinates. 135 136 Coordinates must be "list of cartesians", i.e. a Nx3 array. 137 """ 138 return np.mean(coordinates, axis=0) 139 140def vecnorm(a): 141 """Return a/|a|""" 142 return a / mdamath.norm(a) 143 144def wrapangle(angle): 145 """Wrap angle (in radians) to be within -pi < angle =< pi""" 146 if angle > np.pi: 147 angle -= 2 * np.pi 148 elif angle <= -np.pi: 149 angle += 2 * np.pi 150 return angle 151 152def sample_sd(a, dummy): 153 return np.std(a, ddof=1) 154 155def mean_abs_dev(a, mean_a=None): 156 if mean_a is None: 157 mean_a = np.mean(a) 158 return np.mean(np.fabs(a - mean_a)) 159 160def helanal_trajectory(universe, selection="name CA", 161 begin=None, finish=None, 162 matrix_filename="bending_matrix.dat", 163 origin_pdbfile="origin.pdb", 164 summary_filename="summary.txt", 165 screw_filename="screw.xvg", 166 tilt_filename="local_tilt.xvg", 167 fitted_tilt_filename="fit_tilt.xvg", 168 bend_filename="local_bend.xvg", 169 twist_filename="unit_twist.xvg", 170 prefix="helanal_", ref_axis=None, 171 verbose=False): 172 """Perform HELANAL helix analysis on all frames in `universe`. 173 174 Parameters 175 ---------- 176 universe : Universe 177 selection : str (optional) 178 selection string that selects Calpha atoms [``"name CA"``] 179 begin : float (optional) 180 start analysing for time (ps) >= *begin*; ``None`` starts from the 181 beginning [``None``] 182 finish : float (optional) 183 stop analysis for time (ps) =< *finish*; ``None`` goes to the 184 end of the trajectory [``None``] 185 matrix_filename : str (optional) 186 Output file- bending matrix [``"bending_matrix.dat"``] 187 origin_pdbfile : str (optional) 188 Output file- origin pdb file [``"origin.pdb"``] 189 summary_filename : str (optional) 190 Output file- all of the basic data [``"summary.txt"``] 191 screw_filename : str (optional) 192 Output file- local tilts of individual residues from 2 to n-1 193 [``"screw.xvg"``] 194 tilt_filename : str (optional) 195 Output file- tilt of line of best fit applied to origin axes 196 [``"local_tilt.xvg"``] 197 bend_filename : str (optional) 198 Output file- local bend angles between successive local helix axes 199 [``"local_bend.xvg"``] 200 twist_filename : str (optional) 201 Output file- local unit twist between successive helix turns 202 [``"unit_twist.xvg"``] 203 prefix : str (optional) 204 Prefix to add to all output file names; set to ``None`` to disable 205 [``"helanal__"``] 206 ref_axis : array_like (optional) 207 Calculate tilt angle relative to the axis; if ``None`` then ``[0,0,1]`` 208 is chosen [``None``] 209 verbose : bool (optional) 210 Toggle diagnostic outputs. [``True``] 211 212 Raises 213 ------ 214 FinishTimeException 215 If the specified finish time precedes the specified start time or 216 current time stamp of trajectory object. 217 218 Notes 219 ----- 220 Only a single helix is analyzed. Use the selection to specify the helix, 221 e.g. with "name CA and resid 1:20" or use start=1, stop=20. 222 223 224 .. versionchanged:: 0.13.0 225 New `quiet` keyword to silence frame progress output and most of the 226 output that used to be printed to stdout is now logged to the logger 227 *MDAnalysis.analysis.helanal* (at logelevel *INFO*). 228 229 .. versionchanged:: 0.16.0 230 Removed the `start` and `end` keywords for selecting residues because this can 231 be accomplished more transparently with `selection`. The first and last resid 232 are directly obtained from the selection. 233 234 .. deprecated:: 0.16.0 235 The `quiet` keyword argument is deprecated in favor of the new 236 `verbose` one. 237 238 """ 239 if ref_axis is None: 240 ref_axis = np.array([0., 0., 1.]) 241 else: 242 # enable MDA API so that one can use a tuple of atoms or AtomGroup with 243 # two atoms 244 ref_axis = np.asarray(ref_axis) 245 246 ca = universe.select_atoms(selection) 247 start, end = ca.resids[[0, -1]] 248 trajectory = universe.trajectory 249 250 if finish is not None: 251 if trajectory.ts.time > finish: 252 # you'd be starting with a finish time (in ps) that has already passed or not 253 # available 254 raise FinishTimeException( 255 'The input finish time ({finish} ps) precedes the current trajectory time of {traj_time} ps.'.format( 256 finish=finish, traj_time=trajectory.time)) 257 258 if start is not None and end is not None: 259 logger.info("Analysing from residue %d to %d", start, end) 260 elif start is not None and end is None: 261 logger.info("Analysing from residue %d to the C termini", start) 262 elif start is None and end is not None: 263 logger.info("Analysing from the N termini to %d", end) 264 logger.info("Analysing %d/%d residues", ca.n_atoms, universe.atoms.n_residues) 265 266 if prefix is not None: 267 prefix = str(prefix) 268 matrix_filename = prefix + matrix_filename 269 origin_pdbfile = prefix + origin_pdbfile 270 summary_filename = prefix + summary_filename 271 screw_filename = prefix + screw_filename 272 tilt_filename = prefix + tilt_filename 273 fitted_tilt_filename = prefix + fitted_tilt_filename 274 bend_filename = prefix + bend_filename 275 twist_filename = prefix + twist_filename 276 backup_file(matrix_filename) 277 backup_file(origin_pdbfile) 278 backup_file(summary_filename) 279 backup_file(screw_filename) 280 backup_file(tilt_filename) 281 backup_file(fitted_tilt_filename) 282 backup_file(bend_filename) 283 backup_file(twist_filename) 284 285 global_height = [] 286 global_twist = [] 287 global_rnou = [] 288 global_bending = [] 289 global_bending_matrix = [] 290 global_tilt = [] 291 global_fitted_tilts = [] 292 global_screw = [] 293 294 pm = ProgressMeter(trajectory.n_frames, verbose=verbose, 295 format="Frame %(step)10d: %(time)20.1f ps\r") 296 for ts in trajectory: 297 pm.echo(ts.frame, time=ts.time) 298 frame = ts.frame 299 if begin is not None: 300 if trajectory.time < begin: 301 continue 302 if finish is not None: 303 if trajectory.time > finish: 304 break 305 306 ca_positions = ca.positions 307 twist, bending_angles, height, rnou, origins, local_helix_axes, local_screw_angles = \ 308 main_loop(ca_positions, ref_axis=ref_axis) 309 310 origin_pdb(origins, origin_pdbfile) 311 312 #calculate local bending matrix( it is looking at all i, j combinations) 313 if len(global_bending_matrix) == 0: 314 global_bending_matrix = [[[] for item in local_helix_axes] for item in local_helix_axes] 315 316 for i in range(len(local_helix_axes)): 317 for j in range(i + 1, len(local_helix_axes)): 318 angle = np.rad2deg(np.arccos(np.dot(local_helix_axes[i], local_helix_axes[j]))) 319 global_bending_matrix[i][j].append(angle) 320 #global_bending_matrix[j][i].append(angle) 321 #global_bending_matrix[i][i].append(0.) 322 323 fit_vector, fit_tilt = vector_of_best_fit(origins) 324 global_height += height 325 global_twist += twist 326 global_rnou += rnou 327 #global_screw.append(local_screw_angles) 328 global_fitted_tilts.append(np.rad2deg(fit_tilt)) 329 330 #print out rotations across the helix to a file 331 with open(twist_filename, "a") as twist_output: 332 print(frame, end='', file=twist_output) 333 for loc_twist in twist: 334 print(loc_twist, end='', file=twist_output) 335 print("", file=twist_output) 336 337 with open(bend_filename, "a") as bend_output: 338 print(frame, end='', file=bend_output) 339 for loc_bend in bending_angles: 340 print(loc_bend, end='', file=bend_output) 341 print("", file=bend_output) 342 343 with open(screw_filename, "a") as rot_output: 344 print(frame, end='', file=rot_output) 345 for rotation in local_screw_angles: 346 print(rotation, end='', file=rot_output) 347 print("", file=rot_output) 348 349 with open(tilt_filename, "a") as tilt_output: 350 print(frame, end='', file=tilt_output) 351 for tilt in local_helix_axes: 352 print(np.rad2deg(mdamath.angle(tilt, ref_axis)), 353 end='', file=tilt_output) 354 print("", file=tilt_output) 355 356 with open(fitted_tilt_filename, "a") as tilt_output: 357 print(frame, np.rad2deg(fit_tilt), file=tilt_output) 358 359 if len(global_bending) == 0: 360 global_bending = [[] for item in bending_angles] 361 #global_tilt = [ [] for item in local_helix_axes ] 362 for store, tmp in zip(global_bending, bending_angles): 363 store.append(tmp) 364 #for store,tmp in zip(global_tilt,local_helix_axes): store.append(mdamath.angle(tmp,ref_axis)) 365 366 367 twist_mean, twist_sd, twist_abdev = stats(global_twist) 368 height_mean, height_sd, height_abdev = stats(global_height) 369 rnou_mean, rnou_sd, rnou_abdev = stats(global_rnou) 370 ftilt_mean, ftilt_sd, ftilt_abdev = stats(global_fitted_tilts) 371 372 bending_statistics = [stats(item) for item in global_bending] 373 #tilt_statistics = [ stats(item) for item in global_tilt] 374 375 bending_statistics_matrix = [[stats(col) for col in row] for row in global_bending_matrix] 376 with open(matrix_filename, 'w') as mat_output: 377 print("Mean", file=mat_output) 378 for row in bending_statistics_matrix: 379 for col in row: 380 formatted_angle = "{0:6.1f}".format(col[0]) 381 print(formatted_angle, end='', file=mat_output) 382 print('', file=mat_output) 383 384 print('\nSD', file=mat_output) 385 for row in bending_statistics_matrix: 386 for col in row: 387 formatted_angle = "{0:6.1f}".format(col[1]) 388 print(formatted_angle, end='', file=mat_output) 389 print('', file=mat_output) 390 391 print("\nABDEV", file=mat_output) 392 for row in bending_statistics_matrix: 393 for col in row: 394 formatted_angle = "{0:6.1f}".format(col[2]) 395 print(formatted_angle, end='', file=mat_output) 396 print('', file=mat_output) 397 398 logger.info("Height: %g SD: %g ABDEV: %g (Angstroem)", height_mean, height_sd, height_abdev) 399 logger.info("Twist: %g SD: %g ABDEV: %g", twist_mean, twist_sd, twist_abdev) 400 logger.info("Residues/turn: %g SD: %g ABDEV: %g", rnou_mean, rnou_sd, rnou_abdev) 401 logger.info("Fitted tilt: %g SD: %g ABDEV: %g", ftilt_mean, ftilt_sd, ftilt_abdev) 402 logger.info("Local bending angles:") 403 residue_statistics = list(zip(*bending_statistics)) 404 measure_names = ["Mean ", "SD ", "ABDEV"] 405 if start is None: 406 output = " ".join(["{0:8d}".format(item) 407 for item in range(4, len(residue_statistics[0]) + 4)]) 408 else: 409 output = " ".join(["{0:8d}".format(item) 410 for item in range(start + 3, len(residue_statistics[0]) + start + 3)]) 411 logger.info("ResID %s", output) 412 for measure, name in zip(residue_statistics, measure_names): 413 output = str(name) + " " 414 output += " ".join(["{0:8.1f}".format(residue) for residue in measure]) 415 logger.info(output) 416 417 with open(summary_filename, 'w') as summary_output: 418 print("Height:", height_mean, "SD", height_sd, "ABDEV", height_abdev, '(nm)', file=summary_output) 419 print("Twist:", twist_mean, "SD", twist_sd, "ABDEV", twist_abdev, 420 file=summary_output) 421 print("Residues/turn:", rnou_mean, "SD", rnou_sd, "ABDEV", rnou_abdev, 422 file=summary_output) 423 print("Local bending angles:", file=summary_output) 424 residue_statistics = list(zip(*bending_statistics)) 425 measure_names = ["Mean ", "SD ", "ABDEV"] 426 print("ResID", end='', file=summary_output) 427 if start is None: 428 for item in range(4, len(residue_statistics[0]) + 4): 429 output = "{0:8d}".format(item) 430 print(output, end='', file=summary_output) 431 else: 432 for item in range(start + 3, len(residue_statistics[0]) + start + 3): 433 output = "{0:8d}".format(item) 434 print(output, end='', file=summary_output) 435 print('', file=summary_output) 436 437 for measure, name in zip(residue_statistics, measure_names): 438 print(name, end='', file=summary_output) 439 for residue in measure: 440 output = "{0:8.1f}".format(residue) 441 print(output, end='', file=summary_output) 442 print('', file=summary_output) 443 444 445def tilt_correct(number): 446 """Changes an angle (in degrees) so that it is between 0º and 90º""" 447 if number < 90.: 448 return number 449 else: 450 return 180. - number 451 452 453def backup_file(filename): 454 if os.path.exists(filename): 455 target_name = "#" + filename 456 failure = True 457 if not os.path.exists(target_name): 458 os.rename(filename, target_name) 459 failure = False 460 else: 461 for i in range(20): 462 alt_target_name = target_name + "." + str(i) 463 if os.path.exists(alt_target_name): 464 continue 465 else: 466 os.rename(filename, alt_target_name) 467 failure = False 468 break 469 if failure: 470 raise IOError("Too many backups. Clean up and try again") 471 472 473def stats(some_list): 474 if len(some_list) == 0: 475 return [0, 0, 0] 476 list_mean = np.mean(some_list) 477 list_sd = sample_sd(some_list, list_mean) 478 list_abdev = mean_abs_dev(some_list, list_mean) 479 return [list_mean, list_sd, list_abdev] 480 481 482def helanal_main(pdbfile, selection="name CA", ref_axis=None): 483 """Simple HELANAL_ run on a single frame PDB/GRO. 484 485 Computed data are returned as a dict and also logged at level INFO to the 486 logger *MDAnalysis.analysis.helanal*. A simple way to enable a logger is to 487 use :func:`~MDAnalysis.lib.log.start_logging`. 488 489 Parameters 490 ---------- 491 pdbfile : str 492 filename of the single-frame input file 493 selection : str (optional) 494 selection string, default is "name CA" to select all C-alpha atoms. 495 ref_axis : array_like (optional) 496 Calculate tilt angle relative to the axis; if ``None`` then ``[0,0,1]`` 497 is chosen [``None``] 498 499 Returns 500 ------- 501 result : dict 502 The `result` contains keys 503 * Height: mean, stdev, abs dev 504 * Twist: mean, stdev, abs dev 505 * Residues/turn: mean, stdev, abs dev 506 * Local bending angles: array for computed angles (per residue) 507 * Unit twist angles: array for computed angles (per residue) 508 * Best fit tilt 509 * Rotation angles: local screw angles (per residue) 510 511 512 Notes 513 ----- 514 Only a single helix is analyzed. Use the selection to specify the 515 helix, e.g. with "name CA and resid 1:20". 516 517 518 Example 519 ------- 520 Analyze helix 8 in AdK (PDB 4AKE); the standard logger is started and 521 writes output to the file ``MDAnalysis.log``:: 522 523 MDAnalysis.start_logging() 524 data = MDAnalysis.analysis.helanal_main("4ake_A.pdb", selection="name CA and resnum 161-187") 525 526 527 .. versionchanged:: 0.13.0 528 All output is returned as a dict and logged to the logger 529 *MDAnalysis.analysis.helanal* instead of being printed to stdout. 530 531 .. versionchanged:: 0.16.0 532 Removed the `start` and `end` keywords for selecting residues because this can 533 be accomplished more transparently with `selection`. 534 535 """ 536 537 universe = MDAnalysis.Universe(pdbfile) 538 ca = universe.select_atoms(selection) 539 540 logger.info("Analysing %d/%d residues", ca.n_atoms, universe.atoms.n_residues) 541 542 twist, bending_angles, height, rnou, origins, local_helix_axes, local_screw_angles = \ 543 main_loop(ca.positions, ref_axis=ref_axis) 544 545 #TESTED- origins are correct 546 #print current_origin 547 #print origins 548 549 max_angle = np.max(bending_angles) 550 mean_angle = np.mean(bending_angles) 551 #sd calculated using n-1 to replicate original fortran- assumes a limited sample so uses the sample standard 552 # deviation 553 sd_angle = sample_sd(bending_angles, mean_angle) 554 mean_absolute_deviation_angle = mean_abs_dev(bending_angles, mean_angle) 555 #TESTED- stats correct 556 #print max_angle, mean_angle, sd_angle, mean_absolute_deviation_angle 557 558 #calculate local bending matrix(now it is looking at all i, j combinations) 559 # (not used for helanal_main()) 560# for i in local_helix_axes: 561# for j in local_helix_axes: 562# if (i == j).all(): 563# angle = 0. 564# else: 565# angle = np.rad2deg(np.arccos(np.dot(i, j))) 566# #string_angle = "%6.0f\t" % angle 567# #print string_angle, 568# #print '' 569# #TESTED- local bending matrix! 570 571 #Average helical parameters 572 mean_twist = np.mean(twist) 573 sd_twist = sample_sd(twist, mean_twist) 574 abdev_twist = mean_abs_dev(twist, mean_twist) 575 #TESTED-average twists 576 #print mean_twist, sd_twist, abdev_twist 577 mean_rnou = np.mean(rnou) 578 sd_rnou = sample_sd(rnou, mean_rnou) 579 abdev_rnou = mean_abs_dev(rnou, mean_rnou) 580 #TESTED-average residues per turn 581 #print mean_rnou, sd_rnou, abdev_rnou 582 mean_height = np.mean(height) 583 sd_height = sample_sd(height, mean_height) 584 abdev_height = mean_abs_dev(height, mean_height) 585 #TESTED- average rises 586 587 #calculate best fit vector and tilt of said vector 588 fit_vector, fit_tilt = vector_of_best_fit(origins) 589 590 data = { 591 'Height': np.array([mean_height, sd_height, abdev_height]), 592 'Twist': np.array([mean_twist, sd_twist, abdev_twist]), 593 'Residues/turn': np.array([mean_rnou, sd_rnou, abdev_rnou]), 594 'Local bending angles': np.asarray(bending_angles), 595 'Unit twist angles': np.asarray(twist), 596 'Best fit tilt': fit_tilt, 597 'Rotation Angles': np.asarray(local_screw_angles), 598 } 599 600 logger.info("Height: %g SD: %g ABDEV: %g (Angstroem)", mean_height, sd_height, abdev_height) 601 logger.info("Twist: %g SD: %g ABDEV: %g", mean_twist, sd_twist, abdev_twist) 602 logger.info("Residues/turn: %g SD: %g ABDEV: %g", mean_rnou, sd_rnou, abdev_rnou) 603 604 output = " ".join(["{0:8.1f}\t".format(angle) for angle in bending_angles]) 605 logger.info("Local bending angles: %s", output) 606 607 output = " ".join(["{0:8.1f}\t".format(twist_ang) for twist_ang in twist]) 608 logger.info("Unit twist angles: %s", output) 609 610 logger.info("Best fit tilt: %g", fit_tilt) 611 612 output = " ".join(["{0:.1f}".format(item) for item in local_screw_angles]) 613 logger.info("Rotation Angles from 1 to n-1 (local screw angles): %s", output) 614 615 return data 616 617def origin_pdb(origins, pdbfile): 618 """Write origins to PDB (multi-frame). 619 620 This PDB can be loaded together with the trajectory into, e.g. VMD_, to 621 view the helix axis together with all the atoms. 622 """ 623 with open(pdbfile, 'a') as output: 624 i = 1 625 for xyz in origins: 626 tmp = "ATOM {0:3d} CA ALA {1:3d} {2:8.3f}{3:8.3f}{4:8.3f} 1.00 0.00".format(i, i, xyz[0], xyz[1], xyz[2]) 627 print(tmp, file=output) 628 i += 1 629 print("TER\nENDMDL", file=output) 630 631 632def main_loop(positions, ref_axis=None): 633 # rewrite in cython? 634 635 if ref_axis is None: 636 ref_axis = np.array([0., 0., 1.]) 637 else: 638 ref_axis = np.asarray(ref_axis) 639 twist = [] 640 rnou = [] 641 height = [] 642 origins = [[0., 0., 0.] for item in positions[:-2]] 643 local_helix_axes = [] 644 location_rotation_vectors = [] 645 for i in range(len(positions) - 3): 646 vec12 = positions[i + 1] - positions[i] 647 vec23 = positions[i + 2] - positions[i + 1] 648 vec34 = positions[i + 3] - positions[i + 2] 649 650 dv13 = vec12 - vec23 651 dv24 = vec23 - vec34 652 653 #direction of the local helix axis 654 current_uloc = vecnorm(np.cross(dv13, dv24)) 655 local_helix_axes.append(current_uloc) 656 657 #TESTED- Axes correct 658 #print current_uloc 659 660 dmag = mdamath.norm(dv13) 661 emag = mdamath.norm(dv24) 662 663 costheta = np.dot(dv13, dv24) / (dmag * emag) 664 #rnou is the number of residues per turn 665 current_twist = np.arccos(costheta) 666 twist.append(np.rad2deg(current_twist)) 667 rnou.append(2 * np.pi / current_twist) 668 #radius of local helix cylinder radmag 669 670 costheta1 = 1.0 - costheta 671 radmag = (dmag * emag) ** 0.5 / (2 * costheta1) 672 673 #Height of local helix cylinder 674 current_height = np.dot(vec23, current_uloc) 675 height.append(current_height) 676 #TESTED- Twists etc correct 677 #print current_twist*180/np.pi, 2*np.pi/current_twist, height 678 679 dv13 = vecnorm(dv13) 680 dv24 = vecnorm(dv24) 681 682 #record local rotation 683 location_rotation_vectors.append(dv13) 684 685 rad = [radmag * item for item in dv13] 686 current_origin = [(item[0] - item[1]) for item in zip(positions[i + 1], rad)] 687 origins[i] = current_origin 688 689 #TESTED- origins are correct 690 #print current_origin 691 692 rad = [radmag * item for item in dv24] 693 current_origin = [(item[0] - item[1]) for item in zip(positions[i + 2], rad)] 694 origins[i + 1] = current_origin 695 #Record final rotation vector 696 location_rotation_vectors.append(dv24) 697 698 #local bending angles (eg i > i+3, i+3 > i+6) 699 700 bending_angles = [0 for item in range(len(local_helix_axes) - 3)] 701 for axis in range(len(local_helix_axes) - 3): 702 angle = np.arccos(np.dot(local_helix_axes[axis], local_helix_axes[axis + 3])) 703 bending_angles[axis] = np.rad2deg(angle) 704 #TESTED- angles are correct 705 #print np.rad2deg(angle) 706 707 local_screw_angles = [] 708 #Calculate rotation angles for (+1) to (n-1) 709 fit_vector, fit_tilt = vector_of_best_fit(origins) 710 for item in location_rotation_vectors: 711 local_screw_tmp = np.rad2deg(rotation_angle(fit_vector, ref_axis, item)) 712 #print local_screw_tmp 713 local_screw_angles.append(local_screw_tmp) 714 715 return twist, bending_angles, height, rnou, origins, local_helix_axes, local_screw_angles 716 717 718def rotation_angle(helix_vector, axis_vector, rotation_vector): 719 reference_vector = np.cross(np.cross(helix_vector, axis_vector), helix_vector) 720 second_reference_vector = np.cross(axis_vector, helix_vector) 721 screw_angle = mdamath.angle(reference_vector, rotation_vector) 722 alt_screw_angle = mdamath.angle(second_reference_vector, rotation_vector) 723 updown = np.cross(reference_vector, rotation_vector) 724 725 if not (np.pi < screw_angle < 3 * np.pi / 4): 726 if screw_angle < np.pi / 4 and alt_screw_angle < np.pi / 2: 727 screw_angle = np.pi / 2 - alt_screw_angle 728 elif screw_angle < np.pi / 4 and alt_screw_angle > np.pi / 2: 729 screw_angle = alt_screw_angle - np.pi / 2 730 elif screw_angle > 3 * np.pi / 4 and alt_screw_angle < np.pi / 2: 731 screw_angle = np.pi / 2 + alt_screw_angle 732 elif screw_angle > 3 * np.pi / 4 and alt_screw_angle > np.pi / 2: 733 screw_angle = 3 * np.pi / 2 - alt_screw_angle 734 else: 735 logger.debug("Big Screw Up: screw_angle=%g degrees", np.rad2deg(screw_angle)) 736 737 if mdamath.norm(updown) == 0: 738 logger.warning("PROBLEM (vector is at 0 or 180)") 739 740 helix_dot_rehelix = mdamath.angle(updown, helix_vector) 741 742 #if ( helix_dot_rehelix < np.pi/2 and helix_dot_rehelix >= 0 )or helix_dot_rehelix <-np.pi/2: 743 if (-np.pi / 2 < helix_dot_rehelix < np.pi / 2) or (helix_dot_rehelix > 3 * np.pi / 2): 744 screw_angle = -screw_angle 745 746 return screw_angle 747 748 749def vector_of_best_fit(origins): 750 origins = np.asarray(origins) 751 centroids = center(origins) 752 M = origins - centroids 753 A = np.dot(M.transpose(), M) 754 u, s, vh = np.linalg.linalg.svd(A) 755 vector = vh[0] 756 #Correct vector to face towards first residues 757 rough_helix = origins[0] - centroids 758 agreement = mdamath.angle(rough_helix, vector) 759 if not (-np.pi / 2 < agreement < np.pi / 2): 760 vector *= -1 761 best_fit_tilt = mdamath.angle(vector, [0, 0, 1]) 762 return vector, best_fit_tilt 763