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 get_symm_ops() converts a transform record into a 3×3 rotation matrix M and a translation vector T as numpy.ndarrays:

>>> 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:

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 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)