Generic Software Frameworks for Gc-ms Based Metabolomics

Metabolomics has seen a rapid development of new technologies, methodologies, and data analysis procedures during the past decade. The development of fast gasand liquid-chromatography devices coupled to sensitive mass-spectrometers, supplemented by the unprecedented precision of nuclear magnetic resonance for structure elucidation of small molecules, together with the public availability of database resources associated to metabolites and metabolic pathways, has enabled researchers to approach the metabolome of organisms in a high-throughput fashion. Other "omics" technologies have a longer history in high-throughput, such as next generation sequencing for genomics, RNA microarrays for transcriptomics, and mass spectrometry methods for proteomics. All of these together give researchers a unique opportunity to study and combine multi-omics aspects, forming the discipline of "Systems Biology" in order to study organisms at multiple scales simultaneously.


Introduction
Metabolomics has seen a rapid development of new technologies, methodologies, and data analysis procedures during the past decade.The development of fast gas-and liquid-chromatography devices coupled to sensitive mass-spectrometers, supplemented by the unprecedented precision of nuclear magnetic resonance for structure elucidation of small molecules, together with the public availability of database resources associated to metabolites and metabolic pathways, has enabled researchers to approach the metabolome of organisms in a high-throughput fashion.Other "omics" technologies have a longer history in high-throughput, such as next generation sequencing for genomics, RNA microarrays for transcriptomics, and mass spectrometry methods for proteomics.All of these together give researchers a unique opportunity to study and combine multi-omics aspects, forming the discipline of "Systems Biology" in order to study organisms at multiple scales simultaneously.Like all other "omics" technologies, metabolomics data acquisition is becoming more reliable and less costly, while at the same time throughput is increased.Modern time-of-flight (TOF) mass spectrometers are capable of acquiring full scan mass spectra at a rate of 500Hz from 50 to 750 m/z and with a mass accuracy <5 ppm with external calibration (Neumann & Böcker, 2010).At the opposite extreme of machinery, Fourier-transform ion-cyclotron-resonance (FTICR) mass spectrometers coupled to liquid chromatography for sample separation reach an unprecedented mass accuracy of <1 ppm m/z and very high mass resolution (Miura et al., 2010).These features are key requirements for successful and unique identification of metabolites.Coupled to chromatographic separation devices, these machines create datasets ranging in size from a few hundred megabytes to several gigabytes per run.While this is not a severe limitation for small scale experiments, it may pose a significant burden on projects that aim at studying the metabolome or specific metabolites of a large number of specimens and replicates, for example in medical research studies or in routine diagnostics applications tailored to the metabolome of a specific species (Wishart et al., 2009).
Thus, there is a need for sophisticated methods that can treat these datasets efficiently in terms of computational resources and which are able to extract, process, and compare the relevant information from these datasets.Many such methods have been published, however there is a high degree of fragmentation concerning the availability and accessibility of these methods, which makes it hard to integrate them into a lab's workflow.
The aim of this work is to discuss the necessary and desirable features of a software framework for metabolomics data preprocessing based on gas-chromatography (GC) and comprehensive two-dimensional gas-chromatography (GCxGC) coupled to single-dimension detectors (flame/photo ionization, FID/PID) or multi-dimension detectors (mass spectrometry, MS).We compare the features of publicly available Open Source frameworks that usually have a steep learning curve for end-users and bioinformaticians alike, owing to their inherent complexity.Many users will thus be appaled by the effort it takes to get used to a framework.Thus, the main audience of this work are bioinformaticians and users willing to invest some time in learning to use and/or program in these frameworks in order to set up a lab specific analytical platform.For a review of LC-MS based metabolomics data preprocessing consider (Castillo, Mattila, Miettinen, Orešič & Hyötyläinen, 2011).
Before we actually compare the capabilities of these different frameworks, we will first define a typical workflow for automatic data processing of metabolomics experiments and will discuss available methods within each of the workflow's steps.
We will concentrate on frameworks available under an Open Source license, thus allowing researchers to examine their actual implementation details.This distinguishes these frameworks from applications that are only provided on explicit request, under limited terms of use, or that are not published together with their source code (Lommen, 2009;Stein, 1999), which is still often the case in metabolomics and may hamper comparability and reuse of existing solutions.Additionally, all frameworks compared in this work are available for all major operating systems such as Microsoft Windows, Linux, and Apple Mac OSx as standalone applications or libraries.
Web-based methods are not compared within this work as they most often require a complex infrastructure to be set up and maintained.However, we will give a short overview of recent publications on this topic and provide short links to the parts of the metabolomics pipeline that we discuss in the following section.A survey of web-based methods is provided by Tohge & Fernie (2009).More recent web-based applications for metabolomics include the retention time alignment methods Warp2D (Ahmad et al., 2011) and ChromA (Hoffmann & Stoye, 2009), which are applicable to GC-MS or LC-MS data, and Chromaligner (Wang et al., 2010), which aligns GC and LC data with single-dimension detectors like FID.
In the Application section, we will exemplarily describe two pipelines for metabolomics analyses based on our own Open Source framework Maltcms: ChromA, which is applicable to GC-MS, and ChromA4D, which is applicable to data from comprehensive GCxGC-MS experiments.We show how to set up, configure and execute each pipeline using instructional datasets.These two workflows include the typical steps of raw-data preprocessing in metabolomics, including peak-finding and integration, peak-matching among multiple replicate groups and tentative identification using mass-spectral databases, as well as visualizations of raw and processed data.We will describe the individual steps of the workflows of the two application pipelines to give the reader a thorough understanding of the methods used by ChromA and ChromA4D.Finally, we discuss the current state of the presented Open Source frameworks and give an outlook into the future of software frameworks and data standards for metabolomics.

A typical workflow for a metabolomics experiment
Metabolomics can be defined as the study of the metabolic state of an organism or its response to direct or indirect perturbation.In order to find differences between two or more states, for example before treatment with a drug and after, and among one or multiple specimens, the actual hypothesis for the experiment needs to be defined.Based on this hypothesis, a design for the structure of the experiments and their subsequent analysis can be derived.This involves, among many necessary biological or medical considerations, the choice of sample extraction procedures and preparation methods, as well as the choice of the analytical methods used for downstream sample analysis.
Preprocessing of the data from those experiments begins after the samples have been acquired using the chosen analytical method, such as GC-MS or LC-MS.Owing to the increasing amount of data produced by high-throughput metabolomics experiments, with large sample numbers and high-accuracy/high-speed analytical devices, it is a key requirement that the resulting data is processed with very high level of automation.It is then that the following typical workflow is applied in some variation, as illustrated in Figure 1.

Data acquisition and conversion
The most common formats exported from GC-MS and LC-MS machines today are NetCDF (Rew & Davis, 1990), based on the specifications in the ASTM/AIA standard ANDI-MS (Matthews, 2000), mzXML (Oliver et al., 2004), mzData (Orchard et al., 2005) recently as the successor to the latter two, mzML (Deutsch, 2008;Martens et al., 2010).All of these formats include well-defined data structures for meta-information necessary to interpret data in the right context, such as detector type, chromatographic protocol, detector potential and other details about the separation and acquisition of the data.Furthermore, they explicitly model chromatograms and mass spectra, with varying degrees of detail.
NetCDF is the oldest and probably most widely used format today.It is routinely exported even by older machinery, which offers backwards compatibility to those.It is a general-purpose binary format, with a header that describes the structure of the data contained in the file, grouped into variables and indexed by dimensions.In recent years, efforts were made to establish open formats for data exchange based on a defined grammar in extensible markup language (XML) with extendable controlled vocabularies, to allow new technologies to be easily incorporated into the file format without breaking backwards compatibility.Additionally, XML formats are human readable which narrows the technology gap.mzXML was the first approach to establish such a format.It has been superseded by mzData and, more recently, mzML was designed as a super-set of both, incorporating extensibility through the use of an indexed controlled vocabulary.This allows mzML to be adapted to technologies like GCxGC-MS without having to change its definition, although its origins are in the proteomics domain.One drawback of XML-based formats is often claimed to be their considerably larger space requirements when compared to the supposedly more compact binary data representations.Recent advances in mzML approach this issue by compressing spectral data using gzip compression.
The data is continuously stored in a vendor-dependent native format during sample processing on a GC-MS machine.Along with the mass spectral information, like ion mass (or equivalents) and abundance, the acquisition time of each mass spectrum is recorded.Usually, the vendor software includes methods for data conversion into one of the aforementioned formats.However, especially when a high degreee of automation is desired, it may be beneficial to directly access the data in their native format.This avoids the need to run the vendor's proprietary software manually for every data conversion task.Both the ProteoWizard framework (Kessner et al., 2008) and the Trans Proteomic Pipeline (Deutsch et al., 2010) include multiple vendor-specific libraries for that use case.

Preprocessing
Raw mass specrometry data is usually represented in sparse formats, only recording those masses whose intensities exceed a user-defined threshold.This thresholding is usually applied within the vendor's proprietary software and may lead to artificial gaps within the data.Thus, the first step in preprocessing involves the binning of mass spectra over time into bins of defined size in the m/z dimension, followed by interpolation of missing values.After binning, the data is stored as a rectangular array of values, with the first dimension representing time, the second dimension representing the approximate bin mass values, and the third dimension representing the intensity corresponding to each measured ion.This process is also often described as resampling (Lange et al., 2007).
Depending on various instrumental parameters, the raw exported data may require additional processing.The most commonly reported methods for smoothing are the Savitzky-Golay filter (Savitzky & Golay, 1964), LOESS regression (Smith et al., 2006) and variants of local averaging, for example by a windowed moving average filter.These methods can also be applied to interpolate values where gaps are present in the original data.The top-hat filter (Bertsch et al., 2008;Lange et al., 2007) is used to remove a varying baseline from the signal.More refined methods use signal decomposition and reconstruction methods, such as Fourier transform and continuous wavelet transform (CWT) (Du et al., 2006;Fredriksson et al., 2009;Tautenhahn et al., 2008) in order to remove noise and baseline contributions from the signal and simultaneously find peaks.

Peak detection
Often the process of peak detection is decoupled from the actual preprocessing of the data.XCMS (Smith et al., 2006), for example, uses a Gaussian second derivative peak model with a fixed kernel width and signal-to-noise threshold to find peaks along the chromatographic domain of each ion bin.Other methods extend this approach to use a multi-scale continuous wavelet transform using such a kernel over various widths, tracking the response of the transformed signal in order to locate peak apex positions in scale-space before estimating the true peak widths based on the kernel scale with maximum response (Fredriksson et al., 2009;Tautenhahn et al., 2008).However, these methods usually allow only a small number of co-eluting peaks in different mass-bins, since they were initially designed to work with LC-MS data mainly, where only one parent ion and a limited number of accompanying adduct ions are expected.In GC-MS, electron-ionization creates rich fragmentation mass spectra, which pose additional challenges to deconvolution of co-eluting ions and subsequent association to peak groups.Even though its source code is not publicly available, the method used by AMDIS (Stein, 1999) has seen wide practical application and is well accepted as a reference by the metabolomics and analytical chemistry communities.

Alignment
The alignment problem in metabolomics and proteomics stems from the analytical methods used.These produce sampled sensor readings acquired over time in fixed or programmed intervals, usually called chromatograms.The sensor readings can be one-or multidimensional.In the first case, detectors like ultra violet and visible light absorbance detectors (UV/VIS) or flame ionization detectors (FID) measure the signal response as one-dimensional features, e.g. as the absorbance spectrum or electrical potential, respectively.Multi-dimensional detectors like mass spectrometers record a large number of features simultaneously, e.g.mass and ion count.The task is then to find corresponding and non-corresponding features between different sample acquisitions.This correspondence problem is a term used by Åberg et al. (2009) which describes the actual purpose of alignment, namely to find true correspondences between related analytical signals over a number of sample acquisitions.For GC-MS-and LC-MS-based data, a number of different methods have been developed, some of which are described in more detail by Castillo, Gopalacharyulu, Yetukuri &Orešič (2011) andÅberg et al. (2009).Here, we will concentrate on those methods that have been reported to be applicable to GC-MS.In principle, alignment algorithms can be classified into two main categories: peak-and signal-based methods.Methods of the first type start with a defined set of peaks, which are present in most or all samples that are to be aligned before determining the best correspondences of the peaks between samples in order to then derive a time correction function.Krebs et al. (2006) locate landmark peaks in the TIC and then select pairs of those peaks with a high correlation between their mass spectra in order to fit an interpolating spline between a reference chromatogram and the to-be-aligned one.The method of Robinson et al. (2007) is inspired by multiple sequence alignment algorithms and uses dynamic programming to progressively align peak lists without requiring an explicit reference chromatogram.Other methods, like that of Chae et al. (2008) perform piecewise, block-oriented matching of peaks, either on the TIC, on selected masses, or on the complete mass spectra.Time correction is applied after the peak assignments between the reference chromatogram and the others have been calculated.Signal-based methods include recent variants of correlation optimized warping (Smilde & Horvatovich, 2008), parametric time warping (Christin et al., 2010) and dynamic time warping (Christin et al., 2010;Clifford et al., 2009;Hoffmann & Stoye, 2009;Prince & Marcotte, 2006) and usually consider the complete chromatogram for comparison.However, attempts are made to reduce the computational burden associated with a complete pairwise comparison of mass spectra by partitioning the chromatograms into similar regions (Hoffmann & Stoye, 2009), or by selecting a representative subset of mass traces (Christin et al., 2010).Another distinction in alignment algorithms is the requirement of an explicit reference for alignment.Some methods apply clustering techniques to select one chromatogram that is most similar to all others (Hoffmann & Stoye, 2009;Smilde & Horvatovich, 2008), while other methods choose such a reference based on the number of features contained in a chromatogram (Lange et al., 2007) or by manual user choice (Chae et al., 2008;Clifford et al., 2009).For high-throughput applications, alignments should be fast to calculate and reference selection should be automatic.Thus, a sampling method for time correction has recently been reported by Pluskal et al. (2010) for LC-MS.A comparison of these methods is given in the same publication.

Statistical evaluation
After peaks have been located and integrated for all samples, and their correspondence has been established, peak report tables can be generated, containing peak information for each sample and peak, with associated corrected retention times and peak areas.Additionally, peaks may have been putatively identified by searching against a database, such as the GMD (Hummel et al., 2007) or the NIST mass-spectral database (Babushok et al., 2007).
These peak tables can then be analyzed with further methods, in order to detect e.g.systematic differences between different sample groups.Prior to such an analysis, the peak areas need to be normalized.This is usually done by using a spiked-in compound which is not expected to occur naturally as a reference.The normalization compound is supposed to have the same concentration in all samples.The compound's peak area can then be used to normalize all peak areas of a sample with respect to it (Doebbe et al., 2010).
Different experimental designs allow to analyze correlations of metabolite levels for the same subjects under different conditions (paired), or within and between groups of subjects.For simple paired settings, multiple t-tests with corrections for multiple testing can be applied (Berk et al., 2011), while for comparisons between groups of subjects, Fisher's F-Statistic (Pierce et al., 2006) and various analysis of variance (ANOVA), principal component analysis (PCA) and partial least squares (PLS) methods are applied (Kastenmüller et al., 2011;Wiklund et al., 2008;Xia et al., 2011).

Evaluation of hypothesis
Finally, after peak areas have been normalized and differences have been found between sample groups, the actual results need to be put into context and be interpreted in their biological context.This task is usually not handled by the frameworks described in this chapter.Many web-based analysis tools allow to put the data into a larger context, by providing name-or id-based mapping of the experimentally determined metabolite concentrations onto biochemical pathways like MetaboAnalyst (Xia & Wishart, 2011), MetabolomeExpress (Carroll et al., 2010), or MeltDB (Neuweger et al., 2008).The latter allows association of the metabolomics data with other results for the same subjects under study or with results from other "omics" experiments on the same target subjects, but this is beyond the scope of the frameworks presented herein.

Frameworks for GC-MS analysis
A number of Open Source frameworks have been developed for LC-MS based proteomics frameworks like OpenMS (Bertsch et al., 2008), ProteoWizard (Kessner et al., 2008), and most notably the TransProteomicPipeline (Deutsch et al., 2010).Even though many of the steps required for proteomics apply similarily to metabolomics applications, there are still some essential differences due to the different analytical setups and technologies (e.g.matrix assisted laser desorption ionization mass spectrometry, MALDI-MS) used in the two fields.XCMS (Smith et al., 2006) was among the first frameworks to offer support for data preprocessing in LC-MS based metabolomics.Later, MZmine2 (Pluskal et al., 2010) offered an alternative with a user-friendly interface and easy extendability.Lately, Scheltema et al. (2011) published their PeakML format and mzMatch framework also for LC-MS applications.As of now, there seem to be only a few frameworks available for GC-MS based metabolomics that offer similar methods, namely PyMS (Callaghan et al., 2010;Isaac et al., 2009) and Maltcms/ChromA (Hoffmann & Stoye, 2009;Maltcms, 2011) .These will be presented in more detail in this section.A compact overview of the Open Source frameworks discussed herein is given in Table 1.A detailed feature comparison can be found in Table 2.

XCMS
XCMS (Smith et al., 2006) is a very mature framework and has seen constant development during the last five years.It is mainly designed for LC-MS applications, however its binning, peak finding and alignment are also applicable to GC-MS data.XCMS is implemented in the GNU R programming language, the de-facto standard for Open Source statistics.Since GNU R is an interpreted scripting language, it is easy to write custom scripts that realize additional functionality of the typical GC-MS workflow described above.XCMS is part of the Bioconductor package collection, which offers many computational methods for various "omics" technologies.Further statistical methods are available from GNU R.
XCMS supports input in NetCDF, mzXML, mzData and, more recently, mzML format.This allows XCMS to be used with virtually any chromatography-mass spectrometry data, since vendor software supports conversion to at least one of those formats.XCMS uses the xcmsRaw object as its primary tabular data structure for each binned data file.The xcmsSet object is then used to represent peaks and peak groups and is used by its peak alignment and diffreport features.
The peak finding methods in XCMS are quite different from each other.For data with normal or low mass resolution and accuracy, the matched filter peak finder (Smith et al., 2006) is usually sensitive enough.It uses a Gaussian peak template function with user defined width and signal-to-noise critera to locate peaks on individual binned extracted ion current (EIC) traces over the complete time range of the binned chromatogram.The other method, CentWave (Tautenhahn et al., 2008) is based on a continuous wavelet transform on areas of interest within the raw data matrix.Both peak finding methods report peak boundaries and integrated areas for raw data and for the data reconstructed from the peak finder's signal response values.
Initially designed for LC-MS, XCMS does not have a method to group co-eluting peaks into peak groups, as is a requirement in GC-MS methods using electron ionization.However, CAMERA (Tautenhahn et al., 2007) shows how XCMS can be used as a basis in order to create a derived application, in this case for ion annotation between samples.
Peak alignment in XCMS is performed using local LOESS regression between peak groups with very similar m/z and retention time behaviour and good support within each sample group.This allows a simultaneous alignment and retention time correction of all peaks.The other available method is based on the Obi-Warp dynamic time warping (Prince & Marcotte, 2006) algorithm and is capable of correcting large non-linear retention time distortions.It uses the peak set with the highest number of features as alignment reference, which is comparable to the approach used by Lange et al. (2007).However, it is much more computationally demanding then the LOESS-based alignment.
XCMS's diffreport generates a summary report of significant analyte differences between two sample sets.It uses Welch's two-sample t-statistic to calculate p-values for each analyte group.ANOVA may be used for more than two sample sets.
A number of different visualizations are also available, both for raw and processed data.These include TIC plots, EIC plots, analyte group plots for grouped features, and chromatogram (rt, m/z, intensity) surface plots.
XCMS can use GNU R's Rmpi infrastructure to execute arbitary function calls, such as profile generation and peak finding, in parallel on a local cluster of computers.

PyMS
PyMS (Callaghan et al., 2010;Isaac et al., 2009) is a programming framework for GC-MS metabolomics based on the Python programming language.It can therefore use a large number of scientific libraries which are accessible via the SciPy and NumPy packages (SciPy, 2011).Since Python is a scripting language, it allows to do rapid prototyping, comparable to GNU R.However, Python's syntax may be more familiar for programmers with a background in object-oriented programming languages.
The downloadable version of PyMS currently only supports NetCDF among the more recent open data exchange formats.Nonetheless, it is the only framework in this comparison with support for the JCAMP GC-MS file format.
PyMS provides dedicated data structures for chromatograms, allowing efficient access to EICs, mass spectra, and peak data.
In order to find peaks, PyMS also builds a rectangular profile matrix with the dimensions time, m/z and intensity.Through the use of slightly shifted binning boundaries, they reduce the chance of false assignments of ion signals to neighboring bins, when binning is performed with unit precision (bin width of 1 m/z).PyMS offers the moving average and the Savitzky-Golay (Savitzky & Golay, 1964) filters for signal smoothing of EICs within the profile matrix.Baseline correction can be performed by the top-hat filter (Lange et al., 2007).The actual peak finding is based on the method described by Biller & Biemann (1974) and involves the matching of local peak maxima co-eluting within a defined window.Peaks are integrated for all co-eluting masses, starting from a peak apex to both sides and ending if the increase in area falls below a given threshold.
Peak alignment in PyMS is realized by the method introduced by Robinson et al. (2007).It is related to progressive multiple sequence alignment methods and is based on a generic dynamic programming algorithm for peak lists.It proceeds by first aligning peak lists within sample groups, before aligning the aligned peak lists of different groups, until all groups have been aligned.
Visualizations of chromatogram TICs, EICs, peaks and mass spectra are available and are displayed to the user in an interactive plot panel.
For high-throughput applications, PyMS can be used together with MPI to parallelize tasks within a local cluster of computers.

Maltcms
The framework Maltcms allows to set up and configure individual processing components for various types of computational analyses of metabolomics data.The framework is implemented in JAVA and is modular using the service provider pattern for maximal decoupling of interface and implementation, so that it can be extended in functionality at runtime.
Maltcms can read data from files in NetCDF, mzXML, mzData or mzML format.It uses a pipeline paradigm to model the typical preprocessing workflow in metabolomics, where each processing step can define dependencies on previous steps.This allows automatic pipeline validation and ensures that a user can not define an invalid pipeline.The workflow itself is serialized to XML format, keeping track of all resources created during pipeline execution.Using a custom post-processor, users can define which results of the pipeline should be archived.
Maltcms uses a generalization of the ANDI-MS data schema internally and a data provider interface with corresponding implementations to perform the mapping from any proprietary data format to an internal data object model.This allows efficient access to individual mass spectra and other data available in the raw-data files.Additionally, developers need no special knowledge of any supported file format, since all data can be accessed generically.Results from previous processing steps are referenced in the data model to allow both shadowing of data, e.g.creating a processing result variable with the same name as an already existing variable, and aggregation of processing results.Thus, all previous processing results are transparently accessible for downstream elements of a processing pipeline, unless they have been shadowed.Primary storage of processing results is performed on a per-chromatogram basis in the binary NetCDF file format.Since metabolomics experiments create large amounts of data, a focus is put on efficient data structures, data access, and scalability of the framework.
Embedding Maltcms in existing workflows or interfacing with other software is also possible, as alignments, peak-lists and other feature data can be exported as comma separated value files or in specific xml-based formats, which are well-defined by custom schemas.
To exploit the potential of modern multi-core CPUs and distributed computing networks, Maltcms supports multi-threaded execution on a local machine or within a grid of connected computers using an OpenGrid infrastructure (e.g.Oracle Grid Engine or Globus Toolkit (Foster, 2005)) or a manually connected network of machines via remote method invocation (RMI).
The framework is accompanied by many libraries for different purposes, such as the JFreeChart library for 2D-plotting or, for BLAS compatible linear algebra, math and statistics implementations, the Colt and commons-math libraries.Building upon the base library Cross, which defines the commonly available interfaces and default implementations, Maltcms provides the domain dependent data structures and specializations for processing of chromatographic data.

Name
Version

ChromA
ChromA is a configuration of Maltcms that includes preprocessing, in the form of mass binning, time-scale alignment and annotation of signal peaks found within the data, as well as visualizations of unaligned and aligned data from GC-MS and LC-MS experiments.
The user may supply mandatory alignment anchors as CSV files to the pipeline and a database location for tentative metabolite identification.Further downstream processing can be performed either on the retention time-corrected chromatograms in NetCDF format, or on the corresponding peak tables in either CSV format or XML format.
Peaks can either be imported from other tools, by providing them in CSV format to ChromA, giving at least the scan index of each peak in a file per row.Alternatively, ChromA has a fast peak finder that locates peaks based on derivatives of the smoothed and baseline-corrected TIC, using a moving average filter followed by top-hat filter baseline-substraction, with a predefined minimum peak-width.Peak alignment is based on a star-wise or tree-based application of an enhanced variant of pairwise dynamic time warping (DTW) (Hoffmann & Stoye, 2009).To reduce both runtime and space requirements, conserved signals throughout the data are identified, constraining the search space of DTW to a precomputed closed polygon.The alignment anchors can be augmented or overwritten by user-defined anchors, such as previously identified compounds, characteristic mass or MS/MS identifications.
Then, the candidates are paired by means of a bidirectional best-hits (BBH) criterion, which can compare different aspects of the candidates for similarity.Paired anchors are extended to k-cliques with configurable k, which help to determine the conservation or absence of signals across measurements, especially with respect to replicate groups.Tentative identification of peaks against a database using their mass spectra is possible using the MetaboliteDB module.This module provides access to mass-spectral databases in msp-compatible format, for example the Golm Metabolite Database or the NIST EI-MS database.
ChromA visualizes alignment results including paired anchors in birds-eye view or as a simultaneous overlay plot of the TIC.Additionally, absolute and relative differential charts are provided, which allow easy spotting of quantitative differences.
Peak tables are exported in CSV format, including peak apex positions, area under curve, peak intensity and possibly tentative database identifications.Additionally, information about the matched and aligned peak groups is saved in CSV format.

Frameworks for GCxGC-MS analysis
The automatic and routine analysis of comprehensive GCxGC-MS data is yet to be established.GCxGC-MS couples a second chromatographic column to the first one, thereby achieving a much higher peak capacity and thus a better separation of closely co-eluting analytes (Castillo, Mattila, Miettinen, Orešič & Hyötyläinen, 2011).Usually, for a one-hour run, the raw data file size exceeds a few Gigabytes.Quite a number of algorithms have been published on alignment of peaks in such four-dimensional (first column retention time, second column retention time, mass, and intensity values) data (Kim et al., 2011;Oh et al., 2008;Pierce et al., 2005;Vial et al., 2009;Zhang, 2010), however only a few methods are available for a more complete typical preprocessing workflow.A compact overview of the available frameworks, their licenses and programming languages is given in Table 3. Table 4 gives a more detailed feature matrix of these frameworks.The remainder of this section gives a concise overview of the frameworks Guineu (Castillo, Mattila, Miettinen, Orešič &Hyötyläinen, 2011) andChromA4D (Maltcms, 2011).

Name
Version Supported methods Software license Programming language Guineu 0.8.2 GCxGC-MS (LC-MS) GNU GPL v2 JAVA 6 Maltcms/ChromA4D 1.1 GCxGC-MS GNU L-GPL v3 JAVA 6 Table 3. Feature comparison of Open Source software frameworks for GCxGC-MS based metabolomics

Guineu
Guineu is a recently published graphical user interface and application for the comparative analysis of GCxGC-MS data (Castillo, Mattila, Miettinen, Orešič & Hyötyläinen, 2011).It currently reads LECO ChromaTOF software's peak list output after smoothing, baseline correction, peak finding, deconvolution, database search and retention index (RI) calculation have been performed within ChromaTOF.
The peak lists are aligned pairwise using the score alignment algorithm, which requires user-defined retention time windows for both separation dimensions.Additionally, the one-dimensional retention index (RI) of each peak is used within the score calculation.a threshold for mass spectral similarity is needed in order to create putative peak groups.Additional peak lists are added incrementally to an already aligned path, based on the individual peaks' score against those peaks that are already contained within the path.
Guineu provides different filters to remove peaks by name, group occurrence count, or other features from the ChromaTOF peak table.In order to identify compound classes, the Golm metabolite database (GMD) substructure search is used.Peak areas can be extracted from ChromaTOF using the TIC, or using extracted, informative or unique masses.Peak area normalization is available relative to multiple user-defined standard compounds.
After peak list processing, Guineu produces an output table containing information for all aligned peaks, containing information on the original analyte annotation as given by ChromaTOF, peak areas, average retention times in both dimensions together with the average RI and further chemical information on the functional group and substructure prediction as given by the GMD.It is also possible to link the peak data to KEGG and Pubchem via the CAS annotation, if it is available for the reported analyte.
For statistical analysis of the peak data, Guineu provides fold change-and t-tests, principal component analysis (PCA), analysis of variance (ANOVA) and other methods.
Guineu's statistical analysis methods provide different plots of the data sets, e.g. for showing the principal components of variation within the data sets after analysis with PCA.

ChromA4D
For the comparison of comprehensive two-dimensional gas chromatography-mass spectrometry (GCxGC-MS) data, ChromA4D accepts NetCDF files as input.Additionally, the user needs to provide the total runtime on the second orthogonal column (modulation time) to calculate the second retention dimension information from the raw data files.For tentative metabolite identification, the location of a database can be given by the user.ChromA4D reports the located peaks, their respective integrated TIC areas, their best matching corresponding peaks in other chromatograms, as well as a tentative identification for each peak.Furthermore, all peaks are exported together with their mass spectra to MSP format, which allows for downstream processing and re-analysis with AMDIS and other tools.The exported MSP files may be used to define a custom database of reference spectra for subsequent analyses.
Peak areas are found by a modified seeded region growing algorithm.All local maxima of the TIC representation that exceed a threshold are selected as initial seeds.Then, the peak area is determined by using the distance of the seed mass spectrum to all neighbor mass spectra as a measure of the peak's coherence.The area is extended until the distance exceeds a given threshold.No information about the expected peak shape is needed.The peak integration is based on the sum of TICs of the peak area.An identification of the area's average or apex mass spectrum or the seed mass spectrum is again possible using the MetaboliteDB module.
To represent the similarities and differences between different chromatograms, bidirectional best hits are used to find co-occurring peaks.These are located by using a distance that exponentially penalizes differences in the first and second retention times of the peaks to be compared.To avoid a full computation of all pairs of peaks, only those peaks within a defined window of retention times based on the standard deviation of the exponential time penalty function are evaluated.
ChromA4D's visualizations represent aligned chromatograms as color overlay images, similar to those used in differential proteomics.This allows a direct visual comparison of signals present in one sample, but not present in another sample.
ChromA4D creates peak report tables in CSV format, which include peak apex positions in both chromatographic dimensions, area under curve, peak intensity and possibly tentative database identifications.Additionally, information about the matched and aligned peak groups is saved in CSV format.

Application examples
The following examples for GC-MS and GCxGC-MS are based on the Maltcms framework, using the ChromA and ChromA4D configurations described in the previous sections.In order to run them, the recent version of Maltcms needs to be downloaded and unzipped to a local folder on a computer.Additionally, Maltcms requires a JAVA runtime environment version 6 or newer to be installed.If these requirements are met, one needs to start a command prompt and change to the folder containing the unzipped Maltcms.

An example workflow for GC-MS
The experiment used to illustrate an example workflow for one-dimensional GC-MS consists of two samples of standard compounds, which contain mainly sugars, amino acids, other organic acids and nucleosides, measured after manual (MD) and after automatic derivatization (AD) with the derivatization protocol and substances given below.Group AD consists of a sample of n-alkanes standard and two replicates of mix1, namely mix1-1 and mix1-2.We will show how ChromA can be used to find and integrate peaks, as well as compare and align the peaks between the samples, and finally how the alignment results can be used for quality control.

Acquisition and data processing
The samples were acquired on an Agilent GC 7890N with MSD 5975C triple axis detector.An Agilent HP5ms column with a length of 30 m, a diameter of 0.25 mm, and a film thickness of 0.25 μm (Agilent, Santa Clara CA, USA) was used for the gas-chromatographic separation, followed by a deactivated restriction capillary with 50 cm length and a diameter of 0.18 mm.Per sample, 1 μL was injected onto the column in pulsed splitless mode (30 psi for 2 min).The flow rate was set to 1.5 mL/min of Helium.The linear temperature ramp started at 50 • C for 2 min until it reached its maximum of 325 • C at a rate of 10 • C/min.The raw data were exported to NetCDF format using the Agilent ChemStation software v.B.04.01 (Agilent, Santa Clara CA, USA) with default parameters and without additional preprocessing applied.
A sample containing n-alkanes was measured as an external standard for manual (MD) and automatic derivatization (AD) in order to be able to later determine retention indices for  the other samples.The acquired data were exported to ANDI-MS (NetCDF) format before ChromA was applied.The default ChromA pipeline chroma.propertieswas run from the unzipped Maltcms directory with the following command (issued on a single line of input): > java -Xmx1G -jar maltcms.jar-i ../data/ -o ../output/ -f * .CDF \ -c cfg/chroma.properties-i points to the directory containing the input data, -o points to the directory where output should be placed, -f can be a comma separated list of filenames or, as in this case, a wildcard expression, matching all files in the input directory having a file name ending with .CDF.
The final argument indicated by -c is the path to the configuration file used for definition of the pipeline and its commands.An overlay of the raw TICs of the samples is depicted in Figure 2(a).The default ChromA pipeline configuration creates a profile matrix with nominal mass bin width.Then, the TIC peaks are located separately within each sample data file and are integrated (Figure 2(b)).The peak apex mass spectra are then used in the next step in order to build a multiple peak alignment between all peaks of all samples by finding large cliques, or clusters of peaks exhibiting similar retention time behaviour and having highly similar mass spectra.This coarse alignment could already be used to calculate a polynomial fit, correcting retention time shift for all peaks.However, the ChromA pipeline uses the peak clusters in order to constrain a dynamic time warping (DTW) alignment in the next step, which is calculated between all pairs of samples.The resulting distances are used to determine the reference sample with the lowest sum of distances to all remaining samples.Those are then aligned to the reference using the warp map obtained from the pairwise DTW calculations.The pairwise DTW distances can easily be used for a hierarchical cluster analysis.Similar samples should be grouped into the same cluster, while dissimilar samples should be grouped into different clusters.Figure 3 shows the results of applying a complete linkage clustering algorithm provided by GNU R to the pairwise distance matrix.It is clearly visible that the samples are grouped correctly, without incorporation of any external group assignment.Thus, this method can be used for quality control of multiple sample acquisitions, when the clustering results are compared against a pre-defined number of sample groups.

An example workflow for GCxGC-MS
The instructional samples presented in this section were preprocessed according to the protocol given by Doebbe et al. (2010).The description of the protocol has been adapted from that reference where necessary.

Acquisition and data processing
The sample acquisition was performed on a LECO Pegasus 4D TOF-MS (LECO, St. Joseph, MI, USA).The Pegasus 4D system was equipped with an Agilent 6890 gas chromatograph (Agilent, Santa Clara, CA, USA).The inlet temperature was set to 275 • C.An Rtx-5ms (Restek, Bellefonte, PA, USA) capillary column was used with a length of 30 m, 0.25 mm diameter and 0.25 μm film thickness as the primary column.The secondary column was a BPX-50 (SGE, Ringwood, Victoria, Australia) capillary column with a length of 2 m, a diameter of 0.1 mm and 0.1 μm film thickness.The temperature program of the primary oven was set to the following conditions: 70 • C for 2 min, 4 • C/min to 180 • C, 2 • C/min to 230 • C, 4 • C/min to 325 • C hold 3 min.This program resulted in a total runtime of about 70 min for each sample.The secondary oven was programmed with an offset of 15 • C to the primary oven temperature.The thermal modulator was set 30 • C relative to the primary oven and to a modulation time of 5 seconds with a hot pulse time of 0.4 seconds.The mass spectrometer ion source temperature was set to 200 • C and the ionization was performed at -70eV.The detector voltage was set to 1600V and the stored mass range was 50-750 m/z with an acquisition rate of 200 spectra/second.The raw acquired samples in LECO's proprietary ELU format were exported to NetCDF format using the LECO ChromaTOF software v.4.22 (LECO, St. Joseph, MI, USA).Initial attempts to export the full, raw data failed with a crash beyond a NetCDF file size of 4GBytes.Thus, we resampled the data with ChromaTOF to 100 Hz (resampling factor 2) and exported with automatic signal smoothing and baseline offset correction value of 1 which resulted in file sizes around 3GBytes per sample.The samples presented in this section are named "Standard-Mix1-1" and "Standard-Mix1-2" and were measured on different days (Nov. 29th, 2008 andDec. 12th, 2008).
The default ChromA4D pipeline for peak finding was called from within the unzipped Maltcms directory (issued on a single line of input): > java -Xmx2G -jar maltcms.jar-i ../data/ -o ../output/ \ -f * .cdf-c cfg/4Dpeakfinding.properties The pipeline first preprocesses the data by applying a median filter followed by a top hat filter in order to remove high-and low-frequency noise contributions (Figures 4(a) and 4(b)).ChromA4D then uses a variant of seeded region growing in order to extend peak seeds, which are found as local maxima of the 2D-TIC.These initial seeds are then extended until the mass spectral similarity of the seed and the next evaluated candidate drops below a user-defined threshold, or until the peak area reaches its maximum, pre-defined size (Figure 5(a)).After peak area integration, the pipeline clusters peaks between samples based on their mass spectral similarity and retention time behaviour in both dimensions to form peak cliques (not shown) as multiple peak alignments, which are then exported into CSV format for further downstream processing.Another possible application shown in Figure 5(b) is the visualization of pairwise GCxGC-MS alignments using DTW on the vertical 2D-TIC slices, which can be useful for qualitative comparisons.

Summary and outlook
The present state of Open Source frameworks for metabolomics is very diverse.A number of tools have seen steady development and improvement over the last years, such as XCMS, MZmine, and PyMS, while others are still being developed, such as mzMatch, Guineu, and Maltcms.There is currently no framework available that covers every aspect of metabolomics data preprocessing.Most of the frameworks concentrate on one or a few analytical technologies with the largest distinction being between GC-MS and LC-MS.GCxGC-MS raw data processing is currently only handled by Maltcms' ChromA4D pipeline, while Guineu processes peak lists exported from LECO's ChromaTOF software and offers statistical methods for sample comparison together with a user-friendly graphical interface.
We showed two instructive examples on setting up and running the basic processing pipelines ChromA and ChromA4D for GC-MS and GCxGC-MS raw data.The general structure of these pipelines would be slightly different for each of the Open Source frameworks presented in this chapter, however, the basic concepts behind the processing steps are the same for all tools.Since metabolomics is an evolving field of research, no framework captures all possible use-cases, but it will be interesting to see which frameworks will be flexible and extendable enough to be adapted to new requirements in the near future.
In order to combine experiments from multiple "omics" experiments, another level of abstraction on top of local or web-service based tools for data processing, fusion, and integration of metabolomics experiments is a necessary future requirement.Generic workflow systems like Taverna (Hull et al., 2006) or Conveyor (Linke et al., 2011) offer integration of such resources, augmented with graphical editors for point-and-click user interaction.However, due to their generic nature these systems are far away from being as user-friendly as applications designed for a specific data analysis task and require some expert knowledge when assembling task-specific processing graphs.
One point that requires further attention is the definition and controlled evolution of peak data formats for metabolomics, along with other formats for easier exchange of secondary data between applications and frameworks.A first step in this direction has been taken by Scheltema et al. (2011) by defining the PeakML format.However, it is important that such formats are curated and evolved, possibly by a larger non-profit organization like the HUPO within its proteomics standards initiative HUPO PSI.Primary data is already acessible in a variety of different, defined formats, the most recent addition being mzML (Martens et al., 2010) which is curated by the PSI.Such standardization attempts can however only be successful and gain the required momentum if also the manufacturers of analytical machinery support the formats with their proprietary software within a short time frame after the specification and see a benefit in offering such functionality due to the expressed demand of scientists working in the field as in case of NetCDF, mzData, or mzML.
Fig. 1.A typical workflow for a metabolomics experiment.Steps shown in orange (solid border) are usually handled within the bioinformatics domain, while the steps shown in green (dashed border) often involve co-work with scientists from other disciplines.
(a) Overlay of unaligned data sets, extracted from middle section within a time range of 1100 to 1700 seconds.(b)Overlay with highlighted peak areas (without n-alkanes) after peak finding and integration.Zoomed in to provide more detail.

Fig. 2 .
Fig. 2. TIC overlay plots of the raw GC-MS data sets.
Fig. 3. Clustering of GC-MS samples based on pairwise DTW similarities transformed to distances.The samples are clearly separated into two clusters, one containing the n-alkane standard samples, the other one containing the mix1 samples.
Fig. 4. Visualizations of Standard-Mix1-1 before and after signal filtering with the ChromA4D processing pipeline.

Table 4
Peak detection MAX-SRG: TIC local maxima, seeded region growing based on ms similarity.
. Feature comparison of Open Source software frameworks for preprocessing of GCxGC-MS based metabolomics data.Key to abbreviations: Data formats A: NetCDF, G: ChromaTOF peak lists, H: CSV peak lists.Signal preprocessing MA: moving average, MM: moving median, TH: top-hat filter, CV: coefficient of variation threshold.ANOVA: analysis of variance, FT: F-test, between group vs. within group variance.