diff --git a/docs/source/api.rst b/docs/source/api.rst
index 00a15bc..e5b1456 100644
--- a/docs/source/api.rst
+++ b/docs/source/api.rst
@@ -2,3 +2,6 @@ API Documentation
=================
.. automodule:: helanal.helanal
+ :show-inheritance:
+ :members:
+ :inherited-members:
diff --git a/docs/source/conf.py b/docs/source/conf.py
index b93b4e5..f8fbc9b 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -24,11 +24,11 @@
project = "helanal"
copyright = (
- "2024, Your name (or your organization/company/team). "
+ "2024, Fiona Naughton. "
"Project structure based on the "
"MDAnalysis Cookiecutter version 0.1"
)
-author = "Your name (or your organization/company/team)"
+author = "Fiona Naughton"
# The short X.Y version
version = ""
@@ -54,7 +54,6 @@
"sphinx.ext.intersphinx",
"sphinx.ext.extlinks",
"mdanalysis_sphinx_theme",
- "nbsphinx",
]
autosummary_generate = True
diff --git a/docs/source/getting_started.rst b/docs/source/getting_started.rst
index 7c955b8..cd05103 100644
--- a/docs/source/getting_started.rst
+++ b/docs/source/getting_started.rst
@@ -6,33 +6,93 @@ Installation
*TBA*
-Example use
+Basic Usage
-----------
-Import MDAnalysis and helanal::
+helanal is used to calculate helix properties for helices with at least
+9 residues.
+
+Here, we demonstrate the basic usage and options of helanal. First, import
+MDAnalysis and helanal::
import MDAnalysis as mda
- from MDAnalysis.tests.datafiles import PSF, DCD
import helanal
-You can pass in a single selection::
+For this example, we will use datafiles from the MDAnalysis tests::
+ from MDAnalysis.tests.datafiles import PSF, DCD
u = mda.Universe(PSF, DCD)
- hel = helanl.HELANAL(u, select='name CA and resnum 161-187')
+
+To analyse a single helix, pass a selection with one atom per residue
+(normally this will be the Cα atoms)::
+
+ hel = helanal.HELANAL(u, select='name CA and resnum 161-187')
hel.run()
-All computed properties are available in ``.results``::
+All computed properties for each simulation frame will be available as
+arrays in ``.results``, e.g.::
+
+ helix_tilt = hel.results.global_tilts
+
+For more details about the properties calculated and method used,
+see :doc:`properties_computed`.
+
+
+Further Options
+---------------
+
+- The **simulation frames** over which analysis is performed can be specified
+ using ``start``, ``stop``, and/or ``step``, or by providing a list
+ ``frames`` with which to slice the trajectory::
+
+ hel.run(start=5, step=10)
- print(hel.results.summary)
+- The **reference axis** to which helix tilts and screw angles are calculated
+ can be specified using ``ref_axis``::
-Alternatively, you can analyse several helices at once by passing
-in multiple selection strings::
+ hel_xaxis = helanal.HELANAL(u, select='name CA and resnum 161-187',
+ ref_axis=[1,0,0])
- hel2 = helanal.HELANAL(u, select=('name CA and resnum 100-160',
- 'name CA and resnum 200-230'))
+- To analyse **multiple helices** at once, pass in a list of selection
+ strings::
-The :func:`helix_analysis` function will carry out helix analysis on
-atom positions, treating each row of coordinates as an alpha-carbon
-equivalent::
+ hel_multi = helanal.HELANAL(u, select=('name CA and resnum 100-160',
+ 'name CA and resnum 200-230'))
+
+ Each property in ``.results`` will now be a list of arrays, one for
+ each helix.
+
+- By default, the results of a single-helix analysis are flattened, such
+ that each result is a single array, rather than a list-of-arrays of
+ length 1. This behaviour can be turned off (i.e. to be consistent
+ with the behaviour for multi-helix analysis) using
+ ``flatten_single_helix``::
+
+ hel_noflat = helanal.HELANAL(u, select='name CA and resnum 161-187',
+ flatten_single_helix=False)
+
+- Analysis can be carried out directly on a single set of positions using
+ the :func:`helix_analysis` function::
hel_xyz = helanal.helix_analysis(u.atoms.positions, ref_axis=[0, 0, 1])
+ Results are returned as a dictionary.
+
+
+
+Citations
+---------
+
+helanal is based on the `HELANAL algorithm`_ from [Bansal2000]_, which itself
+uses the method of [Sugeta1967]_ to characterise each local axis. Please cite
+them when using this module in published work.
+
+.. [Sugeta1967] Sugeta, H. and Miyazawa, T. 1967. General method for
+ calculating helical parameters of polymer chains from bond lengths, bond
+ angles and internal rotation angles. *Biopolymers* 5 673 - 679
+
+.. [Bansal2000] Bansal M, Kumar S, Velavan R. 2000.
+ HELANAL - A program to characterise helix geometry in proteins.
+ *J Biomol Struct Dyn.* 17(5):811-819.
+
+.. _`HELANAL algorithm`:
+ https://web.archive.org/web/20090226192455/http://www.ccrnp.ncifcrf.gov/users/kumarsan/HELANAL/helanal.f
diff --git a/docs/source/images/axes.svg b/docs/source/images/axes.svg
new file mode 100755
index 0000000..811d34a
--- /dev/null
+++ b/docs/source/images/axes.svg
@@ -0,0 +1 @@
+
\ No newline at end of file
diff --git a/docs/source/images/bends.svg b/docs/source/images/bends.svg
new file mode 100755
index 0000000..f01cff0
--- /dev/null
+++ b/docs/source/images/bends.svg
@@ -0,0 +1 @@
+
\ No newline at end of file
diff --git a/docs/source/images/global_axis.svg b/docs/source/images/global_axis.svg
new file mode 100755
index 0000000..d858e36
--- /dev/null
+++ b/docs/source/images/global_axis.svg
@@ -0,0 +1 @@
+
\ No newline at end of file
diff --git a/docs/source/images/global_tilts.svg b/docs/source/images/global_tilts.svg
new file mode 100755
index 0000000..a8d4ba6
--- /dev/null
+++ b/docs/source/images/global_tilts.svg
@@ -0,0 +1 @@
+
\ No newline at end of file
diff --git a/docs/source/images/heights.svg b/docs/source/images/heights.svg
new file mode 100755
index 0000000..014ed27
--- /dev/null
+++ b/docs/source/images/heights.svg
@@ -0,0 +1 @@
+
\ No newline at end of file
diff --git a/docs/source/images/helix_directions.svg b/docs/source/images/helix_directions.svg
new file mode 100755
index 0000000..1184794
--- /dev/null
+++ b/docs/source/images/helix_directions.svg
@@ -0,0 +1 @@
+
\ No newline at end of file
diff --git a/docs/source/images/origins.svg b/docs/source/images/origins.svg
new file mode 100755
index 0000000..9419744
--- /dev/null
+++ b/docs/source/images/origins.svg
@@ -0,0 +1 @@
+
\ No newline at end of file
diff --git a/docs/source/images/screw_angles.svg b/docs/source/images/screw_angles.svg
new file mode 100755
index 0000000..21b7a92
--- /dev/null
+++ b/docs/source/images/screw_angles.svg
@@ -0,0 +1 @@
+
\ No newline at end of file
diff --git a/docs/source/images/twists.svg b/docs/source/images/twists.svg
new file mode 100755
index 0000000..312cbb9
--- /dev/null
+++ b/docs/source/images/twists.svg
@@ -0,0 +1 @@
+
\ No newline at end of file
diff --git a/docs/source/images/window.svg b/docs/source/images/window.svg
new file mode 100755
index 0000000..d828c7c
--- /dev/null
+++ b/docs/source/images/window.svg
@@ -0,0 +1 @@
+
\ No newline at end of file
diff --git a/docs/source/index.rst b/docs/source/index.rst
index 0e46402..421ad77 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -7,13 +7,24 @@ Welcome to helanal's documentation!
=========================================================
This module contains code to analyse protein helices using the
-HELANAL_ algorithm
-([Bansal2000]_ , [Sugeta1967]_ ).
+HELANAL_ algorithm ([Bansal2000]_, [Sugeta1967]_).
HELANAL_ quantifies the geometry of helices in proteins on the basis of their
Cα atoms. It can determine local structural features such as the local
-helical twist and rise, virtual torsion angle, local helix origins and
-bending angles between successive local helix axes.
+helical twist and rise, local helix origins and bending angles along the
+helix.
+
+.. toctree::
+ :maxdepth: 1
+ :caption: Contents
+
+ getting_started
+ properties_computed
+ api
+
+
+References
+^^^^^^^^^^
.. _HELANAL: https://pubmed.ncbi.nlm.nih.gov/10798526/
@@ -25,20 +36,3 @@ bending angles between successive local helix axes.
HELANAL - A program to characterise helix geometry in proteins.
*J Biomol Struct Dyn.* 17(5):811-819.
-
-.. toctree::
- :maxdepth: 1
- :caption: Contents:
-
- getting_started
- example
- api
-
-
-
-Indices and tables
-==================
-
-* :ref:`genindex`
-* :ref:`modindex`
-* :ref:`search`
diff --git a/docs/source/properties_computed.rst b/docs/source/properties_computed.rst
new file mode 100644
index 0000000..8ad0d69
--- /dev/null
+++ b/docs/source/properties_computed.rst
@@ -0,0 +1,117 @@
+Computed Properties
+===================
+
+Helanal follows the procedure of Sugeta and Miyazawa [Sugeta1967]_ and is based
+on the `Fortran HELANAL implementation`_ by [Bansal2000]_.
+Properties are computed for a 'window' of four consecutive :math:`C_α` atoms,
+and this window is then slid along the length of the helix in one-residue steps.
+
+For each window consisting of atoms :math:`c_i`, :math:`c_{i+1}`,
+:math:`c_{i+2}`, :math:`c_{i+3}`, the vectors :math:`\mathbf{B_1}`,
+:math:`\mathbf{B_2}`, and :math:`\mathbf{B_3}` joining (respectively) atoms
+:math:`c_i` → :math:`c_{i+1}`, :math:`c_{i+1}` → :math:`c_{i+2}` and
+:math:`c_{i+2}` → :math:`c_{i+3}` are calculated, along with the vectors
+:math:`\mathbf{D_1} = \mathbf{B_1} - \mathbf{B_2}` and
+:math:`\mathbf{D_2} = \mathbf{B_2} - \mathbf{B_3}`.
+
+.. image:: images/window.svg
+ :width: 60%
+ :align: center
+
+From these, the helix properties below are computed in each simulation frame.
+These properties are available in ``.results`` as arrays, the shape of which
+depends on the number of residues :math:`n_{res}` (and the property being
+calculated). **Note that each helix must contain at least 9 residues.**
+
+If multiple helices are being analysed, the results are returned as lists (of
+length :math:`n_{helices}`) of arrays of the indicated shape.
+
+All angles are in degrees.
+
+.. list-table::
+ :widths: 40 20 20
+ :header-rows: 1
+
+ * - Description
+ -
+ - Shape
+ * - ``local_helix_directions``: the normalised vector :math:`\mathbf{D_1}`
+ (or :math:`\mathbf{D_2}`) for each atom :math:`c_{i+1}` (or
+ :math:`c_{i+2}`).
+
+ Assuming ~even spacing of the atoms, this vector will bisect the angle
+ formed by (:math:`c_i,c_{i+1},c_{i+2}`), lie approximately in the plane
+ perpendicular to the helix axis, and point from the projected local
+ helix centre to the atom :math:`c_{i+1}`.
+ - .. image:: images/helix_directions.svg
+ :width: 100%
+ - :math:`(n_{frames},` :math:`n_{res}-2, 3)`
+ * - ``local_twists``: the approximate 'twist' of the helix between atoms
+ :math:`c_{i+1}` and :math:`c_{i+2}`, calculated as the angle :math:`θ`
+ between :math:`\mathbf{D_1}` and :math:`\mathbf{D_2}`.
+ - .. image:: images/twists.svg
+ :width: 100%
+ - :math:`(n_{frames},` :math:`n_{res}-3)`
+ * - ``local_nres_per_turn``: the number of residues that fit in one complete
+ turn of the helix, based on `local_twist`.
+ -
+ - :math:`(n_{frames},` :math:`n_{res}-3)`
+ * - ``local_origins``: the projected centre of each 4-atom window, in line
+ with atom :math:`c_{i+1}`.
+
+ Calculated as the approximate intersection of :math:`\mathbf{D_1}` and
+ :math:`\mathbf{D_2}` projected on the perpendicular plane, assuming ~even
+ spacing of atoms.
+ - .. image:: images/origins.svg
+ :width: 100%
+ - :math:`(n_{frames},` :math:`n_{res}-2,` :math:`3)`
+ * - ``local_axes``: the (normalised) central axis :math:`\mathbf{A}` of the
+ 4-atom window, calculated as the normal to the two vectors
+ :math:`\mathbf{D_1}` and :math:`\mathbf{D_2}`.
+ - .. image:: images/axes.svg
+ :width: 100%
+ - :math:`(n_{frames},` :math:`n_{res}-3,` :math:`3)`
+ * - ``local_heights``: the 'rise' :math:`h` of the helix (in the direction
+ of `local_axes`) between atoms :math:`c_{i+1}` and :math:`c_{i+2}`.
+ - .. image:: images/heights.svg
+ :width: 100%
+ - :math:`(n_{frames},` :math:`n_{res}-3)`
+ * - ``local_bends``: the angle of bending of the helix between adjacent
+ 4-atom windows, i.e. the angle :math:`β` between the `local_axes`
+ :math:`\mathbf{A_i}` (of atoms :math:`c_i,c_{i+1},c_{i+2},c_{i+3}`) and
+ :math:`\mathbf{A_{i+3}}` (of atoms
+ :math:`c_{i+3},c_{i+4},c_{i+5},c_{i+6}`).
+ - .. image:: images/bends.svg
+ :width: 100%
+ - :math:`(n_{frames},` :math:`n_{res}-6)`
+ * - ``all_bends``: pair-wise matrix of angles between all pairs of
+ `local_axes`.
+ -
+ - :math:`(n_{frames},` :math:`n_{res}-3,` :math:`n_{res}-3)`
+ * - ``global_axis``: the length-wise axis :math:`\mathbf{G}` for the overall
+ helix, pointing from the end of the helix to the start. Calculated as the
+ vector of best fit through all `local_origins`.
+ - .. image:: images/global_axis.svg
+ :width: 100%
+ - :math:`(n_{frames},` :math:`3)`
+ * - ``global_tilts``: the angle :math:`γ` between the `global_axis`
+ :math:`\mathbf{G}` and the reference axis (specified by the ``ref_axis``
+ option). If no axis is specified, the z-axis is used.
+ - .. image:: images/global_tilts.svg
+ :width: 100%
+ - :math:`(n_{frames},)`
+ * - ``local_screw_angles``: The cylindrical azimuthal angle :math:`α` of
+ atom :math:`c_{i+1}` (in the range -pi to pi).
+
+ This is measured as the angle between the `ref_axis` to the
+ `local_helix_directions` vector :math:`\mathbf{D}`, when both are
+ projected on a plane perpendicular to `global_axis`.
+ - .. image:: images/screw_angles.svg
+ :width: 100%
+ - :math:`(n_{frames},` :math:`n_{res}-2)`
+
+A summary of the results, including mean, sample standard deviation and mean
+absolute deviation is also provided in ``results.summary``.
+
+.. _`Fortran HELANAL implementation`:
+ https://web.archive.org/web/20090226192455/http://www.ccrnp.ncifcrf.gov/users/kumarsan/HELANAL/helanal.f
diff --git a/helanal/helanal.py b/helanal/helanal.py
index 2a97468..3977093 100644
--- a/helanal/helanal.py
+++ b/helanal/helanal.py
@@ -22,47 +22,6 @@
# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
#
-"""
-HELANAL --- analysis of protein helices
-=======================================
-
-This module contains code to analyse protein helices using the
-HELANAL_ algorithm
-([Bansal2000]_ , [Sugeta1967]_ ).
-
-HELANAL_ quantifies the geometry of helices in proteins on the basis of their
-Cα atoms. It can determine local structural features such as the local
-helical twist and rise, virtual torsion angle, local helix origins and
-bending angles between successive local helix axes.
-
-.. _HELANAL: https://pubmed.ncbi.nlm.nih.gov/10798526/
-
-.. [Sugeta1967] Sugeta, H. and Miyazawa, T. 1967. General method for
- calculating helical parameters of polymer chains from bond lengths, bond
- angles and internal rotation angles. *Biopolymers* 5 673 - 679
-
-.. [Bansal2000] Bansal M, Kumar S, Velavan R. 2000.
- HELANAL - A program to characterise helix geometry in proteins.
- *J Biomol Struct Dyn.* 17(5):811-819.
-
-
-Classes
--------
-
-.. autoclass:: HELANAL
-
-
-Functions
----------
-
-.. autofunction:: helix_analysis
-
-.. autofunction:: vector_of_best_fit
-
-.. autofunction:: local_screw_angles
-
-"""
-
import warnings
import numpy as np
@@ -194,11 +153,11 @@ def helix_analysis(positions, ref_axis=[0, 0, 1]):
# ^ ^
# \ / bi
- # \ /
+ # \ Vi+1 /
# CA_i+2 <----- CA_i+1
# / \ / ^
- # / r \ / \
- # V / \ θ / \
+ # / r \ / \Vi
+ # Vi+2/ \ θ / \
# / \ / CA_i
# v origin
# CA_i+3
@@ -206,6 +165,8 @@ def helix_analysis(positions, ref_axis=[0, 0, 1]):
# V: vectors
# bi: approximate "bisectors" in plane of screen
# Note: not real bisectors, as the vectors aren't normalised
+ # but assuming ~equal spacing of CA atoms, i.e. |Vi| ~ |Vi+1|
+ # should be approximately true.
# θ: local_twists
# origin: origins
# local_axes: perpendicular to plane of screen. Orthogonal to "bisectors"
@@ -214,6 +175,7 @@ def helix_analysis(positions, ref_axis=[0, 0, 1]):
bisectors = vectors[:-1] - vectors[1:] # (n_res-2, 3)
bimags = mdamath.pnorm(bisectors) # (n_res-2,)
adjacent_mag = bimags[:-1] * bimags[1:] # (n_res-3,)
+ local_helix_directions = (bisectors.T/bimags).T # (n_res-2, 3)
# find angle between bisectors for twist and n_residue/turn
cos_theta = mdamath.pdot(bisectors[:-1], bisectors[1:])/adjacent_mag
@@ -238,18 +200,37 @@ def helix_analysis(positions, ref_axis=[0, 0, 1]):
local_bends = np.diagonal(bend_matrix, offset=3) # (n_res-6,)
# radius of local cylinder
+ ## This appears to be calculated based on the following.
+ ##
+ ## CA_i+2 CA_i+2
+ ## r / | \ / ^
+ ## / | \Vi+1 / \
+ ## / θ | \ / bi \
+ ## origin-----------CA_i+1 ----->CA_i+1
+ ## x | y / \ ^
+ ## | /Vi \ /
+ ## | / \ /
+ ## CA_i CA_i
+ ##
+ ## If we assume atoms are ~evenly spaced around a circle, bi are ~bisectors
+ ## and their projections should meet at the centre of the circle
+ ## i.e. r ~= x + y
+ ## bi = Vi-Vi+1 is the digonal of a parallelogram formed by Vi and Vi+1.
+ ## A vector from CAi to CAi+2 is the other diangonal, so is perpendular
+ ## to and bisects bi - i.e. x = rcos0 and y = |bi|/2
+ ## Solve for r and substitute |bi| -> adjacent_mag**0.5 to take into account
+ ## some difference.
radii = (adjacent_mag**0.5) / (2*(1.0-cos_theta)) # (n_res-3,)
+
# special case: angle b/w bisectors is 0 (should virtually never happen)
# guesstimate radius = half bisector magnitude
radii = np.where(cos_theta != 1, radii, (adjacent_mag**0.5)/2)
# height of local cylinder
heights = np.abs(mdamath.pdot(vectors[1:-1], local_axes)) # (n_res-3,)
- local_helix_directions = (bisectors.T/bimags).T # (n_res-2, 3)
-
# get origins by subtracting radius from atom i+1
origins = positions[1:-1].copy() # (n_res-2, 3)
- origins[:-1] -= (radii*local_helix_directions[:-1].T).T
+ origins[:-1] -= radii[:,None]*local_helix_directions[:-1]
# subtract radius from atom i+2 in last one
origins[-1] -= radii[-1]*local_helix_directions[-1]