Structure analysis

Reduced unit cells

The reduction finds special bases of lattices. In practice, these bases are found by removing lattice centering (i.e. obtaining a primitive cell) and using a prescribed interative procedure. As it is worded in the International Tables for Crystallography A 3.1.1.4 (2016), “the reduction procedures employ metrical properties to develop a sequence of basis transformations which lead to a reduced basis and reduced cell”.

There are three popular unit cell reductions:

  • the Minkowski-Buerger reduction, which minimizes a+b+c (in special cases multiple, up to 6 different bases have the same minimal sum a+b+c),

  • the Eisenstein-Niggli reduction, which adds extra conditions to the previous one and makes the result unique,

  • the Selling-Delaunay reduction (the second name is alternatively transliterated as Delone), which minimizes a2+b2+c2+(a+b+c)2.

First names here (Minkowski, Eisenstein, Selling) belong to mathematicians working on the reduction of quadratic forms. The second names – to people applying this math to crystallography. Usually, we use only the second name. The Niggli reduction is the most popular of the three.

Niggli and Buerger reductions

Gemmi implements separately the Niggli and Buerger reductions. The procedures are iterative. Most of the unit cells from the PDB need only 1-2 iterations to get reduced (1.3 on average, not counting the normalization steps as separate iterations). On the other hand, one can always construct a primitive cell with extremely long basis vectors that would require hundreds of iterations. The Buerger reduction is simpler and faster than Niggli, but Niggli is also fast – one iteration takes less than 1μs.

Gemmi implementation is based on the algorithms published by B. Gruber in the 1970’s: Gruber, Acta Cryst. A29, 433 (1973) for the Buerger reduction, and Křivý & Gruber, Acta Cryst. A32, 297 (1976) for the Niggli reduction. Additionally, the Niggli reduction is using ε to compare numbers, as proposed by Grosse-Kunstleve et al, Acta Cryst. A60, 1 (2004).

Gruber’s algorithms use vector named G6, which is somewhat similar to the metric tensor. G6 has six elements named: A, B, C, ξ (xi), η (eta) and ζ (zeta), which correspond to:

(a2, b2, c2, 2b⋅c, 2a⋅c, 2a⋅b)

Gemmi has a class named GruberVector that contains these six numbers and reduction algorithms implemented as methods. This class can be initialized with UnitCell and SpaceGroup:

>>> cell = gemmi.UnitCell(63.78, 63.86, 124.40, 90.0, 90.0, 90.0)
>>> sg = gemmi.SpaceGroup('I 2 2 2')
>>> gv = gemmi.GruberVector(cell, sg)

or with 6-tuple corresponding to G6 of a primitive cell:

>>> g6_param = gv.parameters  # obtain such a tuple
>>> gemmi.GruberVector(g6_param)
<gemmi.GruberVector((5905.34, 5905.34, 5905.34, -7742.79, -7732.57, 3664.69))>

We can check if G6 already corresponds to a Buerger and Niggli cell:

>>> gv.is_niggli()
False
>>> gv.is_buerger()
False

We can access the G6 parameters as a tuple:

>>> gv.parameters
(5905.337, 5905.337, 5905.337, -7742.7856, -7732.5744, 3664.686)

and obtain the corresponding cell parameters (with angles in degrees):

>>> gv.cell_parameters()  # primitive cell
(76.84619053668177, 76.84619053668177, 76.84619053668177, 130.96328311485175, 130.89771578326727, 71.92353702711762)

And most importantly, we can reduce the cell. niggli_reduce() performs the Niggli reduction on G6, returning the number of iterations it took:

>>> gv.niggli_reduce()
3

Now G6 contains smaller numbers:

>>> gv
<gemmi.GruberVector((4067.89, 4078.10, 5905.34, -4078.10, -4067.89, -0.00))>

To create a new UnitCell with reduced parameters do:

>>> gemmi.UnitCell(* gv.cell_parameters())
<gemmi.UnitCell(63.78, 63.86, 76.8462, 114.551, 114.518, 90)>

or use a helper method:

>>> gv.get_cell()
<gemmi.UnitCell(63.78, 63.86, 76.8462, 114.551, 114.518, 90)>

Similarly, we can perform the Buerger reduction:

>>> gv = gemmi.GruberVector(g6_param)
>>> gv.buerger_reduce()
3

In this case both functions gave the same result.

>>> gv.is_niggli()
True
>>> gv.get_cell()
<gemmi.UnitCell(63.78, 63.86, 76.8462, 114.551, 114.518, 90)>

Functions niggli_reduce, is_niggli and is_buerger can take optional parameter epsilon (default: 1e-9) that is used for comparing numbers. Additionally, niggli_reduce can take iteration_limit (default: 100). To check how the computations would work without ε we can set it to 0:

>>> gv.is_buerger(epsilon=0)
True
>>> gv.is_niggli(epsilon=0)
False
>>> gv.niggli_reduce(epsilon=0, iteration_limit=100)
6
>>> gv.get_cell()
<gemmi.UnitCell(63.78, 63.86, 76.8462, 114.551, 114.518, 90)>

Here, the Niggli conditions were initially found not fulfilled, because one expression that should be non-negative was about -5e-13. A few extra iterations sorted it out (without any real changes), but it’s not always the case – that’s why we have iteration_limit to prevent infinite loop.

The original Křivý-Gruber algorithm doesn’t calculate the change-of-basis transformation that leads to the reduced cell. In gemmi, this transformation can be obtained as proposed in the 2004 paper of Grosse-Kunstleve et al: the change-of-basis matrix is updated in each step together with the Gruber vector. Updating this matrix makes the reduction twice slower (but it’s still in tens of ns, so it’s fast enough for any purpose). To track the change of basis, pass the following option:

>>> gv = gemmi.GruberVector(cell, sg, track_change_of_basis=True)

After the Niggli reduction, the transformation will be available in the change_of_basis property:

>>> gv.niggli_reduce()
3
>>> cob = gv.change_of_basis
>>> cob
<gemmi.Op("x-z/2,y-z/2,z/2")>

This operator transforms Niggli cell to the original cell (so it’s actually the inverse of the reduction change-of-basis):

>>> gv.get_cell().changed_basis_forward(cob, set_images=False)
<gemmi.UnitCell(63.78, 63.86, 124.4, 90, 90, 90)>

and the other way around:

>>> cell.changed_basis_backward(cob, set_images=False)
<gemmi.UnitCell(63.78, 63.86, 76.8462, 114.551, 114.518, 90)>

Currently, tracking is implemented only for the Niggli reduction, not for the Buerger reduction.

Selling-Delaunay reduction

Gemmi implementation is based on

  • section 3.1.2.3 “Delaunay reduction and standardization” in the Tables vol. A (2016),

  • Patterson & Love (1957), “Remarks on the Delaunay reduction”, Acta Cryst. 10, 111,

  • Andrews et al (2019), “Selling reduction versus Niggli reduction for crystallographic lattices”, Acta Cryst. A75, 115.

Similarly to the GruberVector, here we have a class named SellingVector that contains the six elements of S6 – the inner products among the four vectors a, b, c and d=–(a+b+c):

s23=b⋅c, s13=a⋅c, s12=a⋅b, s14=a⋅d, s24=b⋅d, s34=c⋅d.

SellingVector can be initialized with UnitCell and SpaceGroup:

>>> sv = gemmi.SellingVector(cell, sg)

or with a tuple of six numbers S6:

>>> sv.parameters
(-3871.3928, -3866.2872, 1832.343, -3871.3928, -3866.2872, 1832.343)
>>> gemmi.SellingVector(_)
<gemmi.SellingVector((-3871.39, -3866.29, 1832.34, -3871.39, -3866.29, 1832.34))>

Similarly as in the previous section, we can check if S6 already corresponds to a Delaunay cell:

>>> sv.is_reduced()
False

Each reduction step decreases Σbi2 (b1, b2, b3 and b4 are alternative symbols for a, b, c and d). The sum Σbi2 can be calculated with:

>>> sv.sum_b_squared()
23621.348

Similarly to niggli_reduce(), the Selling reduction procedure takes optional arguments epsilon and iteration_limit and returns the iteration count:

>>> sv.reduce()
2

Now we can check the result:

>>> sv
<gemmi.SellingVector((-2033.94, -2033.94, -1832.34, -2039.05, -2039.05, 0.00))>
>>> sv.is_reduced()
True
>>> sv.sum_b_squared()
19956.662

Now, the corresponding four vectors can be in any order. We may sort them so that a≤b≤c≤d:

>>> sv.sort()
>>> sv
<gemmi.SellingVector((-2039.05, -2033.94, 0.00, -2033.94, -2039.05, -1832.34))>

Finally, we can get the corresponding UnitCell:

>>> gemmi.UnitCell(* sv.cell_parameters())
<gemmi.UnitCell(63.78, 63.86, 76.8462, 114.551, 114.518, 90)>
>>> sv.get_cell()  # helper function that does the same
<gemmi.UnitCell(63.78, 63.86, 76.8462, 114.551, 114.518, 90)>

S6 can be used to calculate G6, and the other way around:

>>> sv.gruber()
<gemmi.GruberVector((4067.89, 4078.10, 5905.34, -4078.10, -4067.89, 0.00))>
>>> _.selling()
<gemmi.SellingVector((-2039.05, -2033.94, 0.00, -2033.94, -2039.05, -1832.34))>

TBC

Superposition

Gemmi includes the QCP method (Liu P, Agrafiotis DK, & Theobald DL, 2010) for superposing two lists of points in 3D. The C++ function superpose_positions() takes two arrays of positions and an optional array of weights. Before applying this function to chains it is necessary to determine pairs of corresponding atoms. Here, as a minimal example, we superpose backbone of the third residue:

>>> model = gemmi.read_structure('../tests/4oz7.pdb')[0]
>>> res1 = model['A'][2]
>>> res2 = model['B'][2]
>>> atoms = ['N', 'CA', 'C', 'O']
>>> gemmi.superpose_positions([res1.sole_atom(a).pos for a in atoms],
...                           [res2.sole_atom(a).pos for a in atoms])  
<gemmi.SupResult object at 0x...>
>>> _.rmsd
0.006558389527590043

To make it easier, we also have a higher-level function calculate_superposition() that operates on ResidueSpans. This function first performs the sequence alignment. Then the matching residues are superposed, using either all atoms in both residues, or only Cα atoms (for peptides) and P atoms (for nucleotides). Atoms that don’t have counterparts in the other span are skipped. The returned object (SupResult) contains RMSD and the transformation (rotation matrix + translation vector) that superposes the second span onto the first one.

Note that RMSD can be defined in two ways: the sum of squared deviations is divided either by 3N (PyMOL) or by N (SciPy). QCP (and gemmi) returns the former. To get the latter multiply it by √3.

Here is a usage example:

>>> model = gemmi.read_structure('../tests/4oz7.pdb')[0]
>>> polymer1 = model['A'].get_polymer()
>>> polymer2 = model['B'].get_polymer()
>>> ptype = polymer1.check_polymer_type()
>>> sup = gemmi.calculate_superposition(polymer1, polymer2, ptype, gemmi.SupSelect.CaP)
>>> sup.count  # number of atoms used
10
>>> sup.rmsd
0.1462689168993659
>>> sup.transform.mat
<gemmi.Mat33 [-0.0271652, 0.995789, 0.0875545]
             [0.996396, 0.034014, -0.0777057]
             [-0.0803566, 0.085128, -0.993124]>
>>> sup.transform.vec
<gemmi.Vec3(-17.764, 16.9915, -1.77262)>

The arguments to calculate_superposition() are:

  • two ResidueSpans,

  • polymer type (to avoid determining it when it’s already known). The information whether it’s protein or nucleic acid is used during sequence alignment (to detect gaps between residues in the polymer – it helps in rare cases when the sequence alignment alone is ambiguous), and it decides whether to use Cα or P atoms (see the next point),

  • atom selection: one of SupSelect.CaP (only Cα or P atoms), SupSelect.MainChain or SupSelect.All (all atoms),

  • (optionally) altloc – the conformer choice. By default, atoms with non-blank altloc are ignored. With altloc=’A’, only the A conformer is considered (atoms with altloc either blank or A). Etc.

  • (optionally) trim_cycles (default: 0) – maximum number of outlier rejection cycles. This option was inspired by the align command in PyMOL.

    >>> sr = gemmi.calculate_superposition(polymer1, polymer2, ptype,
...                                    gemmi.SupSelect.All, trim_cycles=5)
>>> sr.count
73
>>> sr.rmsd
0.18315488879658484

  • (optionally) trim_cutoff (default: 2.0) – outlier rejection cutoff in RMSD,

To calculate current RMSD between atoms (without superposition) use function calculate_current_rmsd() that takes the same arguments except the ones for trimming:

>>> gemmi.calculate_current_rmsd(polymer1, polymer2, ptype, gemmi.SupSelect.CaP).rmsd
19.660883858565462

The calculated superposition can be applied to a span of residues, changing the atomic positions in-place:

>>> polymer2[2].sole_atom('CB')  # before
<gemmi.Atom CB at (-30.3, -10.6, -11.6)>
>>> polymer2.transform_pos_and_adp(sup.transform)
>>> polymer2[2].sole_atom('CB')  # after
<gemmi.Atom CB at (-28.5, -12.6, 11.2)>
>>> # it is now nearby the corresponding atom in chain A:
>>> polymer1[2].sole_atom('CB')
<gemmi.Atom CB at (-28.6, -12.7, 11.3)>

Selections

Gemmi selection syntax is based on the selection syntax from MMDB, which is sometimes called CID (Coordinate ID). The MMDB syntax is described at the bottom of the pdbcur documentation.

The selection has a form of /-separated parts: /models/chains/residues/atoms. Empty parts can be omitted when it’s not ambiguous. Gemmi (but not MMDB) can take additional properties added at the end after a semicolon (;).

Let us go through the individual filters first:

  • /1 – selects model 1 (if the PDB file doesn’t have MODEL records, it is assumed that the only model is model 1).

  • //D (or just D) – selects chain D.

  • ///10-30 (or 10-30) – residues with sequence IDs from 10 to 30.

  • ///10A-30A (or 10A-30A or ///10.A-30.A or 10.A-30.A) – sequence ID can include insertion code. The MMDB syntax has dot between sequence sequence number and insertion code. In Gemmi the dot is optional.

  • ///(ALA) (or (ALA)) – selects residues with a given name.

  • ////CB (or CB:* or CB[*]) – selects atoms with a given name.

  • ////[P] (or just [P]) – selects phosphorus atoms.

  • ////:B (or :B) – selects atoms with altloc B.

  • ////;q<0.5 (or ;q<0.5) – selects atoms with occupancy below 0.5 (inspired by PyMOL, where it’d be q<0.5).

  • ////;b>40 (or ;b>40) – selects atoms with the isotropic B-factor above a given value.

  • * – selects all atoms.

Note that the chain name and altloc can be an empty. The syntax supports also comma-separated lists and negations with !:

  • (!ALA) – all residues but alanine,

  • [C,N,O] – all C, N and O atoms,

  • [!C,N,O] – all atoms except C, N and O,

  • :,A – altloc either empty or A (which makes one conformation),

  • /1/A,B/20-40/CA[C]:,A – multiple selection criteria, all of them must be fulfilled.

Incompatibility with MMDB. In MMDB, if the chemical element is specified (e.g. [C] or [*]), the alternative location indicator defaults to “” (no altloc), rather than to “*” (any altloc). This might be surprising. In Gemmi, if ‘:’ is absent the altloc is not checked (“*”).

Note: the selections in Gemmi are not widely used yet and the API may evolve.

A selection is is a standalone object with a list of filters that can be applied to any Structure, Model or its part. Empty selection matches all atoms:

>>> sel = gemmi.Selection()  # empty - no filters

Selection initialized with a string parses the string and creates corresponding filters:

>>> # select all Cl atoms
>>> sel = gemmi.Selection('[CL]')

The selection can then be used on any structure. A helper function first() returns the first matching atom:

>>> st = gemmi.read_structure('../tests/1pfe.cif.gz')
>>> # get the first result as pointer to model and CRA (chain, residue, atom)
>>> sel.first(st)
(<gemmi.Model 1 with 2 chain(s)>, <gemmi.CRA A/CL 20/CL>)

Function str() creates a string from the selection:

>>> sel = gemmi.Selection('A/1-4/N9')
>>> sel.str()
'//A/1.-4./N9'

The Selection objects has methods for iterating over the selected items in the hierarchy:

>>> for model in sel.models(st):
...     print('Model', model.name)
...     for chain in sel.chains(model):
...         print('-', chain.name)
...         for residue in sel.residues(chain):
...             print('   -', str(residue))
...             for atom in sel.atoms(residue):
...                 print('          -', atom.name)
...
Model 1
- A
   - 1(DG)
          - N9
   - 2(DC)
   - 3(DG)
          - N9
   - 4(DT)

Selection can be used to create a new structure (or model) with a copy of the selection. In this example, we copy alpha-carbon atoms:

>>> st = gemmi.read_structure('../tests/1orc.pdb')
>>> st[0].count_atom_sites()
559
>>> selection = gemmi.Selection('CA[C]')

>>> # create a new structure
>>> ca_st = selection.copy_structure_selection(st)
>>> ca_st[0].count_atom_sites()
64

>>> # create a new model
>>> ca_model = selection.copy_model_selection(st[0])
>>> ca_model.count_atom_sites()
64

Selection can also be used to remove atoms. In this example we remove atoms with B-factor above 50:

>>> sel = gemmi.Selection(';b>50')
>>> sel.remove_selected(ca_st)
>>> ca_st[0].count_atom_sites()
61
>>> sel.remove_selected(ca_model)
>>> ca_model.count_atom_sites()
61

We can also do the opposite and remove atoms that are not selected:

>>> sel.remove_not_selected(ca_model)
>>> ca_model.count_atom_sites()
0

Each residue and atom has a flag that can be set manually and used to create a selection. In this example we select residues in the radius of 8Å from a selected point:

>>> selected_point = gemmi.Position(20, 40, 30)
>>> ns = gemmi.NeighborSearch(st[0], st.cell, 8.0).populate()
>>> # First, a flag is set for neigbouring residues.
>>> for mark in ns.find_atoms(selected_point):
...     mark.to_cra(st[0]).residue.flag = 's'
>>> # Then, we select residues with this flag.
>>> selection = gemmi.Selection().set_residue_flags('s')
>>> # Next, we can use this selection.
>>> selection.copy_model_selection(st[0]).count_atom_sites()
121

Note: NeighborSearch searches for atoms in all symmetry images. This is why it takes UnitCell as a parameter. To search only in atoms directly listed in the file pass empty cell (gemmi.UnitCell()).

Instead of the whole residues, we can select atoms. Here, we select atoms in the radius of 8Å from a selected point:

>>> # selected_point and ns are reused from the previous example
>>> # First, a flag is set for neigbouring atoms.
>>> for mark in ns.find_atoms(selected_point):
...     mark.to_cra(st[0]).atom.flag = 's'
>>> # Then, we select atoms with this flag.
>>> selection = gemmi.Selection().set_atom_flags('s')
>>> # Next, we can use this selection.
>>> selection.copy_model_selection(st[0]).count_atom_sites()
59

Graph analysis

The graph algorithms in Gemmi are limited to finding the shortest path between atoms (bonds = graph edges). This part of the library is not documented yet.

The rest of this section shows how to use Gemmi together with external graph analysis libraries to analyse the similarity of chemical molecules. To do this, first we set up a graph corresponding to the molecule.

Here we show how it can be done in the Boost Graph Library.

#include <boost/graph/graph_traits.hpp>
#include <boost/graph/adjacency_list.hpp>
#include <gemmi/cif.hpp>             // for cif::read_file
#include <gemmi/chemcomp.hpp>        // for ChemComp, make_chemcomp_from_block

struct AtomVertex {
  gemmi::Element el = gemmi::El::X;
  std::string name;
};

struct BondEdge {
  gemmi::BondType type;
};

using Graph = boost::adjacency_list<boost::vecS, boost::vecS,
                                    boost::undirectedS,
                                    AtomVertex, BondEdge>;

Graph make_graph(const gemmi::ChemComp& cc) {
  Graph g(cc.atoms.size());
  for (size_t i = 0; i != cc.atoms.size(); ++i) {
    g[i].el = cc.atoms[i].el;
    g[i].name = cc.atoms[i].id;
  }
  for (const gemmi::Restraints::Bond& bond : cc.rt.bonds) {
    int n1 = cc.get_atom_index(bond.id1.atom);
    int n2 = cc.get_atom_index(bond.id2.atom);
    boost::add_edge(n1, n2, BondEdge{bond.type}, g);
  }
  return g;
}

And here we use NetworkX in Python:

>>> import networkx

>>> G = networkx.Graph()
>>> block = gemmi.cif.read('../tests/SO3.cif')[-1]
>>> so3 = gemmi.make_chemcomp_from_block(block)
>>> for atom in so3.atoms:
...     G.add_node(atom.id, Z=atom.el.atomic_number)
...
>>> for bond in so3.rt.bonds:
...     G.add_edge(bond.id1.atom, bond.id2.atom)  # ignoring bond type
...

To show a quick example, let us count automorphisms of SO3:

>>> import networkx.algorithms.isomorphism as iso
>>> GM = iso.GraphMatcher(G, G, node_match=iso.categorical_node_match('Z', 0))
>>> # expecting 3! automorphisms (permutations of the three oxygens)
>>> sum(1 for _ in GM.isomorphisms_iter())
6

With a bit more of code we could perform a real cheminformatics task.

Graph isomorphism

In this example we use Python NetworkX to compare molecules from the Refmac monomer library with Chemical Component Dictionary (CCD) from PDB. The same could be done with other graph analysis libraries, such as Boost Graph Library, igraph, etc.

The program below takes compares specified monomer cif files with corresponding CCD entries. Hydrogens and bond types are ignored. It takes less than half a minute to go through the 25,000 monomer files distributed with CCP4 (as of Oct 2018), so we do not try to optimize the program.

# Compares graphs of molecules from cif files (Refmac dictionary or similar)
# with CCD entries.

import sys
import networkx
from networkx.algorithms import isomorphism
import gemmi

CCD_PATH = 'components.cif.gz'

def graph_from_chemcomp(cc):
    G = networkx.Graph()
    for atom in cc.atoms:
        G.add_node(atom.id, Z=atom.el.atomic_number)
    for bond in cc.rt.bonds:
        G.add_edge(bond.id1.atom, bond.id2.atom)
    return G

def compare(cc1, cc2):
    s1 = {a.id for a in cc1.atoms}
    s2 = {a.id for a in cc2.atoms}
    b1 = {b.lexicographic_str() for b in cc1.rt.bonds}
    b2 = {b.lexicographic_str() for b in cc2.rt.bonds}
    if s1 == s2 and b1 == b2:
        #print(cc1.name, "the same")
        return
    G1 = graph_from_chemcomp(cc1)
    G2 = graph_from_chemcomp(cc2)
    node_match = isomorphism.categorical_node_match('Z', 0)
    GM = isomorphism.GraphMatcher(G1, G2, node_match=node_match)
    if GM.is_isomorphic():
        print(cc1.name, 'is isomorphic')
        # we could use GM.match(), but here we try to find the shortest diff
        short_diff = None
        for n, mapping in enumerate(GM.isomorphisms_iter()):
            diff = {k: v for k, v in mapping.items() if k != v}
            if short_diff is None or len(diff) < len(short_diff):
                short_diff = diff
            if n == 10000:  # don't spend too much here
                print(' (it may not be the simplest isomorphism)')
                break
        for id1, id2 in short_diff.items():
            print('\t', id1, '->', id2)
    else:
        print(cc1.name, 'differs')
        if s2 - s1:
            print('\tmissing:', ' '.join(s2 - s1))
        if s1 - s2:
            print('\textra:  ', ' '.join(s1 - s2))

def main():
    ccd = gemmi.cif.read(CCD_PATH)
    absent = 0
    for f in sys.argv[1:]:
        block = gemmi.cif.read(f)[-1]
        cc1 = gemmi.make_chemcomp_from_block(block)
        try:
            block2 = ccd[cc1.name]
        except KeyError:
            absent += 1
            #print(cc1.name, 'not in CCD')
            continue
        cc2 = gemmi.make_chemcomp_from_block(block2)
        cc1.remove_hydrogens()
        cc2.remove_hydrogens()
        compare(cc1, cc2)
    if absent != 0:
        print(absent, 'of', len(sys.argv) - 1, 'monomers not found in CCD')

main()

If we run it on monomers that start with M we get:

$ examples/ccd_gi.py $CLIBD_MON/m/*.cif
M10 is isomorphic
       O9 -> O4
       O4 -> O9
MK8 is isomorphic
       O2 -> OXT
MMR differs
      missing: O12 O4
2 of 821 monomers not found in CCD

So in M10 the two atoms marked green are swapped:

_images/M10-isomorphism.png

(The image was generated in NGL and compressed with Compress-Or-Die.)

Substructure matching

Now a little script to illustrate subgraph isomorphism. The script takes a (three-letter-)code of a molecule that is to be used as a pattern and finds CCD entries that contain such a substructure. As in the previous example, hydrogens and bond types are ignored.

# List CCD entries that contain the specified entry as a substructure.
# Ignoring hydrogens and bond types.

import sys
import networkx
from networkx.algorithms import isomorphism
import gemmi

CCD_PATH = 'components.cif.gz'

def graph_from_block(block):
    cc = gemmi.make_chemcomp_from_block(block)
    cc.remove_hydrogens()
    G = networkx.Graph()
    for atom in cc.atoms:
        G.add_node(atom.id, Z=atom.el.atomic_number)
    for bond in cc.rt.bonds:
        G.add_edge(bond.id1.atom, bond.id2.atom)
    return G

def main():
    assert len(sys.argv) == 2, "Usage: ccd_subgraph.py three-letter-code"
    ccd = gemmi.cif.read(CCD_PATH)
    entry = sys.argv[1]
    pattern = graph_from_block(ccd[entry])
    pattern_nodes = networkx.number_of_nodes(pattern)
    pattern_edges = networkx.number_of_edges(pattern)
    node_match = isomorphism.categorical_node_match('Z', 0)
    for block in ccd:
        G = graph_from_block(block)
        GM = isomorphism.GraphMatcher(G, pattern, node_match=node_match)
        if GM.subgraph_is_isomorphic():
            print(block.name, '\t +%d nodes, +%d edges' % (
                networkx.number_of_nodes(G) - pattern_nodes,
                networkx.number_of_edges(G) - pattern_edges))

main()

Let us check what entries have HEM as a substructure:

$ examples/ccd_subgraph.py HEM
1FH    +6 nodes, +7 edges
2FH    +6 nodes, +7 edges
4HE    +7 nodes, +8 edges
522    +2 nodes, +2 edges
6CO    +6 nodes, +7 edges
6CQ    +7 nodes, +8 edges
89R    +3 nodes, +3 edges
CLN    +1 nodes, +2 edges
DDH    +2 nodes, +2 edges
FEC    +6 nodes, +6 edges
HAS    +22 nodes, +22 edges
HCO    +1 nodes, +1 edges
HDM    +2 nodes, +2 edges
HEA    +17 nodes, +17 edges
HEB    +0 nodes, +0 edges
HEC    +0 nodes, +0 edges
HEM    +0 nodes, +0 edges
HEO    +16 nodes, +16 edges
HEV    +2 nodes, +2 edges
HP5    +2 nodes, +2 edges
ISW    +0 nodes, +0 edges
MH0    +0 nodes, +0 edges
MHM    +0 nodes, +0 edges
N7H    +3 nodes, +3 edges
NTE    +3 nodes, +3 edges
OBV    +14 nodes, +14 edges
SRM    +20 nodes, +20 edges
UFE    +18 nodes, +18 edges

Maximum common subgraph

In this example we use McGregor’s algorithm implemented in the Boost Graph Library to find maximum common induced subgraph. We call the MCS searching function with option only_connected_subgraphs=true, which has obvious meaning and can be changed if needed.

To illustrate this example, we compare ligands AUD and LSA:

_images/aud_lsa.png

The whole code is in examples/with_bgl.cpp. The same file has also examples of using the BGL implementation of VF2 to check graph and subgraph isomorphisms.

// We ignore hydrogens here.
// Example output:
//   $ time with_bgl -c monomers/a/AUD.cif monomers/l/LSA.cif
//   Searching largest subgraphs of AUD and LSA (10 and 12 vertices)...
//   Maximum connected common subgraph has 7 vertices:
//     SAA -> S10
//     OAG -> O12
//     OAH -> O11
//     NAF -> N9
//     CAE -> C7
//     OAI -> O8
//     CAD -> C1
//   Maximum connected common subgraph has 7 vertices:
//     SAA -> S10
//     OAG -> O11
//     OAH -> O12
//     NAF -> N9
//     CAE -> C7
//     OAI -> O8
//     CAD -> C1
//   
//   real	0m0.012s
//   user	0m0.008s
//   sys	0m0.004s

struct PrintCommonSubgraphCallback {
  Graph g1, g2;
  template <typename CorrespondenceMap1To2, typename CorrespondenceMap2To1>
  bool operator()(CorrespondenceMap1To2 map1, CorrespondenceMap2To1,
                  typename GraphTraits::vertices_size_type subgraph_size) {
    std::cout << "Maximum connected common subgraph has " << subgraph_size
              << " vertices:\n";
    for (auto vp = boost::vertices(g1); vp.first != vp.second; ++vp.first) {
      Vertex v1 = *vp.first;
      Vertex v2 = boost::get(map1, v1);
      if (v2 != GraphTraits::null_vertex())
        std::cout << "  " << g1[v1].name << " -> " << g2[v2].name << std::endl;
    }
    return true;
  }
};

void find_common_subgraph(const char* cif1, const char* cif2) {
  gemmi::ChemComp cc1 = make_chemcomp(cif1).remove_hydrogens();
  Graph graph1 = make_graph(cc1);
  gemmi::ChemComp cc2 = make_chemcomp(cif2).remove_hydrogens();
  Graph graph2 = make_graph(cc2);
  std::cout << "Searching largest subgraphs of " << cc1.name << " and "
      << cc2.name << " (" << cc1.atoms.size() << " and " << cc2.atoms.size()
      << " vertices)..." << std::endl;
  mcgregor_common_subgraphs_maximum_unique(
      graph1, graph2,
      get(boost::vertex_index, graph1), get(boost::vertex_index, graph2),
      [&](Edge a, Edge b) { return graph1[a].type == graph2[b].type; },
      [&](Vertex a, Vertex b) { return graph1[a].el == graph2[b].el; },
      /*only_connected_subgraphs=*/ true,
      PrintCommonSubgraphCallback{graph1, graph2});
}

Torsion angles

This section presents functions dedicated to calculation of the dihedral angles φ (phi), ψ (psi) and ω (omega) of the protein backbone. These functions are built upon the more general calculate_dihedral function, introduced in the section about coordinates, which takes four points in the space as arguments.

calculate_omega() calculates the ω angle, which is usually around 180°:

>>> from math import degrees
>>> chain = gemmi.read_structure('../tests/5cvz_final.pdb')[0]['A']
>>> degrees(gemmi.calculate_omega(chain[0], chain[1]))
159.9092215006572
>>> for res in chain[:5]:
...     next_res = chain.next_residue(res)
...     if next_res:
...         omega = gemmi.calculate_omega(res, next_res)
...         print(res.name, degrees(omega))
...
ALA 159.9092215006572
ALA -165.26874513591127
ALA -165.85686681169696
THR -172.99968385093572
SER 176.74223937657652

The φ and ψ angles are often used together, so they are calculated in one function calculate_phi_psi():

>>> for res in chain[:5]:
...     prev_res = chain.previous_residue(res)
...     next_res = chain.next_residue(res)
...     phi, psi = gemmi.calculate_phi_psi(prev_res, res, next_res)
...     print('%s %8.2f %8.2f' % (res.name, degrees(phi), degrees(psi)))
...
ALA      nan   106.94
ALA  -116.64    84.57
ALA   -45.57   127.40
THR   -62.01   147.45
SER   -92.85   161.53

In C++ these functions can be found in gemmi/calculate.hpp.

The torsion angles φ and ψ can be visualized on the Ramachandran plot. Let us plot angles from all PDB entries with the resolution higher than 1.5A. Usually, glycine, proline and the residue preceding proline (pre-proline) are plotted separately. Here, we will exclude pre-proline and make separate plot for each amino acid. So first, we calculate angles and save φ,ψ pairs in a set of files – one file per residue.

import sys
from math import degrees
import gemmi

ramas = {aa: [] for aa in [
    'LEU', 'ALA', 'GLY', 'VAL', 'GLU', 'SER', 'LYS', 'ASP', 'THR', 'ILE',
    'ARG', 'PRO', 'ASN', 'PHE', 'GLN', 'TYR', 'HIS', 'MET', 'CYS', 'TRP']}

for path in gemmi.CoorFileWalk(sys.argv[1]):
    st = gemmi.read_structure(path)
    if 0.1 < st.resolution < 1.5:
        model = st[0]
        for chain in model:
            for res in chain.get_polymer():
                # previous_residue() and next_residue() return previous/next
                # residue only if the residues are bonded. Otherwise -- None.
                prev_res = chain.previous_residue(res)
                next_res = chain.next_residue(res)
                if prev_res and next_res and next_res.name != 'PRO':
                    v = gemmi.calculate_phi_psi(prev_res, res, next_res)
                    try:
                        ramas[res.name].append(v)
                    except KeyError:
                        pass

# Write data to files
for aa, data in ramas.items():
    with open('ramas/' + aa + '.tsv', 'w') as f:
        for phi, psi in data:
            f.write('%.4f\t%.4f\n' % (degrees(phi), degrees(psi)))

The script above works with coordinate files in any of the formats supported by gemmi (PDB, mmCIF, mmJSON). As of 2019, processing a local copy of the PDB archive in the PDB format takes about 20 minutes.

In the second step we plot the data points with Matplotlib. We use a script that can be found in examples/rama_plot.py. Six of the resulting plots are shown here (click to enlarge):

_images/ramachandran-per-aa.png

Topology

A macromolecular refinement program typically starts from reading a coordinate file and a monomer library. The monomer library specifies restraints (bond distances, angles, …) in monomers as well as modifications introduced by links between monomers.

The coordinates and restraints are combined into what we call here a topology. It contains restraints applied to the model. A monomer library may specify angle CD2-CE2-NE1 in TRP. In contrast, the topology specifies angles between concrete atoms (say, angle #721-#720-#719).

Together with preparing a topology, macromolecular programs (in particular, Refmac) may also add or shift hydrogens (to the riding positions) and reorder atoms. In Python we have one function that does it all:

gemmi.prepare_topology(st: gemmi.Structure,
                       monlib: gemmi.MonLib,
                       model_index: int = 0,
                       h_change: gemmi.HydrogenChange = HydrogenChange.NoChange,
                       reorder: bool = False,
                       warnings: object = None,
                       ignore_unknown_links: bool = False) -> gemmi.Topo

where

  • monlib is an instance of an undocumented MonLib class. For now, here is an example how to read the CCP4 monomer library (a.k.a Refmac dictionary):

    monlib_path = os.environ['CCP4'] + '/lib/data/monomers'
resnames = st[0].get_all_residue_names()
monlib = gemmi.read_monomer_lib(monlib_path, resnames)

  • h_change is one of:

    • HydrogenChange.NoChange – no change,

    • HydrogenChange.Shift – shift existing hydrogens to ideal (riding) positions,

    • HydrogenChange.Remove – remove all H and D atoms,

    • HydrogenChange.ReAdd – discard and re-create hydrogens in ideal positions,

    • HydrogenChange.ReAddButWater – the same, but doesn’t add H in waters,

  • reorder – changes the order of atoms inside each residue to match the order in the corresponding monomer cif file,

  • warnings – by default, exception is raised when a chemical component is missing in the monomer library, or when link is missing, or the hydrogen adding procedure comes across an unexpected configuration. You can set warnings=sys.stderr to only print a warning to stderr and continue. sys.stderr can be replaced with any object that has methods write(str) and flush().

If hydrogen position is not uniquely determined its occupancy is set to zero.

TBC

Local copy of the PDB archive

Some of the examples in this documentation work with a local copy of the Protein Data Bank archive. This subsection describes the assumed setup.

Like in BioJava, we assume that the $PDB_DIR environment variable points to a directory that contains structures/divided/mmCIF – the same arrangement as on the PDB’s FTP server.

$ cd $PDB_DIR
$ du -sh structures/*/*  # as of Jun 2017
34G    structures/divided/mmCIF
25G    structures/divided/pdb
101G   structures/divided/structure_factors
2.6G   structures/obsolete/mmCIF

A traditional way to keep an up-to-date local archive is to rsync it once a week:

#!/bin/sh -x
set -u  # PDB_DIR must be defined
rsync_subdir() {
  mkdir -p "$PDB_DIR/$1"
  # Using PDBe (UK) here, can be replaced with RCSB (USA) or PDBj (Japan),
  # see https://www.wwpdb.org/download/downloads
  rsync -rlpt -v -z --delete \
      rsync.ebi.ac.uk::pub/databases/pdb/data/$1/ "$PDB_DIR/$1/"
}
rsync_subdir structures/divided/mmCIF
#rsync_subdir structures/obsolete/mmCIF
#rsync_subdir structures/divided/pdb
#rsync_subdir structures/divided/structure_factors

Gemmi has a helper function for using the local archive copy. It takes a PDB code (case insensitive) and a symbol denoting what file is requested: P for PDB, M for mmCIF, S for SF-mmCIF.

>>> os.environ['PDB_DIR'] = '/copy'
>>> gemmi.expand_if_pdb_code('1ABC', 'P') # PDB file
'/copy/structures/divided/pdb/ab/pdb1abc.ent.gz'
>>> gemmi.expand_if_pdb_code('1abc', 'M') # mmCIF file
'/copy/structures/divided/mmCIF/ab/1abc.cif.gz'
>>> gemmi.expand_if_pdb_code('1abc', 'S') # SF-mmCIF file
'/copy/structures/divided/structure_factors/ab/r1abcsf.ent.gz'

If the first argument is not in the PDB code format (4 characters for now) the function returns the first argument.

>>> arg = 'file.cif'
>>> gemmi.is_pdb_code(arg)
False
>>> gemmi.expand_if_pdb_code(arg, 'M')
'file.cif'

Multiprocessing

(Python-specific)

Most of the gemmi objects cannot be pickled. Therefore, they cannot be passed between processes when using the multiprocessing module. Currently, the only picklable classes (with protocol >= 2) are: UnitCell and SpaceGroup.

Usually, it is possible to organize multiprocessing in such a way that gemmi objects are not passed between processes. The example script below traverses subdirectories and asynchronously analyses coordinate files. It uses 4 worker processes in parallel. The processes get file path and return a tuple.

import multiprocessing
import sys
import gemmi

def f(path):
    st = gemmi.read_structure(path)
    weight = st[0].calculate_mass()
    return (st.name, weight)

def main(top_dir):
    with multiprocessing.Pool(processes=4) as pool:
        it = pool.imap_unordered(f, gemmi.CoorFileWalk(top_dir))
        for (name, weight) in it:
            print('%s %.1f kDa' % (name, weight / 1000))

if __name__ == '__main__':
    main(sys.argv[1])