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.7.1-dev
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 Global Phasing Ltd.
https://github.com/project-gemmi/gemmi
Usage: gemmi [--version] [--help] <command> [<args>]
Commands:
align sequence alignment (global, pairwise, affine gap penalty)
blobs list unmodelled electron density blobs
cif2json translate (mm)CIF to (mm)JSON
cif2mtz convert structure factor mmCIF to MTZ
cifdiff compare tags in two (mm)CIF files
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
crd prepare topology file (.crd) for Refmac
ecalc calculate normalized amplitudes E
fprime calculate anomalous scattering factors f' and f"
grep search for tags in CIF file(s)
h add or remove hydrogen atoms
json2cif translate mmJSON to mmCIF
map print info or modify a CCP4 map
map2sf transform CCP4 map to map coefficients (in MTZ or mmCIF)
mask make a bulk-solvent mask in the CCP4 format
merge merge intensities from multi-record reflection file
mondiff compare two monomer CIF files
mtz print info about MTZ reflection file
mtz2cif convert MTZ to structure factor mmCIF
reindex reindex MTZ file
residues list residues from a coordinate file
rmsz validate geometry using monomer library
set modify coordinate files (think pdbset)
sf2map transform map coefficients (from MTZ or mmCIF) to map
sfcalc calculate structure factors from a model
sg info about space groups
tags list tags from CIF file(s)
validate validate CIF 1.1 syntax
wcn calculate local density / contact numbers (WCN, CN, ACN, LDM)
xds2mtz convert XDS_ASCII to MTZ
Hint. To process multiple files in parallel, in multiple threads, use a tool such as GNU parallel:
$ find $PDB_DIR/structures/divided/mmCIF/ -name '*.cif.gz' | parallel gemmi grep _exptl.method
validate¶
A CIF validator. Apart from checking the syntax it can check most of the rules imposed by DDL1 and DDL2 dictionaries.
If you want to validate mmCIF files, the current version of the PDBx/mmCIF specification, maintained by the PDB, is distributed as one file (mmcif_pdbx_v50.dic), which can be used to validate all kinds of mmCIF files: coordinate files, reflection files, and CCD monomers. Note that such validation can spot only certain types of mistakes. It won’t tell you if the file is appropriate for deposition to the PDB. Dictionary-based validation can’t even tell if the file contains all the necessary tables; it is unaware of what the file represents – coordinates, reflection data or something else. On the other hand, the mmCIF files deposited to the PDB do not need to strictly conform to the PDBx/mmCIF spec. Not even the files distributed by the PDB are fully compliant (partly because not everything can be expressed in DDL2 syntax; usually it’s about child-parent relationships; PDB’s own validator, program CifCheck from mmcif-dict-suite, has a few exceptions hardcoded in C++, so that non-conformance is not accidental).
$ gemmi validate -h
Usage: gemmi validate [options] FILE [...]
Options:
-h, --help Print usage and exit.
-V, --version Print version and exit.
-v, --verbose Verbose output.
-q, --quiet Show only errors.
-f, --fast Syntax-only check.
-s, --stat Show token statistics
-d, --ddl=PATH DDL for validation.
-c, --context Check _pdbx_{category|item}_context.type.
--no-regex Skip regex checking (when using DDL2)
--no-mandatory Skip checking if mandatory tags are present.
--no-unique Skip checking if category keys (DDL2) are unique.
-p Check if parent items (DDL2) are present.
-r, --recursive Recurse directories and process all CIF files.
Validation specific to CCP4 monomer files:
-m, --monomer Run checks specific to monomer dictionary.
--z-score=Z Use Z for validating _chem_comp_atom.[xyz] (default: 2.0).
--ccd=PATH CCD file for comparison.
--audit-on=DATE Check only if CCD component was remediated on DATE.
grep¶
Searches for values in CIF files that correspond to the given tag(s).
Option -h
shows the usage.
Some of the command-line options (-c
, -l
, -H
, -n
)
correspond to the options of GNU grep.
$ 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.
By default, recursive directory search checks only *.cif(.gz) files.
To change it, specify --name=* or --name=*.hkl.
By default (without -E), TAG it a tag or a glob pattern. Must start with _.
Options:
-h, --help display this help and exit
-V, --version display version information and exit
-v, --verbose print additional messages to stderr
-f, --file=FILE obtain file (or PDB ID) list from FILE
--name=PATTERN filename glob pattern used in recursive grep; by
default, *.cif and *.cif.gz files are searched
-S, --pdb-sf if PDB ID is given, search structure factor file
-m, --max-count=NUM print max NUM values per file
-O, --one-block optimize assuming one block per file
-E, --extended-regexp interpret TAG as egrep-like regular expression
-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
-T, --only-tags print only matching tags, not values
-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
When used without any options, it prints the values associated with the specified tag. Each line contains the block name and the value:
$ 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
The output can be easily processed with well-known Unix utilities, therefore gemmi-grep does not have internal options for sorting and filtering.
For text-based filtering, use
grep
:$ # list mmCIF files with links to EMDB $ gemmi grep _pdbx_database_related.db_name /pdb/mmCIF/aa/* | grep EMDB 4AAS:EMDB 5AA0:EMDB
For numeric filtering, use
awk
(see the example with cell angles below).For sorting the output, use
sort
:$ # the heaviest monomers in CCD $ gemmi grep _chem_comp.formula_weight components.cif.gz | sort -k2 -t: -nr | head -3 JGH:5496.513 WO2:4363.030 HFW:4343.697
For counting different enumeration values, use
sort | uniq -c
(see the example with _entity_poly.type below).If you add
cut
to the mix, you could even count the number of distinct values (we won’t go further than this):$ # chemical components with the greatest variety of elements $ # LC_COLLATE=C avoids ignoring ':', which could result in order HOH:H HO:HO HOH:O $ gemmi grep _chem_comp_atom.type_symbol components.cif.gz | LC_COLLATE=C sort -u | cut -d: -f1 | uniq -c | sort -nr | head -5 8 YL6 8 YF7 8 VM5
Gemmi-grep parses the CIF syntax and prints values from both name-value
pairs and loops. In particular, gemmi grep _one
will give the same
output for both _one 1
and loop_ _one _two 1 2
.
It makes it more robust than using Unix grep. For example,
if one was using Unix grep to get R-free from mmCIF file,
it would fail in special cases such as the PDB entry 5MOO, which has two
Rfree values in a loop (see above, the second example in this section).
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 occurrence 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 greatest 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
or which chemical components in the CCD have the biggest number of atoms:
$ gemmi grep -c _chem_comp_atom.atom_id components.cif.gz | sort -t: -k2 -nr | head -3
JSG:440
KDL:355
TQN:348
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
gemmi-grep operates on the syntax level and cannot match values from
different tables, as this would require consulting a DDL dictionary.
In the example above we have two values from the same table (_refine
)
and a deposition date (single value). This works well. However,
would not be able to add the corresponding wavelengths from _diffrn_source
.
If an extra tag (specified with -a
) is not in the same table
as the main tag, gemmi-grep will use the first value for this tag,
not the corresponding one.
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/xray.html
If -a
is used together with -c
, the values are counted
independently for each tag:
$ 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.
To count also null values ?
and .
, add --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 analyze PDB data deposited in 2016, we may first create 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 year2016.txt 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
If we make a list of filenames only:
$ awk -F: '{print $1}' year2016.txt > year2016_.txt
we can then use -f
to grep files from this list:
$ gemmi grep -f year2016_.txt _diffrn.ambient_temp
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
cif2json¶
Syntax-level conversion from CIF 1.1 to JSON. 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).
$ gemmi cif2json -h
Usage:
gemmi cif2json [options] INPUT_FILE OUTPUT_FILE
Convert CIF file (any CIF files, including mmCIF) to JSON.
The output can be COMCIFS CIF-JSON (-c), mmJSON (-m),
or a custom JSON flavor (default).
General options:
-h, --help Print usage and exit.
-V, --version Print version and exit.
-v, --verbose Verbose output.
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).
Modifications:
--skip-category=CAT Do not output tags starting with _CAT
--sort Sort tags in alphabetical order.
When output file is -, write to standard output.
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
json2cif¶
The opposite of cif2json, but currently the only supported input is mmJSON.
$ gemmi json2cif -h
Usage:
gemmi json2cif [options] INPUT_FILE OUTPUT_FILE
Convert mmJSON to mmCIF.
Options:
-h, --help Print usage and exit.
-V, --version Print version and exit.
-v, --verbose Verbose output.
--style=STYLE One of: default, pdbx (categories separated with #),
aligned (left-aligned columns).
--cif2cif Read CIF not JSON.
--skip-category=CAT Do not output tags starting with _CAT
--sort Sort tags in alphabetical order.
When output file is -, write to standard output.
cifdiff¶
Compares categories and tags that are present in two (mm)CIF files. It does not compare values.
$ gemmi cifdiff -h
Usage:
gemmi cifdiff [options] FILE1.cif FILE2.cif
gemmi cifdiff [options] -n FILE.cif
Compares (or just prints) categories and tags in CIF files.
First block only.
-h, --help Print usage and exit.
-V, --version Print version and exit.
-v, --verbose Verbose output.
-t TAG Compare values of TAG.
-q Print only categories.
-n No comparison, just list categories and tags.
The output resembles the output of diff, so it can be colored by editors that highlight “syntax” or by the colordiff script. Here is an example:
$ gemmi diff deposit.cif output.cif
--- Reading deposit.cif
+++ Reading output.cif
_entity. rows: 8 -> 11
id
type
+ _chem_comp. rows: 0 -> 25
+ id
+ type
_exptl. rows: 1
entry_id
method
+ crystals_number
_reflns. rows: 1
+ entry_id
pdbx_ordinal
+ pdbx_diffrn_id
number_obs
d_resolution_high
d_resolution_low
percent_possible_obs
pdbx_redundancy
pdbx_Rmerge_I_obs
- pdbx_Rrim_I_all
- pdbx_Rpim_I_all
+ pdbx_Rsym_value
pdbx_netI_over_sigmaI
- pdbx_CC_half
[...]
convert¶
Conversion between macromolecular coordinate formats: PDB, mmCIF and mmJSON.
$ gemmi convert -h
Usage:
gemmi convert [options] INPUT_FILE OUTPUT_FILE
Allows conversion between PDB, mmCIF, and mmJSON formats.
General options:
-h, --help Print usage and exit.
-V, --version Print version and exit.
-v, --verbose Verbose output.
--from=FORMAT Input format (default: inferred from file extension).
--to=FORMAT Output format (default: inferred from file extension).
mmCIF output options:
--style=STYLE one of: default, pdbx (categories separated with #),
aligned (left-aligned columns).
--all-auth Write _atom_site.auth_atom_id (same as label_atom_id)
and auth_comp_id (same as label_comp_id).
-b NAME, --block=NAME Set block name and default _entry.id
--sort Sort tags in alphabetical order.
--skip-category=CAT Do not output tags starting with _CAT
PDB input options:
--segment-as-chain Append segment id to label_asym_id (chain name).
--old-pdb Read only the first 72 characters of each line.
-L, --force-label Add label_seq_id even if SEQRES is missing
PDB output options:
--short-ter Write PDB TER records without numbers (iotbx compat.).
--linkr Write LINKR record (for Refmac) if link_id is known.
--copy-remarks (pdb->pdb only) Copy REMARK records.
General output options:
--minimal Write only the most essential records.
--shorten Shorten chain names to 1 (if # < 63) or 2 characters.
--rename-chain=OLD:NEW Rename chain OLD to NEW (--rename-chain=:A adds
missing chain IDs).
--shorten-tlc Change 5-character monomer names to 3-char. aliases.
--monomer=OLD:NEW Change monomer name (CCD code) OLD to NEW.
-s FILE Use sequence(s) from FILE in PIR or FASTA format. Each
chain is assigned the best matching sequence, if any.
--sifts-num[=AC,...] Set sequence ID to SIFTS-mapped UniProt positions, add
5000+ to non-mapped seqnums. See docs for details.
-B MIN[:MAX] Set isotropic B-factors to a single value or constrain
them to a range.
--scale-biso=MULT Multiply isotropic B-factors by MULT.
--anisou=yes|no|heavy Add or remove ANISOU records.
--assign-records[=A|H] Re-assign ATOM/HETATM (w/o argument: auto).
Macromolecular operations:
--select=SEL Output only the specified selection.
--remove=SEL Remove the specified selection.
--apply-symop=OP Apply operation, e.g. '-x,y+1/2,-z' or 'x,y,z+0.1'.
--reframe Standardize the coordinate system (frame).
--expand-ncs=dup|num|x Expand strict NCS from MTRIXn or _struct_ncs_oper.
Argument controls naming of new chains; see docs.
--assembly=ID Output bioassembly with specified ID (1, 2, ...).
--remove-h Remove hydrogens.
--remove-waters Remove waters.
--remove-lig-wat Remove ligands and waters.
--trim-to-ala Trim aminoacids to alanine.
FORMAT can be specified as one of: mmcif, mmjson, pdb. chemcomp (read-only).
chemcomp = coordinates of a component from CCD or monomer library (see docs).
When output file is -, write to standard output (default format: pdb).
Gemmi formats PDB records in the same way as wwPDB (including
trailing spaces) to enable file comparison using diff
.
CCD files include two sets of coordinates: example model and ideal.
By default, when converting a CCD component to another format,
both sets are written as separate models.
To output only a single model, specify --from=chemcomp:m
(for model)
or --from=chemcomp:i
(for ideal).
--expand-ncs
¶
This option expands strict NCS, as defined in the MTRIX
records (PDB)
or the _struct_ncs_oper
table (mmCIF).
The difficult part is naming the added chains; three options are provided:
=dup
– the chain names remain unchanged, with distinct segment IDs added,=num
– the names are appended with numbers to avoid duplication,=x
– unused 1- or 2-character names are assigned, supporting up to 62 + 62×62 = 3906 chains and allowing output in the PDB format.
--sifts-num
¶
This option changes sequence IDs to the corresponding sequence positions from UniProt. Residues that don’t have UniProt correspondence have their sequence numbers increased by an offset of 5000 (like in PDBrenum). In rare cases where UniProt positions are around 5000 or higher, the offset is increased to 6000 or a larger round number.
Note that the mapping between PDB and UniProt is based on
SIFTS (not DBREF) annotations in mmCIF files.
Currently, these annotations are present only in the PDB NextGen Archive
and in PDBe “updated” files – --sifts-num
works only with these files!
If accession codes (ACs) are specified in this option, only the matching UniProt ACs are used, and non-matching chains have their sequence numbers increased by an offset.
For chimeric chains that correspond to 2+ UniProt sequences, we use the sequence with the most corresponding residues. This choice can be overridden by explicitly specifying the preferred UniProtKB identifiers (usually just one AC, but it’s possible to specify multiple comma-separated ACs, in the order of preference).
*
in the argument list (e.g. --sifts-num=P01234,*
)
means that all chains with corresponding UniProt entries are renumbered
to match UniProt positions (but P01234 is preferred for chimeric chains).
An absent argument (--sifts-num
) is equivalent to *
(--sifts-num=*
).
Chains that don’t match the UniProt ACs have their sequence numbers
bumped by 5000+, similarly to ligands and waters in matching chains.
To leave the non-matching chains unchanged, add =
at the end,
e.g. --sifts-num=*,=
.
set¶
Modifies atom positions, isotropic B-factors, and/or occupancies in a PDB or mmCIF file. It serves as a partial replacement for CCP4 PDBSET.
Unlike most other gemmi tools, it doesn’t parse the entire input file. Instead, it reads only what is necessary to locate the relevant numbers and replaces them, leaving the rest of the file unchanged.
$ gemmi set -h
Usage:
gemmi set [options] INPUT_FILE OUTPUT_FILE
Modify atom attributes in a coordinate file.
Options:
-h, --help Print usage and exit.
-V, --version Print version and exit.
-v, --verbose Verbose output.
-B MIN[:MAX] Set isotropic B-factors to a single value MIN
or clamp them to MIN:MAX.
-O MIN[:MAX] Set occupancies to a single value or clamp them to MIN:MAX.
--noise M Add random shifts, uniform in (-M,M), to x,y,z.
--shift='DX DY DZ' Translate the model coordinates (units: Angstroms).
--select=SEL Apply transformations only to selected atoms (MMDB syntax).
map¶
Shows a summary of a CCP4 map file, optionally performing simple transformations.
It plots a histogram of values in the console.
All histograms plotted by gemmi have a default width of 80 characters,
which can be overridden by the environment variable COLUMNS.
In Bash shell this variable is set automatically, but you need
export COLUMNS
to pass it to the program.
$ gemmi map -h
Usage:
gemmi map [options] CCP4_MAP[...]
-h, --help Print usage and exit.
-V, --version Print version and exit.
-v, --verbose Verbose output.
-d, --dump Print a map summary (default action).
--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.
--write-mask=FILE Make a mask by thresholding the map.
Options for making a mask:
--threshold Explicit threshold value for 0/1 mask.
--fraction Threshold is selected to have this fraction of 1's.
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 bulk-solvent mask in the CCP4 format.
INPUT is a coordinate file (mmCIF, PDB, etc).
-h, --help Print usage and exit.
-V, --version Print version and exit.
-v, --verbose Verbose output.
--timing Print how long individual steps take.
-s, --spacing=D Max. sampling for the grid (default: 1A).
-g, --grid=NX,NY,NZ Grid sampling.
-r, --radius=R Use constant radius of atom spheres.
--r-probe=Rp Use VdW radius + Rp (default: 1.0A).
--r-shrink=Rs Finally, remove a shell of thickness Rs (default: 1.1A).
--island-limit=VOL Remove "islands" up to VOL A^3.
--hydrogens Don't ignore hydrogens.
--any-occupancy Don't ignore zero-occupancy atoms.
--cctbx-compat Use vdW, Rprobe, Rshrink radii from cctbx.
--refmac-compat Use radii compatible with Refmac.
-I, --invert 0 for solvent, 1 for molecule.
mtz¶
$ gemmi mtz -h
Usage:
gemmi mtz [options] MTZ_FILE[...]
Print information from an mtz file.
-h, --help Print usage and exit.
-V, --version Print version and exit.
-v, --verbose Verbose output.
-H, --headers Print raw headers, until the END record.
-d, --dump Print a subset of CCP4 mtzdmp information.
-B N, --batch=N Print data from batch header N.
-b, --batches Print data from all batch headers.
-e (with -B or -b) expanded info from batch headers.
-A, --appendix Print appended text.
--tsv Print all the data as tab-separated values.
-s, --stats Print column statistics (completeness, mean, etc).
--histogram=LABEL Print histogram of values in column LABEL.
--cells Print cell parameters only.
--check-asu=ccp4|tnt Check if reflections are in ASU.
--compare=FILE Compare two MTZ files.
--toggle-endian Toggle assumed endianness (little <-> big).
--no-isym Do not apply symmetry from M/ISYM column.
mtz2cif¶
Converts reflection data from MTZ to mmCIF.
$ gemmi mtz2cif -h
Usage:
gemmi mtz2cif [options] MTZ_FILE [MTZ_FILE] CIF_FILE
Options:
-h, --help Print usage and exit.
-V, --version Print version and exit.
-v, --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).
--id=ID value for _entry.id (default: xxxx).
--skip-empty[=COLS] Skip reflections with no values. If COLS are given, eg.
'I(+),I(-)', only values in those columns are checked.
--skip-negative-sigi Skip reflections with sigma(I)<0 in unmerged data.
--no-comments Do not write comments in the mmCIF file.
--no-history Do not write MTZ history in the mmCIF file.
--no-staraniso-tensor Do not write _reflns.pdbx_aniso_B_tensor_*.
--run-from=WRAPPER Add note in _software.description.
--wavelength=LAMBDA Set wavelength (default: from input file).
--validate For two MTZ files: validate the intensities match.
--less-ano Skip anomalous columns (even if they are in the spec).
Used once, skips I(+)/I(-) if <I> and F(+)/F(-) are
present. Used twice, skips all anomalous columns.
--separate Write merged and unmerged data in separate blocks.
--depo Prepare merged+unmerged mmCIF file for deposition.
--nfree=N Flag value used for the free set (default: auto)
--trim=N (for testing) output only reflections -N <= h,k,l <=N.
One or two MTZ files are taken as the input. If two files are given,
one must be merged and the other unmerged.
Alternatively, the two input files can be reflection mmCIF file and
unmerged MTZ - this adds unmerged data block to the existing mmCIF file.
XDS_ASCII.HKL file can be used instead of unmerged MTZ file (also when
converting a single file), but then the spec file is ignored.
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]
or
$SPECIAL TAG
for example:
SIGF_native * SIGF_meas_au 12.5e
FREE I pdbx_r_free_flag 3.0f
$counter id
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 column 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.
$SPECIAL is $counter, $dataset, $image, $. or $?, as appropriate for
_diffrn_refln.id, _diffrn_refln.diffrn_id, image number or null (./?)
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.
-v, --verbose Verbose output.
-b NAME, --block=NAME mmCIF block to convert, by name.
-B INDEX mmCIF block to convert, by index (default: 1).
--add CIF_FILE Use additional input mmCIF file, first block.
-l, --list Dry run and list blocks in mmCIF file.
-d DIR, --dir=NAME Output directory.
--spec=FILE Conversion spec.
--print-spec Print default spec and exit.
--title MTZ title.
-H LINE, --history=LINE Add a history line.
--wavelength=LAMBDA Set wavelength (default: from input file).
-u, --unmerged Write unmerged MTZ file(s).
--refln-to=MODE Read refln category as: normal, friedel (converts Fs
to F+/F-), unmerged or auto (default).
--sort Order reflections according to Miller indices.
--asu=ccp4|tnt Move merged reflections into CCP4 or TNT ASU.
--skip-negative-sigma Skip reflections with sigma<0 (in any MTZ Q column).
--zero-to-mnf If value and sigma are 0, set both to MNF.
--local Take file from local copy of the PDB archive in
$PDB_DIR/structures/divided/structure_factors/
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.
To convert data from multiple CIF files into one MTZ file use:
gemmi cif2mtz [options] CIF1 --add CIF2 --add CIF3 MTZ_FILE
It's for special cases, e.g. when map coefficients are in separate files.
Similarly to mtz2cif, this converter can also be customized using spec files.
The default spec files can be generated using options
--print-spec
for merged data and --print-spec --unmerged
for unmerged.
$ gemmi cif2mtz --print-spec
# Each line in the spec contains 4-5 words:
# - tag (without category) from _refln or _diffrn_refln
# - MTZ column label
# - MTZ column type
# - MTZ dataset for the column (must be 0 or 1)
# - (optional) how to map mmCIF symbols to MTZ numbers
# For MERGED data only. Use --print-spec --unmerged for unmerged.
# Conversion of the first MTZ columns (H K L) is hardcoded - not in the spec.
pdbx_r_free_flag FreeR_flag I 0
status FreeR_flag I 0 o=1,f=0
intensity_meas IMEAN J 1
F_squared_meas IMEAN J 1
intensity_sigma SIGIMEAN Q 1
F_squared_sigma SIGIMEAN Q 1
pdbx_I_plus I(+) K 1
pdbx_I_plus_sigma SIGI(+) M 1
pdbx_I_minus I(-) K 1
pdbx_I_minus_sigma SIGI(-) M 1
F_meas FP F 1
F_meas_au FP F 1
F_meas_sigma SIGFP Q 1
F_meas_sigma_au SIGFP Q 1
pdbx_F_plus F(+) G 1
pdbx_F_plus_sigma SIGF(+) L 1
pdbx_F_minus F(-) G 1
pdbx_F_minus_sigma SIGF(-) L 1
pdbx_anom_difference DP D 1
pdbx_anom_difference_sigma SIGDP Q 1
F_calc FC F 1
F_calc_au FC F 1
phase_calc PHIC P 1
pdbx_F_calc_with_solvent F-model F 1
pdbx_phase_calc_with_solvent PHIF-model P 1
fom FOM W 1
weight FOM W 1
pdbx_HL_A_iso HLA A 1
pdbx_HL_B_iso HLB A 1
pdbx_HL_C_iso HLC A 1
pdbx_HL_D_iso HLD A 1
pdbx_FWT FWT F 1
pdbx_PHWT PHWT P 1
pdbx_DELFWT DELFWT F 1
pdbx_DELPHWT PHDELWT P 1
$ gemmi cif2mtz --print-spec --unmerged
# Each line in the spec contains 4-5 words:
# - tag (without category) from _refln or _diffrn_refln
# - MTZ column label
# - MTZ column type
# - MTZ dataset for the column (must be 0 or 1)
# - (optional) how to map mmCIF symbols to MTZ numbers
# For UNMERGED data only. Conversion of the first MTZ columns
# (H K L M/ISYM BATCH) is hardcoded - not in the spec.
intensity_meas I J 0
intensity_net I J 0
intensity_sigma SIGI Q 0
pdbx_detector_x XDET R 0
pdbx_detector_y YDET R 0
pdbx_scan_angle ROT R 0
sf2map¶
Transforms map coefficients from either MTZ or SF mmCIF to CCP4 map.
$ gemmi sf2map -h
Usage:
gemmi sf2map [options] INPUT_FILE MAP_FILE
gemmi sf2map --check INPUT_FILE
INPUT_FILE must be an MTZ or mmCIF file 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.
Option --check can take column name mapping as arg, for example:
--check=dfc:FCALC,fosc:FOBS (typically, use --check without an arg).
Possible keys are fwt, phwt, delfwt, phdelwt, dfc, fosc, fom, free.
Options:
-h, --help Print usage and exit.
-V, --version Print version and exit.
-v, --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).
--weight=COLUMN (normally not needed) weighting for F.
-g, --grid=NX,NY,NZ Minimum grid size.
--exact Use the exact grid size specified with --grid.
-s, --sample=NUMBER Set spacing to d_min/NUMBER (3 is common).
--zyx Invert axis order in output map: Z=fast and X=slow.
-G Print size of the grid that would be used and exit.
--timing Print calculation times.
--normalize Scale the map to standard deviation 1 and mean 0.
--mapmask=FILE Output only map covering the structure from FILE,
similarly to CCP4 MAPMASK with XYZIN.
--margin=N (w/ --mapmask) Border in Angstroms (default: 5).
--select=SEL (w/ --mapmask) Atom selection for mask, MMDB syntax.
--check[=cols] Analyze map coefficient columns in MTZ file.
The --sample
option is named after the GRID SAMPLE
keyword of
the venerable CCP4 FFT program; its value has the same meaning.
--check
¶
This option is inspired by Ian Tickle’s MTZFIX program. It checks how the map coefficients were calculated by determining relationships between MTZ columns, separately for acentric, centric and missing reflections. It works well with MTZ files from BUSTER, Refmac and Servalcat, which include map coefficients along with scaled Fobs, D.Fc and the figure-of-merit m. Here is an example output:
Columns used in checking map coefficients:
for FM (normal map): FWT PHWT
for FD (difference map): DELFWT PHDELWT
for D.Fc: FC_ALL PHIC_ALL
for scaled Fo: FP
for figure-of-merit m: FOM
for free flags: FREE
Is FD (DELFWT) set to 0 or NaN for any of 20 FREE flags ... yes, for 0
-> free reflections are NOT used for maps
Phases PHWT and PHDELWT ... match (mod 180)
Phases PHIC_ALL and PHWT ... match (mod 180)
For all 127888 acentric reflections (excl. missing/unused):
FM = 2m.Fo - D.Fc
FD = 2(m.Fo - D.Fc) = 2(FM - m.Fo) = FM - D.Fc
For all 2486 centric reflections (excl. missing/unused):
FM = 2m.Fo - D.Fc
FD = 2(m.Fo - D.Fc) = 2(FM - m.Fo) = FM - D.Fc
For all 16359 missing/unused reflections:
FM = D.Fc
FD = 0
If the MTZ file lacks some of these columns (for instance, if only unscaled Fobs are present), the output is less informative.
Why is this useful at all?
Typically, an MTZ file from refinement has coefficients for two maps:
a density map, 2mFo–DFc , sometimes referred to as 2Fo–Fc even if m and D are used,
and a difference map, mFo–DFc , sometimes referred to as Fo–Fc.
However, the exact formulas vary across implementations:
For acentric reflections, the difference map may or may not include a factor of 2. In the literature, the formula is actually 2(mFo–DFc), not the one given above. In an MTZ file it may not be 2×. (This only affects the map’s absolute values – not that important.)
The formulas for centric reflections have been derived as: mFo for the density map and mFo–DFc (without 2×) for the difference map; however, not all programs use these formulas.
For reflections without observed data (missing Fo), the density map coefficients are either filled with DFc or not.
For reflections in a free set, Fo is either used (to improve map quality) or not (to ensure that the free set is used solely for validation).
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.
-v, --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).
--spacegroup=SG Overwrite space group from map header.
merge¶
Merge intensities from multi-record reflection file.
$ gemmi merge -h
Usage:
gemmi merge [options] INPUT_FILE OUTPUT_FILE
gemmi merge --compare [options] UNMERGED_FILE MERGED_FILE
gemmi merge --compare [options] MMCIF_FILE_WITH_BOTH
Options:
-h, --help Print usage and exit.
-V, --version Print version and exit.
-v, --verbose Verbose output.
--anom output/compare I(+) and I(-) instead of IMEAN.
--no-sysabs Do not output systematic absences.
--nobs Add MTZ column NOBS with the number of merged
reflections.
-b NAME, --block=NAME output mmCIF block name: data_NAME (default: merged).
--compare compare unmerged and merged data (no output file).
--print-all print all compared reflections.
The input file can be SF-mmCIF with _diffrn_refln, MTZ or XDS_ASCII.HKL.
The output file can be either SF-mmCIF or MTZ.
ecalc¶
Calculates normalized amplitudes E from amplitudes F. Uses “Karle” approach, similar to CCP4 ECALC.
$ gemmi ecalc -h
Usage:
gemmi ecalc [options] INPUT.mtz OUTPUT.mtz
Calculates E (normalised structure amplitudes) from F.
-h, --help Print usage and exit.
-V, --version Print version and exit.
-v, --verbose Verbose output.
--F=LABEL Label of the input column (default: F).
--E=LABEL Label of the output column (default: E).
--no-sigma Do not use sigma columns (SIGF->SIGE).
--method=METHOD Method of calculating <F^2>, one of:
ma - moving average (not implemented yet),
2 - resolution bins equispaced in 1/d^2,
3 - resolution bins in 1/d^3 (default),
ec - bins with equal count of reflections.
--binsize=N Number of reflections per bin or in moving window
(default: 200).
--maxbins=N Maximum number of bins (default: 100).
Column for SIGF is always the next column after F (the type must be Q).
If INPUT.mtz has column E, it is replaced in OUTPUT.mtz.
Otherwise, a new column is appended.
The name for SIGE, if it is added as a column, is 'SIG' + label of E.
sfcalc¶
Calculates structure factors from a model.
$ gemmi sfcalc -h
Usage:
gemmi sfcalc [options] INPUT_FILE
Calculates structure factors of a model (PDB, mmCIF or SMX CIF file).
Uses FFT to calculate all reflections up to requested resolution for MX
files. Otherwise, for SMX and --hkl, F's are calculated directly.
This program can also compare F's calculated directly with values
calculated through FFT or with values read from a reflection file.
Options:
-h, --help Print usage and exit.
-V, --version Print version and exit.
-v, --verbose Verbose output.
--hkl=H,K,L Calculate structure factor F_hkl.
--dmin=NUM Calculate structure factors up to given resolution.
--for=TYPE TYPE is xray (default), electron, neutron or mott-bethe.
--normalize-it92 Normalize X-ray form factors (a tiny change).
--use-charge Use X-ray form factors with charges when available.
--ciffp Read f' from _atom_type_scat_dispersion_real in CIF.
--wavelength=NUM Wavelength [A] for calculation of f' (use --wavelength=0
or -w0 to ignore anomalous scattering).
--unknown=SYMBOL Use form factor of SYMBOL for unknown atoms.
--noaniso Ignore anisotropic ADPs.
--margin=NUM For non-crystal use bounding box w/ margin (default: 10).
Options for density and FFT calculations (with --dmin):
--rate=NUM Shannon rate used for grid spacing (default: 1.5).
--blur=NUM B added for Gaussian blurring (default: auto).
--rcut=Y Use atomic radius r such that rho(r) < Y (default: 1e-5).
--test[=CACHE] Calculate exact values and report differences (slow).
--write-map=FILE Write density (excl. bulk solvent) as CCP4 map.
--to-mtz=FILE Write Fcalc to a new MTZ file.
Options for anisotropic scaling (only w/ FFT):
--scale-to=FILE:COL Anisotropic scaling to F from MTZ file.
Argument: FILE[:FCOL[:SIGFCOL]] (defaults: F and SIGF).
--sigma-cutoff=NUM Use only data with F/SIGF > NUM (default: 0).
Options for bulk solvent correction and scaling (only w/ FFT):
--d-mask=NUM Spacing of mask grid (default: same as for model).
--radii-set=SET Set of per-element radii, one of: vdw, cctbx, refmac.
--r-probe=NUM Value added to VdW radius (default: 1.0A).
--r-shrink=NUM Value for shrinking the solvent area (default: 1.1A).
--solvent-mask=FILE Use mask from file instead of calculating it.
--ksolv=NUM Value (if optimizing: initial value) of k_solv.
--bsolv=NUM Value (if optimizing: initial value) of B_solv.
--kov=NUM Value (if optimizing: initial value) of k_overall.
--baniso=B11:...:B23 Anisotropic scale matrix (6 colon-separated numbers: B11,
B22, B33, B12, B13, B23).
Options for comparing calculated values with values from a file:
--compare=FILE Re-calculate Fcalc and report differences.
--f=LABEL MTZ column label (default: FC) or small molecule cif tag
(default: F_calc or F_squared_calc).
--phi=LABEL MTZ column label (default: PHIC).
In general, structure factors can be calculated
either directly, by summing contributions from each atom to each reflection,
or by calculating an electron density on a grid and using discrete Fourier transform.
This program can measure the errors resulting from the latter method (in addition to its main function – calculation of the structure factors). The errors depend on
the grid spacing – controlled by the oversampling
--rate=R
; the maximum spacing is dmin/2R,atomic radius – we neglect electron density of the atom beyond this radius; only density contributions above the (absolute) value specified with
--rcut
are taken into account,Gaussian dampening (blurring) factor – artificial temperature factor Bextra added to all atomic B-factors (the structure factors are later corrected to cancel it out); either specified with
--blur
or picked automatically.
Choosing these parameters is a trade-off between efficiency and accuracy,
as described elsewhere.
The option --test
can be used to see how accuracy and efficiency
depends on the choice of parameters. For example, this shell script
performs a series of calculations with differing Bextra:
model=1mru.pdb
dmin=2.5
gemmi sfcalc --dmin=$dmin --test $model >cache.tsv
for i in `seq -20 5 20`; do
printf -- "$i\t" >&2
gemmi sfcalc --dmin=$dmin --rate=1.5 --rcut=1e-4 --blur=$i --test=cache.tsv $model
done >/dev/null
Running it prints:
-20 RMSE=0.93304 0.5495% max|dF|=38.80 R=0.301% 0.27671s
-15 RMSE=0.37007 0.2179% max|dF|=41.26 R=0.094% 0.28366s
-10 RMSE=0.27075 0.1595% max|dF|=44.35 R=0.041% 0.29322s
-5 RMSE=0.27228 0.1604% max|dF|=47.59 R=0.029% 0.30459s
0 RMSE=0.28903 0.1702% max|dF|=50.95 R=0.029% 0.31399s
5 RMSE=0.30806 0.1814% max|dF|=54.35 R=0.032% 0.32527s
10 RMSE=0.32847 0.1934% max|dF|=57.92 R=0.036% 0.33360s
15 RMSE=0.35028 0.2063% max|dF|=61.66 R=0.041% 0.34181s
20 RMSE=0.37283 0.2196% max|dF|=65.44 R=0.047% 0.35380s
The error used in RMSE is the magnitude of the difference of two vectors: |Fapprox – Fexact|. The next column is RMSE normalized by the sum of |Fcalc|. Then we have maximum error for a single reflection, and the wall time of computations. We can see that in this case negative “dampening” (subtracting about 10A2 from all B-factors) improves both accuracy and performance.
fprime¶
Calculate anomalous scattering factors (f' and f"). Uses Cromer-Liberman algorithm with corrections from Kissel and Pratt. This and different approaches are discussed in the documentation of the underlying functions.
$ gemmi fprime -h
Usage:
gemmi fprime [options] ELEMENT[...]
Prints anomalous scattering factors f' and f".
Options:
-h, --help Print usage and exit.
-V, --version Print version and exit.
-e, --energy=ENERGY Energy [eV]
-w, --wavelength=LAMBDA Wavelength [A]
Here is an example how to print f’ and f” using gemmi, XrayDB, CCP4 crossec and cctbx (pyFprime is not included because it is a GUI-only program). The Chantler’s data from XrayDB is probably the most reliable one:
$ gemmi fprime --wavelength=1.2 Se
Element E[eV] Wavelength[A] f' f"
Se 10332.0 1.2 -1.4186 0.72389
$ python3 -c "import xraydb; print(xraydb.f1_chantler('Se', 10332.0), xraydb.f2_chantler('Se', 10332.0))"
-1.4202028957329489 0.7100533627953146
$ echo -e "atom SE\n cwav 1 1.2 0\n END" | crossec | grep ^SE
SE 1.2000 -1.5173 0.7240
$ cctbx.eltbx.show_fp_fdp --wavelength=1.2 --elements=Se
Wavelength: 1.2 Angstrom
Element: Se
Henke et al. : f'=-1.44568 , f''=0.757958
Sasaki et al. : f'=-1.5104 , f''=0.724000
diff f''=-2.29 %
reindex¶
Experimental. See description of the reindexing function for more details.
Reindex reflections in MTZ file.
$ gemmi reindex -h
Usage:
gemmi reindex [options] INPUT_MTZ OUTPUT_MTZ
Options:
-h, --help Print usage and exit.
-V, --version Print version and exit.
-v, --verbose Verbose output.
--hkl=OP Reindexing transform as triplet (e.g. k,h,-l).
--no-history Do not add 'Reindexed with...' line to mtz HISTORY.
--no-sort Do not reorder reflections.
--asu=ccp4|tnt Write merged data in CCP4 (default) or TNT ASU.
Input file can be gzipped.
residues¶
List residues from a coordinate file, one per line.
$ gemmi residues -h
Usage:
gemmi residues [options] INPUT[...]
Prints chains, residues and atoms. Or without atoms. Or only entities.
-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.
-l, --label Print 'label' chain ID and seq ID in brackets.
--check-seqid Check if sequence IDs are unique and exit.
--no-alt Do not print altlocs.
-s, --short Shorter output (no atom info). Can be given 2x or 3x.
-e, --entities List (so-called, in mmCIF speak) entities.
-c, --chains List chain IDs.
INPUT is a coordinate file (mmCIF, PDB, etc).
The optional 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
Note the options -s
(short) that can be used up to 3 times,
making the output more concise:
$ echo $PDB_DIR
$ gemmi residues -ss 1mru
/data/structures/divided/mmCIF/mr/1mru.cif.gz
A polymer THR PRO SER HIS LEU SER ASP ARG TYR GLU... (269 residues)
B polymer MET THR THR PRO SER HIS LEU SER ASP ARG... (271 residues)
A non-polymers MG MG AGS
B non-polymers MG MG AG
Option -e
lists so-called entities:
$ gemmi residues -e 3oov
/data/structures/divided/mmCIF/oo/3oov.cif.gz
Polymers
entity 1, polypeptide(L), length 169, subchains:
- A from strand A, 164 residues: 6-169
- B from strand B, 162 residues: 1-166 except 138-141
Others
entity 2, non-polymer (GOL), subchains: C D E F
entity 3, water (HOH), subchains: G H
align¶
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 align -h
Pairwise sequence alignment with scoring matrix and affine gap penalty.
Usage:
gemmi align [options] FILE[...]
Aligns sequence from the model to the full sequence (SEQRES).
Both are from the same FILE - either in the PDB or mmCIF format.
If the mmCIF format is used, option --check-mmcif can be used.
gemmi align [options] --query=CHAIN1 --target=CHAIN2 FILE [FILE2]
Aligns CHAIN1 from FILE to CHAIN2 from FILE2 (if given) or FILE.
By default, the sequence of residues in the model is used.
To use SEQRES prepend '+' to the chain name (e.g. --query=+A),
or, when using mmCIF, prepend '@' to entity name (--query=@1).
gemmi align [options] --text-align STRING1 STRING2
Aligns two ASCII strings (used for testing).
Options:
-h, --help Print usage and exit.
-V, --version Print version and exit.
--check-mmcif checks alignment against _atom_site.label_seq_id
--query=[+|@]CHAIN Align CHAIN from file INPUT1.
--target=[+|@]CHAIN Align CHAIN from file INPUT2.
--text-align Align characters in two strings (for testing).
Scoring (absolute values):
--blosum62 Use BLOSUM62 score matrix.
--partial=y|n Use scoring meant to align partially-modelled polymer to
its full sequence (default in 1st mode).
--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).
Output options:
-p Print formatted alignment with one-letter codes.
--rmsd In addition to aligning two CHAINs (--query and
--target), superpose structures and print RMSD.
-v, --verbose Verbose output.
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 align -p --match=0 --gapo=0 --text-align Saturday Sunday
Score: -3 CIGAR: 1M2I5M
Saturday
| |.|||
S--unday
This tool uses modified code from ksw2. See the Sequence alignment section for more details.
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.
-v, --verbose Verbose output.
--asu=N Draw ASU in NxNxN map grid and exit. Uses N(N+1) columns.
contents¶
Analyzes and summarizes the content of a coordinate file.
Inspired by the CCP4 program rwcontents
.
By default, it prints the atom count, estimated number of hydrogens in the protein, molecular weight of the protein, ASU volume, Matthews coefficient, and the fractional volume of solvent in the crystal.
It has options to print other information – see the help message below.
$ 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.
-v, --verbose Verbose output.
--select=SEL Use only the selection.
-b Print statistics of isotropic ADPs (B-factors).
--dihedrals Print peptide dihedral angles.
-n Do not print content (for use with other options).
To print a list of chains and residues, or entities, use residues.
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.
-v, --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].
--minocc=MIN Ignore atoms with occupancy < MIN.
--ignore=N Ignores atom pairs from the same: 0=none, 1=residue, 2=same or
adjacent residue, 3=chain, 4=asu.
--nosym Ignore contacts between symmetry mates.
--assembly=ID Output bioassembly with given ID (1, 2, ...).
--noh Ignore hydrogen (and deuterium) atoms.
--nowater Ignore water.
--noligand Ignore ligands and water.
--count Print only a count of atom pairs.
--twice Print each atom pair A-B twice (A-B and B-A).
--sort Sort output by distance.
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.
-v, --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).
--weight=COLUMN (normally not needed) weighting for F.
-g, --grid=NX,NY,NZ Minimum grid size.
--exact Use the exact grid size specified with --grid.
-s, --sample=NUMBER Set spacing to d_min/NUMBER (3 is common).
-G Print size of the grid that would be used and exit.
--timing Print calculation times.
h¶
Adds or removes hydrogens.
Hydrogen positions are determined from restraints from a monomer library. If the restraints allow various directions to a hydrogen atom, that atom is assigned one of possible positions and its occupancy is zeroed.
$ 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.
-v, --verbose Verbose output.
--monomers=DIR Monomer library directory (default: $CLIBD_MON).
-L CIF, --lib=CIF User's restraint file(s). See more info below.
--format=FORMAT Input format (default: from the file extension).
--sort Order atoms in residues according to _chem_comp_atom.
--update If deprecated atom names (from _chem_comp_atom.alt_atom_id)
are used in the model, change them.
Hydrogen options, mutually exclusive. Default: add hydrogens, but not to water.
--water Add hydrogens also to water.
--unique Add only hydrogens with uniquely determined positions.
--keep Do not add/remove hydrogens, only change positions.
--no-change Do not change hydrogens, not even positions.
--remove Only remove hydrogens.
Option -L/--lib can be used multiple times, in order of priority.
Its argument is either a file path or one of the two special values:
'+' = monomer blocks in mmCIF INPUT_FILE (ignored by default)
'@' = the priority of the monomer library (ML, default: lowest)
Example 1: -L file.cif -L+ order: file.cif, input file, ML
Example 2: -L@ -L file.cif order: ML, file.cif
To test the addition of hydrogen to a single chemical component from a monomer library, run:
gemmi h --format=chemcomp -L+ XYZ.cif output.pdb
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.
-v, --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).
crd¶
Prepares an intermediate (a.k.a. prepared, crd or topology) file for Refmac.
It reads a coordinate file, a monomer library and, optionally,
additional monomer CIF files. It writes out a topology
file in a format understood by CCP4 Refmac.
Using it as XYZIN requires adding make coordinates prepared
:
$ gemmi crd in.pdb in.crd
$ refmac5 xyzin in.crd hklin in.mtz xyzout out.pdb hklout out.mtz << eof
make coordinates prepared
# other refmac keywords ("refinement type restrained", "ncycle 10", etc)
end
eof
Such intermediate step is always present in Refmac, even if it is not visible
to the user. In normal operation, Refmac writes, reads and deletes temporary
files with extensions crd and rst. gemmi crd
writes a single file
with the content of both crd and rst.
$ gemmi crd -h
Usage:
gemmi crd [options] INPUT_FILE OUTPUT_FILE
Prepare intermediate Refmac files.
INPUT_FILE can be in PDB, mmCIF or mmJSON format.
Options:
-h, --help Print usage and exit.
-V, --version Print version and exit.
-v, --verbose Verbose output.
--monomers=DIR Monomer library directory (default: $CLIBD_MON).
-L CIF, --lib=CIF User's restraint file(s). See more info below.
--auto-cis=Y|N Assign cis/trans ignoring CISPEP record (default: Y).
--auto-link=Y|N Find links not included in LINK/SSBOND (default: N).
--auto-ligand=Y|N If ligand has no definition make ad-hoc restraints (N).
--no-aliases Ignore _chem_comp_alias.
Hydrogen options (default: remove and add on riding positions):
-H, --no-hydrogens Remove (and do not add) hydrogens.
--keep-hydrogens Preserve hydrogens from the input file.
Option -L/--lib can be used multiple times, in order of priority.
Its argument is either a file path or one of the two special values:
'+' = monomer blocks in mmCIF INPUT_FILE (ignored by default)
'@' = the priority of the monomer library (ML, default: lowest)
Example 1: -L file.cif -L+ order: file.cif, input file, ML
Example 2: -L@ -L file.cif order: ML, file.cif
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 altogether. 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’s CC as well as the 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 experimental 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’s 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.
-v, --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).
xds2mtz¶
Converts XDS ASCII file to MTZ format. Optionally, filters the reflections and applies polarization correction.
$ gemmi xds2mtz -h
Usage:
gemmi xds2mtz [options] XDS_FILE MTZ_FILE
Options:
-h, --help Print usage and exit.
-V, --version Print version and exit.
-v, --verbose Verbose output.
--title MTZ title.
-H LINE, --history=LINE Add a history line.
--project=PROJECT Project in MTZ hierarchy (default: 'XDSproject')
--crystal=CRYSTAL Crystal in MTZ hierarchy (default: 'XDScrystal')
--dataset=DATASET Dataset in MTZ hierarchy (default: 'XDSdataset')
--batchmin=BATCHMIN Delete reflections with BATCH<BATCHMIN (default: 1)
Polarization correction and overload elimination options for INTEGRATE.HKL
files:
--polarization=VALUE XDS parameter FRACTION_OF_POLARIZATION
--normal='Pnx Pny Pnz' XDS POLARIZATION_PLANE_NORMAL (default: '0 1 0')
--overload=OVERLOAD XDS parameter OVERLOAD to eliminate reflections with
MAXC>OVERLOAD
If XDS_FILE is -, the input is read from stdin.