Gemmi program

The library comes with a command-line program which is also named gemmi; running a program is easier than calling a library function.

This program is actually a set of small programs, each of them corresponding to a subcommand:

$ gemmi -h
gemmi 0.3.2dev
Command-line utility that accompanies the GEMMI library,
which is a joint project of CCP4 and Global Phasing Ltd.
Licence: Mozilla Public License 2.0.
Copyright 2017-2019 Global Phasing Ltd.
https://github.com/project-gemmi/gemmi

Usage: gemmi [--version] [--help] <command> [<args>]

Commands:
 blobs         list unmodelled electron density blobs
 cif2mtz       convert structure factor mmCIF to MTZ
 contact       searches for contacts (neighbouring atoms)
 contents      info about content of a coordinate file (pdb, mmCIF, ...)
 convert       convert file (CIF - JSON, mmCIF - PDB) or modify structure
 grep          search for tags in CIF file(s)
 h             add or remove hydrogen atoms
 map           print info or modify a CCP4 map
 map2sf        transform CCP4 map to map coefficients (in MTZ or mmCIF)
 mask          make mask in the CCP4 format
 mondiff       compare two monomer CIF files
 mtz           print info about MTZ reflection file
 mtz2cif       convert MTZ to structure factor mmCIF
 residues      list residues from a coordinate file
 rmsz          validate geometry using monomer library
 seq           sequence alignment (global, pairwise, affine gap penalty)
 sf2map        transform map coefficients (from MTZ or mmCIF) to map
 sg            info about space groups
 validate      validate CIF 1.1 syntax
 wcn           calculate local density / contact numbers (WCN, CN, ACN, LDM)

validate

A CIF validator. Apart from checking the syntax it can check most of the rules imposed by DDL1 and DDL2 dictionaries.

$ gemmi validate -h
Usage: gemmi validate [options] FILE [...]

Options:
  -h, --help      Print usage and exit.
  -V, --version   Print version and exit.
  -f, --fast      Syntax-only check.
  -s, --stat      Show token statistics
  --verbose       Verbose output.
  -q, --quiet     Show only errors.
  -d, --ddl=PATH  DDL for validation.
  -m, --monomer   Extra checks for Refmac dictionary files.

grep

Searches for a specified tag in CIF files and prints the associated values, one value per line:

$ gemmi grep _refine.ls_R_factor_R_free 5fyi.cif.gz
5FYI:0.2358
$ gemmi grep _refine.ls_R_factor_R_free mmCIF/mo/?moo.cif.gz
1MOO:0.177
3MOO:0.21283
4MOO:0.22371
5MOO:0.1596
5MOO:0.1848
$ gemmi grep -b _software.name 5fyi.cif.gz
DIMPLE
PHENIX

Some of the command-line options correspond to the options of GNU grep (-c, -l, -H, -n). As with other utilities, option --help shows the usage:

$ gemmi grep -h
Usage: gemmi grep [options] TAG FILE_OR_DIR_OR_PDBID[...]
       gemmi grep -f FILE [options] TAG
Search for TAG in CIF files.

Options:
  -h, --help               display this help and exit
  -V, --version            display version information and exit
  -f, --file=FILE          obtain file (or PDB ID) list from FILE
  -m, --max-count=NUM      print max NUM values per file
  -O, --one-block          optimize assuming one block per file
  -a, --and=tag            Append delimiter (default ';') and the tag value
  -d, --delimiter=DELIM    CSV-like output with specified delimiter
  -n, --line-number        print line number with output lines
  -H, --with-filename      print the file name for each match
  -b, --no-blockname       suppress the block name on output
  -t, --with-tag           print the tag name for each match
  -l, --files-with-tag     print only names of files with the tag
  -L, --files-without-tag  print only names of files without the tag
  -c, --count              print only a count of values per block or file
  -r, --recursive          ignored (directories are always recursed)
  -w, --raw                include '?', '.', and string quotes
  -s, --summarize          display joint statistics for all files

This is a minimalistic program designed to be used together with Unix text-processing utilities. For example, it cannot filter values itself, but one may use grep:

$ gemmi grep _pdbx_database_related.db_name /pdb/mmCIF/aa/* | grep EMDB
4AAS:EMDB
5AA0:EMDB

Gemmi-grep tries to be simple to use like Unix grep, but at the same time it is aware of the CIF syntax rules. In particular, gemmi grep _one will give the same output for both _one 1 and loop_ _one _two 1 2. This is helpful in surprising corner cases. For example, when a PDB entry has two Rfree values (see the 5MOO example above).

Gemmi-grep does not support regular expression, only globbing (wildcards): ? represents any single character, * represents any number of characters (including zero). When using wildcards you may also want to use the -t option which prints the tag:

$ gemmi grep -t _*free 3gem.cif
3GEM:[_refine.ls_R_factor_R_free] 0.182
3GEM:[_refine.ls_percent_reflns_R_free] 5.000
3GEM:[_refine.ls_number_reflns_R_free] 3951
3GEM:[_refine.correlation_coeff_Fo_to_Fc_free] 0.952
3GEM:[_refine_ls_shell.R_factor_R_free] 0.272
3GEM:[_refine_ls_shell.number_reflns_R_free] 253

Let say we want to find extreme unit cell angles in the PDB. _cell.angle_*a will match _cell.angle_alpha as well as beta and gamma, but not _cell.angle_alpha_esd etc.

$ gemmi grep -d' ' _cell.angle_*a /pdb/mmCIF/ | awk '$2 < 50 || $2 > 140 { print $0; }'
4AL2 144.28
2EX3 45.40
2GMV 145.09
4NX1 140.060
4OVP 140.070
1SPG 141.90
2W1I 146.58

The option -O is used to make gemmi-grep faster. With this option the program finds only the first occurence of the tag in file. Note that if the file has only one block (like mmCIF coordinate files) and the tag is specified without wildcards then we cannot have more than one match anyway.

Searching the whole compressed mmCIF archive from the PDB (35GB of gzipped files) should take on an average computer between 10 and 30 minutes, depending where the searched tag is located. This is much faster than with other CIF parsers (to my best knowledge) and it makes the program useful for ad-hoc PDB statistics:

$ gemmi grep -O -b _entity_poly.type /pdb/mmCIF | sort | uniq -c
      1 cyclic-pseudo-peptide
      4 other
      2 peptide nucleic acid
   9905 polydeoxyribonucleotide
    156 polydeoxyribonucleotide/polyribonucleotide hybrid
     57 polypeptide(D)
 168923 polypeptide(L)
   4559 polyribonucleotide
     18 polysaccharide(D)

Option -c counts the values in each block or file. As an example we may check which entries have the biggest variety of chemical components (spoiler: ribosomes):

$ gemmi grep -O -c _chem_comp.id /pdb/mmCIF | sort -t: -k2 -nr | head
5J91:58
5J8A:58
5J7L:58
5J5B:58
4YBB:58
5JC9:57
5J88:57
5IT8:57
5IQR:50
5AFI:50

Going back to moo, we may want to know to what experimental method the Rfree values correspond:

$ gemmi grep _refine.ls_R_factor_R_free -a _refine.pdbx_refine_id mmCIF/mo/?moo.cif.gz
1MOO:0.177;X-RAY DIFFRACTION
3MOO:0.21283;X-RAY DIFFRACTION
4MOO:0.22371;X-RAY DIFFRACTION
5MOO:0.1596;X-RAY DIFFRACTION
5MOO:0.1848;NEUTRON DIFFRACTION

Option -a (--and) can be specified many times. If we would add -a _pdbx_database_status.recvd_initial_deposition_date we would get the deposition date in each line. In this case it would be repeated for 5MOO:

5MOO:0.1596;X-RAY DIFFRACTION;2016-12-14
5MOO:0.1848;NEUTRON DIFFRACTION;2016-12-14

To output TSV (tab-separated values) add --delimiter='\t'. What are the heaviest chains?

$ gemmi grep --delimiter='\t' _entity.formula_weight -a _entity.pdbx_description /hdd/mmCIF/ | sort -nrk2 | head -3
6EK0    1641906.750     28S ribosomal RNA
5T2C    1640238.125     28S rRNA
5LKS    1640238.125     28S ribosomal RNA

With some further processing the option -a can be used to generate quite sophisticated reports. Here is a little demo: https://project-gemmi.github.io/pdb-stats/

The major limitation here is that gemmi-grep cannot match corresponding values from different tables (it is not possible on the syntax level). In the example above we have two values from the same table (_refine) and a deposition date (single value). This works well. But we are not able to add corresponding wavelengths from _diffrn_source. If an extra tag (specified with -a) is not in the same table as the main tag, gemmi-grep uses only the first value for this tag.

Unless we just count the number of value. Counting works for any combination of tags:

$ gemmi grep -c _refln.intensity_meas -a _diffrn_refln.intensity_net r5paysf.ent.gz
r5paysf:63611;0
r5payAsf:0;356684

(The file used in this example is structure factor (SF) mmCIF. Strangely these files in the PDB have extension ent not cif.)

The first number in the output above is the number of specified intensities. If you would like to count in also values ? and . specify the option --raw:

$ gemmi grep --raw -c _refln.intensity_meas r5paysf.ent.gz
r5paysf:63954
r5payAsf:0

Gemmi-grep can work with any CIF files but it has one feature specific to the PDB data. When $PDB_DIR is set one may use PDB codes: just 5moo or 5MOO instead of the path to 5moo.cif.gz. And for convenience, using a PDB code implies option -O.

The file paths or PDB codes can be read from a file. For example, if we want to analyse PDB data deposited in 2016 we may first make a file that lists all such files:

$ gemmi grep -H -O _pdbx_database_status.recvd_initial_deposition_date $PDB_DIR/structures/divided/mmCIF | \
        grep 2016 >year2016.txt

The 2016.txt file file has lines that start with the filename:

/hdd/structures/divided/mmCIF/ww/5ww9.cif.gz:5WW9:2016-12-31
/hdd/structures/divided/mmCIF/ww/5wwc.cif.gz:5WWC:2016-12-31

and a command such as:

$ gemmi grep -f year2016.out _diffrn.ambient_temp

will grep only the listed cif files.

Exit status of gemmi-grep has the same meaning as in GNU grep: 0 if a line is selected, 1 if no lines were selected, and 2 if an error occurred.

Examples

comp_id check

The monomer library (Refmac dictionary) has tags such as _chem_comp_atom.comp_id, _chem_comp_bond.comp_id that are expected to be consistent with the block name:

$ gemmi grep _*.comp_id $CLIBD_MON/a/ASN.cif
comp_ASN:ASN
[repeated 106 times]

We can quickly check if the names are always consistent by filtering the output above with awk, for all monomer files, to print only lines where the block name and comp_id differ:

$ gemmi grep _*.comp_id $CLIBD_MON/? | awk -F: 'substr($1, 6) != $2'
comp_M43:N09
...

planarity

The monomer library includes planarity restraints. Each row in the _chem_comp_plane_atom table with the same plane_id represents atom belonging to the same plane. What is the maximum number of atoms in one plane?

$ gemmi grep _chem_comp_plane_atom.plane_id $CLIBD_MON/? | uniq -c | sort -nr | head -3
 38 comp_LG8:plan-1
 36 comp_UCM:plan-1
 36 comp_SA3:plan-1

convert

$ gemmi convert -h
Usage:
 gemmi convert [options] INPUT_FILE OUTPUT_FILE

with possible conversions CIF-JSON, and mmCIF-PDB-mmJSON.
FORMAT can be specified as one of: cif, json, pdb.

General options:
  -h, --help             Print usage and exit.
  -V, --version          Print version and exit.
  --verbose              Verbose output.
  --from=FORMAT          Input format (default: from the file extension).
  --to=FORMAT            Output format (default: from the file extension).

CIF output options:
  --pdbx-style           Similar styling (formatting) as in wwPDB.
  --skip-category=CAT    Do not output tags starting with _CAT
  -b NAME, --block=NAME  Set block name and default _entry.id
  --sort                 Sort tags in alphabetical order.

JSON output options:
  -c, --comcifs          Conform to the COMCIFS CIF-JSON standard draft.
  -m, --mmjson           Compatible with mmJSON from PDBj.
  --bare-tags            Output tags without the first underscore.
  --numb=quote|nosu|mix  Convert the CIF numb type to one of:
                           quote - string in quotes,
                           nosu - number without s.u.,
                           mix (default) - quote only numbs with s.u.
  --dot=STRING           JSON representation of CIF's '.' (default: null).

PDB input/output options:
  --segment-as-chain     Append segment id to label_asym_id (chain name).
  --short-ter            Write PDB TER records without numbers (iotbx compat.).

Macromolecular operations:
  --expand-ncs=dup|addn  Expand strict NCS specified in MTRIXn or equivalent.
                         New chain names are the same or have added numbers.
  --remove-h             Remove hydrogens.
  --remove-waters        Remove waters.
  --remove-lig-wat       Remove ligands and waters.
  --trim-to-ala          Trim aminoacids to alanine.

When output file is -, write to standard output.

This programs combines a few functions.

CIF – JSON

Syntax-level conversion. The JSON representation of the CIF data can be customized. In particular we support CIF-JSON standard from COMCIFS and mmJSON standard from PDBj (the latter is specific to mmCIF files).

The major difference between the two is that CIF-JSON is dictionary-agnostic: it cannot recognize categories (mmJSON groups by categories), and it cannot recognize numbers (so it quotes the numbers). CIF-JSON adds also two extra objects: “CIF-JSON” and “Metadata”. The minor differences are:

CIF CIF-JSON mmJSON
data_a a data_a
_tag _tag tag
_CasE _case CasE
. false null
? null null

mmCIF – PDB – mmJSON

Conversion between macromolecular coordinate formats.

We made an effort to format the PDB files we write in the same way as the software used internally by the PDB, apart from writing fewer records. Thanks to this, in some scenarios a diff tool can be used to compare a PDB file written by Gemmi with an official PDB file from PDB.

The library and the converter also have an option (--iotbx-compat) that formats PDB files similarly to iotbx from cctbx (for example, in this mode TER records have no numbers). We do not aim to be fully compatible with CCTBX, but in many cases the difference will be only in whitespace.

NCS expansion

The option --expand-ncs expands strict NCS, defined in the MTRIX record (PDB) or in the _struct_ncs_oper table (mmCIF). By default, new chains have different names than the original ones. But when used together with --iotbx-compat, the program mimicks iotbx.pdb.expand_ncs and leaves the same chain names while adding distinct segment IDs.

map

Shows a summary of a CCP4 map file, optionally performing simple transformations.

$ gemmi map -h
Usage:
 gemmi map [options] CCP4_MAP[...]

  -h, --help         Print usage and exit.
  -V, --version      Print version and exit.
  --verbose          Verbose output.
  --deltas           Statistics of dx, dy and dz.
  --check-symmetry   Compare the values of symmetric points.
  --write-xyz=FILE   Write transposed map with fast X axis and slow Z.
  --write-full=FILE  Write map extended to cover whole unit cell.

mask

Makes a mask in the CCP4 format. It has two functions:

  • masking atoms if the input file is a coordinate file,
  • using a threshold to convert a CCP4 map file to a mask file.
$ gemmi mask -h
Usage:
 gemmi mask [options] INPUT output.msk

Makes a mask in the CCP4 format.
If INPUT is a CCP4 map the mask is created by thresholding the map.
If INPUT is a coordinate file (mmCIF, PDB, etc) the atoms are masked.
  -h, --help           Print usage and exit.
  -V, --version        Print version and exit.
  --verbose            Verbose output.
  --from=coor|map      Input type (default: from file extension).

Options for making a mask from a map:
  -t, --threshold      The density cutoff value.
  -f, --fraction       The volume fraction to be above the threshold.

Options for masking a model:
  -s, --spacing=D      Max. sampling for the grid (default: 1A).
  -g, --grid=NX,NY,NZ  Grid sampling.
  -r, --radius         Radius of atom spheres (default: 3.0A).

mtz

$ gemmi mtz -h
Usage:
 gemmi mtz [options] MTZ_FILE[...]
Print informations from an mtz file.
  -h, --help       Print usage and exit.
  -V, --version    Print version and exit.
  --verbose        Verbose output.
  -H, --headers    Print raw headers, until the END record.
  -d, --dump       Print a subset of CCP4 mtzdmp informations.
  --tsv            Print all the data as tab-separated values.
  --stats          Print column statistics (completeness, mean, etc).
  --check-asu      Check if reflections are in conventional ASU.
  --toggle-endian  Toggle assumed endiannes (little <-> big).

mtz2cif

Converts reflection data from MTZ to mmCIF.

$ gemmi mtz2cif -h
Usage:
  gemmi mtz2cif [options] MTZ_FILE CIF_FILE
Options:
  -h, --help             Print usage and exit.
  -V, --version          Print version and exit.
  --verbose              Verbose output.
  --spec=FILE            Column and format specification.
  --print-spec           Print default spec and exit.
  -b NAME, --block=NAME  mmCIF block name: data_NAME (default: mtz).
  --skip-empty           Skip reflections with no values.
  --no-comments          Do not write comments in the mmCIF file.

If CIF_FILE is -, the output is printed to stdout.
If spec is -, it is read from stdin.

Lines in the spec file have format:
  [FLAG] COLUMN TYPE TAG [FORMAT]
for example:
  SIGF_native * SIGF_meas_au 12.5e
  FREE I pdbx_r_free_flag 3.0f
FLAG (optional) is either ? or &:
  ? = ignored if no column in the MTZ file has this name.
  & = ignored if the previous line was ignored.
  Example:
      ? I    J intensity_meas
      & SIGI Q intensity_sigma
COLUMN is MTZ column label. Columns H K L are added if not specified.
  Alternative labels can be separated with | (e.g. FREE|FreeR_flag).
TYPE is used for checking the columm type, unless it is '*'.
TAG does not include category name, it is only the part after _refln.
FORMAT (optional) is printf-like floating-point format:
 - one of e, f, g with optional flag, width and precision
 - flag is one of + - # _; '_' stands for ' ', for example '_.4f'
 - since all numbers in MTZ are stored as float, the integer columns use
   the same format as float. The format of _refln.status is ignored.

cif2mtz

Converts reflection data from mmCIF to MTZ.

$ gemmi cif2mtz -h
Usage:
  gemmi cif2mtz [options] CIF_FILE MTZ_FILE
  gemmi cif2mtz [options] CIF_FILE --dir=DIRECTORY
Options:
  -h, --help               Print usage and exit.
  -V, --version            Print version and exit.
  --verbose                Verbose output.
  -b NAME, --block=NAME    mmCIF block to convert.
  -d DIR, --dir=NAME       Output directory.
  --title                  MTZ title.
  -H LINE, --history=LINE  Add a history line.

First variant: converts the first block of CIF_FILE, or the block
specified with --block=NAME, to MTZ file with given name.

Second variant: converts each block of CIF_FILE to one MTZ file
(block-name.mtz) in the specified DIRECTORY.

If CIF_FILE is -, the input is read from stdin.

sf2map

Transforms map coefficients from either MTZ or SF mmCIF to CCP4 map.

$ gemmi sf2map -h
Usage:
  gemmi sf2map [options] INPUT_FILE MAP_FILE

INPUT_FILE must be either MTZ or mmCIF with map coefficients.

By default, the program searches for 2mFo-DFc map coefficients in:
  - MTZ columns FWT/PHWT or 2FOFCWT/PH2FOFCWT,
  - mmCIF tags _refln.pdbx_FWT/pdbx_PHWT.
If option "-d" is given, mFo-DFc map coefficients are searched in:
  - MTZ columns DELFWT/PHDELWT or FOFCWT/PHFOFCWT,
  - mmCIF tags _refln.pdbx_DELFWT/pdbx_DELPHWT.


Options:
  -h, --help           Print usage and exit.
  -V, --version        Print version and exit.
  --verbose            Verbose output.
  -d, --diff           Use difference map coefficients.
  --section=NAME       MTZ dataset name or CIF block name
  -f COLUMN            F column (MTZ label or mmCIF tag).
  -p COLUMN            Phase column (MTZ label or mmCIF tag).
  -g, --grid=NX,NY,NZ  Grid size (user-specified minimum).
  --exact              Use the exact grid size specified by --grid.
  -s, --sample=NUMBER  Set spacing to d_min/NUMBER (3 is usual).
  --zyx                Output axes Z Y X as fast, medium, slow (default is X Y
                       Z).
  -G                   Print size of the grid that would be used and exit.
  --normalize          Scale the map to standard deviation 1 and mean 0.

The --sample option is named after the GRID SAMPLE keyword of the venerable CCP4 FFT program; its value has the same meaning.

map2sf

Transforms CCP4 map into map coefficients.

$ gemmi map2sf -h
Usage:
  gemmi map2sf [options] MAP_FILE OUTPUT_FILE COL_F COL_PH

Writes map coefficients (amplitude and phase) of a map to OUTPUT_FILE.
The output is MTZ if it has mtz extension, otherwise it is mmCIF.

Options:
  -h, --help       Print usage and exit.
  -V, --version    Print version and exit.
  --verbose        Verbose output.
  -b, --base=PATH  Add new columns to the data from this file.
  --section=NAME   Add new columns to this MTZ dataset or CIF block.
  --dmin=D_MIN     Resolution limit.
  --ftype=TYPE     MTZ amplitude column type (default: F).
  --phitype=TYPE   MTZ phase column type (default: P).

residues

List residues from a coordinate file, one per line.

$ gemmi residues -h
Usage:
 gemmi residues [options] INPUT[...]
Prints one residue per line, with atom names.
  -h, --help          Print usage and exit.
  -V, --version       Print version and exit.
  --format=FORMAT     Input format (default: from the file extension).
  -mSEL, --match=SEL  Print residues/atoms matching the selection.
INPUT is a coordinate file (mmCIF, PDB, etc).
The selection SEL has MMDB syntax:
/mdl/chn/s1.i1(res)-s2.i2/at[el]:aloc (all fields are optional)

Example:

$ gemmi residues -m '/3/*/(CYS,CSD)' 4pth.pdb
Model 3
A   85  CYS: N CA C O CB SG H HA HB2 HB3 HG
A  152  CSD: N CA CB SG C O OD1 OD2 HA HB2 HB3

seq

Sequence alignment (global, pairwise, affine gap penalty). Used primarily for aligning the residues in the model’s chains to the full sequence from the SEQRES record.

$ gemmi seq -h
Usage:
 gemmi seq [options] INPUT[...]
Compares sequence (SEQRES) with model from a PDB or mmCIF file.
For testing, it can also compare strings with option --text-align.
Performs global alignment with scoring matrix and affine gap penalty.
Currently, we only use match/mismatch scoring matrix.

  -h, --help     Print usage and exit.
  -V, --version  Print version and exit.
  --verbose      Verbose output.
  --check-mmcif  Compare alignment with mmCIF _atom_site.label_seq_id
  --text-align   Align characters in two strings (for testing).

Scoring (absolute values):
  --match=INT    Match score (default: 1).
  --mism=INT     Mismatch penalty (default: 1).
  --gapo=INT     Gap opening penalty (default: 1).
  --gape=INT     Gap extension penalty (default: 1).

For the testing purpose, it can align text strings. For example, the Levenshtein distance can be calculated by setting the gap opening penalty to zero:

$ ./gemmi seq --match=0 --gapo=0 --text-align Saturday Sunday
Score: -3   CIGAR: 1M2I5M
=II=X===
Saturday
S--unday

This tool uses modified code from ksw2.

sg

Prints information about given space group.

$ gemmi sg -h
Usage:
 gemmi sg [options] SPACEGROUP[...]
Prints information about the space group.
  -h, --help     Print usage and exit.
  -V, --version  Print version and exit.

contents

Analyses and summarizes content of a coordinate file. Inspired by CCP4 program rwcontents.

$ gemmi contents -h
Usage:
 gemmi contents [options] INPUT[...]
Analyses content of a PDB or mmCIF.
  -h, --help     Print usage and exit.
  -V, --version  Print version and exit.
  --verbose      Verbose output.
  --dihedrals    Print peptide dihedral angles.

contact

Searches for contacts in a model.

$ gemmi contact -h
Usage:
 gemmi contact [options] INPUT[...]
Searches for contacts in a model (PDB or mmCIF).
  -h, --help     Print usage and exit.
  -V, --version  Print version and exit.
  --verbose      Verbose output.
  -d, --maxdist=D  Maximal distance in A (default 3.0)
  --cov=TOL      Use max distance = covalent radii sum + TOL [A].
  --covmult=M    Use max distance = M * covalent radii sum + TOL [A].
  --occsum=MIN   Ignore atom pairs with summed occupancies < MIN.
  --any          Output any atom pair, even from the same residue.
  --noh          Ignore hydrogen (and deuterium) atoms.
  --count        Print only a count of atom pairs.

blobs

Searches for unmodelled blobs in electron density. Similar to “Validate > Unmodelled blobs…” in Coot. For use in Dimple.

$ gemmi blobs -h
Usage:
 gemmi blobs [options] MTZ_OR_MMCIF PDB_OR_MMCIF

Search for umodelled blobs of electron density.

Options:
  -h, --help            Print usage and exit.
  -V, --version         Print version and exit.
  --verbose             Verbose output.

The area around model is masked to search only unmodelled density.
  --mask-radius=NUMBER  Mask radius (default: 2.0 A).
  --mask-water          Mask water (water is not masked by default).

Searching blobs of density above:
  --sigma=NUMBER        Sigma (RMSD) level (default: 1.0).
  --abs=NUMBER          Absolute level in electrons/A^3.

Blob criteria:
  --min-volume=NUMBER   Minimal volume (default: 10.0 A^3).
  --min-score=NUMBER    Min. this electrons in blob (default: 15.0).
  --min-sigma=NUMBER    Min. peak rmsd (default: 0.0).
  --min-peak=NUMBER     Min. peak density (default: 0.0 el/A^3).

Options for map calculation:
  -d, --diff            Use difference map coefficients.
  --section=NAME        MTZ dataset name or CIF block name
  -f COLUMN             F column (MTZ label or mmCIF tag).
  -p COLUMN             Phase column (MTZ label or mmCIF tag).
  -g, --grid=NX,NY,NZ   Grid size (user-specified minimum).
  --exact               Use the exact grid size specified by --grid.
  -s, --sample=NUMBER   Set spacing to d_min/NUMBER (3 is usual).
  -G                    Print size of the grid that would be used and exit.

h

Adds or removes hydrogens. Hydrogen are put in positions based only on restraints from a monomer library.

$ gemmi h -h
Usage:
 gemmi h [options] INPUT_FILE OUTPUT_FILE

Add hydrogens in positions specified by the monomer library.
By default, it removes and re-adds all hydrogens.
By default, hydrogens are not added to water.

Options:
  -h, --help      Print usage and exit.
  -V, --version   Print version and exit.
  --verbose       Verbose output.
  --monomers=DIR  Monomer library dir (default: $CLIBD_MON).
  --remove        Only remove hydrogens.
  --keep          Do not add/remove hydrogens, only change positions.
  --water         Add hydrogens also to waters.
  --sort          Order atoms in residues according to _chem_comp_atom.

mondiff

Compares restraints from two monomer CIF files. It is intended for comparing restraints for the same monomer, but generated with different programs (or different versions of the same program).

The files should have format used by the CCP4/Refmac monomer library. This format is supported by all major macromolecular refinement programs.

$ gemmi mondiff -h
Usage:
  gemmi mondiff [options] FILE1 FILE2
Options:
  -h, --help         Print usage and exit.
  -V, --version      Print version and exit.
  --verbose          Verbose output.

Minimal reported differences:
  --bond=DELTA       difference in distance value (default: 0.01).
  --bond-esd=DELTA   difference in distance esd (default: 0.1).
  --angle=DELTA      difference in angle value (default: 0.1).
  --angle-esd=DELTA  difference in angle esd (default: 1.0).
  --rel=SIGMA         abs(value difference) / esd > SIGMA (default: 0.0).

wcn

Calculates Weighted Contact Number (WCN) and a few other similar metrics.

WCN can be used to predicts B-factors (ADPs) from coordinates, and to compare this prediction with the values from refinement.

Background

Protein flexibility and dynamic properties can be to some degree inferred from the atomic coordinates of the structure. Various approaches are used in the literature: molecular dynamics, Gaussian or elastic network models, normal mode analysis, calculation of solvent accessibility or local packing density, and so on.

Here we apply the simplest approach, which is pretty effective. It originates from the 2002 PNAS paper in which Bertil Halle concluded that B-factors are more accurately predicted by counting nearby atoms than by Gaussian network models. This claim was based on the analysis of only 38 high resolution structures (and a neat theory), but later on the method was validated on many other structures.

In particular, in 2007 Manfred Weiss brought this method to the attention of crystallographers by analysing in Acta Cryst D different variants of the methods on a wider set of more representative crystals. Recently, the parameters fine-tuned by Weiss have been used for guessing which high B-factors (high comparing with the predicted value) result from the radiation damage.

Only a few months later, in 2008, Chih-Peng Lin et al. devised a simple yet significant improvement to the original Halle’s method: weighting the counted atoms by 1/d2, the inverse of squared distance. It nicely counters the increasing average number of atoms with the distance (~ d2). This method was named WCN – weighted contact number (hmm.. “contact”).

These two methods are so simple that it seems easy to find a better one. But according to my quick literature search, no better method of this kind has been found yet. In 2009 Li and Bruschweiler proposed weighting that decreases exponentially (that model was named LCBM), but in my hands it does not give better results than WCN.

In 2016 Shahmoradi and Wilke did a data analysis aiming to disentangle the effects of local and longer-range packing in the above methods. They were not concerned with B-factors, though, but with the rate of protein sequence evolution. Because the “contact” methods predict many things. Interestingly, if the exponent in WCN is treated as a parameter (equal -2 in the canonical version), the value -2.3 gives the best results in predicting evolution.

TLS

We also need to note that TLS-like methods that model B-factors as rigid-body motion of molecules are reported to give better correlation with experimental B-factors than other methods. But because such models use experimental B-factors on the input and employ more parameters, they are not directly comparable with WCN.

Unlike the TLS that is routinely used in the refinement of diffraction data, the TLS modelling described here is isotropic. It uses 10 parameters (anisotropic TLS requires 20) as described in a paper by Kuriyan and Weis (1991). Soheilifard et al (2008) got even better results by increasing B-factors at the protein ends, using 13 parameters all together. This model was named eTLS (e = extended).

The high effectiveness of the TLS model does not mean that B-factors are dominated by the rigid-body motion. As noted by Kuriyan and Weis, the TLS model captures also the fact that atoms in the interior of a protein molecule generally have smaller displacements than those on the exterior. Additionally, authors of the LCBM paper find that the TLS model fitted to only half of the protein poorly fits the other half, which suggests overfitting.

We may revisit rigid-body modelling in the future, but now we get back to the contact numbers.

Details

The overview above skipped a few details.

  • While the WCN method is consistently called WCN, the Halle’s method was named LDM (local density model) in the original paper, and is called CN (contact number) in some other papers. CN is memorable when comparing with WCN (which adds ‘W’ – weighting) and with ACN (which adds ‘A’ – atomic).
  • These method are used either per-atom (for predicting B-factors, etc.) or per-residue (for evolutionary rate, etc.). So having “A” in ACN clarifies how it is used. To calculate the contact number per-residue one needs to pick a reference point in the residue (Cβ, the center of mass or something else), but here we do only per-atom calculations.
  • The CN method requires a cut-off, and the cut-off values vary widely, from about 5 to 18Å. In the original paper it was 7.35Å, Weiss got 7.0Å as the optimal value, Shahmoradi 14.3Å.
  • The CN can be seen as weighted by Heaviside step function, and smoothing it helps a little bit (as reported by both Halle and Weiss).
  • Similarly to eTLS, the LCBM method has eLCBM variant that adds “end effects” – special treatment of the termini.
  • Finally, these methods may or may not consider the symmetry mates in the crystal. Halle checked that including symmetric images improves the prediction. Weiss (ACN) and Li and Bruschweiler (LCBM) are also taking symmetry into account. But I think other papers don’t.

Metrics for comparison

To compare a number of nearby atoms with B-factor we either rescale the former, or we use a metric that does not require rescaling. The Pearson correlation coefficient (CC) is invariant under linear transformation, so it can be calculated directly unless we would like to apply non-linear scaling. Which was tried only in the Manfred Weiss’ paper: scaling function with three parameters improved CC by 0.012 comparing with linear function (that has two parameters). Here, to keep it simple, we only do linear scaling.

As noted by Halle, Pearson CC as well as mean-square deviation can be dominated by a few outliers. Therefore Halle used relative mean absolute deviation (RMAD): sum of absolute differences divided by the average absolute deviation in the experimantal values. Halle justifies this normalization writing that it allows to compare structures determined at different temperatures. This is debatable as can be seen from ccp4bb discussions on how to compare B-factors between two structures. But for sure RMAD is a more robust metric, so we also use it. It adds another complication, though. To minimize the absolute deviation we cannot use least-squares fitting, but rather quantile regression with q=0.5.

Another metric is the rank correlation. It is interesting because it is invariant under any monotonic scaling. But it is not guaranteed to be a good measure of similarity.

Results

To be wrapped up and published. But in the meantime here are some thoughts:

  • The optimal exponent is slightly larger than 2; the difference is small, so we prefer to use 2 (i.e. w=1/r2).
  • Accounting for all symmetry mates (i.e. for intermolecular contacts in the crystal) improves the results – and then the cut-off is necessary.
  • The optimal cut-off is around 15A – let’s use 15A.
  • Averaging predicted B-factors of nearby atoms helps; we use Gaussian smoothing (blurring) with σ around 2A.
  • Pearson CC around 0.8 may seem high, but it corresponds to R2=0.64, i.e. it we explain only 64% of the B-factor variance. Even less of the absolute deviation – below 50%.
  • Minimizing absolute deviation (with quantile regression) gives similar results as the ordinary least squares (OLS). The difference in terms of RMAS is only ~0.03.
  • Combining WCN with CN is helping only a tiny bit (i.e. both are highly correlated) at the cost of additional parameter that is fitted. Combining WCN with rotation-only model (squared distance from the center of mass) increases CC slightly more, but still not much.
  • Accounting for symmetry mates worsens prediction of evolutionary rates. I used data from Shahmoradi and Wilke to check this.

Program

gemmi-wcn implements combination of the CN and WCN methods above.

Being based on a general-purpose crystallographic library it handles corner cases that are often ignored. A good example is searching for contacts. For most of the structures, considering only the same and neighbouring unit cells (1+26) is enough. But some structures have contacts between molecules several unit cells apart, even with only a single chain in the asu.

TBC

$ gemmi wcn -h
Usage:
 gemmi wcn [options] INPUT[...]
Calculation of local density / contact numbers: WCN, CN, ACN, LDM, etc.
  -h, --help       Print usage and exit.
  -V, --version    Print version and exit.
  --verbose        Verbose output.
  -f, --file=FILE  Obtain paths or PDB IDs from FILE, one per line.
  -l, --list       List per-residue values.
  --min-dist=DIST  Minimum distance for "contacts" (default: 0.8).
  --cutoff=DIST    Maximum distance for "contacts" (default: 15).
  --pow=P          Exponent in the weighting (default: 2).
  --blur=SIGMA     Apply Gaussian smoothing of predicted B-factors.
  --rom            Rotation only model: |pos-ctr_of_chain|^P instead of WCN.
  --chain=CHAIN    Use only one chain from the INPUT file.
  --sanity         Run sanity checks first.
  --sidechains=X   One of: include, exclude, only (default: include).
  --no-crystal     Ignore crystal symmetry and intermolecular contacts.
  --omit-ends=N    Ignore N terminal residues from each chain end.
  --print-res      Print also resolution and R-free.
  --xy-out=DIR     Write DIR/name.xy files with WCN and B(exper).