Biological assemblies and symmetry
==================================

The biologically relevant molecule is often larger than the contents of the
asymmetric unit — it is generated by applying a set of rotation/translation
operators to some subset of chains.  The PDB records this under ``REMARK 350``,
and pidibble parses each operator into its own indexed sub-record.

>>> from pidibble.pdbparse import PDBParser
>>> p = PDBParser(source_db='rcsb', source_id='4zmj').parse()

Assembly transforms
-------------------

Each transform for biological assembly 1 is a record keyed
``REMARK.350.BIOMOLECULE1.TRANSFORM<n>``:

>>> t1 = p.parsed['REMARK.350.BIOMOLECULE1.TRANSFORM1']

The ``header`` attribute of the *first* transform of an assembly lists the chains
the operators apply to:

>>> t1.header
['G', 'B', 'A', 'C', 'D']

Turning a transform into matrices
---------------------------------

The helper :func:`~pidibble.pdbparse.get_symm_ops` converts a transform record
into a 3×3 rotation matrix ``M`` and a translation vector ``T`` as
:class:`numpy.ndarray`\ s:

>>> from pidibble.pdbparse import get_symm_ops
>>> M, T = get_symm_ops(t1)
>>> M
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])
>>> T
array([0., 0., 0.])

The first transform is the identity (the asymmetric unit itself); the others
build out the assembly.  For 4ZMJ — a trimer — transforms 2 and 3 are the ±120°
rotations of a three-fold axis:

>>> M, T = get_symm_ops(p.parsed['REMARK.350.BIOMOLECULE1.TRANSFORM2'])
>>> M
array([[-0.5     , -0.866025,  0.      ],
       [ 0.866025, -0.5     ,  0.      ],
       [ 0.      ,  0.      ,  1.      ]])
>>> T
array([107.18   , 185.64121,   0.     ])

Applying an operator to coordinates is then just ``M @ xyz + T``:

.. code-block:: python

   import numpy as np

   assembly = []
   for key in sorted(k for k in p.parsed if k.startswith('REMARK.350.BIOMOLECULE1.TRANSFORM')):
       M, T = get_symm_ops(p.parsed[key])
       for atom in p.parsed['ATOM']:
           if atom.residue.chainID in p.parsed['REMARK.350.BIOMOLECULE1.TRANSFORM1'].header:
               xyz = np.array([atom.x, atom.y, atom.z])
               assembly.append(M @ xyz + T)

Crystallographic symmetry
-------------------------

Crystallographic symmetry operators (the space-group operations under
``REMARK 290``) are parsed the same way, into
``REMARK.290.CRYSTSYMMTRANS.<n>`` sub-records that also work with
:func:`~pidibble.pdbparse.get_symm_ops`:

>>> ops = sorted(k for k in p.parsed if k.startswith('REMARK.290.CRYSTSYMMTRANS'))
>>> len(ops)
6
>>> M, T = get_symm_ops(p.parsed[ops[0]])    # the first operator is the identity

The unit cell and space group themselves are in ``CRYST1`` (with fields ``a``,
``b``, ``c``, ``alpha``, ``beta``, ``gamma``, ``sGroup`` and ``z``):

>>> c = p.parsed['CRYST1']
>>> c.sGroup
'P 63'
>>> c.a, c.b, c.c
(107.18, 107.18, 103.06)