An algorithm for linear metabolic pathway alignment

Abstract

Metabolic pathway alignment represents one of the most powerful tools for comparative analysis of metabolism. It involves recognition of metabolites common to a set of functionally-related metabolic pathways, interpretation of biological evolution processes and determination of alternative metabolic pathways. Moreover, it is of assistance in function prediction and metabolism modeling. Although research on genomic sequence alignment is extensive, the problem of aligning metabolic pathways has received less attention. We are motivated to develop an algorithm of metabolic pathway alignment to reveal the similarities between metabolic pathways. A new definition of the metabolic pathway is introduced. The algorithm has been implemented into the PathAligner system; its web-based interface is available at http://bibiserv.techfak.uni-bielefeld.de/pathaligner/.

Keywords: metabolic pathways, alignment, algorithm, PathAligner

Introduction

During the past decade bioinformatics has undergone rapid development. A vast amount of information on different organisms, both on genetic and metabolic levels, has been accumulated and systematically stored in specific databases that are available via the Internet [Collado-Vides and Hofestädt, 2002]. This rapid accumulation of biological data provides the possibility of studying metabolic pathways systematically. Analysis of metabolic pathways is an essential topic in understanding the relationship between genotype to phenotype. Metabolic pathway alignment is of importance to study biology evolution, pharmacological targets and other biotechnological applications [Dandekar et al., 1999].

Research on genomic sequence alignment has, so far, been intensively conducted [Needleman and Wunsch, 1970; Smith and Waterman, 1981]. Applications and tools, such as FASTA [Pearson and Lipman, 1988] (http://www.ebi.ac.uk/fasta3) and BLAST [Altschul et al., 1990; Altschul et al., 1997] (http://www.ncbi.nlm.nih.gov/BLAST) have been developed to further understand the biological homology and estimate evolutionary distance. Although the information provided by sequenced genomes can yield insights into their evolution and cellular metabolism, knowledge of the genome sequence alone is really only the starting point of the real work. The true functional units of metabolic systems are metabolic pathways. In fact, the emphasis of research efforts begins to turn back from gene sequences to cell functions as the completion of a huge number of genomes and the accumulation of knowledge about metabolism have made the comparison of metabolic pathways possible.

Several approaches to metabolic pathway alignment have been made in the past years. Forst and Schulten, 1999, 2001, extended the DNA sequence alignment methods to define distances between metabolic pathways by combining sequence information of involved genes. They combined sequence information of involved genes with information of the corresponding pathway. For each functional role of the pathway, all genes in the genomes that code for this functional role were used. Dandekar et al., 1999, compared glycolysis, Entner-Doudoroff pathway and pyruvate processing in 17 organisms based on the genomic and metabolic pathway data by aligning specific pathway-related enzyme-encoding genes on the genomes. Liao et al., 2002, developed a computational method to compare organisms based on whole metabolic pathway analysis. The presence and absence of metabolic pathways in organisms was presented as a Boolean vector. Based on this methodology and by using some specific distance measures on these profiles, pairwise comparison of a set of completed genomes were performed. Tohsato et al., 2000, proposed a multiple (local) alignment algorithm by utilizing information content that was extended to symbols having hierarchically structured EC numbers. They considered that reaction similarities can be expressed by the similarities between EC numbers of the respective enzymes and applied their method to pathway analyses of sugar, DNA and amino acid metabolisms. Heymans and Singh, 2003, presented a technique for the phylogenetic analysis of metabolic pathways based on the topology of the underlying graphs. A distance measure between graphs was defined using the similarity between nodes (enzymes) of the graphs and the structural relationship between them. The distance measure was applied to enzyme-enzyme relational graphs (two enzymes are related if they activate reactions which share at least one chemical compound) derived from metabolic pathways. Using this approach, pathways and groups of pathways of different organisms were compared to each other and the resulting distance matrix was used to obtain a phylogenetic tree. Kelley et al., 2003, devised a method called PathBLAST to align two protein interaction networks. The alignments were scored by the degree of protein sequence similarity at each pathway position and by the quality of the protein interactions they contain.

We are going to present a new pathway alignment strategy to analyze and fully characterize metabolic pathways in the cell. The alignment algorithm differs from classical string alignments as well as previous approaches. It considers functional similarities of enzymatic reactions involved in two metabolic pathways. This paper is structured as follows: the next section explains definitions of metabolic pathways. In the main section "Metabolic Pathway Alignment" section, firstly, theory basics are described. Then, similarity function, strip and index function are introduced. Finally, the alignment algorithm, its properties and its implementation are presented. A discussion about the algorithm is addressed in the form of conclusions.

Metabolic pathway definitions

Cellular metabolism is most often described and interpreted in terms of the enzyme-catalyzed biochemical reactions that make up the metabolic network. Traditionally, metabolic pathways have been explained in the context of their historical discovery, often named after key molecules (e. g. "glycolysis", "urea cycle", "pentose phosphate pathway", "citric acid cycle" and so on). The classification of molecules into pathways has historical reasons and is not explicitly based on qualitative differences of the interactions. One reason for this is that it is easier to think about the network as pathways which are densely connected.

So far, metabolic pathways are defined in literature [Mavrovouniotis, 1995; Voet and Voet, 1995; Hartwell, 1997; Selkov et al., 1998; Forst and Schulten, 1999; Schuster et al., 2000] in various ways with varying degrees of formality. Mavrovouniotis, 1995, defined a biochemical pathway as an abstraction of a subset of intricate networks in the soup of interacting biomolecules. A prevailing definition is that a metabolic pathway is a special case of a metabolic network with distinct start and end points, initial and terminal vertices, respectively, and a unique path between them, i. e. a directed reaction graph with substrates as vertices and arcs denoting enzymatic reactions [Forst and Schulten, 1999]. Typical metabolic pathways are given by the wall chart of Boehringer Mannheim [Michal, 1999] and KEGG [Kanehisa and Goto, 2000] (http://www.genome.ad.jp/kegg/metabolism.html), which have been verified with a number of printed and on-line sources. Some databases such as KEGG or WIT [Overbeek et al., 2000] (http://wit.mcs.anl.gov/WIT2/) represent metabolic pathway graphs with labeled arcs indicating the involved enzymes. Schuster et al., 2000, provided a general definition of metabolic pathways based on the concept of elementary flux modes (or minimal T-invariants in Petri nets [Voss et al., 2003]). It allows to test whether sets of enzymes form a consistent pathway allowing mass balancing for each intermediate and complying with the directionality of reactions (irreversibility). However, it represents modes under idealized situations without regulation and feedback. We should emphasize that the definition of a metabolic pathway is not exact, there are always interactions among pathways.

Moreover, the basic strategy to represent and compute pathways is the reactant-product binary relation. Properties of the pathway that rely upon the integration of two or more input molecules and unrelated output molecules and feedback effects are ignored (Figure 1). In that context, pathway refers to a more or less linear path from substrate to product, even though usually some of the molecules involved are also found in other pathways.

Figure 1: Abstract pathway with binary relations, compared to pathways in reality.

As a result, the traditional well-known pathways such as glycolysis or TCA will not be considered as a well-defined pathway, because they often contain some branches or alternative pathways. In fact there are several pathways inside them. We emphasize that a metabolic pathway is a subset of reactions that describe the biochemical conversion of a given reactant to its desired end product. In other words, several biochemical reactions act together in a pathway to transform a set of initial substrates into products with very different structures.

Let M = {m₁...m_n} be a set of metabolites in cells. Let e_i : M M be a function for enzymatic reactions taking place in the cells. For all m₁, m₂, m₃ M, the following property holds:

e₁(m₁) = m₂ and e₂(m₂) = m₃ e₂(e₁(m₁)) = m₃

Let e₁(m₁) = m₂, e₂(m₂) = m₃,..., e_k(m_k) = m_k+1, we define e₁e₂...e_k (m₁) = e_k(e_k-1...e₁(m₁)) = m_k+1. A new proposed definition of a metabolic pathway is presented and discussed in the following paragraphs.

Definition 1

Given e_i : M M, a metabolic pathway is defined as a subset of successive enzymatic reactions P = e₁e₂...e_k.

For each e_i (1 i k), there exist a pair of metabolites (m_i, m_i+1), m_i M, m_i+1 M, involved in the enzymatic reaction as substrate and product. P is a composition function of the reactions e_i, and we have P(m) = e_ke_k-1...e₁(m). Given an initial substrate m₁, the ending product m_k+1 is predictable. If m_k+1 = m_i (1 i k), i. e. the product of e_k is one of the metabolites involved in the previous steps, then the pathway is called as a cyclic pathway. Obviously, if m_k+1 m_i (1 i k), then it is a linear (non-cyclic) pathway.

Each enzymatic reaction e_i (1 i k) is catalyzed by a certain enzyme that is denoted as a unique EC number. The EC number is expressed with a 4-level hierarchical scheme that has been developed by the International Union of Biochemistry and Molecular Biology (IUBMB) [Webb, 1992]. The 4-digit EC number, d₁.d₂.d₃.d₄ represents a sub-sub-subclass indication of biochemical reaction. For instance, arginase is numbered by EC 3.5.3.1, which indicates that the enzyme is a hydrolase (EC 3.*.*.*), and acts on the "carbon-nitrogen bonds, other than peptide bonds" (sub-class EC 3.5.*.*) in linear amidines (sub-sub-class EC 3.5.3.*). These enzymes are normally separated-enzymes. Those enzymes may form a multienzyme complex (noncovalent aggregates of enzymes) or a membrane-bound system; a representative enzyme can be chosen unless there is a unique term for it. Thus, we can adapt the EC number as a unique name for the responding enzyme catalyzed reaction.

In brief, pathways are abstractions of sets of enzymatic reactions. They are substructures that are partitions of the metabolism.

Example 1
e_3.5.3.1 means the biochemical reaction that is catalyzed by the enzyme 3.5.3.1., which catalyzes arginine into urea.

and

indicates that a metabolic pathway e_2.1.3.3e_6.3.4.5e_4.3.2.1e_3.5.3.1 starts from the enzymatic reaction 2.1.3.3 with carbamoyl-P as reactant and after a series of reactions (2.1.3.3, 6.3.4.5, 3.5.3.1) results in urea as the product.

Metabolic pathway alignment

Theory basics

In order to score the similarity (identity percentage) between two metabolic pathways, we define the similarity function that is adapted from standard alignment definitions on strings.

Definition 2

Let E be a finite set of e functions, an edit operation is an ordered pair (α, β)  (E {ε}) (E {ε})\{(ε, ε)}.

α and β denote 4-digit EC strings of enzymatic reaction function, e. g. α = e_1.1.1.1 β = e_2.3.4.5, ε denotes the empty string for null function. However, if α ε and β ε, then the edit operation (α, β) is identified with a pair of enzymatic reaction function.

An edit operation (α, β) is written as α β (we simply write α, β as EC numbers). There are three kinds of edit operations:

α ε denotes the deletion of the enzymatic reaction function α

ε β denotes the insertion of the enzymatic reaction function β

α β denotes the replacement of the enzymatic reaction function α by the enzymatic reaction function β

notice that ε ε never happens.

Definition 3

Let E₁ = e₁e₂ ... e_m and E₂ = e'₁e'₂...e'_n be two metabolic pathways, an alignment of E₁ and E₂ is a pair sequence

(α₁ β₁,..., α_h β_h)

of edit operations such that E'₁ = α₁,..., α_h and E'₂ = β₁,..., β_h.

Note that the unique alignment of ε and ε is the empty alignment, that is, the empty sequence of edit operations. Empty element ε can be inserted at any position, i. e. at the beginning or end. An alignment is usually written by placing the EC numbers of the two aligned pathways on different lines.

Example 2
The alignment A = (2.4.2.3 2.4.2.4, 3.5.4.5 ε, 3.1.3.5 3.1.3.5, ε 2.7.4.9) of the pathways e_2.4.2.3e_3.5.4.5e_3.1.3.5 and e_2.4.2.4e_3.1.3.5e_2.7.4.9 is written as follows, one above the other:

Example 3
Five alignments of E₁ = e_2.7.4.14e_3.1.3.5e_3.5.4.5e_2.4.2.3 and E₂ =  e_2.7.4.14e_3.2.2.10e_3.5.4.1e _2.4.2.3e_3.5.4.5

Lemma 1

Let A = (α₁  β₁,..., α_h  β_h) be an alignment of E₁ = e₁e₂... e_m and E₂ = e'₁e'₂... e'_n. Then m + n  h  max{m, n}.

Proof
1.) The alignment

of E₁ and E₂ is of maximal length. Its length is m + n, hence m + n  h.

2.) Let m   n. Then

3.) The case m < n is similar to case 2. Hence h  max{m, n}.

To estimate the number of pathway alignments, we define that Aligns(m,n) is the number of alignments of one pathway E₁ of m EC numbers with another pathway E₂ of n EC numbers.

Lemma 2

For all m, n  0, Aligns(0,0) = 1, Aligns(m,0) =1 and Aligns(0,n) = 1, then

Aligns(m,n) = Aligns(m-1,n) + Aligns(m, n-1) + Aligns(m-1, n-1)

Proof
The idea is to focus on the end of the alignment. If e_m is deleted, then there exist Aligns(m-1, n) alignments of the earlier part of the pathway. If e'_n is deleted, then Aligns(m, n-1) alignments result. If e_m and e'_n are aligned, Aligns(m-1, n-1) alignments result. Therefore,

If not to count and as distinct, the new way of counting alignments is to identify aligned pairs and to ignore permutations of . The notation of index alignment is introduced.

Definition 4

An index alignment of E₁ = e₁e₂...e_m and E₂ = e'₁e'₂...e'_n is a set of index pairs, (i₁, j₁),(i₂, j₂),...,(i_r, j_r) satisfying: 1 i₁ < i₂ < ... < i_r m and 1 j₁ < j₂ < ... < j_r n.

For 1 h r, the index pair (i_h, j_h) stands for the replacement e_ih e'_jh. We say that e_ih is matched or aligned with e'_jh. All EC numbers in E₁ and E₂ not occurring in an index alignment are considered to be deleted in E₁ or E₂. In a graphical representation, the index pairs of the index alignment appear as lines connecting the EC numbers (Example 4).

Example 4
The following index alignment of E₁ = e_2.7.4.14e_3.1.3.5e_3.5.4.5e_2.4.2.3 and E₂ = e_2.7.4.14e_3.2.2.10e_3.5.4.1e _2.4.2.3e_3.5.4.5 represents the alignments of Example 3. In particular, P₁ represents A₁, A₂ and A₃, while P₂ represents A₄, and P₃ represents A₅.

Lemma 3

Let Indexaligns(m, n) be the number of index alignment of two fixed pathways of length m and n. Then

Proof
1.) For each r [0, min{m, n}] we have: for the ordered selection of the indices i₁,...,i_r there are possibilities; for the ordered selection of the indices j₁,...,j_r there are possibilities.

All these possibilities have to be combined: .

2.) The key to the last equality is to consider the binomial expansion (x + y)n. For details see Appendix 1.

Obviously, Indexaligns(m, n) is a special case of Aligns(m, n), and it is possible to further reduce the number of alignments by requiring conditional matches. There are at least different alignments between E₁ and E₂.

Similarity function

The notion of alignment requires some scoring or optimization criterion. A scoring scheme must account for replacements, insertions and deletions. Scores are measures of sequence similarity (similar sequences have high scores), which is given by a similarity function. The similarity function is measured by the following definition:

Definition 5

A similarity function σ assigns to each edit operation (α β) a nonnegative real number. The similarity σ (α ε) and σ (ε β) of the deletion operation (α ε) and insertion operation (ε β) is 0. For all replacement operations (α β) α ε, β ε, say, α = d₁.d₂.d₃.d₄ and β = d₁'.d₂'.d₃'.d₄', then the similarity function σ (α β) is defined by:

The definition is a heuristic device to measure the difference of enzymatic reactions that are catalyzed by known (classified) enzymes. It does not exclude the possibility that d₄, d₃.d₄, and d₂.d₃.d₄ can be respectively expressed as wide card symbols *, *.* and *.*.* which means no clear classification of the enzyme.

However, single pair of EC string comparison is simply measuring the difference between EC strings. Often it is additionally of interest to analyze the total difference between two strings into a collection of individual elementary differences. The most important mode of such analyses is an alignment of the pathways. The function σ can be extended to alignments in a straightforward way: the similarity σ (A) of an alignment A = (α₁  β₁,..., α_h  β_h) is the sum of the similarities of the edit operations A consists of.

Example 5
The similarity of the alignment A₅ in Example 3 is:

σ (A5) =	σ (2.7.4.14  2.7.4.14) + σ (3.1.3.5  3.2.2.10) + σ (3.5.4.5  3.5.4.1) + σ (2.4.2.3  2.4.2.3) + σ (ε  3.5.4.5)
=	1 + 0.25 + 0.75 + 1 + 0
=	3.0

When considering the lengths of pathways, an alignment scoring scheme is given.

Definition 6

An alignment scoring scheme, Score(E₁, E₂) of two metabolic pathways is the minimal mean similarity of their alignment

Lemma 4

0 Score(E₁, E₂) 1. Score(E₁, E₂) = 1 if and only if E₁ = E₂.

Proof
The worst case is that no index pair exists in the alignment, i. e. all edit operations are deletions and/or insertions. Hence similarity σ (A) is zero. The best case is that the similarities of all edit operations are 1, i. e. α = β, there is neither deletion nor insertion. Hence similarity σ (A) is n, and E₁ and E₂ share the same length, m = n. Therefore

The scoring scheme defined by Definition 6 is generic. Algorithms for optimal alignment can seek either to minimize a dissimilarity measure or maximize a scoring function. In order to achieve a maximum possible score of the alignment, the edit distance could be adapted to measure the similarity between two pathways by calculating the minimal cost of the edit operations [Levenshtein, 1966; Wagner and Fisher, 1974; Sellers, 1980]. However, our scheme is based on the index pair matching with strip functions (Definition 7) under the following consideration.

A classical dynamic alignment algorithm is to achieve the optimal matching by minimizing the edit distance. However, when taking the biological aspects of metabolic pathways into account, especially when considering that two evolutionarily related pathways are diverged for some certain biological purpose, the alignment with the edit distance is arbitrary. For example, the same metabolic pathway from two organisms may have diverged since the organisms evolved from a common ancestor, and individual metabolites and enzymes may have been changed or added or lost in one pathway. Two theories exist. The "retrograde evolution" theory [Horowitz, 1945] states that sequential disappearance of key intermediary metabolites induces the recruitment of similar available substrates via new enzymes. The "substrate ambiguity" theory [Jensen, 1976] indicates enzyme recruitment from a pool of ancestral enzymes with basic functions and substrate ambiguity. Therefore, metabolic pathways exist for certain cellular functions that transform substrates into (intermediate) products. Based on these considerations, we restrict the alignment algorithm and define a new function to perform the removal of unmatched elements from both ends of the pathway.

Strip and index function

Also set

when e_i e'_j (1 i m and 1 j n) or e₁e₂...e_m = Ψ or e'₁e'₂...e'_n = Ψ or e₁e₂...e_m = e'₁e'₂...e'_n = Ψ. The result of δ (E₁, E₂) can be further performed with δ until δ (E₁, E₂) = (E₁, E₂) (Figure 2).

Figure 2: A schematic alignment with strip function.

The set of all matched points will be denoted as λ* (E₁, E₂). Such that

where λ (δ ⁰ (E₁, E₂))= λ (E₁, E₂); 1 i < i' <...< k' < k m, 1 j < j' < ... < l' < l n. According to Definition 4, we know that (i, j),(i', j'),...(k', l'),(k, l) is an index alignment of E₁ and E₂.

Example 6
Let two pathways E₁ = e_2.7.7.6e_3.6.1.5e_2.7.4.14e_3.1.3.5e_3.5.4.5e_2.4.2.3e_1.3.1.1e_3.5.2.2e_3.5.1.6 and  E₂ = e_2.7.7.6e_2.7.4.6e_2.7.4.14e_3.2.2.10e_3.5.4.1e_1.3.1.2e_3.5.2.2e_3.5.1.6, then

Lemma 5

Let Indexaligns(δ ^r(E₁, E₂)) be the number of index alignment of two fixed pathways of length m and n with r matched points. Then

Proof
See Appendix 2.

Observation
Let E₁ and E₂ be two pathways of length m and n, there exists a minimum r that enables , then .

Algorithms and implementation

Before describing our algorithms, we introduce four different subscript notations of the Strip function δ and the Index function λ. From its definition we know that the Strip function δ is able to strip two pathways if there are two pairs of elements that are matched, e. g. e_i = e'_j and e_k = e'_l. Here, e_i, e'_j, e_k and e'_l are 4-hierarchical numbers. We count all four numbers by default. Now if we count only the first three numbers for matching, then the Strip function δ can be written as δ₃ in order to distinguish the default setting δ that can also be written as δ₄ in this case. The same goes for δ₂ and δ₁. Accordingly, λ₄ is defined as the default 4-number indexing, and λ₃, and so on.

Pairwise alignment

In general, we align metabolic pathways one above the other. The alignment algorithm is based on likelihood calculations of index pairs. Given two metabolic pathways P₁ and P₂ (Figure 3), the implemented algorithm is given by Figure 4 and the following pseudo-code:

Figure 3: Two aligned metabolic pathways.

Figure 4: Overview of algorithm used to align metabolic pathways based on a 4-step strip function.

Line 1-3: Initialize the set of unaligned EC number sequences, their lengths and score value.

Line 4-29: Starting from both ends towards the middle, align one sequence to another and attempt to find all EC numbers with same 4-level hierarchical numbers. Score the similarities. Recall the alignment positions where EC numbers are identical and cut the sequences into more sub-sequences by removing the identical EC numbers.

Line 7-27: Each pair of sub-sequences is initialized to begin a new round of 3-level hierarchical EC number matching until all pairs of sub-sequences are aligned. A similarity score is calculated afterwards.

Line 12-25: Applying the same rule again, find the similarities of rest unaligned sub-sub-sequences based on 2-level hierarchical EC number matching.

Line 17-23: Then sub-sub-sub-sequences on 1-level matching if any.

Properties

Definition 8

Given two pathways E₁ and E₂, a mapping M(E₁, E₂) is defined as a set of position correspondences (i, j) satisfying 1 i m and 1 j n such that e_i :::: e'_j. The notation "::::" denotes that e_i and e_j' are compared to be identical in turn according to their 4-digit hierarchical patterns. A mapping is maximal if no other pair (l, k) exists such that e_l :::: e'_k.

Obviously, each mapping site has two characteristics, site position (location) i and site feature name e_i. The map E₁ = e₁e₂...e_m is defined as a sequence of pairs e_i = (a_i, r_i), where a_i = the location of the i-th site in number of pathways and r_i = the feature name at the i-th site. Similarly, E₂ = e'₁e'₂...e'_n is a map where e'_j = (b_j, s_j). For further analysis, let the symbol denote a natural order sequence of a_i, and denote the sequence of corresponding position b_j. The number of E₁, E₂ maps is defined as and . Clearly, = m, n.

Example 7
The maximal mapping of two pathways E₁= e_3.6.1.5e_2.7.4.14e_3.1.3.5e_2.7.1.48 and E₂ = e_3.6.1.8e_3.1.3.5e_2.7.1.48e_2.7.4.14 is {(1,1),(2,4),(3,2),(4,3)}. = 1, 2, 3, 4 and = 1, 4, 2, 3. The mapping can be represented as line over line Suppose a mapping (Figure 5), the following properties hold:

1. = = d,

2. a_i1 < a_i2 <...< a_id-1 < a_id

3. b_j1 < b_j2 <...< b_jd-1 < b_jd is not always true.

Figure 5: An illustrative map of two pathways.

As emphasized above, locations of is not necessary in natural order. We define a map alignment as a mapping where is a sequence in natural order.

Lemma 6

Given the maximal mapping M of two pathways E₁ = e₁e₂...e_m and E₂ = e'₁e'₂...e'_n, and let = a_i1a_i2...a_id-1a_id and = b_j1b_j2...b_jd-1b_jd be two maps of M. Then we have:

1.) A sub-sequence of in natural order with the longest length is the maximal map alignment of E₁ and E₂.

2.) A sub-sequence of in natural order between b_j1 and b_jd is the index alignment.

Proof
1.) Suppose b_jk0b_jk1...b_jkt is a sub-sequence of in natural order with the longest length, we can obtain the responding positions of this sub-sequence: a_ik0a_ik1...a_ikt, so that they are a map alignment of E₁ and E₂. For contradiction, let us assume that there exists another mapping pair (a'_ik, b'_jk) excluded from the maximal map alignment, then b_jk0b_jk1...b_jkt is not the longest subsequence, which is a contradiction.

2.) According to the definition of index alignment, the first step is to map from both ends of E₁, b_j1 and b_jd are found. Next step is to map a_i2 and a_id-1, b_j2 is one map of only when b_j2 is greater than b_j1, and b_jd-1 is one map only when b_jd-1 is less than b_jd and greater than b_j2. Repeat the mapping till no more b_jk is satisfied.

Example 8
The map of the two pathways in Example 6 is:

The maximal map alignment

The index alignment is

Proof
Given two pathways E₁ = e_3.6.1.5e_2.7.4.14e_3.1.3.5e_2.7.1.48 and E₂ = e_3.6.1.8e_3.1.3.5e_2.7.1.48e_2.7.4.14, then

Score(E₁, E₂) =  ¼*σ(2.7.4.142.7.4.14) + σ(3.6.1.53.6.1.8) = ¼*(1 + 0.75) = 0.43

While the maximum alignment will be

The maximal matching is

Time complexity analysis

The best case to strip two pathways and to obtain their index pairs is that they are identical, which costs O(m) computing time. While the worst case is that they are not matchable, it needs O(mn) to cover all elements. In general, we define O(m n) to show the "average" time complexity. Therefore, to create λ₄*(E₁, E₂), it takes O(m n). λ₃*(E₁', E₂') takes , λ²*(E₁", E₂") takes and λ₁*(E₁"', E₂"') takes , where K denotes numbers of index alignment of (E₁, E₂), R_l denotes each stripped sub-pathway of (E₁', E₂'), S_i denotes each stripped sub-sub-pathway of (E₁", E₂"), and _jm_i^l and _jn_i^l represent the length of each stripped sub-sub-sub-pathway (E₁"', E₂"').
Therefore, the total time complexity is .

The best-case complexity of the algorithm is the minimum number of steps taken on any instance of size m and n. It represents the alignment of two identical pathways, which takes O(n).

The worst-case complexity of the algorithm is the maximum number of steps taken on any instance of size m and n. It represents the alignment of two pathways with no first 2-number is same, which takes at least O(3mn).

Multiple alignment

The multiple metabolic pathway alignment allows us to extract and represent biologically important but faintly / widely dispersed pathway similarities which, for instance, makes it possible to identify pathways preserved by evolution that play an important role in the cellular function and can give us hints about the evolutionary history of certain pathways.

By allowing the alignment of more than two metabolic pathways, the pairwise alignment algorithm can be extended. A multiple alignment of metabolic pathways E₁, E₂, ..., E_k, can be seen as a generalization of pairwise metabolic pathway alignment - instead of aligning two pathways, k pathways are aligned simultaneously, where k is any number greater than two. A heuristic algorithm is used to perform the multiple metabolic pathway alignment. The general idea of the method is to construct a succession of pairwise pathway alignments:

Step 1: Choose one pathway E₁ from E = {E₁, E₂,...,E_k} and align with {E₂,...,E_k} one after another to find the most similar pathway E_i (2 i k).

Step 2: Choose the pathway E_i and align with E\{E₁, E_i} to get the most similar pathway E'_i.

Step 3: Iterate step 2 until all pathways are aligned.

The time complexity of multiple alignment of k pathways is , where p  (1, k-1), the lengths of pathway p and p + 1 are ^pm and ^p+1n.

Implementation

The algorithm has been implemented in the PathAligner system [Chen and Hofestädt, 2004] (http://bibiserv.techfak.uni-bielefeld.de/pathaligner). Three web-based alignment interfaces are implemented: "E-E Alignment", "M-E-M Alignment" and "Multiple Alignment". "E-E Alignment" uses the basic algorithm to align two linear metabolic pathways (represented as EC number sequences). The user can also align any such a metabolic pathway against our pool database to find a list of hits. "M-E-M Alignment" considers the differences of metabolites in two pathways, which are presented as "Metabolite-EC number-Metabolite" patterns of sequence. It is possible to pick up two such pathways and align them to identify whether they are alternative pathways. "Multiple Alignment" allows the alignment of more than two metabolic pathways.

Conclusion

Identification and analysis of metabolic networks is a complex task due to the complexity of the metabolic system. Abstract pathway defined as a linear reaction sequence is practical for our alignment algorithm. We have presented an algorithm to study the problem of metabolic pathway alignment. Our algorithm calculates the hierarchical similarities of EC numbers mapping from both ends of the sequences. The entire processing for pairwise alignment takes time , where m, n are the lengths of two (stripped) metabolic pathways. The algorithm described here has been successfully implemented into the PathAligner system.

Our algorithm calculates the similarity of reactions; it can be adopted to comparatively analyze other kinds of biopathways, such as signaling pathways. Although there is no nomenclature system for signal transductions at the moment, Chen et al., 2004, presented a classification of signal transduction events which can be used for this purpose. Since the ST number (each signal transduction is assigned a unique ST number) is also a 4-level hierarchical structure, it makes it possible to exploit our metabolic pathway alignment algorithm to perform signaling pathway alignment.

Currently a metabolic pathway is practically defined as a linear molecule sequence. However, when the topology of network is concerned, more information related to the components of the pathways like their length, size, number of feedback cycles, number of crosstalks between pathways and area reachable from any point in the network should be considered as much as possible. In this case, the metabolic pathway comparison will rather be the comparison of a network instead of a single linear pathway alignment. As a result, this leads to another type of pathway computation, which can be categorized as the comparison of biological networks. By comparing such types of networks to different biological systems, it is possible to identify similarities and variations among different species.

Acknowledgements

The work was supported by the German Research Foundation (DFG) graduate program "Bioinformatics" in the Uni-Bielefeld.