# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*-
# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
#
# MDAnalysis --- https://www.mdanalysis.org
# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors
# (see the file AUTHORS for the full list of names)
#
# Released under the GNU Public Licence, v2 or any higher version
#
# Please cite your use of MDAnalysis in published work:
#
# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler,
# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein.
# MDAnalysis: A Python package for the rapid analysis of molecular dynamics
# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th
# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy.
# doi: 10.25080/majora-629e541a-00e
#
# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein.
# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations.
# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
#
"""
Nucleic acid analysis --- :mod:`MDAnalysis.analysis.nuclinfo`
=============================================================
:Author: Elizabeth Denning
:Year: 2011
:Copyright: GNU Public License v3
The module provides functions to analyze nucleic acid structures, in
particular
- backbone dihedrals,
- chi dihedrals,
- AS or CP phase angles,
- Watson-Crick N1-N3 distances, C2-O2 distances, N6-O4 distances, O6-N4 distances.
For applications of this kind of analysis see
:footcite:p:`Denning2011,Denning2012`.
All functions take a :class:`~MDAnalysis.core.universe.Universe` as an
argument together with further parameters that specify the base or bases in
question. Angles are in degrees. The functions use standard CHARMM names for
nucleic acids and atom names.
.. rubric:: References
.. footbibliography::
Distances
---------
.. autofunction:: wc_pair
.. autofunction:: minor_pair
.. autofunction:: major_pair
Phases
------
.. autofunction:: phase_cp
.. autofunction:: phase_as
Dihedral angles
---------------
.. autofunction:: tors
.. autofunction:: tors_alpha
.. autofunction:: tors_beta
.. autofunction:: tors_gamma
.. autofunction:: tors_delta
.. autofunction:: tors_eps
.. autofunction:: tors_zeta
.. autofunction:: tors_chi
.. autofunction:: hydroxyl
.. autofunction:: pseudo_dihe_baseflip
"""
import numpy as np
from math import pi, sin, cos, atan2, sqrt, pow
from MDAnalysis.lib import mdamath
[docs]
def wc_pair(universe, i, bp, seg1="SYSTEM", seg2="SYSTEM"):
"""Watson-Crick basepair distance for residue `i` with residue `bp`.
The distance of the nitrogen atoms in a Watson-Crick hydrogen bond is
computed.
Parameters
----------
universe : Universe
:class:`~MDAnalysis.core.universe.Universe` containing the trajectory
i : int
resid of the first base
bp : int
resid of the second base
seg1 : str (optional)
segment id for first base ["SYSTEM"]
seg2 : str (optional)
segment id for second base ["SYSTEM"]
Returns
-------
float
Watson-Crick base pair distance
Notes
-----
If failure occurs be sure to check the segment identification.
.. versionadded:: 0.7.6
"""
if universe.select_atoms(" resid {0!s} ".format(i)).resnames[0] in ["DC", "DT", "U", "C", "T", "CYT", "THY", "URA"]:
a1, a2 = "N3", "N1"
if universe.select_atoms(" resid {0!s} ".format(i)).resnames[0] in ["DG", "DA", "A", "G", "ADE", "GUA"]:
a1, a2 = "N1", "N3"
wc_dist = universe.select_atoms("(segid {0!s} and resid {1!s} and name {2!s}) "
"or (segid {3!s} and resid {4!s} and name {5!s}) "
.format(seg1, i, a1, seg2, bp, a2))
wc = mdamath.norm(wc_dist[0].position - wc_dist[1].position)
return wc
[docs]
def minor_pair(universe, i, bp, seg1="SYSTEM", seg2="SYSTEM"):
"""Minor-Groove basepair distance for residue `i` with residue `bp`.
The distance of the nitrogen and oxygen atoms in a Minor-groove hydrogen
bond is computed.
Parameters
----------
universe : Universe
:class:`~MDAnalysis.core.universe.Universe` containing the trajectory
i : int
resid of the first base
bp : int
resid of the second base
seg1 : str (optional)
segment id for first base ["SYSTEM"]
seg2 : str (optional)
segment id for second base ["SYSTEM"]
Returns
-------
float
Minor groove base pair distance
Notes
-----
If failure occurs be sure to check the segment identification.
.. versionadded:: 0.7.6
"""
if universe.select_atoms(" resid {0!s} ".format(i)).resnames[0] in ["DC", "DT", "U", "C", "T", "CYT", "THY", "URA"]:
a1, a2 = "O2", "C2"
if universe.select_atoms(" resid {0!s} ".format(i)).resnames[0] in ["DG", "DA", "A", "G", "ADE", "GUA"]:
a1, a2 = "C2", "O2"
c2o2_dist = universe.select_atoms("(segid {0!s} and resid {1!s} and name {2!s}) "
"or (segid {3!s} and resid {4!s} and name {5!s})"
.format(seg1, i, a1, seg2, bp, a2))
c2o2 = mdamath.norm(c2o2_dist[0].position - c2o2_dist[1].position)
return c2o2
[docs]
def major_pair(universe, i, bp, seg1="SYSTEM", seg2="SYSTEM"):
"""Major-Groove basepair distance for residue `i` with residue `bp`.
The distance of the nitrogen and oxygen atoms in a Major-groove hydrogen
bond is computed.
Parameters
----------
universe : Universe
:class:`~MDAnalysis.core.universe.Universe` containing the trajectory
i : int
resid of the first base
bp : int
resid of the second base
seg1 : str (optional)
segment id for first base ["SYSTEM"]
seg2 : str (optional)
segment id for second base ["SYSTEM"]
Returns
-------
float
Major groove base pair distance
Notes
-----
If failure occurs be sure to check the segment identification.
.. versionadded:: 0.7.6
"""
if universe.select_atoms(" resid {0!s} ".format(i)).resnames[0] in ["DC", "DG", "C", "G", "CYT", "GUA"]:
if universe.select_atoms(" resid {0!s} ".format(i)).resnames[0] in ["DC", "C", "CYT"]:
a1, a2 = "N4", "O6"
else:
a1, a2 = "O6", "N4"
if universe.select_atoms(" resid {0!s} ".format(i)).resnames[0] in ["DT", "DA", "A", "T", "U", "ADE", "THY", "URA"]:
if universe.select_atoms(" resid {0!s} ".format(i)).resnames[0] in ["DT", "T", "THY", "U", "URA"]:
a1, a2 = "O4", "N6"
else:
a1, a2 = "N6", "O4"
no_dist = universe.select_atoms("(segid {0!s} and resid {1!s} and name {2!s}) "
"or (segid {3!s} and resid {4!s} and name {5!s}) "
.format(seg1, i, a1, seg2, bp, a2))
major = mdamath.norm(no_dist[0].position - no_dist[1].position)
return major
[docs]
def phase_cp(universe, seg, i):
"""Pseudo-angle describing the phase of the ribose pucker for residue `i` using the CP method.
The angle is computed by the positions of atoms in the ribose ring.
Parameters
----------
universe : Universe
:class:`~MDAnalysis.core.universe.Universe` containing the trajectory
seg : str
segment id for base
i : int
resid of the first base
Returns
-------
float
phase angle in degrees
.. versionadded:: 0.7.6
"""
atom1 = universe.select_atoms(" atom {0!s} {1!s} O4\' ".format(seg, i))
atom2 = universe.select_atoms(" atom {0!s} {1!s} C1\' ".format(seg, i))
atom3 = universe.select_atoms(" atom {0!s} {1!s} C2\' ".format(seg, i))
atom4 = universe.select_atoms(" atom {0!s} {1!s} C3\' ".format(seg, i))
atom5 = universe.select_atoms(" atom {0!s} {1!s} C4\' ".format(seg, i))
data1 = atom1.positions
data2 = atom2.positions
data3 = atom3.positions
data4 = atom4.positions
data5 = atom5.positions
r0 = (data1 + data2 + data3 + data4 + data5) * (1.0 / 5.0)
r1 = data1 - r0
r2 = data2 - r0
r3 = data3 - r0
r4 = data4 - r0
r5 = data5 - r0
R1 = ((r1 * sin(2 * pi * 0.0 / 5.0)) + (r2 * sin(2 * pi * 1.0 / 5.0)) +
(r3 * sin(2 * pi * 2.0 / 5.0)) + (r4 * sin(2 * pi * 3.0 / 5.0)) +
(r5 * sin(2 * pi * 4.0 / 5.0)))
R2 = ((r1 * cos(2 * pi * 0.0 / 5.0)) + (r2 * cos(2 * pi * 1.0 / 5.0)) +
(r3 * cos(2 * pi * 2.0 / 5.0)) + (r4 * cos(2 * pi * 3.0 / 5.0)) +
(r5 * cos(2 * pi * 4.0 / 5.0)))
x = np.cross(R1[0], R2[0])
n = x / sqrt(pow(x[0], 2) + pow(x[1], 2) + pow(x[2], 2))
r1_d = np.dot(r1, n)
r2_d = np.dot(r2, n)
r3_d = np.dot(r3, n)
r4_d = np.dot(r4, n)
r5_d = np.dot(r5, n)
D = ((r1_d * sin(4 * pi * 0.0 / 5.0)) + (r2_d * sin(4 * pi * 1.0 / 5.0))
+ (r3_d * sin(4 * pi * 2.0 / 5.0)) + (r4_d * sin(4 * pi * 3.0 / 5.0))
+ (r5_d * sin(4 * pi * 4.0 / 5.0))) * -1 * sqrt(2.0 / 5.0)
C = ((r1_d * cos(4 * pi * 0.0 / 5.0)) + (r2_d * cos(4 * pi * 1.0 / 5.0))
+ (r3_d * cos(4 * pi * 2.0 / 5.0)) + (r4_d * cos(4 * pi * 3.0 / 5.0))
+ (r5_d * cos(4 * pi * 4.0 / 5.0))) * sqrt(2.0 / 5.0)
phase_ang = (atan2(D, C) + (pi / 2.)) * 180. / pi
return phase_ang % 360
[docs]
def phase_as(universe, seg, i):
"""Pseudo-angle describing the phase of the ribose pucker for residue `i` using the AS method
The angle is computed by the position vector of atoms in the ribose ring.
Parameters
----------
universe : Universe
:class:`~MDAnalysis.core.universe.Universe` containing the trajectory
seg : str
segment id for base
i : int
resid of the first base
Returns
-------
float
phase angle in degrees
.. versionadded:: 0.7.6
"""
angle1 = universe.select_atoms(" atom {0!s} {1!s} C1\' ".format(seg, i),
" atom {0!s} {1!s} C2\' ".format(seg, i),
" atom {0!s} {1!s} C3\' ".format(seg, i),
" atom {0!s} {1!s} C4\' ".format(seg, i))
angle2 = universe.select_atoms(" atom {0!s} {1!s} C2\' ".format(seg, i),
" atom {0!s} {1!s} C3\' ".format(seg, i),
" atom {0!s} {1!s} C4\' ".format(seg, i),
" atom {0!s} {1!s} O4\' ".format(seg, i))
angle3 = universe.select_atoms(" atom {0!s} {1!s} C3\' ".format(seg, i),
" atom {0!s} {1!s} C4\' ".format(seg, i),
" atom {0!s} {1!s} O4\' ".format(seg, i),
" atom {0!s} {1!s} C1\' ".format(seg, i))
angle4 = universe.select_atoms(" atom {0!s} {1!s} C4\' ".format(seg, i),
" atom {0!s} {1!s} O4\' ".format(seg, i),
" atom {0!s} {1!s} C1\' ".format(seg, i),
" atom {0!s} {1!s} C2\' ".format(seg, i))
angle5 = universe.select_atoms(" atom {0!s} {1!s} O4\' ".format(seg, i),
" atom {0!s} {1!s} C1\' ".format(seg, i),
" atom {0!s} {1!s} C2\' ".format(seg, i),
" atom {0!s} {1!s} C3\' ".format(seg, i))
data1 = angle1.dihedral.value()
data2 = angle2.dihedral.value()
data3 = angle3.dihedral.value()
data4 = angle4.dihedral.value()
data5 = angle5.dihedral.value()
B = ((data1 * sin(2 * 2 * pi * (1 - 1.) / 5.))
+ (data2 * sin(2 * 2 * pi * (2 - 1.) / 5.))
+ (data3 * sin(2 * 2 * pi * (3 - 1.) / 5.))
+ (data4 * sin(2 * 2 * pi * (4 - 1.) / 5.))
+ (data5 * sin(2 * 2 * pi * (5 - 1.) / 5.))) * -2. / 5.
A = ((data1 * cos(2 * 2 * pi * (1 - 1.) / 5.))
+ (data2 * cos(2 * 2 * pi * (2 - 1.) / 5.))
+ (data3 * cos(2 * 2 * pi * (3 - 1.) / 5.))
+ (data4 * cos(2 * 2 * pi * (4 - 1.) / 5.))
+ (data5 * cos(2 * 2 * pi * (5 - 1.) / 5.))) * 2. / 5.
phase_ang = atan2(B, A) * 180. / pi
return phase_ang % 360
[docs]
def tors(universe, seg, i):
"""Calculation of nucleic backbone dihedral angles.
The dihedral angles are alpha, beta, gamma, delta, epsilon, zeta, chi.
The dihedral is computed based on position of atoms for resid `i`.
Parameters
----------
universe : Universe
:class:`~MDAnalysis.core.universe.Universe` containing the trajectory
seg : str
segment id for base
i : int
resid of the first base
Returns
-------
[alpha, beta, gamma, delta, epsilon, zeta, chi] : list of floats
torsion angles in degrees
Notes
-----
If failure occurs be sure to check the segment identification.
.. versionadded:: 0.7.6
"""
a = universe.select_atoms(" atom {0!s} {1!s} O3\' ".format(seg, i - 1),
" atom {0!s} {1!s} P ".format(seg, i),
" atom {0!s} {1!s} O5\' ".format(seg, i),
" atom {0!s} {1!s} C5\' ".format(seg, i))
b = universe.select_atoms(" atom {0!s} {1!s} P ".format(seg, i),
" atom {0!s} {1!s} O5\' ".format(seg, i),
" atom {0!s} {1!s} C5\' ".format(seg, i),
" atom {0!s} {1!s} C4\' ".format(seg, i))
g = universe.select_atoms(" atom {0!s} {1!s} O5\' ".format(seg, i),
" atom {0!s} {1!s} C5\' ".format(seg, i),
" atom {0!s} {1!s} C4\' ".format(seg, i),
" atom {0!s} {1!s} C3\' ".format(seg, i))
d = universe.select_atoms(" atom {0!s} {1!s} C5\' ".format(seg, i),
" atom {0!s} {1!s} C4\' ".format(seg, i),
" atom {0!s} {1!s} C3\' ".format(seg, i),
" atom {0!s} {1!s} O3\' ".format(seg, i))
e = universe.select_atoms(" atom {0!s} {1!s} C4\' ".format(seg, i),
" atom {0!s} {1!s} C3\' ".format(seg, i),
" atom {0!s} {1!s} O3\' ".format(seg, i),
" atom {0!s} {1!s} P ".format(seg, i + 1))
z = universe.select_atoms(" atom {0!s} {1!s} C3\' ".format(seg, i),
" atom {0!s} {1!s} O3\' ".format(seg, i),
" atom {0!s} {1!s} P ".format(seg, i + 1),
" atom {0!s} {1!s} O5\' ".format(seg, i + 1))
c = universe.select_atoms(" atom {0!s} {1!s} O4\' ".format(seg, i),
" atom {0!s} {1!s} C1\' ".format(seg, i),
" atom {0!s} {1!s} N9 ".format(seg, i),
" atom {0!s} {1!s} C4 ".format(seg, i))
if len(c) < 4:
c = universe.select_atoms(" atom {0!s} {1!s} O4\' ".format(seg, i),
" atom {0!s} {1!s} C1\' ".format(seg, i),
" atom {0!s} {1!s} N1 ".format(seg, i),
" atom {0!s} {1!s} C2 ".format(seg, i))
alpha = a.dihedral.value() % 360
beta = b.dihedral.value() % 360
gamma = g.dihedral.value() % 360
delta = d.dihedral.value() % 360
epsilon = e.dihedral.value() % 360
zeta = z.dihedral.value() % 360
chi = c.dihedral.value() % 360
return [alpha, beta, gamma, delta, epsilon, zeta, chi]
[docs]
def tors_alpha(universe, seg, i):
"""alpha backbone dihedral
The dihedral is computed based on position atoms for resid `i`.
Parameters
----------
universe : Universe
:class:`~MDAnalysis.core.universe.Universe` containing the trajectory
seg : str
segment id for base
i : int
resid of the first base
Returns
-------
alpha : float
torsion angle in degrees
.. versionadded:: 0.7.6
"""
a = universe.select_atoms(" atom {0!s} {1!s} O3\' ".format(seg, i - 1),
" atom {0!s} {1!s} P ".format(seg, i),
" atom {0!s} {1!s} O5\' ".format(seg, i),
" atom {0!s} {1!s} C5\' ".format(seg, i))
alpha = a.dihedral.value() % 360
return alpha
[docs]
def tors_beta(universe, seg, i):
"""beta backbone dihedral
The dihedral is computed based on position atoms for resid `i`.
Parameters
----------
universe : Universe
:class:`~MDAnalysis.core.universe.Universe` containing the trajectory
seg : str
segment id for base
i : int
resid of the first base
Returns
-------
beta : float
torsion angle in degrees
.. versionadded:: 0.7.6
"""
b = universe.select_atoms(" atom {0!s} {1!s} P ".format(seg, i),
" atom {0!s} {1!s} O5\' ".format(seg, i),
" atom {0!s} {1!s} C5\' ".format(seg, i),
" atom {0!s} {1!s} C4\' ".format(seg, i))
beta = b.dihedral.value() % 360
return beta
[docs]
def tors_gamma(universe, seg, i):
""" Gamma backbone dihedral
The dihedral is computed based on position atoms for resid `i`.
Parameters
----------
universe : Universe
:class:`~MDAnalysis.core.universe.Universe` containing the trajectory
seg : str
segment id for base
i : int
resid of the first base
Returns
-------
gamma : float
torsion angle in degrees
.. versionadded:: 0.7.6
"""
g = universe.select_atoms(" atom {0!s} {1!s} O5\' ".format(seg, i),
" atom {0!s} {1!s} C5\' ".format(seg, i),
" atom {0!s} {1!s} C4\' ".format(seg, i),
" atom {0!s} {1!s} C3\' ".format(seg, i))
gamma = g.dihedral.value() % 360
return gamma
[docs]
def tors_delta(universe, seg, i):
"""delta backbone dihedral
The dihedral is computed based on position atoms for resid `i`.
Parameters
----------
universe : Universe
:class:`~MDAnalysis.core.universe.Universe` containing the trajectory
seg : str
segment id for base
i : int
resid of the first base
Returns
-------
delta : float
torsion angle in degrees
.. versionadded:: 0.7.6
"""
d = universe.select_atoms(" atom {0!s} {1!s} C5\' ".format(seg, i),
" atom {0!s} {1!s} C4\' ".format(seg, i),
" atom {0!s} {1!s} C3\' ".format(seg, i),
" atom {0!s} {1!s} O3\' ".format(seg, i))
delta = d.dihedral.value() % 360
return delta
[docs]
def tors_eps(universe, seg, i):
"""Epsilon backbone dihedral
The dihedral is computed based on position atoms for resid `i`.
Parameters
----------
universe : Universe
:class:`~MDAnalysis.core.universe.Universe` containing the trajectory
seg : str
segment id for base
i : int
resid of the first base
Returns
-------
epsilon : float
torsion angle in degrees
.. versionadded:: 0.7.6
"""
e = universe.select_atoms(" atom {0!s} {1!s} C4\' ".format(seg, i),
" atom {0!s} {1!s} C3\' ".format(seg, i),
" atom {0!s} {1!s} O3\' ".format(seg, i),
" atom {0!s} {1!s} P ".format(seg, i + 1))
epsilon = e.dihedral.value() % 360
return epsilon
[docs]
def tors_zeta(universe, seg, i):
"""Zeta backbone dihedral
The dihedral is computed based on position atoms for resid `i`.
Parameters
----------
universe : Universe
:class:`~MDAnalysis.core.universe.Universe` containing the trajectory
seg : str
segment id for base
i : int
resid of the first base
Returns
-------
zeta : float
torsion angle in degrees
.. versionadded:: 0.7.6
"""
z = universe.select_atoms(" atom {0!s} {1!s} C3\' ".format(seg, i),
" atom {0!s} {1!s} O3\' ".format(seg, i),
" atom {0!s} {1!s} P ".format(seg, i + 1),
" atom {0!s} {1!s} O5\' ".format(seg, i + 1))
zeta = z.dihedral.value() % 360
return zeta
[docs]
def tors_chi(universe, seg, i):
"""chi nucleic acid dihedral
The dihedral is computed based on position atoms for resid `i`.
Parameters
----------
universe : Universe
:class:`~MDAnalysis.core.universe.Universe` containing the trajectory
seg : str
segment id for base
i : int
resid of the first base
Returns
-------
chi : float
torsion angle in degrees
.. versionadded:: 0.7.6
"""
c = universe.select_atoms(" atom {0!s} {1!s} O4\' ".format(seg, i),
" atom {0!s} {1!s} C1\' ".format(seg, i),
" atom {0!s} {1!s} N9 ".format(seg, i),
" atom {0!s} {1!s} C4 ".format(seg, i))
if len(c) < 4:
c = universe.select_atoms(" atom {0!s} {1!s} O4\' ".format(seg, i),
" atom {0!s} {1!s} C1\' ".format(seg, i),
" atom {0!s} {1!s} N1 ".format(seg, i),
" atom {0!s} {1!s} C2 ".format(seg, i))
chi = c.dihedral.value() % 360
return chi
[docs]
def hydroxyl(universe, seg, i):
"""2-hydroxyl dihedral. Useful only for RNA calculations.
.. Note:: This dihedral calculation will only work if using atom
names as documented by charmm force field parameters,
namely "C1', C2', O2', H2'".
Parameters
----------
universe : Universe
:class:`~MDAnalysis.core.universe.Universe` containing the trajectory
seg : str
segment id for base
i : int
resid of the first base
Returns
-------
hydroxyl_angle : float
torsion angle in degrees
.. versionadded:: 0.7.6
"""
h = universe.select_atoms("atom {0!s} {1!s} C1'".format(seg, i),
"atom {0!s} {1!s} C2'".format(seg, i),
"atom {0!s} {1!s} O2'".format(seg, i),
"atom {0!s} {1!s} H2'".format(seg, i))
try:
hydr = h.dihedral.value() % 360
except ValueError:
errmsg = (f"Resid {i} does not contain atoms C1', C2', O2', H2' but "
f"atoms {list(h.atoms)}")
raise ValueError(errmsg) from None
return hydr
[docs]
def pseudo_dihe_baseflip(universe, bp1, bp2, i,
seg1="SYSTEM", seg2="SYSTEM", seg3="SYSTEM"):
"""pseudo dihedral for flipped bases. Useful only for nucleic acid base flipping
The dihedral is computed based on position atoms for resid `i`
.. Note:: This dihedral calculation will only work if using atom names as
documented by charmm force field parameters.
Parameters
----------
universe : Universe
:class:`~MDAnalysis.core.universe.Universe` containing the
trajectory
bp1 : int
resid that base pairs with `bp2`
bp2 : int
resid below the base that flips
i : int
resid of the base that flips
segid1 : str (optional)
segid of resid base pairing with `bp2`
segid2 : str (optional)
segid, same as that of segid of flipping resid `i`
segid3 : str (optional)
segid of resid `i` that flips
Returns
-------
float
pseudo dihedral angle in degrees
.. versionadded:: 0.8.0
"""
bf1 = universe.select_atoms(
" ( segid {0!s} and resid {1!s} and nucleicbase ) "
"or ( segid {2!s} and resid {3!s} and nucleicbase ) "
.format( seg1, bp1, seg2, bp2))
bf4 = universe.select_atoms("(segid {0!s} and resid {1!s} and nucleicbase) ".format(seg3, i))
bf2 = universe.select_atoms("(segid {0!s} and resid {1!s} and nucleicsugar) ".format(seg2, bp2))
bf3 = universe.select_atoms("(segid {0!s} and resid {1!s} and nucleicsugar) ".format(seg3, i))
x = [bf1.center_of_mass(), bf2.center_of_mass(),
bf3.center_of_mass(), bf4.center_of_mass()]
pseudo = mdamath.dihedral(x[0] - x[1], x[1] - x[2], x[2] - x[3])
pseudo = np.rad2deg(pseudo) % 360
return pseudo