BACKGROUND
the "genomic era" is being characterized by the massive determination of the molecular components of living systems.
genome sequencing projects are yielding complete genome sequences for hundreds of organisms  <cit> .
although not as massive as envisaged, structural genomics projects are speeding up the rate of protein structure determinations  <cit> .
many other initiatives also address large-scale repertories of molecular components and their relationships, seeking to decipher the "transcriptome"  <cit> , the "interactome"  <cit> , etc.
these massive data contain much information about living organisms.
studied as a whole, from a systemic perspective, they provide global pictures of different aspects of biology, which can help to answer very basic questions about how life evolved and how organisms do what they do with their "molecular toolkits".
for example, the repertory of protein sequences  and their evolutionary relationships  can be represented in a "sequence space"  <cit> .
studying this global landscape of protein sequences as a whole has produced important information on the estimated total number of protein families , functional groups, evolutionary relationships, etc.
<cit> .
similarly, the wealth of information on protein three-dimensional structures has allowed comprehensive hierarchical classifications such as scop  <cit>  and cath  <cit>  to be generated.
these classifications can subsequently be mined to perform global studies on the estimated total number of folds , the evolution of protein structures, the relationship between sequence and functional spaces, etc.
<cit> .
metabolism was one of the last aspects to be studied from an "-omic" perspective  <cit> .
the "metabolome" can be defined as the total complement of small-molecule metabolites within a cell, while the "reactome" is the set of chemical reactions that transforms these metabolites.
in many cases, "metabolone" refers not only to the complement of small molecules but to quantitative measurements of their concentrations as well.
an important difference from the "proteome" and other "-omics" is that most metabolic information has not been obtained by a high-throughput process but by the accumulation of knowledge over many years of detailed biochemical characterization of compounds, reactions and enzymes, knowledge that has been stored in databases such as kegg  <cit>  and brenda  <cit> .
only recently, high-througput techniques have been applied to the massive determination of metabolomes  <cit>  and, even more recently, reactomes  <cit> .
this has advantages and drawbacks.
on the one hand, metabolic data are more reliable and do not suffer from the high degree of error common to high-throughput-omics data, such as the "interactome".
on the other, biochemical knowledge can be biased to certain pathways and processes that have been more intensely studied than others.
one of the main problems in performing global studies of the reactome or the metabolome , equivalent to the studies discussed above for the proteome, is that there is no good way of representing and quantifying metabolic information.
it is quite straightforward to quantify the similarity between the sequences or three-dimensional structures of two proteins.
it is even possible to discretize these continuous scales of similarity to define intervals of evolutionary and functional relationship  <cit> , expected structural similarity for a given sequence similarity  <cit> , etc.
the possibility of representing protein sequences and structures, as well as quantifying their corresponding similarities, is crucial for generating global classifications and subsequent global studies such as those cited above.
however, quantifying similarities among chemical compounds is not as straightforward as it is for protein sequences and structures.
there are many different ways of quantifying similarities between small molecules  <cit> , which can be used to study the global characteristics of the "metabolic space"  <cit>  in a similar way to sequence and structure spaces.
the main problem arises with chemical reactions : there is no easy way of quantifying the "similarity between two chemical reactions".
the few existing approaches  <cit>  are not intended to generate detailed global classifications of biochemical reactions in an unsupervised way or to perform global studies of metabolomes and reactomes aimed at quantifying their chemical diversities.
indeed, the most widely-used classification of chemical reactions  schema) is neither quantitative nor based on existing data.
rather, ec is based on an a priori imposed schema that organizes all existing enzymatic activities into a hierarchical classification of four levels with six main classes in the first level.
it was designed mainly for nomenclature purposes and lacks a clear reflection at the molecular or evolutionary levels.
although very useful for many tasks, this classification was not designed for quantitative global classifications or comparisons of biochemical activities.
in this work, we propose a very simple way of representing chemical reactions in a vectorial form, which allows them to be compared quantitatively.
this representation is based on an equivalent vectorial representation of chemical compounds, which imposes no a priori definition of functional groups or other "important" parts of molecules.
it allows any chemical species to be represented in the same vectorial space.
with this description, chemical compounds, reactions and pathways can be concomitantly represented in the same vectorial space.
we demonstrate that the vectorial representation of chemical compounds, in spite of its simplicity, efficiently captures their molecular properties and can be used for predictive purposes.
we used this method to generate an unsupervised global clustering of metabolic activities.
this clustering makes more biological sense than equivalent groupings defined by the ec hierarchy, an improvement that is evaluated statistically using different criteria.
moreover, it results in an optimal set of  <dig> groups of reactions for which we look for their biological meaning.
methods
vectorial representation of chemical compounds
we looked for a metabolic representation that:  allows all chemical species to be represented in the same space regardless of their composition, molecular size, etc.
;  imposes no a priori definition of important chemical features , since a representation based on a pre-defined set of functional groups would necessarily restrict the results obtained to those that can be given in terms of those groups; and  allows chemical reactions and pathways to be represented concomitantly in the same space .
we represent a chemical structure with a vector in which each component represents a triplet of neighboring atoms, and its value indicates the frequency of that triplet in the molecule .
using pairs of atoms instead of triplets would lead to shorter vectors , but would not capture many aspects of the molecule's functionality; the chemical properties of an atom depend mostly on its neighbors up to distance  <dig>  for example, the different properties of the carbonyl group c = o depend on the neighbors of the c: oh-cr = o , h-cr = o , r-cr = o , etc.
on the other hand, groups of  <dig>   <dig>  ... neighboring atoms would produce vectors with too many components .
we previously used a similar representation to train a machine learning system for predicting the biodegradability features of chemical compounds  <cit>  .
this representation is not perfectly univocal, since in some cases different molecules can have the same vector .
nevertheless, it provides a good balance between information content and utility for the goal of this work.
the representation contains more structural information than is first apparent.
for example the bond types, although not explicitly considered in the representation, are implicitly taken into account to some extent: as we go from hydrocarbons with few double and triple bonds  to hydrocarbons with many multiple bonds , we gain ccc triplets at the expense of hch and hcc.
in the "results" section, we demonstrate that this representation captures important physico-chemical properties of compounds and can be used for predictive purposes.
vectorial representation of biochemical reactions
the vector for a reaction that transforms one chemical compound into another is just the difference between the vectors of the two compounds  .
the stoichiometry is taken into account by multiplying the compound vectors by the corresponding stoichiometric coefficients: e.g.
for a reaction a -> 2b the reaction vector would be va->b = 2·vb - va.
since this vectorial representation of reactions is based on that of the compounds, it inheres many of its characteristics:  any chemical reaction can be represented in the same space regardless of its characteristics ;  there is no a priori definition of important reaction features.
moreover, the components of the reaction vector represent the triplets that are changed in the chemical reaction so they neatly encode the characteristics of the transformation .
another intuitive feature is that the vector of a given reaction is the opposite of the vector of the inverse reaction  .
in this representation, all chemical reactions have to be simplified as transformations of one single compound into another.
for a reaction with more than one reactant and/or product, a "main" pair substrate ->product has to be chosen and information about the other substrates/products becomes part of the vector components.
this is not a problem since in metabolism we can usually define a "main transformation" and as such is annotated in metabolic databases .
so, in this multidimensional "metabolic space", chemical compounds are points , chemical reactions can be seen as arrows  going from one point  to another , and biochemical pathways can be seen as sets of consecutive reaction vectors  .
datasets
we took from kegg  <cit>  all the metabolic reactions annotated for the model organism e. coli.
for each reaction, the substrate→product transformations annotated as "main" in the "rpair" field in kegg are taken.
the molecular structures for the two compounds involved in this transformation are also retrieved from kegg in .mol format.
hydrogens are added with the openbabel package http://www.openbabel.org.
open babel is also used to transform the .mol files into .pdb since the bond information required to calculate the triplets is easier to obtain from .pdb files.
finally, we calculate the two triplet vectors for the compounds, and the corresponding triplet vector for the reaction as the difference .
we end up with  <dig> compound vectors and  <dig> reaction vectors.
we can have more than one vector for a given ec code, depending on whether that enzyme activity is associated with  more than one reaction.
both reaction and compound vectors are defined by  <dig> components .
these triplets were not imposed a priori.
they are the ones with the highest frequencies in our dataset.
we excluded some triplets with very low frequency, involving rare atoms and appearing in very few compounds.
some triplets involve "unknown" atoms .
these come from molecules in kegg for which some parts are not detailed at the atomic level but indicated by a general group "-r".
these "r" groups  usually represent biological polymers .
in order to study the relationships between our vectorial representation of reactions and their functional characteristics, additional information on the enzymes catalyzing these reactions was obtained from other resources.
from brenda  <cit>  we retrieved the "reaction type", a set of one or more keywords describing the transformation carried out by the enzyme.
the total number of unique keywords in this dataset is  <dig>  from brenda we also retrieved the prosite motifs associated with the enzymes.
prosite  <cit>  is a database of short sequence motifs.
we retrieved  <dig> unique prosite motifs associated with the enzymes in our dataset.
interpro  <cit>  is a resource that integrates information from diverse databases defining protein domains according to different structural and sequence criteria, such as pfam, prodom, smart, scop, etc.
from interpro, we retrieved  <dig> unique signatures  associated with the enzymes in our dataset.
gene ontology  <cit>  defines three sets of terms for representing three different and independent aspects of the complex phenomena of protein function .
these terms are related by parenthood relationships defining three hierarchies which go from very general to more specific terms.
a protein  becomes annotated by associating it to one or more of these terms.
from gene ontology we retrieved the go terms in the "molecular function"  and "biological process"  sub-ontologies for our enzymes:  <dig> go terms .
in the following, we use the general term "keyword" referring to brenda keywords, prosite motifs, interpro domains and go terms.
additionally, to study in more detail the relationships between the compound vectors and the physico-chemical properties of the molecules they represent, a supplementary set of chemical compounds with associated physico-chemical information was obtained from http://cheminformatics.org/ that dataset was previously used by karthikeyan et al.
to develop a melting-point prediction method  <cit> .
the molecules in it are provided in smiles format http://www.daylight.com/smiles.
we convert from smiles to .pdb with openbabel, and then to triplet vectors as explained above.
for this dataset, owing to the more diverse atomic compositions of the molecules, we needed  <dig> atom triplets.
from the original  <dig> molecules in the dataset,  <dig> could not be represented in vectorial form for several reasons , resulting in a final dataset of  <dig> compounds.
hence, our final dataset contains  <dig> molecules with very diverse chemical properties; e.g.
the range of molecular weights is 84- <dig> and the range of logp is - <dig>  to + <dig> .
we used three physicochemical properties in this work: logp ; mw ; and tpsa .
the program spss was used to perform a multivariate linear regression analysis to relate the compound vectors to these chemical properties.
the weka package  <cit>  was used to perform classification studies aimed at evaluating the predictive value of this vectorial representation.
for that, the dataset of  <dig> compounds was randomly divided into  <dig> subsets.
nine of these  were used to perform the linear regression analysis, and the result was then used to predict the chemical property of the remaining group .
cycling this procedure for the other nine groups allowed predictions for all the compounds to be obtained.
the predicted values for the three properties were compared with the experimental values by linear correlation.
clusters of reaction vectors
a distance tree for all the reaction vectors was generated by the upgma algorithm implemented in the r package http://www.r-project.org, based on the euclidean distances between the vectors.
cutting this tree at different levels  produces different groupings with increasing number of clusters.
the "mclust" function of the "mclust" library of the same package was used to calculate the optimal number of clusters.
this function is based on the bic parameter .
we explored from  <dig> to  <dig> clusters and the optimal number according to bic criteria was  <dig>
we compared this clustering with equivalent classifications defined by the ec hierarchy using the enzyme annotations in the brenda, prosite, interpro and go resources .
for our reaction vectors  we can have any number of clusters depending on the level at which we cut the upgma tree, although the optimal number of clusters is  <dig>  on the contrary, the ec classification defines four possible groupings only, depending on the ec level .
in order to compare ec with reactome for equivalent number of clusters we obtain the groupings of the reactome with the same number of clusters .
for a given set of keywords  and a given clustering we created a contingency table containing the frequencies of each "keyword"  in each cluster, that is, the number of enzymes within that cluster associated with the keyword.
a chi-squared test was applied to each contingency table.
we used the yates correction  <cit> , which is applied to account for the low frequencies of some keywords and results in a more conservative test.
since different classifications have different degrees of freedom, we convert the chi-squared values  to z-scores  using fisher's approximation  <cit> , which approximates the chi-squared distribution by a normal distribution with mean  <dig> and σ =  <dig>
where ni are the degrees of freedom: *, k being the number of keywords in each dataset and n the number of clusters.
a high z-score indicates a relationship between the clustering and the set of keywords: elements  clustering together tend to have the same keywords, and vice versa.
to elucidate the biological meanings and characteristics of the  <dig> groups which form the "optimal" clustering of the reactome , we extracted the keywords and the components of the reaction vectors  specifically associated with each.
we did that by calculating, for each keyword/triplet in each of the  <dig> clusters, the p-value of its frequency within that cluster, taking into account the background frequency of that keyword/triplet.
for that we used the hypergeometric distribution and the poisson distribution for low frequency cases   <cit> .
RESULTS
assessment of the compound vectors
as observed in the discussion, there is a plethora of methods for representing chemical compounds in ways manageable by computers.
it is not the objective of this work to propose another description and compare it with existing ones.
nevertheless, we obviously have to evaluate the ability of this description to capture the chemical properties of compounds.
we do that by multivariate linear regression analysis using the chemical property as dependent variable and the components of the vector  as independent variables.
we apply this analysis to a database of  <dig> chemical compounds and to three different chemical properties: hydrophobicity , molecular weight  and polar surface area .
see methods for a full description of the procedure.
r is the pearson correlation coefficient.
such a clear relationship between the compound vectors and the chemical properties can be used for predictive purposes.
as described in methods, we re-calculate the linear regressions using 9/ <dig> of the vectors and use that regression to predict the chemical properties of the remaining 1/ <dig>  cycling these sets in a 10-fold cross-validation.
figure  <dig> shows the relationship between the predicted and experimental logp for the chemical compounds.
the predicted logp values tend to be similar to the experimental ones .
the results of the predictions for the other two properties are available as additional material .
these results are better than the ones for logp, as shown by the correlation  values .
nevertheless, in spite of their methodological value, these results are not useful for a real-world application .
that is why we focus on the results for logp here in spite of being slightly worse, since these are the more interesting ones.
hierarchical clustering of the reaction vectors for the e. coli metabolome
it can be seen that, in general, there is agreement between the ec and the vector-based classifications; clusters of similar vectors tend to contain enzymes with the same ec 1st number.
nevertheless, vectors with the same ec number can be far apart in this classification.
there are other important discrepancies, highlighting the value of the method presented here for quantifying similarities between enzymes.
for example, there is a group of reactions with identical vectors  carried out by very different enzymatic activities according to the ec classification:  <dig> .
<dig> ,  <dig> .
<dig> ,  <dig> .
<dig> ,  <dig> .
<dig> ,  <dig> .
<dig> ,  <dig> .
<dig> ,  <dig> .
<dig>  and  <dig> .
<dig> .
these enzymes, spanning three of the six main ec classes , ligases  and transferases ), are all nevertheless involved in the same chemical transformation, the conversion of an amide to a carboxylic acid: r-conh2→r-cooh.
so our vectorial representation captures a functional similarity between these enzymes that is not represented by the ec classification.
another interesting example is the group representing enzymes involved in phosphate hydrolysis.
these belong to three main ec classes  but are clustered together in our tree.
the hydrolysis of a pyrophosphate group is represented by another cluster, distant from the first, reflecting the different natures of these two processes.
this enzyme tree, in newick format, is available as supplementary material .
assessment of the biological significance of the clustering in comparison with the ec classification
we compared this clustering with the one defined by the ec classification for different number of clusters and according with different descriptors of biological significance .
table  <dig> shows the z-scores of the fisher test for different clustering schemas and different sets of keywords.
a z-score close to  <dig>  indicates no relationship between the clustering and the set of keywords, these keywords being uniformly distributed across all clusters.
it can be seen that both clustering schemas have very high z-scores for all sets of keywords, meaning that both schemas group similar enzymes together .
nevertheless, the z-scores for the reactome clusters automatically obtained from our vectorial representation are much higher than the equivalent ones defined by ec, and this is the case for all sets of keywords.
this means that the reactome clustering is more meaningful in biological terms.
the dependence between clusters and keywords is higher for this reactome clustering, which is especially striking because some of these sets of keywords, defined by experts, are partly based on the ec classification.
this is the case for brenda keywords, and is reflected in the observation that ec has higher z-scores  with brenda than with any of the other three datasets .
perhaps brenda shows the lowest difference between the two clustering schemas, and is least in agreement with the reactome clustering , because its keywords are partially based on ec.
the optimal reactome clustering with  <dig> groups, for which no equivalent exists in ec, is also included.
some ec groupings results in clusters without any enzyme associated to keywords in that particular dataset.
these are excluded from the analysis and are responsible for the small variations in the number of clusters .
for comparison, the reactome is clustered with the same  number of clusters.
characterization of the  <dig> clusters
according to the bic criteria , the optimal number of enzyme clusters within our dataset was  <dig>  we investigated the characteristics of these groups.
we did that by looking for keywords  and components of the enzyme vectors  that are differentially associated with each of the  <dig> clusters.
for each keyword and each triplet we calculated a p-value that quantifies the statistical significance of their association with a given cluster .
the table with detailed information on the significant keywords/triplets for the  <dig> clusters is available as supplementary material .
the biological meanings of most of the clusters can be elucidated on the basis of their associated keywords/triplets.
for example, cluster # <dig> is related to processes involving transfer/elimination of relatively large groups, as indicated by triplets involving c-c and c-o bonds  and keywords such as "elimination" and "condensation" , "*transferase", "*synthase" .
cluster # <dig> seems to be related to the addition of amino acid groups.
cluster # <dig> is clearly related to redox reactions: triplets involving c-o, c-h and o-h bonds, "oxidation", "redox reaction" , "oxidoreductase", "nad binding", "iron/sulphur cluster" , etc.
cluster # <dig> is related to nucleoside formation via n-glycosyl bonds.
cluster # <dig> seems to be related to sulfur metabolism.
cluster # <dig> is related to acyl group transfers and coa ligation.
cluster # <dig> is related to udp-n-acetylmuramate metabolism , so it could be specific to bacterial metabolism.
clusters # <dig> to # <dig> are all related to different aspects of phosphate group transfer, reflecting the importance of this process for the cell.
but each cluster seems to reflect a different aspect of this general process: while cluster # <dig> seems to be related to the transfer of pyrophosphate  groups, cluster # <dig> is related to hydrolysis of phosphate groups and # <dig> and # <dig> to metabolic kinases .
clusters # <dig>  # <dig> and # <dig> could be methodological "artifacts" arising from the way molecules are represented in kegg and the way we translate them into atom triplets.
in kegg, large groups are sometimes represented by "r-", without specifying their constituent atoms.
these are considered single atoms of unknown character by our system and represented by "x" .
so in our vectors we have some triplets involving "x".
the three clusters mentioned above are associated with these triplets.
nevertheless, since most of the "x" represent large transferable groups  these clusters tend to reflect processes involving transfers of large polymers, so they retain some biological meaning.
thus, these groups of reactions/enzymes detected in the "reactome" in an unsupervised manner provide a global view of the main processes in e. coli metabolism.
discussion
there are two main approaches to constructing biological classifications, which could be termed "top-down" and "bottom-up", or "supervised/unsupervised"  <cit> .
top-down approaches define a set of classes a priori, usually based on expert knowledge, and then assign the objects to them.
in contrast, bottom-up approaches construct the classification  from the data in an unsupervised way, i.e.
hierarchically clustering the data and looking for the resulting groups.
both approaches have advantages and drawbacks.
for protein sequences and structures we can find examples of both approaches, e.g.
scop  <cit>   vs. dali/fssp  <cit>   for classifying protein 3d structures.
for metabolic data, the top-down  approach dominates.
the main schema for classifying biochemical reactions  is very valuable for nomenclature purposes, but in many cases it is not clearly reflected at the molecular level.
while other ways of classifying chemical reactions based more on molecular details  are being developed  <cit> , the ec classification is widely extrapolated beyond its primary goal  on the assumption that it fully reflects molecular details.
the main problem for constructing bottom-up classifications of metabolic phenomena may be that there is no standard way of quantifying similarities among chemical compounds and, more especially, metabolic reactions, in contrast to protein sequences and structures.
there are many methods for representing the structures of small molecules and quantifying their similarities, each intended to a particular purpose  <cit> .
actually this lack of standard shows that the issue more complex than in the case of proteins.
some of these representations have been used to construct bottom-up classifications of the "metabolome"  <cit> , but this pioneering and interesting work has some drawbacks.
for example, it represents metabolic compounds by "imposing" the functional groups considered important.
from the outset, this restricts the kind of results that can be obtained.
the system classifies metabolites on the basis of the functional groups we a priori decide are important.
such "supervision" might hinder the discovery of features that are not related to those functional groups.
the situation is worse for the "reactome".
there are few approaches to quantifying the similarities among biochemical reactions.
these approaches either tackle very specific aspects of chemical transformations and have been applied to small sets of reactions  <cit> , or do not allow chemical compounds to be represented in the same space so that relationships between the "reactome" and the "metabolome" can be studied  <cit> , or are aimed at measuring similarities in mechanism rather than the transformation itself  <cit> .
for these reasons, although these approaches produced very interesting results for specific tasks, they are not easily applicable to the task of generating general global classifications of reactions.
CONCLUSIONS
as far as we know the work presented here is the first that aims to classify  the "reactome" in an unsupervised  way without any imposition of the kind of transformation or functional groups involved.
it is not the objective of this work to propose another way of representing molecular structures and compare it with existing ones.
the main goal is to allow chemical compounds, enzymatic activities and pathways to be represented concomitantly in the same vectorial space.
since we wish to unravel the global properties of the reactome, we want to keep our representation as simple as possible  so that the kind of results we can obtain is not restricted.
we show that, in spite of its simplicity, this representation captures important physicochemical properties of the molecules.
the representation of chemical compounds we use here does not impose any definition of "important groups".
we have previously used a similar representation for predicting the environmental fate of chemical compounds  <cit> .
that representation included information on the bond type in the triplets, as well as the molecular weight and the water solubility of the compounds.
these "heterogeneous" vectors with components representing different things and in different scales , are not the most adequate for representing the chemical reactions in the same space.
moreover, our results show that the information on the water solubility and the mw is implicit in the atom triplets.
the vectorial representation of chemical reactions used in this work does not take into account some characteristics of chemical transformations.
for example, it disregards stereochemistry, and aspect known to be very important in some catalyzed reactions.
on the other hand, it is general enough to represent and study whole reactomes.
we show that the groups of enzymatic activities arising from the unsupervised clustering of the e. coli reactome are more significant than the equivalent ec clusters and have biological meaning.
we can think of these groups as representing the main trends of the metabolic landscape of this organism.
the presence of many groups related to phosphate metabolism is interesting, highlighting the importance of this process for the cell.
the presence of a group devoted to reactions related to peptidoglycan metabolism shows that these results are organism-dependent, and opens interesting possibilities for comparing reactomes between different organisms.
we think that this new representation of chemical reactions opens many possibilities for the global study of nature's metabolic capabilities.
in the same way that global classification of the spaces of protein sequences and structures provided important biological information, studies of metabolism based on global classifications such as the one proposed here may also provide valuable data.
authors' contributions
fp conceived the original idea.
jct implemented the idea and carried out the tests.
fp and jct wrote the paper.
all authors read and approved the final manuscript.
supplementary material
additional file 1
correlation between experimental and predicted molecular weight  and total polar surface area .
click here for file
additional file 2
upgma tree of reaction vectors, in newick format.
this file can be viewed with most phylogenetic tree viewers, including the on-line tool itol .
click here for file
additional file 3
file with the information needed for coloring the tree "add2_upgma.nw" in itol  according to ec number.
click here for file
additional file 4
keywords associated with the  <dig> clusters.
click here for file
acknowledgements
we sincerely thank alfonso valencia , victor de lorenzo  and manuel gomez  for interesting discussions.
we also thank antonio rausell  for statistical advice.
this work was partially funded by projects bio2006- <dig> and pie 200620i <dig> from the spanish ministry for education and science.

 the "genomic era" is being characterized by the massive determination of the molecular components of living systems. genome sequencing projects are yielding complete genome sequences for hundreds of organisms. although not as massive as envisaged, structural genomics projects are speeding up the rate of protein structure determinations. many other initiatives also address large-scale repertories of molecular components and their relationships, seeking to decipher the "transcriptome", the "interactome", etc. these massive data contain much information about living organisms. studied as a whole, from a systemic perspective, they provide global pictures of different aspects of biology, which can help to answer very basic questions about how life evolved and how organisms do what they do with their "molecular toolkits". for example, the repertory of protein sequences and their evolutionary relationships can be represented in a "sequence space". studying this global landscape of protein sequences as a whole has produced important information on the estimated total number of protein families, functional groups, evolutionary relationships, etc.. similarly, the wealth of information on protein three-dimensional structures has allowed comprehensive hierarchical classifications such as scop and cath to be generated. these classifications can subsequently be mined to perform global studies on the estimated total number of folds, the evolution of protein structures, the relationship between sequence and functional spaces, etc.. metabolism was one of the last aspects to be studied from an "-omic" perspective. the "metabolome" can be defined as the total complement of small-molecule metabolites within a cell, while the "reactome" is the set of chemical reactions that transforms these metabolites. in many cases, "metabolone" refers not only to the complement of small molecules but to quantitative measurements of their concentrations as well. an important difference from the "proteome" and other "-omics" is that most metabolic information has not been obtained by a high-throughput process but by the accumulation of knowledge over many years of detailed biochemical characterization of compounds, reactions and enzymes, knowledge that has been stored in databases such as kegg and brenda. only recently, high-througput techniques have been applied to the massive determination of metabolomes and, even more recently, reactomes. this has advantages and drawbacks. on the one hand, metabolic data are more reliable and do not suffer from the high degree of error common to high-throughput-omics data, such as the "interactome". on the other, biochemical knowledge can be biased to certain pathways and processes that have been more intensely studied than others. one of the main problems in performing global studies of the reactome or the metabolome, equivalent to the studies discussed above for the proteome, is that there is no good way of representing and quantifying metabolic information. it is quite straightforward to quantify the similarity between the sequences or three-dimensional structures of two proteins. it is even possible to discretize these continuous scales of similarity to define intervals of evolutionary and functional relationship, expected structural similarity for a given sequence similarity, etc. the possibility of representing protein sequences and structures, as well as quantifying their corresponding similarities, is crucial for generating global classifications and subsequent global studies such as those cited above. however, quantifying similarities among chemical compounds is not as straightforward as it is for protein sequences and structures. there are many different ways of quantifying similarities between small molecules, which can be used to study the global characteristics of the "metabolic space" in a similar way to sequence and structure spaces. the main problem arises with chemical reactions: there is no easy way of quantifying the "similarity between two chemical reactions". the few existing approaches are not intended to generate detailed global classifications of biochemical reactions in an unsupervised way or to perform global studies of metabolomes and reactomes aimed at quantifying their chemical diversities. indeed, the most widely-used classification of chemical reactions schema) is neither quantitative nor based on existing data. rather, ec is based on an a priori imposed schema that organizes all existing enzymatic activities into a hierarchical classification of four levels with six main classes in the first level. it was designed mainly for nomenclature purposes and lacks a clear reflection at the molecular or evolutionary levels. although very useful for many tasks, this classification was not designed for quantitative global classifications or comparisons of biochemical activities. in this work, we propose a very simple way of representing chemical reactions in a vectorial form, which allows them to be compared quantitatively. this representation is based on an equivalent vectorial representation of chemical compounds, which imposes no a priori definition of functional groups or other "important" parts of molecules. it allows any chemical species to be represented in the same vectorial space. with this description, chemical compounds, reactions and pathways can be concomitantly represented in the same vectorial space. we demonstrate that the vectorial representation of chemical compounds, in spite of its simplicity, efficiently captures their molecular properties and can be used for predictive purposes. we used this method to generate an unsupervised global clustering of metabolic activities. this clustering makes more biological sense than equivalent groupings defined by the ec hierarchy, an improvement that is evaluated statistically using different criteria. moreover, it results in an optimal set of groups of reactions for which we look for their biological meaning. methods vectorial representation of chemical compounds we looked for a metabolic representation that: allows all chemical species to be represented in the same space regardless of their composition, molecular size, etc.; imposes no a priori definition of important chemical features, since a representation based on a pre-defined set of functional groups would necessarily restrict the results obtained to those that can be given in terms of those groups; and allows chemical reactions and pathways to be represented concomitantly in the same space. we represent a chemical structure with a vector in which each component represents a triplet of neighboring atoms, and its value indicates the frequency of that triplet in the molecule. using pairs of atoms instead of triplets would lead to shorter vectors, but would not capture many aspects of the molecule's functionality; the chemical properties of an atom depend mostly on its neighbors up to distance for example, the different properties of the carbonyl group c = o depend on the neighbors of the c: oh-cr = o, h-cr = o, r-cr = o, etc. on the other hand, groups of... neighboring atoms would produce vectors with too many components. we previously used a similar representation to train a machine learning system for predicting the biodegradability features of chemical compounds. this representation is not perfectly univocal, since in some cases different molecules can have the same vector. nevertheless, it provides a good balance between information content and utility for the goal of this work. the representation contains more structural information than is first apparent. for example the bond types, although not explicitly considered in the representation, are implicitly taken into account to some extent: as we go from hydrocarbons with few double and triple bonds to hydrocarbons with many multiple bonds, we gain ccc triplets at the expense of hch and hcc. in the "results" section, we demonstrate that this representation captures important physico-chemical properties of compounds and can be used for predictive purposes. vectorial representation of biochemical reactions the vector for a reaction that transforms one chemical compound into another is just the difference between the vectors of the two compounds. the stoichiometry is taken into account by multiplying the compound vectors by the corresponding stoichiometric coefficients: e.g. for a reaction a -> 2b the reaction vector would be va->b = 2·vb - va. since this vectorial representation of reactions is based on that of the compounds, it inheres many of its characteristics: any chemical reaction can be represented in the same space regardless of its characteristics; there is no a priori definition of important reaction features. moreover, the components of the reaction vector represent the triplets that are changed in the chemical reaction so they neatly encode the characteristics of the transformation. another intuitive feature is that the vector of a given reaction is the opposite of the vector of the inverse reaction. in this representation, all chemical reactions have to be simplified as transformations of one single compound into another. for a reaction with more than one reactant and/or product, a "main" pair substrate ->product has to be chosen and information about the other substrates/products becomes part of the vector components. this is not a problem since in metabolism we can usually define a "main transformation" and as such is annotated in metabolic databases. so, in this multidimensional "metabolic space", chemical compounds are points, chemical reactions can be seen as arrows going from one point to another, and biochemical pathways can be seen as sets of consecutive reaction vectors. datasets we took from kegg all the metabolic reactions annotated for the model organism e. coli. for each reaction, the substrate→product transformations annotated as "main" in the "rpair" field in kegg are taken. the molecular structures for the two compounds involved in this transformation are also retrieved from kegg in.mol format. hydrogens are added with the openbabel package http://www.openbabel.org. open babel is also used to transform the.mol files into.pdb since the bond information required to calculate the triplets is easier to obtain from.pdb files. finally, we calculate the two triplet vectors for the compounds, and the corresponding triplet vector for the reaction as the difference. we end up with compound vectors and reaction vectors. we can have more than one vector for a given ec code, depending on whether that enzyme activity is associated with more than one reaction. both reaction and compound vectors are defined by components. these triplets were not imposed a priori. they are the ones with the highest frequencies in our dataset. we excluded some triplets with very low frequency, involving rare atoms and appearing in very few compounds. some triplets involve "unknown" atoms. these come from molecules in kegg for which some parts are not detailed at the atomic level but indicated by a general group "-r". these "r" groups usually represent biological polymers. in order to study the relationships between our vectorial representation of reactions and their functional characteristics, additional information on the enzymes catalyzing these reactions was obtained from other resources. from brenda we retrieved the "reaction type", a set of one or more keywords describing the transformation carried out by the enzyme. the total number of unique keywords in this dataset is from brenda we also retrieved the prosite motifs associated with the enzymes. prosite is a database of short sequence motifs. we retrieved unique prosite motifs associated with the enzymes in our dataset. interpro is a resource that integrates information from diverse databases defining protein domains according to different structural and sequence criteria, such as pfam, prodom, smart, scop, etc. from interpro, we retrieved unique signatures associated with the enzymes in our dataset. gene ontology defines three sets of terms for representing three different and independent aspects of the complex phenomena of protein function. these terms are related by parenthood relationships defining three hierarchies which go from very general to more specific terms. a protein becomes annotated by associating it to one or more of these terms. from gene ontology we retrieved the go terms in the "molecular function" and "biological process" sub-ontologies for our enzymes: go terms. in the following, we use the general term "keyword" referring to brenda keywords, prosite motifs, interpro domains and go terms. additionally, to study in more detail the relationships between the compound vectors and the physico-chemical properties of the molecules they represent, a supplementary set of chemical compounds with associated physico-chemical information was obtained from http://cheminformatics.org/ that dataset was previously used by karthikeyan et al. to develop a melting-point prediction method. the molecules in it are provided in smiles format http://www.daylight.com/smiles. we convert from smiles to.pdb with openbabel, and then to triplet vectors as explained above. for this dataset, owing to the more diverse atomic compositions of the molecules, we needed atom triplets. from the original molecules in the dataset, could not be represented in vectorial form for several reasons, resulting in a final dataset of compounds. hence, our final dataset contains molecules with very diverse chemical properties; e.g. the range of molecular weights is 84- and the range of logp is - to +. we used three physicochemical properties in this work: logp; mw; and tpsa. the program spss was used to perform a multivariate linear regression analysis to relate the compound vectors to these chemical properties. the weka package was used to perform classification studies aimed at evaluating the predictive value of this vectorial representation. for that, the dataset of compounds was randomly divided into subsets. nine of these were used to perform the linear regression analysis, and the result was then used to predict the chemical property of the remaining group. cycling this procedure for the other nine groups allowed predictions for all the compounds to be obtained. the predicted values for the three properties were compared with the experimental values by linear correlation. clusters of reaction vectors a distance tree for all the reaction vectors was generated by the upgma algorithm implemented in the r package http://www.r-project.org, based on the euclidean distances between the vectors. cutting this tree at different levels produces different groupings with increasing number of clusters. the "mclust" function of the "mclust" library of the same package was used to calculate the optimal number of clusters. this function is based on the bic parameter. we explored from to clusters and the optimal number according to bic criteria was we compared this clustering with equivalent classifications defined by the ec hierarchy using the enzyme annotations in the brenda, prosite, interpro and go resources. for our reaction vectors we can have any number of clusters depending on the level at which we cut the upgma tree, although the optimal number of clusters is on the contrary, the ec classification defines four possible groupings only, depending on the ec level. in order to compare ec with reactome for equivalent number of clusters we obtain the groupings of the reactome with the same number of clusters. for a given set of keywords and a given clustering we created a contingency table containing the frequencies of each "keyword" in each cluster, that is, the number of enzymes within that cluster associated with the keyword. a chi-squared test was applied to each contingency table. we used the yates correction, which is applied to account for the low frequencies of some keywords and results in a more conservative test. since different classifications have different degrees of freedom, we convert the chi-squared values to z-scores using fisher's approximation, which approximates the chi-squared distribution by a normal distribution with mean and σ = where ni are the degrees of freedom: *, k being the number of keywords in each dataset and n the number of clusters. a high z-score indicates a relationship between the clustering and the set of keywords: elements clustering together tend to have the same keywords, and vice versa. to elucidate the biological meanings and characteristics of the groups which form the "optimal" clustering of the reactome, we extracted the keywords and the components of the reaction vectors specifically associated with each. we did that by calculating, for each keyword/triplet in each of the clusters, the p-value of its frequency within that cluster, taking into account the background frequency of that keyword/triplet. for that we used the hypergeometric distribution and the poisson distribution for low frequency cases.  assessment of the compound vectors as observed in the discussion, there is a plethora of methods for representing chemical compounds in ways manageable by computers. it is not the objective of this work to propose another description and compare it with existing ones. nevertheless, we obviously have to evaluate the ability of this description to capture the chemical properties of compounds. we do that by multivariate linear regression analysis using the chemical property as dependent variable and the components of the vector as independent variables. we apply this analysis to a database of chemical compounds and to three different chemical properties: hydrophobicity, molecular weight and polar surface area. see methods for a full description of the procedure. r is the pearson correlation coefficient. such a clear relationship between the compound vectors and the chemical properties can be used for predictive purposes. as described in methods, we re-calculate the linear regressions using 9/ of the vectors and use that regression to predict the chemical properties of the remaining 1/ cycling these sets in a 10-fold cross-validation. figure shows the relationship between the predicted and experimental logp for the chemical compounds. the predicted logp values tend to be similar to the experimental ones. the results of the predictions for the other two properties are available as additional material. these results are better than the ones for logp, as shown by the correlation values. nevertheless, in spite of their methodological value, these results are not useful for a real-world application. that is why we focus on the results for logp here in spite of being slightly worse, since these are the more interesting ones. hierarchical clustering of the reaction vectors for the e. coli metabolome it can be seen that, in general, there is agreement between the ec and the vector-based classifications; clusters of similar vectors tend to contain enzymes with the same ec 1st number. nevertheless, vectors with the same ec number can be far apart in this classification. there are other important discrepancies, highlighting the value of the method presented here for quantifying similarities between enzymes. for example, there is a group of reactions with identical vectors carried out by very different enzymatic activities according to the ec classification:.,.,.,.,.,.,. and.. these enzymes, spanning three of the six main ec classes, ligases and transferases ), are all nevertheless involved in the same chemical transformation, the conversion of an amide to a carboxylic acid: r-conh2→r-cooh. so our vectorial representation captures a functional similarity between these enzymes that is not represented by the ec classification. another interesting example is the group representing enzymes involved in phosphate hydrolysis. these belong to three main ec classes but are clustered together in our tree. the hydrolysis of a pyrophosphate group is represented by another cluster, distant from the first, reflecting the different natures of these two processes. this enzyme tree, in newick format, is available as supplementary material. assessment of the biological significance of the clustering in comparison with the ec classification we compared this clustering with the one defined by the ec classification for different number of clusters and according with different descriptors of biological significance. table shows the z-scores of the fisher test for different clustering schemas and different sets of keywords. a z-score close to indicates no relationship between the clustering and the set of keywords, these keywords being uniformly distributed across all clusters. it can be seen that both clustering schemas have very high z-scores for all sets of keywords, meaning that both schemas group similar enzymes together. nevertheless, the z-scores for the reactome clusters automatically obtained from our vectorial representation are much higher than the equivalent ones defined by ec, and this is the case for all sets of keywords. this means that the reactome clustering is more meaningful in biological terms. the dependence between clusters and keywords is higher for this reactome clustering, which is especially striking because some of these sets of keywords, defined by experts, are partly based on the ec classification. this is the case for brenda keywords, and is reflected in the observation that ec has higher z-scores with brenda than with any of the other three datasets. perhaps brenda shows the lowest difference between the two clustering schemas, and is least in agreement with the reactome clustering, because its keywords are partially based on ec. the optimal reactome clustering with groups, for which no equivalent exists in ec, is also included. some ec groupings results in clusters without any enzyme associated to keywords in that particular dataset. these are excluded from the analysis and are responsible for the small variations in the number of clusters. for comparison, the reactome is clustered with the same number of clusters. characterization of the clusters according to the bic criteria, the optimal number of enzyme clusters within our dataset was we investigated the characteristics of these groups. we did that by looking for keywords and components of the enzyme vectors that are differentially associated with each of the clusters. for each keyword and each triplet we calculated a p-value that quantifies the statistical significance of their association with a given cluster. the table with detailed information on the significant keywords/triplets for the clusters is available as supplementary material. the biological meanings of most of the clusters can be elucidated on the basis of their associated keywords/triplets. for example, cluster # is related to processes involving transfer/elimination of relatively large groups, as indicated by triplets involving c-c and c-o bonds and keywords such as "elimination" and "condensation", "*transferase", "*synthase". cluster # seems to be related to the addition of amino acid groups. cluster # is clearly related to redox reactions: triplets involving c-o, c-h and o-h bonds, "oxidation", "redox reaction", "oxidoreductase", "nad binding", "iron/sulphur cluster", etc. cluster # is related to nucleoside formation via n-glycosyl bonds. cluster # seems to be related to sulfur metabolism. cluster # is related to acyl group transfers and coa ligation. cluster # is related to udp-n-acetylmuramate metabolism, so it could be specific to bacterial metabolism. clusters # to # are all related to different aspects of phosphate group transfer, reflecting the importance of this process for the cell. but each cluster seems to reflect a different aspect of this general process: while cluster # seems to be related to the transfer of pyrophosphate groups, cluster # is related to hydrolysis of phosphate groups and # and # to metabolic kinases. clusters # # and # could be methodological "artifacts" arising from the way molecules are represented in kegg and the way we translate them into atom triplets. in kegg, large groups are sometimes represented by "r-", without specifying their constituent atoms. these are considered single atoms of unknown character by our system and represented by "x". so in our vectors we have some triplets involving "x". the three clusters mentioned above are associated with these triplets. nevertheless, since most of the "x" represent large transferable groups these clusters tend to reflect processes involving transfers of large polymers, so they retain some biological meaning. thus, these groups of reactions/enzymes detected in the "reactome" in an unsupervised manner provide a global view of the main processes in e. coli metabolism. discussion there are two main approaches to constructing biological classifications, which could be termed "top-down" and "bottom-up", or "supervised/unsupervised". top-down approaches define a set of classes a priori, usually based on expert knowledge, and then assign the objects to them. in contrast, bottom-up approaches construct the classification from the data in an unsupervised way, i.e. hierarchically clustering the data and looking for the resulting groups. both approaches have advantages and drawbacks. for protein sequences and structures we can find examples of both approaches, e.g. scop vs. dali/fssp for classifying protein 3d structures. for metabolic data, the top-down approach dominates. the main schema for classifying biochemical reactions is very valuable for nomenclature purposes, but in many cases it is not clearly reflected at the molecular level. while other ways of classifying chemical reactions based more on molecular details are being developed, the ec classification is widely extrapolated beyond its primary goal on the assumption that it fully reflects molecular details. the main problem for constructing bottom-up classifications of metabolic phenomena may be that there is no standard way of quantifying similarities among chemical compounds and, more especially, metabolic reactions, in contrast to protein sequences and structures. there are many methods for representing the structures of small molecules and quantifying their similarities, each intended to a particular purpose. actually this lack of standard shows that the issue more complex than in the case of proteins. some of these representations have been used to construct bottom-up classifications of the "metabolome", but this pioneering and interesting work has some drawbacks. for example, it represents metabolic compounds by "imposing" the functional groups considered important. from the outset, this restricts the kind of results that can be obtained. the system classifies metabolites on the basis of the functional groups we a priori decide are important. such "supervision" might hinder the discovery of features that are not related to those functional groups. the situation is worse for the "reactome". there are few approaches to quantifying the similarities among biochemical reactions. these approaches either tackle very specific aspects of chemical transformations and have been applied to small sets of reactions, or do not allow chemical compounds to be represented in the same space so that relationships between the "reactome" and the "metabolome" can be studied, or are aimed at measuring similarities in mechanism rather than the transformation itself. for these reasons, although these approaches produced very interesting results for specific tasks, they are not easily applicable to the task of generating general global classifications of reactions.  as far as we know the work presented here is the first that aims to classify the "reactome" in an unsupervised way without any imposition of the kind of transformation or functional groups involved. it is not the objective of this work to propose another way of representing molecular structures and compare it with existing ones. the main goal is to allow chemical compounds, enzymatic activities and pathways to be represented concomitantly in the same vectorial space. since we wish to unravel the global properties of the reactome, we want to keep our representation as simple as possible so that the kind of results we can obtain is not restricted. we show that, in spite of its simplicity, this representation captures important physicochemical properties of the molecules. the representation of chemical compounds we use here does not impose any definition of "important groups". we have previously used a similar representation for predicting the environmental fate of chemical compounds. that representation included information on the bond type in the triplets, as well as the molecular weight and the water solubility of the compounds. these "heterogeneous" vectors with components representing different things and in different scales, are not the most adequate for representing the chemical reactions in the same space. moreover, our results show that the information on the water solubility and the mw is implicit in the atom triplets. the vectorial representation of chemical reactions used in this work does not take into account some characteristics of chemical transformations. for example, it disregards stereochemistry, and aspect known to be very important in some catalyzed reactions. on the other hand, it is general enough to represent and study whole reactomes. we show that the groups of enzymatic activities arising from the unsupervised clustering of the e. coli reactome are more significant than the equivalent ec clusters and have biological meaning. we can think of these groups as representing the main trends of the metabolic landscape of this organism. the presence of many groups related to phosphate metabolism is interesting, highlighting the importance of this process for the cell. the presence of a group devoted to reactions related to peptidoglycan metabolism shows that these results are organism-dependent, and opens interesting possibilities for comparing reactomes between different organisms. we think that this new representation of chemical reactions opens many possibilities for the global study of nature's metabolic capabilities. in the same way that global classification of the spaces of protein sequences and structures provided important biological information, studies of metabolism based on global classifications such as the one proposed here may also provide valuable data. authors' contributions fp conceived the original idea. jct implemented the idea and carried out the tests. fp and jct wrote the paper. all authors read and approved the final manuscript. supplementary material additional file 1 correlation between experimental and predicted molecular weight and total polar surface area. click here for file additional file 2 upgma tree of reaction vectors, in newick format. this file can be viewed with most phylogenetic tree viewers, including the on-line tool itol. click here for file additional file 3 file with the information needed for coloring the tree "add2_upgma.nw" in itol according to ec number. click here for file additional file 4 keywords associated with the clusters. click here for file acknowledgements we sincerely thank alfonso valencia, victor de lorenzo and manuel gomez for interesting discussions. we also thank antonio rausell for statistical advice. this work was partially funded by projects bio2006- and pie 200620i from the spanish ministry for education and science.
we show that the representation of chemical compounds reflects their physicochemical properties and can be used for predictive purposes. we show that this unsupervised clustering results in groups of enzymes more coherent in biological terms than equivalent groupings obtained from the ec hierarchy. we use the vectorial representations of reactions to perform a global classification of the reactome of the model organism e. coli. in this representation, chemical compounds, reactions and pathways may be represented in the same vectorial space. because there is no good way of measuring similarities between chemical reactions, there is a barrier to the development of global classifications of "metabolic space" and subsequent studies comparable to those done for protein sequences and structures. this hierarchical clustering produces an optimal set of groups which we analyzed for their biological meaning. in this work, we propose a vectorial representation of chemical reactions, which allows them to be compared and classified. global studies of the protein repertories of organisms are providing important information on the characteristics of the protein space.
 global studies of the protein repertories of organisms are providing important information on the characteristics of the protein space. many of these studies entail classification of the protein repertory on the basis of structure and/or sequence similarities. the situation is different for metabolism. because there is no good way of measuring similarities between chemical reactions, there is a barrier to the development of global classifications of "metabolic space" and subsequent studies comparable to those done for protein sequences and structures.  in this work, we propose a vectorial representation of chemical reactions, which allows them to be compared and classified. in this representation, chemical compounds, reactions and pathways may be represented in the same vectorial space. we show that the representation of chemical compounds reflects their physicochemical properties and can be used for predictive purposes. we use the vectorial representations of reactions to perform a global classification of the reactome of the model organism e. coli.  we show that this unsupervised clustering results in groups of enzymes more coherent in biological terms than equivalent groupings obtained from the ec hierarchy. this hierarchical clustering produces an optimal set of groups which we analyzed for their biological meaning. 