chem.cr

GitHub Workflow Status Version License

A modern library written in Crystal primarily designed for manipulating molecular files created by computational chemistry programs. It aims to be both fast and easy to use.

IMPORTANT: this library is in alpha stage, meaning that there is missing functionality, documentation, etc. and there will be breaking changes.

Features | Installation | Usage | Tools | Benchmarks | Roadmap | Testing | Similar software | Contributing | Contributors | License

Features

Installation

Add this to your application's shard.yml:

dependencies:
  chem:
    github: franciscoadasme/chem.cr
    version: ~> 0.1

Requirements

Usage

First require the chem module:

require "chem"
include Chem # avoids typing Chem:: prefix

Let's first read a structure:

st = Strucuture.read "/path/to/file.pdb"
st # => <Structure "1cbn": 644 atoms, 47 residues, periodic>

You can also use a custom read method that accepts specific options:

Structure.from_pdb "/path/to/file.pdb"
Structure.from_pdb "/path/to/file.pdb", chains: ['A'] # read only chain A
Structure.from_pdb "/path/to/file.pdb", het: false # skip HET atoms
Structure.from_pdb "/path/to/file.pdb", alt_loc: 'A' # select alternate location A

You can access PDB header information via the #experiment property:

if expt = st.experiment # check if experiment data is present
  expt.title # => "ATOMIC RESOLUTION (0.83 ANGSTROMS) CRYSTAL STRUCTURE..."
  expt.kind # => XRayDiffraction
  expt.resolution # => 0.83
  expt.deposition_date # => 1991-10-11
  ...
end

You can also read many structures at once:

# read all models
Array(Structure).from_pdb "/path/to/file.pdb"
# read 2th and 5th models
Array(Structure).from_pdb "/path/to/file.pdb", indexes: [1, 4]

Alternatively, you could use an IO iterator to read one by one:

PDB::Parser.new("/path/to/file.pdb").each { |st| ... }
PDB::Parser.new("/path/to/file.pdb").each(indexes: [1, 4]) { |st| ... }

Topology access

You can access topology objects using the bracket syntax (like a hash or associative array or dictionary):

st['A'] # => <Chain A>
st['A'][10] # => <Residue A:ARG10>
st['A'][10]["CA"] # => <Atom A:ARG10:CA(146)>

Alternatively, you can use the #dig and #dig? methods:

st.dig 'A' # => <Chain A>
st.dig 'A', 10 # => <Residue A:ARG10>
st.dig 'A', 10, "CA" # => <Atom A:ARG10:CA(146)>

st.dig 'A', 10, "CJ" # causes an error because "CJ" doesn't exist
st.dig? 'A', 10, "CJ" # => nil

Each topology object have several modifiable properties:

atom = st.dig 'A', 10, "CA"
atom.element.name # => Carbon
atom.coords # => [8.47 4.577 8.764]
atom.occupancy # => 1.0
atom.bonded_atoms.map &.name # => ["N", "C", "HA", "CB"]

Thanks to Crystal's powerful standard library, manipulating topology objects is very easy:

# ramachandran angles
st.residues.map { |r| {r.phi, r.psi} } # => [{129.5, 90.1}, ...]
# renumber residues starting from 1
st.residues.each_with_index { |res, i| res.number = i + 1 }
# constrain Z-axis
st.atoms.each { |atom| atom.constraint = :z }
# total charge
st.atoms.sum_of &.partial_charge
# iterate over secondary structure elements
st.residues.chunk(&.sec).each do |sec, residues|
  sec # => HelixAlpha
  residues # => [<Residue A:ARG1>, <Residue A:LEU2>, ...]
end

Here #residues and #atoms return an array of Residue and Atom instances, respectively. Collections also provide iterator-based access, e.g., #each_atom, that avoids expensive memory allocations:

st.atoms.any? &.constraint # array allocation to just check a condition
st.each_atom.any? &.constraint # faster!

Atom selection

Right now, there is no custom language to select a subset of atoms. However, thanks to Crystal, one can achieve a similar result with an intuitive syntax:

st.atoms.select { |atom| atom.partial_charge > 0 }
# or
st.atoms.select &.partial_charge.>(0)
# compared to a custom language
st.atoms.select "partial_charge > 0"

# select atoms within a cylinder of radius = 4 A and centered at the origin
st.atoms.select { |atom| atom.x**2 + atom.y**2 < 4 }
# compared to a custom language
st.atoms.select "sqrt(x) + sqrt(y) < 4" # or "x**2 + y**2 < 4"

One advantage to using Crystal itself is that it provides type-safety: doing something like atom.name**2 will result in a compilation error, whereas using a custom language will probably produce a confusing error during runtime. Additionally, the code block can be as big and complex as necessary with multiple intermediary computations. Furthermore, a negative condition may be confusing and not be trivial to write, but in Crystal you would simply use #reject instead.

Finally, the above also works for chain and residue collections:

# select protein chains
st.chains.select &.each_residue.any?(&.protein?)
# select only solvent residues
st.residues.select &.solvent?
# select residues with any atom within 5 A of the first CA atom
# (this is equivalent to "same residue as" or "fillres" in other libraries)
ca = st.dig 'A', 1, "CA"
st.residues.select do |res|
  res.each_atom.any? { |atom| Spatial.distance(atom, ca) < 5 }
end
# or
st.atoms.select { |atom| Spatial.distance(atom, ca) < 5 }.residues

Coordinates manipulation

All coordinates manipulation is done using a CoordinatesProxy instance, available for any atom collection (i.e., structure, chain or residue) via #coords:

# geometric center
st.coords.center
# center at origin
st.coords.translate! -st.coords.center
# wraps atoms into the primary unit cell
st.coords.wrap
...

Tools

Additional tools built upon the chem module are available in the repository in the cli directory. These tools can be download individually from the releases page or built from source by running the shards build command on a local copy.

Benchmarks

chem.cr is implemented in pure Crystal, making it as fast or even faster than some C-powered packages.

The benchmark is designed as follows:

IMPORTANT: direct comparison of parsing times should be taken with a grain of salt because each library does something slightly different, e.g., error checking. Some of this functionality is listed below. Nonetheless, these results gives an overall picture in terms of the expected performance.

| | Biopython | chem.cr | Chemfiles | MDAnalysis | MDTraj | schrodinger | VMD | | -------------------- | --------: | ------: | --------: | ---------: | -----: | ----------: | ----: | | Parse 1CRN [ms] | 6.521 | 1.028 | 1.668 | 5.059 | 11.923 | 45.497 | 2.285 | | Parse 3JYV [s] | 0.837 | 0.086 | 0.199 | 0.404 | 1.490 | 0.766 | 0.162 | | Parse 1HTQ [s] | 16.146 | 1.673 | 2.540 | 1.387 | 18.969 | 11.997 | 0.236 | | Count [ms] | 0.210 | 0.009 | 0.322 | 0.041 | 0.079 | 25.997 | 0.165 | | Distance [ms] | 0.172 | 0.000 | 1.016 | 0.382 | 0.990 | 43.101 | 0.379 | | Ramachandran [ms] | 110.450 | 0.607 | - | 690.201 | 4.947 | 68.758 | 1.814 | | | | | | | | | | | License | Biopython | MIT | BSD | GPLv2 | LGPL | Proprietary | VMD | | Parse Header | yes | yes | yes | no | no | no | no | | Parse CONECT | no | yes | yes | no | yes | yes | yes | | Guess bonds | no | no | yes | no | yes | yes | yes | | Hybrid36 | no | yes | no | yes | no | no | no | | Hierarchical parsing | yes | yes | no | no | no | no | no | | Supports disorder | yes | yes | no | no | yes | yes | no |

Latest update: 2019-11-10

Scripts and details are provided at pdb-bench.

Roadmap

Topology manipulation

Input and Output

File formats

Analysis

Other

Testing

Run the tests with crystal spec.

Similar software

There are several libraries providing similar functionality. Here we list some of them and provide one or more reasons for why we didn't use them. However, each one of them is good in their own right so please do check them out if chem.cr does not work for you.

Contributing

  1. Fork it (https://github.com/franciscoadasme/chem.cr/fork)
  2. Create your feature branch (git checkout -b my-new-feature)
  3. Commit your changes (git commit -am 'Add some feature')
  4. Push to the branch (git push origin my-new-feature)
  5. Create a new Pull Request

Contributors

License

Licensed under the MIT license, see the separate LICENSE file.