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 23""" 24Nucleic acid analysis --- :mod:`MDAnalysis.analysis.nuclinfo` 25============================================================= 26 27:Author: Elizabeth Denning 28:Year: 2011 29:Copyright: GNU Public License v3 30 31The module provides functions to analyze nuclic acid structures, in 32particular 33 34- backbone dihedrals, 35- chi dihedrals, 36- AS or CP phase angles, 37- Watson-Crick N1-N3 distances, C2-O2 distances, N6-O4 distances, O6-N4 distances. 38 39For applications of this kind of analysis see [Denning2011]_ and [Denning2012]_. 40 41All functions take a :class:`~MDAnalysis.core.universe.Universe` as an 42argument together with further parameters that specify the base or bases in 43question. Angles are in degrees. The functions use standard CHARMM names for 44nucleic acids and atom names. 45 46 47.. rubric:: References 48 49.. [Denning2011] E.J. Denning, U.D. Priyakumar, L. Nilsson, and A.D. Mackerell, Jr. Impact of 50 2'-hydroxyl sampling on the conformational properties of RNA: update of the 51 CHARMM all-atom additive force field for RNA. *J. Comput. Chem.* 32 (2011), 52 1929--1943. doi: `10.1002/jcc.21777`_ 53 54.. [Denning2012] E.J. Denning and A.D. MacKerell, Jr. Intrinsic Contribution of the 2'-Hydroxyl to 55 RNA Conformational Heterogeneity. *J. Am. Chem. Soc.* 134 (2012), 2800--2806. 56 doi: `10.1021/ja211328g`_ 57 58 59.. _`10.1002/jcc.21777`: http://dx.doi.org/10.1002/jcc.21777 60.. _`10.1021/ja211328g`: http://dx.doi.org/10.1021/ja211328g 61 62 63Distances 64--------- 65 66.. autofunction:: wc_pair 67 68.. autofunction:: minor_pair 69 70.. autofunction:: major_pair 71 72 73Phases 74------ 75 76.. autofunction:: phase_cp 77 78.. autofunction:: phase_as 79 80 81Dihedral angles 82--------------- 83 84.. autofunction:: tors 85 86.. autofunction:: tors_alpha 87 88.. autofunction:: tors_beta 89 90.. autofunction:: tors_gamma 91 92.. autofunction:: tors_delta 93 94.. autofunction:: tors_eps 95 96.. autofunction:: tors_zeta 97 98.. autofunction:: tors_chi 99 100.. autofunction:: hydroxyl 101 102.. autofunction:: pseudo_dihe_baseflip 103 104""" 105from __future__ import division, absolute_import 106 107import numpy as np 108from math import pi, sin, cos, atan2, sqrt, pow 109 110from MDAnalysis.lib import mdamath 111 112 113def wc_pair(universe, i, bp, seg1="SYSTEM", seg2="SYSTEM"): 114 """Watson-Crick basepair distance for residue `i` with residue `bp`. 115 116 The distance of the nitrogen atoms in a Watson-Crick hydrogen bond is 117 computed. 118 119 Parameters 120 ---------- 121 universe : Universe 122 :class:`~MDAnalysis.core.universe.Universe` containing the trajectory 123 i : int 124 resid of the first base 125 bp : int 126 resid of the second base 127 seg1 : str (optional) 128 segment id for first base ["SYSTEM"] 129 seg2 : str (optional) 130 segment id for second base ["SYSTEM"] 131 132 Returns 133 ------- 134 float 135 Watson-Crick base pair distance 136 137 138 Notes 139 ----- 140 If failure occurs be sure to check the segment identification. 141 142 143 .. versionadded:: 0.7.6 144 """ 145 if universe.select_atoms(" resid {0!s} ".format(i)).resnames[0] in ["DC", "DT", "U", "C", "T", "CYT", "THY", "URA"]: 146 a1, a2 = "N3", "N1" 147 if universe.select_atoms(" resid {0!s} ".format(i)).resnames[0] in ["DG", "DA", "A", "G", "ADE", "GUA"]: 148 a1, a2 = "N1", "N3" 149 wc_dist = universe.select_atoms("(segid {0!s} and resid {1!s} and name {2!s}) " 150 "or (segid {3!s} and resid {4!s} and name {5!s}) " 151 .format(seg1, i, a1, seg2, bp, a2)) 152 wc = mdamath.norm(wc_dist[0].position - wc_dist[1].position) 153 return wc 154 155 156def minor_pair(universe, i, bp, seg1="SYSTEM", seg2="SYSTEM"): 157 """Minor-Groove basepair distance for residue `i` with residue `bp`. 158 159 The distance of the nitrogen and oxygen atoms in a Minor-groove hydrogen 160 bond is computed. 161 162 Parameters 163 ---------- 164 universe : Universe 165 :class:`~MDAnalysis.core.universe.Universe` containing the trajectory 166 i : int 167 resid of the first base 168 bp : int 169 resid of the second base 170 seg1 : str (optional) 171 segment id for first base ["SYSTEM"] 172 seg2 : str (optional) 173 segment id for second base ["SYSTEM"] 174 175 Returns 176 ------- 177 float 178 Minor groove base pair distance 179 180 Notes 181 ----- 182 If failure occurs be sure to check the segment identification. 183 184 185 .. versionadded:: 0.7.6 186 """ 187 if universe.select_atoms(" resid {0!s} ".format(i)).resnames[0] in ["DC", "DT", "U", "C", "T", "CYT", "THY", "URA"]: 188 a1, a2 = "O2", "C2" 189 if universe.select_atoms(" resid {0!s} ".format(i)).resnames[0] in ["DG", "DA", "A", "G", "ADE", "GUA"]: 190 a1, a2 = "C2", "O2" 191 c2o2_dist = universe.select_atoms("(segid {0!s} and resid {1!s} and name {2!s}) " 192 "or (segid {3!s} and resid {4!s} and name {5!s})" 193 .format(seg1, i, a1, seg2, bp, a2)) 194 c2o2 = mdamath.norm(c2o2_dist[0].position - c2o2_dist[1].position) 195 return c2o2 196 197 198def major_pair(universe, i, bp, seg1="SYSTEM", seg2="SYSTEM"): 199 """Major-Groove basepair distance for residue `i` with residue `bp`. 200 201 The distance of the nitrogen and oxygen atoms in a Major-groove hydrogen 202 bond is computed. 203 204 205 Parameters 206 ---------- 207 universe : Universe 208 :class:`~MDAnalysis.core.universe.Universe` containing the trajectory 209 i : int 210 resid of the first base 211 bp : int 212 resid of the second base 213 seg1 : str (optional) 214 segment id for first base ["SYSTEM"] 215 seg2 : str (optional) 216 segment id for second base ["SYSTEM"] 217 218 Returns 219 ------- 220 float 221 Major groove base pair distance 222 223 Notes 224 ----- 225 If failure occurs be sure to check the segment identification. 226 227 228 .. versionadded:: 0.7.6 229 """ 230 if universe.select_atoms(" resid {0!s} ".format(i)).resnames[0] in ["DC", "DG", "C", "G", "CYT", "GUA"]: 231 if universe.select_atoms(" resid {0!s} ".format(i)).resnames[0] in ["DC", "C", "CYT"]: 232 a1, a2 = "N4", "O6" 233 else: 234 a1, a2 = "O6", "N4" 235 if universe.select_atoms(" resid {0!s} ".format(i)).resnames[0] in ["DT", "DA", "A", "T", "U", "ADE", "THY", "URA"]: 236 if universe.select_atoms(" resid {0!s} ".format(i)).resnames[0] in ["DT", "T", "THY", "U", "URA"]: 237 a1, a2 = "O4", "N6" 238 else: 239 a1, a2 = "N6", "O4" 240 no_dist = universe.select_atoms("(segid {0!s} and resid {1!s} and name {2!s}) " 241 "or (segid {3!s} and resid {4!s} and name {5!s}) " 242 .format(seg1, i, a1, seg2, bp, a2)) 243 major = mdamath.norm(no_dist[0].position - no_dist[1].position) 244 return major 245 246 247def phase_cp(universe, seg, i): 248 """Pseudo-angle describing the phase of the ribose pucker for residue `i` using the CP method. 249 250 The angle is computed by the positions of atoms in the ribose ring. 251 252 253 Parameters 254 ---------- 255 universe : Universe 256 :class:`~MDAnalysis.core.universe.Universe` containing the trajectory 257 seg : str 258 segment id for base 259 i : int 260 resid of the first base 261 262 Returns 263 ------- 264 float 265 phase angle in degrees 266 267 268 .. versionadded:: 0.7.6 269 """ 270 atom1 = universe.select_atoms(" atom {0!s} {1!s} O4\' ".format(seg, i)) 271 atom2 = universe.select_atoms(" atom {0!s} {1!s} C1\' ".format(seg, i)) 272 atom3 = universe.select_atoms(" atom {0!s} {1!s} C2\' ".format(seg, i)) 273 atom4 = universe.select_atoms(" atom {0!s} {1!s} C3\' ".format(seg, i)) 274 atom5 = universe.select_atoms(" atom {0!s} {1!s} C4\' ".format(seg, i)) 275 276 data1 = atom1.positions 277 data2 = atom2.positions 278 data3 = atom3.positions 279 data4 = atom4.positions 280 data5 = atom5.positions 281 282 r0 = (data1 + data2 + data3 + data4 + data5) * (1.0 / 5.0) 283 r1 = data1 - r0 284 r2 = data2 - r0 285 r3 = data3 - r0 286 r4 = data4 - r0 287 r5 = data5 - r0 288 289 R1 = ((r1 * sin(2 * pi * 0.0 / 5.0)) + (r2 * sin(2 * pi * 1.0 / 5.0)) + 290 (r3 * sin(2 * pi * 2.0 / 5.0)) + (r4 * sin(2 * pi * 3.0 / 5.0)) + 291 (r5 * sin(2 * pi * 4.0 / 5.0))) 292 293 R2 = ((r1 * cos(2 * pi * 0.0 / 5.0)) + (r2 * cos(2 * pi * 1.0 / 5.0)) + 294 (r3 * cos(2 * pi * 2.0 / 5.0)) + (r4 * cos(2 * pi * 3.0 / 5.0)) + 295 (r5 * cos(2 * pi * 4.0 / 5.0))) 296 297 x = np.cross(R1[0], R2[0]) 298 n = x / sqrt(pow(x[0], 2) + pow(x[1], 2) + pow(x[2], 2)) 299 300 r1_d = np.dot(r1, n) 301 r2_d = np.dot(r2, n) 302 r3_d = np.dot(r3, n) 303 r4_d = np.dot(r4, n) 304 r5_d = np.dot(r5, n) 305 306 D = ((r1_d * sin(4 * pi * 0.0 / 5.0)) + (r2_d * sin(4 * pi * 1.0 / 5.0)) 307 + (r3_d * sin(4 * pi * 2.0 / 5.0)) + (r4_d * sin(4 * pi * 3.0 / 5.0)) 308 + (r5_d * sin(4 * pi * 4.0 / 5.0))) * -1 * sqrt(2.0 / 5.0) 309 310 C = ((r1_d * cos(4 * pi * 0.0 / 5.0)) + (r2_d * cos(4 * pi * 1.0 / 5.0)) 311 + (r3_d * cos(4 * pi * 2.0 / 5.0)) + (r4_d * cos(4 * pi * 3.0 / 5.0)) 312 + (r5_d * cos(4 * pi * 4.0 / 5.0))) * sqrt(2.0 / 5.0) 313 314 phase_ang = (atan2(D, C) + (pi / 2.)) * 180. / pi 315 return phase_ang % 360 316 317 318def phase_as(universe, seg, i): 319 """Pseudo-angle describing the phase of the ribose pucker for residue `i` using the AS method 320 321 The angle is computed by the position vector of atoms in the ribose ring. 322 323 Parameters 324 ---------- 325 universe : Universe 326 :class:`~MDAnalysis.core.universe.Universe` containing the trajectory 327 seg : str 328 segment id for base 329 i : int 330 resid of the first base 331 332 Returns 333 ------- 334 float 335 phase angle in degrees 336 337 338 .. versionadded:: 0.7.6 339 """ 340 angle1 = universe.select_atoms(" atom {0!s} {1!s} C1\' ".format(seg, i), 341 " atom {0!s} {1!s} C2\' ".format(seg, i), 342 " atom {0!s} {1!s} C3\' ".format(seg, i), 343 " atom {0!s} {1!s} C4\' ".format(seg, i)) 344 345 angle2 = universe.select_atoms(" atom {0!s} {1!s} C2\' ".format(seg, i), 346 " atom {0!s} {1!s} C3\' ".format(seg, i), 347 " atom {0!s} {1!s} C4\' ".format(seg, i), 348 " atom {0!s} {1!s} O4\' ".format(seg, i)) 349 350 angle3 = universe.select_atoms(" atom {0!s} {1!s} C3\' ".format(seg, i), 351 " atom {0!s} {1!s} C4\' ".format(seg, i), 352 " atom {0!s} {1!s} O4\' ".format(seg, i), 353 " atom {0!s} {1!s} C1\' ".format(seg, i)) 354 355 angle4 = universe.select_atoms(" atom {0!s} {1!s} C4\' ".format(seg, i), 356 " atom {0!s} {1!s} O4\' ".format(seg, i), 357 " atom {0!s} {1!s} C1\' ".format(seg, i), 358 " atom {0!s} {1!s} C2\' ".format(seg, i)) 359 360 angle5 = universe.select_atoms(" atom {0!s} {1!s} O4\' ".format(seg, i), 361 " atom {0!s} {1!s} C1\' ".format(seg, i), 362 " atom {0!s} {1!s} C2\' ".format(seg, i), 363 " atom {0!s} {1!s} C3\' ".format(seg, i)) 364 365 data1 = angle1.dihedral.value() 366 data2 = angle2.dihedral.value() 367 data3 = angle3.dihedral.value() 368 data4 = angle4.dihedral.value() 369 data5 = angle5.dihedral.value() 370 371 B = ((data1 * sin(2 * 2 * pi * (1 - 1.) / 5.)) 372 + (data2 * sin(2 * 2 * pi * (2 - 1.) / 5.)) 373 + (data3 * sin(2 * 2 * pi * (3 - 1.) / 5.)) 374 + (data4 * sin(2 * 2 * pi * (4 - 1.) / 5.)) 375 + (data5 * sin(2 * 2 * pi * (5 - 1.) / 5.))) * -2. / 5. 376 377 A = ((data1 * cos(2 * 2 * pi * (1 - 1.) / 5.)) 378 + (data2 * cos(2 * 2 * pi * (2 - 1.) / 5.)) 379 + (data3 * cos(2 * 2 * pi * (3 - 1.) / 5.)) 380 + (data4 * cos(2 * 2 * pi * (4 - 1.) / 5.)) 381 + (data5 * cos(2 * 2 * pi * (5 - 1.) / 5.))) * 2. / 5. 382 383 phase_ang = atan2(B, A) * 180. / pi 384 return phase_ang % 360 385 386 387def tors(universe, seg, i): 388 """Calculation of nucleic backbone dihedral angles. 389 390 The dihedral angles are alpha, beta, gamma, delta, epsilon, zeta, chi. 391 392 The dihedral is computed based on position of atoms for resid `i`. 393 394 Parameters 395 ---------- 396 universe : Universe 397 :class:`~MDAnalysis.core.universe.Universe` containing the trajectory 398 seg : str 399 segment id for base 400 i : int 401 resid of the first base 402 403 Returns 404 ------- 405 [alpha, beta, gamma, delta, epsilon, zeta, chi] : list of floats 406 torsion angles in degrees 407 408 Notes 409 ----- 410 If failure occurs be sure to check the segment identification. 411 412 413 .. versionadded:: 0.7.6 414 415 """ 416 a = universe.select_atoms(" atom {0!s} {1!s} O3\' ".format(seg, i - 1), 417 " atom {0!s} {1!s} P ".format(seg, i), 418 " atom {0!s} {1!s} O5\' ".format(seg, i), 419 " atom {0!s} {1!s} C5\' ".format(seg, i)) 420 421 b = universe.select_atoms(" atom {0!s} {1!s} P ".format(seg, i), 422 " atom {0!s} {1!s} O5\' ".format(seg, i), 423 " atom {0!s} {1!s} C5\' ".format(seg, i), 424 " atom {0!s} {1!s} C4\' ".format(seg, i)) 425 426 g = universe.select_atoms(" atom {0!s} {1!s} O5\' ".format(seg, i), 427 " atom {0!s} {1!s} C5\' ".format(seg, i), 428 " atom {0!s} {1!s} C4\' ".format(seg, i), 429 " atom {0!s} {1!s} C3\' ".format(seg, i)) 430 431 d = universe.select_atoms(" atom {0!s} {1!s} C5\' ".format(seg, i), 432 " atom {0!s} {1!s} C4\' ".format(seg, i), 433 " atom {0!s} {1!s} C3\' ".format(seg, i), 434 " atom {0!s} {1!s} O3\' ".format(seg, i)) 435 436 e = universe.select_atoms(" atom {0!s} {1!s} C4\' ".format(seg, i), 437 " atom {0!s} {1!s} C3\' ".format(seg, i), 438 " atom {0!s} {1!s} O3\' ".format(seg, i), 439 " atom {0!s} {1!s} P ".format(seg, i + 1)) 440 441 z = universe.select_atoms(" atom {0!s} {1!s} C3\' ".format(seg, i), 442 " atom {0!s} {1!s} O3\' ".format(seg, i), 443 " atom {0!s} {1!s} P ".format(seg, i + 1), 444 " atom {0!s} {1!s} O5\' ".format(seg, i + 1)) 445 c = universe.select_atoms(" atom {0!s} {1!s} O4\' ".format(seg, i), 446 " atom {0!s} {1!s} C1\' ".format(seg, i), 447 " atom {0!s} {1!s} N1 ".format(seg, i), 448 " atom {0!s} {1!s} C2 ".format(seg, i)) 449 if len(c) < 4: 450 c = universe.select_atoms(" atom {0!s} {1!s} O4\' ".format(seg, i), 451 " atom {0!s} {1!s} C1\' ".format(seg, i), 452 " atom {0!s} {1!s} N9 ".format(seg, i), 453 " atom {0!s} {1!s} C4 ".format(seg, i)) 454 455 alpha = a.dihedral.value() % 360 456 beta = b.dihedral.value() % 360 457 gamma = g.dihedral.value() % 360 458 delta = d.dihedral.value() % 360 459 epsilon = e.dihedral.value() % 360 460 zeta = z.dihedral.value() % 360 461 chi = c.dihedral.value() % 360 462 463 return [alpha, beta, gamma, delta, epsilon, zeta, chi] 464 465 466def tors_alpha(universe, seg, i): 467 """alpha backbone dihedral 468 469 The dihedral is computed based on position atoms for resid `i`. 470 471 Parameters 472 ---------- 473 universe : Universe 474 :class:`~MDAnalysis.core.universe.Universe` containing the trajectory 475 seg : str 476 segment id for base 477 i : int 478 resid of the first base 479 480 Returns 481 ------- 482 alpha : float 483 torsion angle in degrees 484 485 486 .. versionadded:: 0.7.6 487 """ 488 a = universe.select_atoms(" atom {0!s} {1!s} O3\' ".format(seg, i - 1), 489 " atom {0!s} {1!s} P ".format(seg, i), 490 " atom {0!s} {1!s} O5\' ".format(seg, i), 491 " atom {0!s} {1!s} C5\' ".format(seg, i)) 492 alpha = a.dihedral.value() % 360 493 return alpha 494 495 496def tors_beta(universe, seg, i): 497 """beta backbone dihedral 498 499 The dihedral is computed based on position atoms for resid `i`. 500 501 Parameters 502 ---------- 503 universe : Universe 504 :class:`~MDAnalysis.core.universe.Universe` containing the trajectory 505 seg : str 506 segment id for base 507 i : int 508 resid of the first base 509 510 Returns 511 ------- 512 beta : float 513 torsion angle in degrees 514 515 516 .. versionadded:: 0.7.6 517 """ 518 b = universe.select_atoms(" atom {0!s} {1!s} P ".format(seg, i), 519 " atom {0!s} {1!s} O5\' ".format(seg, i), 520 " atom {0!s} {1!s} C5\' ".format(seg, i), 521 " atom {0!s} {1!s} C4\' ".format(seg, i)) 522 beta = b.dihedral.value() % 360 523 return beta 524 525 526def tors_gamma(universe, seg, i): 527 """ Gamma backbone dihedral 528 529 The dihedral is computed based on position atoms for resid `i`. 530 531 Parameters 532 ---------- 533 universe : Universe 534 :class:`~MDAnalysis.core.universe.Universe` containing the trajectory 535 seg : str 536 segment id for base 537 i : int 538 resid of the first base 539 540 Returns 541 ------- 542 gamma : float 543 torsion angle in degrees 544 545 546 .. versionadded:: 0.7.6 547 """ 548 g = universe.select_atoms(" atom {0!s} {1!s} O5\' ".format(seg, i), 549 " atom {0!s} {1!s} C5\' ".format(seg, i), 550 " atom {0!s} {1!s} C4\' ".format(seg, i), 551 " atom {0!s} {1!s} C3\' ".format(seg, i)) 552 gamma = g.dihedral.value() % 360 553 return gamma 554 555 556def tors_delta(universe, seg, i): 557 """delta backbone dihedral 558 559 The dihedral is computed based on position atoms for resid `i`. 560 561 Parameters 562 ---------- 563 universe : Universe 564 :class:`~MDAnalysis.core.universe.Universe` containing the trajectory 565 seg : str 566 segment id for base 567 i : int 568 resid of the first base 569 570 Returns 571 ------- 572 delta : float 573 torsion angle in degrees 574 575 576 .. versionadded:: 0.7.6 577 """ 578 d = universe.select_atoms(" atom {0!s} {1!s} C5\' ".format(seg, i), 579 " atom {0!s} {1!s} C4\' ".format(seg, i), 580 " atom {0!s} {1!s} C3\' ".format(seg, i), 581 " atom {0!s} {1!s} O3\' ".format(seg, i)) 582 delta = d.dihedral.value() % 360 583 return delta 584 585 586def tors_eps(universe, seg, i): 587 """Epsilon backbone dihedral 588 589 The dihedral is computed based on position atoms for resid `i`. 590 591 Parameters 592 ---------- 593 universe : Universe 594 :class:`~MDAnalysis.core.universe.Universe` containing the trajectory 595 seg : str 596 segment id for base 597 i : int 598 resid of the first base 599 600 Returns 601 ------- 602 epsilon : float 603 torsion angle in degrees 604 605 606 .. versionadded:: 0.7.6 607 """ 608 e = universe.select_atoms(" atom {0!s} {1!s} C4\' ".format(seg, i), 609 " atom {0!s} {1!s} C3\' ".format(seg, i), 610 " atom {0!s} {1!s} O3\' ".format(seg, i), 611 " atom {0!s} {1!s} P ".format(seg, i + 1)) 612 epsilon = e.dihedral.value() % 360 613 return epsilon 614 615 616def tors_zeta(universe, seg, i): 617 """Zeta backbone dihedral 618 619 The dihedral is computed based on position atoms for resid `i`. 620 621 Parameters 622 ---------- 623 universe : Universe 624 :class:`~MDAnalysis.core.universe.Universe` containing the trajectory 625 seg : str 626 segment id for base 627 i : int 628 resid of the first base 629 630 Returns 631 ------- 632 zeta : float 633 torsion angle in degrees 634 635 636 .. versionadded:: 0.7.6 637 """ 638 z = universe.select_atoms(" atom {0!s} {1!s} C3\' ".format(seg, i), 639 " atom {0!s} {1!s} O3\' ".format(seg, i), 640 " atom {0!s} {1!s} P ".format(seg, i + 1), 641 " atom {0!s} {1!s} O5\' ".format(seg, i + 1)) 642 zeta = z.dihedral.value() % 360 643 return zeta 644 645 646def tors_chi(universe, seg, i): 647 """chi nucleic acid dihedral 648 649 The dihedral is computed based on position atoms for resid `i`. 650 651 Parameters 652 ---------- 653 universe : Universe 654 :class:`~MDAnalysis.core.universe.Universe` containing the trajectory 655 seg : str 656 segment id for base 657 i : int 658 resid of the first base 659 660 Returns 661 ------- 662 chi : float 663 torsion angle in degrees 664 665 666 .. versionadded:: 0.7.6 667 """ 668 c = universe.select_atoms(" atom {0!s} {1!s} O4\' ".format(seg, i), 669 " atom {0!s} {1!s} C1\' ".format(seg, i), 670 " atom {0!s} {1!s} N1 ".format(seg, i), 671 " atom {0!s} {1!s} C2 ".format(seg, i)) 672 if len(c) < 4: 673 c = universe.select_atoms(" atom {0!s} {1!s} O4\' ".format(seg, i), 674 " atom {0!s} {1!s} C1\' ".format(seg, i), 675 " atom {0!s} {1!s} N9 ".format(seg, i), 676 " atom {0!s} {1!s} C4 ".format(seg, i)) 677 chi = c.dihedral.value() % 360 678 return chi 679 680 681def hydroxyl(universe, seg, i): 682 """2-hydroxyl dihedral. Useful only for RNA calculations. 683 684 .. Note:: This dihedral calculation will only work if using atom 685 names as documented by charmm force field parameters, 686 namely "C1', C2', O2', H2'". 687 688 Parameters 689 ---------- 690 universe : Universe 691 :class:`~MDAnalysis.core.universe.Universe` containing the trajectory 692 seg : str 693 segment id for base 694 i : int 695 resid of the first base 696 697 Returns 698 ------- 699 hydroxyl_angle : float 700 torsion angle in degrees 701 702 703 .. versionadded:: 0.7.6 704 705 """ 706 h = universe.select_atoms(" atom {0!s} {1!s} C1\' ".format(seg, i), 707 " atom {0!s} {1!s} C2\' ".format(seg, i), 708 " atom {0!s} {1!s} O2\' ".format(seg, i), 709 " atom {0!s} {1!s} H2\'\' ".format(seg, i)) 710 try: 711 hydr = h.dihedral.value() % 360 712 except ValueError: 713 raise ValueError("Resid {0} does not contain atoms C1', C2', O2', H2' but atoms {1}" 714 .format(i, str(list(h.atoms)))) 715 return hydr 716 717 718def pseudo_dihe_baseflip(universe, bp1, bp2, i, 719 seg1="SYSTEM", seg2="SYSTEM", seg3="SYSTEM"): 720 """pseudo dihedral for flipped bases. Useful only for nucleic acid base flipping 721 722 The dihedral is computed based on position atoms for resid `i` 723 724 .. Note:: This dihedral calculation will only work if using atom names as 725 documented by charmm force field parameters. 726 727 Parameters 728 ---------- 729 universe : Universe 730 :class:`~MDAnalysis.core.universe.Universe` containing the 731 trajectory 732 bp1 : int 733 resid that base pairs with `bp2` 734 bp2 : int 735 resid below the base that flips 736 i : int 737 resid of the base that flips 738 segid1 : str (optional) 739 segid of resid base pairing with `bp2` 740 segid2 : str (optional) 741 segid, same as that of segid of flipping resid `i` 742 segid3 : str (optional) 743 segid of resid `i` that flips 744 745 Returns 746 ------- 747 float 748 pseudo dihedral angle in degrees 749 750 751 .. versionadded:: 0.8.0 752 """ 753 bf1 = universe.select_atoms( 754 " ( segid {0!s} and resid {1!s} and nucleicbase ) " 755 "or ( segid {2!s} and resid {3!s} and nucleicbase ) " 756 .format( seg1, bp1, seg2, bp2)) 757 bf4 = universe.select_atoms("(segid {0!s} and resid {1!s} and nucleicbase) ".format(seg3, i)) 758 bf2 = universe.select_atoms("(segid {0!s} and resid {1!s} and nucleicsugar) ".format(seg2, bp2)) 759 bf3 = universe.select_atoms("(segid {0!s} and resid {1!s} and nucleicsugar) ".format(seg3, i)) 760 x = [bf1.center_of_mass(), bf2.center_of_mass(), 761 bf3.center_of_mass(), bf4.center_of_mass()] 762 pseudo = mdamath.dihedral(x[0] - x[1], x[1] - x[2], x[2] - x[3]) 763 pseudo = np.rad2deg(pseudo) % 360 764 return pseudo 765