hPDB – Haskell library for processing atomic biomolecular structures in protein data bank format

Background Protein DataBank file format is used for the majority of biomolecular data available today. Haskell is a lazy functional language that enjoys a high-level class-based type system, a growing collection of useful libraries and a reputation for efficiency. Findings I present a fast library for processing biomolecular data in the Protein Data Bank format. I present benchmarks indicating that this library is faster than other frequently used Protein Data Bank parsing programs. The proposed library also features a convenient iterator mechanism, and a simple API modeled after BioPython. Conclusion I set a new standard for convenience and efficiency of Protein Data Bank processing in a Haskell library, and release it to open source.


Background
The Protein Data Bank (PDB) is a widely used data repository of atomic resolution, three-dimensional protein and nucleic acid structures [1]. The rapid growth of structural data enables key endeavors to bring knowledge of genomes [2] to the structure and function of large biomolecules. In addition to sequence searches and genome assemblies, efficient and reliable structural data processing are one of the most important and common structural bioinformatics tasks [3].
Haskell is a modern, lazy, pure functional language [4,5] that enjoys fluid syntax, and clarity comparable to Python [6], as well as an efficient compiler that often generates code approaching the speeds of industry standard languages such as C [7] or C++ [8].

Library interface
The library is a comprehensive solution for the parsing, rapid processing and writing of PDB files. I introduce the Correspondence: miga@nmr.mpibpc.mpg.de NMR-2, Max Planck Institute for Biophysical Chemistry, Am Faßberg 11, Göttingen, Germany library by providing examples and describing the underlying data structures a , and finally, I present an evaluation of its efficiency.

Simple use example
A parser example b -a script that reads a structure containing multiple models and splits the structure into single models. A simple PDB.parse action returns a structure, to which I apply splitModels and zip a list of results with model numbers. These results are used to generate names of output files, that are then written using the PDB.write IO action. ByteString is used rather than [Char] within the library for everything except file paths (FilePath), due to efficiency considerations e .

Data structure describing molecules
Different levels of collection hierarchy can be seen on Figure 1: Structure that contains information about whole PDB entry; Model that shows a single model of the molecule; Chain describing a single polymer chain; Residue for a single monomer (aminoacid residue or nucleic acid base) within the polymer; Atom for a single atom location.
Names of these types correspond to the names used by PDB file format definition [1]. Those atoms which may have multiple locations within the model are described by several records, and those residues that have alternative mutants are also described by different records in accord with current practice of PDB [1].

Iterating with Iterable
For the different types of objects, I devised a custom Iterable class that allows iteration over all objects of each of many types contained within an argument object. This class generalizes map, and foldr iteration patterns used over lists, to hierarchical containers, and allows iteration over any type of contained objects, potentially in an indirect manner: The itmap method allows mutation of any of the objects of a given type b contained within type a. To compute a statistic over all contained objects itfoldr or itfoldl can be used.
The class Iterable a b may thus be viewed as a generalization of Functor, Foldable, and Traversable to hierarchical collections that may contain more than one class of objects. All Iterable instances, and method types are shown in Figure 1.
It is often convenient to use a monadic state when computing a statistic, or when renumbering residues. For this purpose, I use monadic variants of the following methods: For efficiency I introduce a rapid itlength method, and a strict itfoldl variant: itfoldl' :: Note that itlength is the only method using a dummy first argument to indicate the type of contained object to be counted. As all other methods use a function argument, automatic type inference finds the proper class instance without requiring a type declaration, as shown in the examples below f .

Structure analysis example
In the following examples I skip the command line interface, assuming that all functions input a parsed Structure object. The most convenient interface for a complex cascade of container types within a PDB structure is composition based on fold, and map analogs.
To compute the center of mass of all Atom objects contained within a structure, I use a two pass algorithm: Here I use itfoldl' instantiated to the following type, automatically inferred from the types of addCoord and s: I will generalize this type to other types within the structure, when showing class Iterable. This generalization allows a function to have a type showing that it can take any object that contains Atom objects within: center :: (Iterable a Atom) => a -> Vector3 Then, I can subtract the computed center from all atomic coordinates with the itmap method, analogous to map for lists. In this example, mapping PDB.Atom objects within a PDB.Structure g , is used: This use of itmap has an instantiated type, that is automatically inferred from its first argument: itmap :: (Atom -> Atom) -> Structure -> Structure

Stateful modification
Simple itmap, and itfoldl' methods are not sufficient to perform a complex stateful operation such as renumbering residues starting from 1.
In this case, I use monadic analogs, such as itmapM h : Renumbering atoms within each model is more involved, because the PDB data format [10] mandates that TER records i are counted along with the atoms, and these records do not have direct representation in hPDB data structures. renumberAtoms = itmap (\m -> runCounterM $ itmapM forEachChain (m :: Model)) where forEachChain :: Chain -> CounterM Chain forEachChain ch = do newCh <-itmapM forEachAtom ch modify (+1) --for TER return newCh forEachAtom :: Atom -> CounterM Atom forEachAtom at = do v <-get modify (+1) return $ at { atSerial = v } In this example renumberAtoms may be labelled with the monomorphic (and useful) type: Although the automatically inferred type allows this function to act not only on the entire Structure but also on any single Model that contains the Chain objects.

Example applications
hPDB's speed and ease of use has allowed for rapid implementation of typical functions such as: orienting structure so that the longest diameter corresponds to the Y axis, http://www.biomedcentral.com/1756-0500/6/483  Table 3). and the second longest cross-sectional dimension corresponds to the X axis (CanonicalAxes in hPDB-examples package), normalizing PDB files for use by applications restrictive with respect to file format (CleanPDB), and examining the sequence of main polymer chain or geometric parameters of small-angle scattering shape reconstructions (Rg example) with minimal code.
In the case of libraries, I used operating system calls or ps program to determine the upper bounds of memory used in Table 2 (including purely virtual allocations).
Haskell memory is reported for the current heap, in addition to the target space for copying garbage collector [17].
Note that Jmol and BioJava may both use more than one thread, which significantly reduces time-to-completion when using a multicore machine as reported in Table 3.
The benchmarks were measured on a quad-core Intel® Core™ i7 2600 processor running at 3.4 GHz j , 16 GB of 1333 MHz memory, and a SAMSUNG 470 Series solid-state disk. The system was running a 64-bit Ubuntu 13.04 with a standard Linux kernel package 3.8.0-31.
While hPDB may be expected to stand out in runtime comparisons to the bytecode-based dynamic language libraries BioRuby and BioPython, surprisingly, serial hPDB is faster than other parsers in compiled languages, with the exception of PyMol. The parallel version of the hPDB parser may be the fastest PDB parser on machines with at least 4 independent processing cores.
It was noted that memory use, even with a necessary overhead (2×) of Haskell's copying garbage collector, compared favorably with memory used by other libraries.
Parsing the entire PDB archive (as of January 6th 2013, compressed, 16 GB) takes approximately 14.5 minutes using 4 cores in parallel, with total CPU and I/O time reported to be 50 minutes. No crashes are reported, but 8k lines (mostly meta data) are reported as erroneous k because they are inconsistent with strict interpretation of PDB format [10].
Benchmarks show that in this specific application, the mildly optimized Haskell parser may provide speeds competitive with compiled languages such as Java and even lower level explicitly allocated languages such as C. Memory usage is also less than any other aforementioned library.
There is another Haskell library parsing PDB files on Hackage [18] called PDBtools, but it was not able to fully parse any of our example files because it does not handle errors in the read routine.

Conclusions
I have shown clear uses of a nice high-level interface for the analysis and modification of molecule descriptions encoded in the PDB file format [10].
While there are many similar parsers written in other languages, this is the first one I am aware of in Haskell, that parses entire coordinate contents within the PDB repository. It is also efficient both in runtime and memory use, and thus, the preferable choice for sophisticated, high volume analyses.  While future work on analysis API extensions would likely further improve utility of this library, I believe that it is ready for production use, as indicated by the many code examples.
I conclude that in this specific application, Haskell has both ease of use and abstraction of high-level dynamic languages, along with a speed competitive with lower level explicit-allocation languages such as C.

Availability and requirements
Source code is available as Additional files 1, 2 and 3 attached to the manuscript or from GitHub repository https://github.com/mgajda/hPDB, and released on Hackage as hPDB. It has been tested with several GHC versions including 7.0.3, 7.2.2, 7.4.2, and the recently released 7.6.2. It has few dependencies, and all are available from Hackage [18].