From b74a80192bd20d202da4f6463f65d65e1c46fc78 Mon Sep 17 00:00:00 2001 From: Fiona Naughton Date: Wed, 6 Nov 2024 14:25:08 +1100 Subject: [PATCH] Updating documentation (#1) Documentation overhaul --- docs/source/api.rst | 3 + docs/source/conf.py | 4 +- docs/source/getting_started.rst | 88 +++++++++++++++--- docs/source/images/axes.svg | 1 + docs/source/images/bends.svg | 1 + docs/source/images/global_axis.svg | 1 + docs/source/images/global_tilts.svg | 1 + docs/source/images/heights.svg | 1 + docs/source/images/helix_directions.svg | 1 + docs/source/images/origins.svg | 1 + docs/source/images/screw_angles.svg | 1 + docs/source/images/twists.svg | 1 + docs/source/images/window.svg | 1 + docs/source/index.rst | 36 +++----- docs/source/properties_computed.rst | 117 ++++++++++++++++++++++++ helanal/helanal.py | 75 ++++++--------- 16 files changed, 249 insertions(+), 84 deletions(-) create mode 100755 docs/source/images/axes.svg create mode 100755 docs/source/images/bends.svg create mode 100755 docs/source/images/global_axis.svg create mode 100755 docs/source/images/global_tilts.svg create mode 100755 docs/source/images/heights.svg create mode 100755 docs/source/images/helix_directions.svg create mode 100755 docs/source/images/origins.svg create mode 100755 docs/source/images/screw_angles.svg create mode 100755 docs/source/images/twists.svg create mode 100755 docs/source/images/window.svg create mode 100644 docs/source/properties_computed.rst 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 120d24a..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 = "" 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 @@ +D2D1ci+1ci+2D1D2A \ 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 @@ +Ai+3Aicici+1ci+2ci+3ci+4ci+5ci+6ci+7ci+8βci+9 \ 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 @@ +originsG \ 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 @@ +ref_axisγG \ 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 @@ +hB2ci+1ci+2 \ 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 @@ +ci+1D1D2 \ 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 @@ +D2D1ci+1ci+2D1D2origin \ 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 @@ +ref_axisGD1αci+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 @@ +θD2D1ci+1ci+2 \ 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 @@ +D1D2B1B2B3cici+1ci+2ci+3D1cici+1ci+2ci+3D2top viewside view \ 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]