Construction of phylogenetic trees by kernel-based comparative analysis of metabolic networks

Background To infer the tree of life requires knowledge of the common characteristics of each species descended from a common ancestor as the measuring criteria and a method to calculate the distance between the resulting values of each measure. Conventional phylogenetic analysis based on genomic sequences provides information about the genetic relationships between different organisms. In contrast, comparative analysis of metabolic pathways in different organisms can yield insights into their functional relationships under different physiological conditions. However, evaluating the similarities or differences between metabolic networks is a computationally challenging problem, and systematic methods of doing this are desirable. Here we introduce a graph-kernel method for computing the similarity between metabolic networks in polynomial time, and use it to profile metabolic pathways and to construct phylogenetic trees. Results To compare the structures of metabolic networks in organisms, we adopted the exponential graph kernel, which is a kernel-based approach with a labeled graph that includes a label matrix and an adjacency matrix. To construct the phylogenetic trees, we used an unweighted pair-group method with arithmetic mean, i.e., a hierarchical clustering algorithm. We applied the kernel-based network profiling method in a comparative analysis of nine carbohydrate metabolic networks from 81 biological species encompassing Archaea, Eukaryota, and Eubacteria. The resulting phylogenetic hierarchies generally support the tripartite scheme of three domains rather than the two domains of prokaryotes and eukaryotes. Conclusion By combining the kernel machines with metabolic information, the method infers the context of biosphere development that covers physiological events required for adaptation by genetic reconstruction. The results show that one may obtain a global view of the tree of life by comparing the metabolic pathway structures using meta-level information rather than sequence information. This method may yield further information about biological evolution, such as the history of horizontal transfer of each gene, by studying the detailed structure of the phylogenetic tree constructed by the kernel-based method.


Background
The availability of pathway databases such as Kyoto Encyclopedia of Genes and Genomes (KEGG), What is there? (WIT3), PathDB, and MetaCyc opens up various new possibilities for comparative analysis. In particular, information about metabolic pathways in different organisms yields important information about their evolution and offers a complementary approach to phylogenetic analysis. Here we present a comparative metabolomic approach to constructing phylogenetic trees that uses physiological functions of the organisms by computing the structural similarity of metabolic networks. The consideration of metabolic components complements the conventional approaches to phylogeny based on genome sequences. Recognizing the similarities and differences in metabolic functions between species may provide insights into other applications in biotechnology, ecology, and evolutionary studies. Several researchers have attempted to rebuild evolutionary history by comparing ribosomal RNA sequences [1], by phylogenomics [2], or by comparing whole genomes to overcome the limitations of the gene-sequence analyses [3][4][5].
Several recent studies have extended conventional phylogenetic analysis to incorporate metabolic pathway information. Forst and Schulten [6,7] presented one of the earliest approaches to extend the conventional sequence comparison and phylogenetic analysis of individual enzymes to metabolic networks. They also presented a method to calculate distances between metabolic networks based on sequence information of the biomolecules involved and information about the corresponding reaction networks. Dandekar et al. (1999) combined strategies in a systematic comparison of the enzymes and corresponding sequence information of the glycolytic pathway [8]. Other approaches involving the reconstructed phylogenies from gene-order data have been based on simulating genome evolution [9], and studying the genome evolution resulting from the metabolic adaptation of the organism to the surrounding environment. Liao et al. (2002) presented a method to group organisms by comparing the profiles of metabolic pathways, where the profiling was based simply on binary attributes (e.g., by denoting the presence or absence of pathways in the organisms) [10].
Whereas the previous approaches incorporated information about the additional metabolic pathways, systematic methods to calculate the similarities between metabolic networks are lacking or contain gaps in some of the biological assumptions. In this paper, we introduce the concept of graph kernels to calculate the similarities between two different network structures. The graph kernel-based approach can compute more efficiently the similarity of two graph structures by the kernel function that can extract important features from the graph. Our approach contrasts with that of Forst and Schulten [6,7] in that the graph kernel calculates the distance based on the network level instead of on its sequence information on the biomolecules involved.
In their comparative analysis of metabolic pathways, Heymans and Singh [11] showed that phylogenetic trees could be made from the graph similarities of metabolic networks. They applied a distance measure between metabolic graphs of the glycolytic pathway and the citric acid cycle from 16 organisms. However, some of their data on phylogenetic inference did not correspond entirely with the conventional taxonomy and did not provide a global view of the specialization of species according to the scale of analyzed species and metabolic pathways.
Several more recent attempts have reconstructed genome trees using different formalisms such as gene ordering [12,13], measuring gene contents [3,14], comparing sequence similarities [15,16], comparing proteome strings [17], and phylogenomics [18,19]. All are based on the principle of genome sequences, but none has applied the concepts of effectible physiology to the phylogenetic analyses. We report on our comparative results and discuss our findings.

Results and discussion
We chose nine pathways of carbohydrate metabolism and 81 species to perform a comparative analysis of metabolic networks that satisfied the most abundant dataset from KEGG (Table 1, 2 and 3). Our sample comprised 13 species of Archaea, eight species of Eukaryota, and 60 species of Eubacteria. The central pathways of metabolism include the glycolytic and pentose phosphate pathways, and the citric acid cycle, which generate biological energy and form the metabolic precursors essential for almost all living cells. To validate the data, we investigated the distribution of each enzyme in the nine pathways. Several enzymes appear at a high frequency in all species, and this frequency decreases rather exponentially as the value of xaxis increases (Figure 1), a phenomenon we observed in all pathways studied. Because the characteristics of the enzyme distribution did not differ, all pathways can be used in the phylogenetic analysis. We observed no obvious tendency of shift or deviation of the distributions from the overall pattern of pathways.
The phylogeny took two directions: conventional taxonomy that focused on the morphological and physiological features to classify species, and a numerical taxonomy that stressed the historical changes in biological sequences. Phylogenies based on the ribosomal RNA molecules led to the proposal of a new tripartite scheme of three domains: Bacteria, Archaea, and Eukarya [20]. Although each approach is feasible on its own, it cannot provide a holistic view of the organism. Current phylogenetic studies indicate that horizontal gene transfer may have played a vital role in the evolution of major lineages [21]. Lake and Moore [22] also noted the pitfalls of comparative genomics based on molecular sequences. Our kernelbased method provides an alternative to the inference of an evolutionary scenario and allows for a higher-level comparison of the phylogenetic trees by measuring the distances between pathways using metabolic network data to infer an evolutionary scenario. The distribution of enzymes in nine reference pathways Figure 1 The distribution of enzymes in nine reference pathways. The x axis is the index of the maximum value of the order sorted by the frequency of enzyme over all pathways. All plots showed similar distributions.  Figure 2 compares the phylogenetic tree constructed by the kernel method using (a) metabolic networks, and (b) by the sequence-based analysis of phosphoglycerate kinase and phosphopyruvate hydratase. The overall tree structure appears to be consistent with the current classification of the three domains, Eubacteria, Archaea, and Eukaryota ( Figure 3). In the domain of Archaea, methanogens such as Methanosarcina, Methanococcus, and Methanobacterium are clustered in the early branch of Pyrococcus.

Consistency with conventional taxonomy
Five thermophiles are also clustered in the same branch with a halophile, Halobacterium. The structure of collective carbohydrate pathways from Arabidopsis thaliana has the smallest distance from yeasts ( Figure 2  parison of the information-transfer pathways and pathway-level organization between two domains [23], whereas eukaryotic metabolic enzymes are primarily of bacterial origin [24].

Inferring hidden order by network clustering
The conventional sequence-based analysis passes over or does not embrace the discordant evolution of each species or the horizontal gene transfer [25]. Our method can cope with this limitation by taking into account the structural features of individual metabolic networks. The disagreement between the molecular sequence data of operational genes and the rRNA tree suggests that different genes have different evolutionary histories [26,27]. To address this problem, Li studied the mitochondrial genomes in relation to the problem of whole-genome phylogeny, where evolutionary events, such as genetic rearrangements that include gene transfer from the exterior, make genome alignments difficult [28].
To compare the kernel-based comparative analysis of metabolic networks to the sequence-based phylogenetic analysis, we analyzed two enzyme sequences that participate in carbohydrate metabolism in all 81 species together using a multiple sequence alignment (Figure 2(b)). In the resulting phylogenetic tree, Archaea and Eukaryota are clustered at each terminal; however, short-distance neighboring node members belong to fairly distant taxonomic groups. The overall structure of the tree eventually becomes remote from not only that of the kernel-based method (Figure 2(a)), but also from that of current taxonomy ( Figure 3). Although the phylogenetic tree constructed from the multiple alignment of two enzyme sequences shows a few unusual characteristics, our approach provides a good solution. The cluster mainly comprised archaeal species including the bacterial members Chlamydia and Chlamydophila, and had long branches at the root of the tree. Three eubacterial members (Ana, Tel, and Syn) are more closely related to Eukaryota than Current classification of biological taxonomy Eubacteria. Moreover, eubacterial groups are separated over the topology, and the Eukaryota are inserted between them (Figure 2(b)). This result shows an example of the sequence-based phylogenetic approach when we intend to perform phylogenetic analyses for as many species as possible. The sequence-based phylogenetic analysis can fail to precisely represent evolutionary history without an analysis using a set of whole sequences. Unfortunately, analyses using whole sequences require massive computing power and are highly complex. A sequence-based phylogenetic analysis can still be limited to cover a number of species, although alternative approaches exist. However, our method can easily solve this problem by utilizing given resources. In order to measure the quality of our constructed trees, we compared the phylogenetic tree based on our graph kernel method with that of Heymans and Singh [11] in terms of their similarity to conventional taxonomy. We used a software tool, 'Cousins' [31], which compares two alternative phylogenetic trees based on common cousin pairs in the trees. The comparison by cousin pairs is said to more focus on local similarity between two trees because it evaluates the similarity based on a cousin pair within a certain degree. Table 4 shows the similarity score of our kernel-based method in comparison with [11] for the glycolysis pathways of 65 organisms. Our method shows a better result in terms of the similarity score with the conventional NCBI taxonomy.
In this paper, we intended to present a meta-level analysis of biological systems to construct a unitary phylogenetic tree that could be used to interpret the context of biological evolution. Our results suggest that the phylogenetic analysis with submetabolic network information might also allow us to infer horizontal or lateral gene transfer. Our results also support the tripartite scheme of the three domains, Bacteria, Archaea and Eukaryota [20].
Comparing the pathogenic bacterial genomes by focusing on the pathways of bacterial and eukaryotic aminoacyl-tRNA synthesis showed that this pathway is uniquely prokaryotic/archaeal and that it is found widely among the pathogenic bacteria. This suggests that members of this pathway can be used as targets for novel antimicrobial drugs [32]. Metabolic analysis of pathogenic organisms may play a critical role in the selective treatment or pre-vention of diseases caused by these organisms by using this innovative concept to develop new drugs.

Conclusion
Biological classification, taxonomy, and systematics are the profound themes in biology. Using phylogeny in evolutionary classification implies functional and morphological innovation, adaptive range, parallelism, and convergence. We have used a method based on the graph kernel to compare information on each metabolic network including cardinality, distance, and topology relating to metabolic networks as a type of undirected graph. Our results showed that our approach has potential in the macroscopic analysis of phylogenetic relationships among organisms in relation to horizontal gene transfer. To obtain information about each causal mechanism in the context of a similar phenotype, one should first analyze the phenomena at the level of a protein network. The analysis of a metabolic network is an example of this type of analysis. Biological entities that interact with the environment and eventually influence adaptation are a function of the activity of proteins and other bioactive molecules rather than gene order or genetic history.
The overall structure of the phylogenetic tree constructed from our experiments supports the tripartite scheme of the three domains Archaea, Eubacteria, and Eukaryota as described in an early report of Woese et al. [20]. The structures of metabolic pathway deduced from Archaea are more similar to those from Eukaryota than to those from Eubacteria. This agrees with the rooted universal tree of life [33,34] and the tree of life [35,36]. The metabolic network structures of organisms reflect their functional relationship with the environment, and the similarity might provide a measure of the organism's physiological functions. The trajectory of an organism's adaptation can be explained using the structure of its metabolic contents. Our approach can be extended to more organisms and applied to other types of biomolecular interactions, such as physical protein interactions in regulatory networks, to provide a basis for understanding the functional relationships between biological networks in different organisms.

Methods
We attempt to cluster organisms by comparing sets of metabolic pathways. Our basic assumption is that different species exhibit overlapping components of metabolic pathways. To construct phylogenetic trees, the features of the organism-specific pathways are automatically extracted by considering the reference pathway. Here, the features represent the connection information between two enzymes. Figure 4 summarizes the procedure for the phylogenetic clustering of organisms by metabolic pathways using four steps:

Method
Similarity score Our method 0.196 [11] 0.154 The procedure for data processing for phylogenetic tree construction from the metabolic networks Figure 4 The procedure for data processing for phylogenetic tree construction from the metabolic networks.
Organism-specific enzyme-enzyme relation lists Step 1: Build the enzyme-enzyme relation lists.
-Step 2: Convert the lists to graph structures.
We evaluated this method using known experimental data on a collection of nine metabolic pathways from 81 representative organisms. Figure 5 shows the simple concept used to compute the distance between two metabolic networks. The resulting phylogenetic trees were cross-compared for consistency with existing methods to analyze phylogenies. The next section describes the datasets collected and the preparation methods, followed by the definition of graph kernels and their use in comparing and clustering metabolic networks for phylogenetic analysis.

Data preparation Dataset
We chose the KEGG database [37] as the resource for previous phylogenetic analysis. KEGG provides both an online map of pathways and the ability to focus on metabolic reactions in specific organisms. Each reaction may be uni-or bidirectional.

Enzyme-enzyme relation lists of organisms
The pathways provided in KEGG are visualized on manually drawn pathway maps or XML-based graphics. To construct enzyme-enzyme relation lists, we used information about chemical compounds and chemical reactions contained in the LIGAND database [38]. The LIGAND database provides detailed molecular information about one type of the generalized protein-protein interaction, namely, the enzyme-enzyme relation. LIGAND is a composite database of ENZYME and COMPOUND. The ENZYME section contains information about enzymatic reactions and enzyme molecules, and the COMPOUND section contains more than 6,000 chemical compounds. The enzyme-enzyme relationship can be extracted from information about enzymes contained in the COM-POUND entries. We automatically extracted information The simple concept for computing the similarity between two metabolic networks using the kernel method Figure 5 The simple concept for computing the similarity between two metabolic networks using the kernel method. about enzymes of a specific organism from the ENZYME section.
The enzyme-enzyme relationships of a specific organism were extracted from the enzyme-enzyme relation list. If two enzymes of r k in an enzyme-enzyme relation list existed in the enzyme list of a specific organism, we inserted r k into R'(O i ). P'(O i ) can be constructed by R'(O i ).

Data analysis Metabolic networks as labeled graphs
Our approach to estimate the distance between two metabolic networks is based on the graph comparison.
; ∑ a L i (l, a) = 0 otherwise. To compare the structures between metabolic networks of two organisms represented in graphs as described above, we adopted a kernel-based approach called the exponential graph kernel [39].

Comparison of metabolic networks: graph kernel
Given two graphs Γ i = (L i , H i ) and Γ j = (L j , H j ), the first simple approach to the graph comparison is to count the common vertices with the same labels in both Γ i and Γ j .
This similarity (or kernel) can be calculated by k(Γ i , Γ j ) = < , >, where the inner product <M i , M j > between two matrices of the same dimension is defined as Based on the definition of the label matrix in Equation where M i (l, l) = 1 only when f(v) = ᐍ l (∃v ∈ i ), and M i (l, l) = 0 otherwise. However, this approach considers only the presence or absence of vertices (enzymes) but does not consider the structure of the graph, such that the successive enzymes or reaction steps cannot be considered when comparing metabolic networks. To capture the structure of the graph, one must also consider vertices that can be reached from a vertex by a subsequent traverse.
In the exponential graph-kernel method, the similarity between two graphs Γ i and Γ j is defined as is a real-valued parameter and its value is chosen by performing many tries. When β = 0, it recovers the simple common vertex-counting measure since exp(0H) = I, the | | × | | identity matrix. Each element H n (a, b) of the matrix H n in Equation (5) represents the number of walks of length n (admitting cycles) from v a to v b , and allows the representation of the global structure of a graph.
Substituting Equation (5) into Equation (4), we can decompose the kernel function k(Γ i , Γ j ) into two meaningful parts, k(Γ i , Γ j ) = k 1 (Γ i , Γ j ) + k 2 (Γ i , Γ j ), where The kernel function k 1 contributes by considering walks of the same length in both graphs, and k 2 can take into account the insertion or deletion of vertices in the graph [39]. As the number of movements in a graph increases, the significance of walks of length n decreases by .
Eventually, the exponential matrix e βH can be interpreted as the product of a continuous process H, from which the identity matrix expands gradually to the matrix of the global structure of Γ [40].
The exponential graph kernel requires the exponentiation of square matrices Hs. This can be performed by matrix diagonalization, with time complexity of about O(| i | 3 ) for H i thus Γ i . [39]. The time complexity of the elementwise product of two matrices in k(Γ i , Γ j ) is O(max (| i | 2 , | j | 2 )). With N graphs, finally, the total time complexity for constructing the kernel matrix K = {k ij } (1 ≤ i, j ≤ N) is From the kernel k(Γ i , Γ j ), the dissimilarity metric is defined in the standard manner, that is, If we use the normalized kernel, then the distance metric is simplified as To summarize, metabolic networks constructed from reference pathways of N organisms were first converted to labeled undirected graphs. Each graph Γ i (1 ≤ i ≤ N) was then represented by two matrices: the vertex-label matrix L i and the adjacency matrix H i . Using these two matrices, we can take into account only the local structure (the direct connectivities between enzymes in pathways) of networks. To compare networks in terms of their global structure, we adopted a kernel-based method, which we named the exponential graph kernel. Finally, the distance matrix acquired from the kernel function was fed into a hierarchical clustering algorithm to construct the phylogenetic trees.

Constructing phylogenetic trees
The distance between two organisms was calculated by comparing their metabolic networks using the measures mentioned earlier. To construct a phylogenetic tree, we used an unweighted pair-group method with arithmetic mean (UPGMA) [41,42], a hierarchical agglomerative clustering algorithm. Given N organisms, the algorithm starts by initializing N clusters, each of which contains exactly one distinct organism, and proceeds by iteratively merging the two nearest clusters until only one cluster (called the root of the tree) remains. The dendrogram derived from UPGMA is a binary tree, which we consider may represent a binary phylogenetic tree.