Multiple sequence alignment system and method

[0002] The present application claims priority to Provisional U.S. Patent Application No. 60/062,306, filed Oct. 17, 1997 by A. Floratos and I. Rigoutsos, herein incorporated by reference in its entirety.

[0003] 1. Field of the Invention

[0004] The present invention relates to the field of aligning data sequences, and, more specifically, to a method and system for aligning multiple data sequences to achieve the best commonality of the sequences.

[0005] 2. Discussion of the Prior Art

[0006] The problem of aligning a large number of sequences appears very often and in many different application areas (e.g. biology, data mining in databases, computer security etc.). It has been observed that most traditional methods of aligning multiple sequences do not work well on large number of sequences on the order of N>10.

[0007] The multiple sequence alignment problem has been studied for at least the last fifteen years, in the context of Computational Biology, it is well known that this is a rather difficult problem. Contributing to the difficulty in solving the problem is that in its most general form, it is difficult to model to the satisfaction of Evolutionary Biologists, Geneticists and other users. The most popular and successful approach to date has been a dynamic programming technique using different mechanisms of scores that is a function of the edit distance along with gap penalties to evaluate the similarity of the different sequences. Dynamic programming relies on identifying good penalty scores for matches and mis-matches. This is difficult to realize in real-world applications. The method is best suited to small sequence numbers on the order of N<6.

[0008] For the situation where there are more than two sequences to be aligned, one of the prior art approaches has been to perform a pairwise alignment where two of the N sequences are analyzed at a time and an N-wise alignment is built from the pairwise alignments. This approach works well for small values of N on the order of N<6, however, for large values of N additional constraints are required to give meaningful alignments.

[0009] It is apparent that the problem of aligning multiple sequences where the number of sequences is large is computationally a very demanding one.

[0010] The problems stated above and the related problems of the prior art are solved with the method and system according to the present invention. A method is provided that is specifically suited to aligning a large number of sequences, on the order of N>10, simultaneously, that brings out the best commonality of the N-sequences. The method is built on top of an existing method for identifying or discovering motifs in sequences. A motif can be informally defined as a repeated pattern that appears in two or more input sequences. The method uses some, if not all, of these identified motifs to give a good alignment. The present invention also provides an additional constraint, K, as a user defined control parameter. The additional parameter constrains the alignment of the N sequences to have at least K of the N sequences agree on a character, whenever possible, in the alignment. The alignment number, K, provides a natural constraint for dealing with a large number of sequences in that a commonality across most, if not all sequences is required to be detected.

[0011] Given a set of N input character sequences, the present invention aligns the N sequences in a way that brings out the best commonality of the N character sequences. At a motif discovery stage, the method of the present invention first identifies any motifs which may exist in the input sequences. At a pruning and alignment stage, the method detects one or more infeasible sets from the motifs identified at the motif discovery stage. Each infeasible set is comprised of a subset of those identified motifs which preclude an alignment which brings out the best commonality. Once identified, one or more motifs from the infeasible sets are removed (“pruned”), from consideration when aligning the sequences at an alignment stage. The remaining motifs comprise a feasible set which is then used to align or render the sequences to maximize a pre-specified cost function.

[0012] The sequence alignment method of the present invention is advantageous for the following reasons; 1) The approach is efficient, particularly on large numbers of long sequences, 2) a direct N-wise alignment is achieved by not composing the alignments from lower order (pairwise) alignments as practiced in the prior art, 3) the approach is adaptable to underlying biological or other models, using the idea of in-exact motif, 4) it allows for easy user intervention and specification in terms of which motif the user must have in the alignments, 5) it gives various alternate solutions, if required, to the user based on different configurations that the user wants to see rather than give a cost value which the user may or may not be able to relate to, 6) the resulting alignment is independent of the order of the input sequences.

[0013] The invention and its operation will become more apparent from the following description of an exemplary embodiment and from the accompanying drawings.

[0014] FIGS. 1a&bshow a directed graph and associated input sequence and corresponding identified motifs from which the graph is derived.

[0015] FIG. 2 is a directed graph that illustrates an infeasible set of the second type,closed cycle.

[0016] FIG. 3 is a directed graph that illustrates an infeasible set of the third type, closed path inconsistency.

[0017] FIG. 4 is a flow diagram for a two-path comparison to identify a closed path inconsistency.

[0018] FIG. 5 Experimental results that illustrate the method of the present invention as applied to biological data.

[0019] The subject invention will be described with reference to numerous details set forth below and the accompanying drawings which illustrate the invention. The following description is illustrative of the invention and is not to be construed as limiting the invention. Numerous specific details are described to provide a thorough understanding of the present invention. It will be apparent, however, to one of ordinary skill in the art that the present invention may be practiced without these specific details. In other instances, certain well known or conventional details are not described to preclude obscuring the present invention in detail.

[0020] The method of the present invention may be implemented as a computer program readable by a general purpose computer for executing computer instructions directed to aligning multiple sequences to bring out the best commonality of the sequences.

[0021] To better appreciate the sequence alignment method of the present invention, an illustrative example is provided:

[0022] Generally, a sequence is a string of characters taken from a finite character set or alphabet. As an arbitrary example, the following sequences could have arose from the English alphabet.

[0023] sequence 1: W A B C G Z T V A Q C B H N G Z T

[0024] sequence 2: A P C G Z E A T C T V B H N W A B C Z A Y

[0025] sequence 3: U I W D V P L H R Q Z S D

[0026] A pattern is a specific set of characters of any length from an alphabet where the pattern has a particular spatial arrangement. An example of such a pattern from sequence 1 above would be “G Z T”, or from sequence 2, “V B H”. Any combination of consecutive letters from any of the sequences would constitute a valid pattern. Patterns may also contain “don't care” characters such as “A . C”, where the dots are used to indicate that the respective position could be occupied by “any” character (don't care character). Any character which is not a don't care character, such as “A” and “C” in the pattern, is referred to as a solid character.

[0027] As previously stated, the present invention is partially directed to the process of motif discovery or identification. As defined in the prior art, a motif is a repeated pattern that appears more than once in a sequence. For the example provided, the pattern “G Z T” qualifies as a motif in sequence 1 above under the prior art definition. In the context of the present application, however, a motif is defined as a pattern that appears in two or more sequences. Under the present definition the pattern “A . C” qualifies as a motif in the example above because it appears in at least two sequences, the first and the second.

[0028] The sequence alignment method of the present invention solves the following problem.

[0029] Problem:

[0030] Given a set N of sequences each of length n₁, n₂, . . . , n_Non any alphabet set, and an alignment number K, where ((2<=K<=N), as input, align the N sequences in such a way that brings out the best commonality of the N sequences, with the possible introduction of gaps.

[0031] The problem formulation as described can be more readily understood by considering each sequence as comprising a row of a matrix. Filler characters are provided wherever necessary to make row lengths equal. What is sought is an alignment of the matrix rows that maximizes a pre-specified cost function. In an exemplary embodiment, the cost function is defined as the number of columns in which a common character appears at least k times. Other cost functions may be realized by the method.

[0032] The sequence alignment method of the present invention proceeds in two stages: 1. motif discovery and 2. pruning and alignment. In the motif discovery stage, the set of input sequences, N, are searched to identify or discover any motifs which may exist in the N input sequences. During the pruning and sequence alignment stage a determination is made regarding whether all or some subset of the motifs discovered in the first stage can be used to align or render the N sequences. More particularly, it must be determined whether one or more identified motifs precludes bringing out the best commonality among the N sequences. If so, one or more offending motifs are pruned or removed from consideration when performing the alignment. Otherwise, if all the motifs are usable then an alignment is formed that utilizes or respects all the identified motifs.

[0033] 1. Motif Discovery

[0034] The method of motif discovery used in the present invention is described in patent application Ser. No. 60/062,306, herein incorporated by reference in its entirety. During the motif discovery stage, motifs are identified or discovered from the N input sequences. A motif is defined in the present context as a pattern that appears in two or more sequences. In most instances a very large number of motifs will be identified. The present invention can accommodate this result at the expense of processing time. To alleviate this potential bottleneck, however, the method, in an alternate embodiment, may also search for only irredundant motifs. Irredundant motifs are a special sub-class of motifs, whose number is only quadratic in input size. A polynomial time algorithm may then be used to extract the irredundant motifs from the input thereby bounding the extraction time.

[0035] 2. Pruning and Sequence Alignment

[0036] A pruning and sequence alignment stage follows the motif discovery stage. At this stage any “offending” motifs are identified and removed from consideration during sequence alignment. Offending motifs are defined as those motifs which preclude performing a sequence alignment which achieves the best commonality among the N sequences. Offending motifs are identified by the method as a consequence of performing three infeasibility tests. The three tests are; 1. a pairwise compatibility test, 2. a smallest cycles test, and 3. a closed path with inconsistencies test. Each will be described in detail. For each of the three tests, one or more offending motifs may be identified. Each group of identified offending motifs is said to form an infeasible set. Therefore, the infeasible sets, when found, comprise one or more offending motifs which violate the respective infeasibility test.

[0037] 2.1 Pairwise Compatibility Test (Determination of the 1st Infeasible Set)

[0038] The first infeasible set is comprised of those offending motifs which fail the test for pairwise compatibility. Those offending motifs are identified in pairs and are therefore defined as pairwise incompatible or forbidden pairs. The pairwise compatibility test will be performed for only those motifs, identified at stage 1, which occur simultaneously in at least one input sequence. The pairwise compatibility test is comprised of two sub-tests 1) a domain crossing mismatch test, and 2) an overlap mismatch test. Each test will be explained by way of example.

[0039] With specific reference to the first test, domain crossing mismatch, assume the following sequences are input to the process;

[0040] Sequence 1: “A B C D E”

[0041] Sequence 2: “D E A B C”

[0042] Two motifs are identified from sequences one and two, they are “A B C”, and “D E”. By virtue of the relative positions of the two motifs in the respective sequences it is apparent that there is no feasible sequence alignment of sequence one and two that will permit the simultaneous alignment of the two motifs. This incompatibilty arises as a consequence of the juxtaposition of the two motifs in the respective sequences. As a consequence, the motif pair fails the domain crossing mismatch test. Failing the first sub-test of the two-part pairwise incompatibility test is sufficient to declare the motif pair incompatible and removes the necessity of performing the second sub-test.

[0043] The second test, overlap mismatch, is also illustrated by example. Given the two sequences,

[0044] Sequence 1: “A B C E E”

[0045] Sequence 2: “A C E E F”,

[0046] two motifs are identified. They are “A . . E”, where the dots indicate don't care characters and “C E E”. It is apparent that no alignment of the two sequences will respect the simultaneous alignment of the two motifs because alignment of one motif requires a linear translation of one sequence such that alignment of the other motif is prevented. When that situation occurs an overlap mismatch between the two motifs exists. From the pairwise incompatibility test a first infeasible set is generated, assuming there is at least one pair of motifs which violate the test.

[0047] A directed graph is then constructed. Directed graphs are well known in the art, example of such graphs may be found in chapter 23 of Introduction to Algorithms by T. H. Cormen, C. E. Leiserson, and R. L. Rivest, herein incorporated by reference. The vertices of the directed graph correspond to the motifs discovered at the first stage, motif discovery. That is, each vertex in the graph will correspond to one of the discovered motifs. Further, a directed edge will connect two vertices (motifs) in the graph only when there exists at least one input sequence in which the corresponding motifs occur together. Stated in the alternative, if two motifs do not appear together in any of the N input sequences, no directed edge will connect their associated vertices in the graph. Directed edges are defined in the graph by an arrow which points to one of the two vertices it connects. For example, for motifs A and B, their connectivity in the graph may be represented as;

[0048] Vertex A (Motif A)→Vertex B (Motif B)

[0049] The directionality of the directed edge indicates that motif A is to the left of motif B in the input sequence list. For those motif pairs which fail the pairwise compatibility test, however, no left/right designation is possible, for reasons stated above, and as a consequence the directed edge will be labeled as a “forbidden pair” with no arrow designation. For all other motif pairs in the graph having an associated directed edge, that directed edge will be labeled as either overlap or non-overlap. An edge will be labeled overlap whenever any instance of the two motifs overlap in any one of the N input sequences. For example, motifs “A B . D” and “C D E” overlap in the sequence 1 below but not in sequence 2.

[0050] Sequence 1=“A B C D E”

[0051] Sequence 2=“A B X D Y G H C D E”

[0052] Note that irrespective of whether the directed edge is labeled overlap or non-overlap, the motifs associated with that edge are pairwise compatible. Only those edges labeled “forbidden pair” have associated motifs which are pairwise incompatible.

[0053] Referring now to FIGS. 1aand 1b, a directed graph 7 is illustrated in FIG. 1a. The graph 7 is constructed from an input sequence list 8, illustrated in FIG. 1b, having three input sequences 8, N=3, and five identified motifs 9.

[0054] Each identified motif 9 corresponds to a respective vertex in the directed graph 7 as illustrated in FIG. 1a. Directed edges DE 20-38 are illustrated and serve to connect the identified motifs 9. The arrow 13 associated with each directed edge indicates that the motif associated with the vertex being pointed to by the arrow 13 occurs to the right of the motif whose vertex is located on the opposite side of the directed edge. For example, directed edge 20 with associated arrow 13 indicates that the directionality of the edge is from vertex 10 to 12. The directionality specifies that motif (vertex) 10 is located to the left of motif (vertex) 12 in the sequence list 8. In addition to the directionality associated with each edge, additional labeling is illustrated. Most edges in the graph have a dual identification with the first identifier characterizing the respective motifs as “forbidden pair”, “overlap”, or “non-overlap”, and a second identifier describing a distance metric, shown in brackets 15 in the directed graph 7. The distance metric is a positive integer which is a measure of the minimum separation distance in the input sequence list 8 between the two motifs (vertices) connected by the directed edge. The separation distance is measured from the first character in each measured motif string. The computation of the distance metric for a directed edge is best illustrated by example.

[0055] Assume, for example, a directed edge connects two vertices (motifs); “A B C” and “D E”. Further assume that the two motifs were identified from the following input sequences;

[0056] As illustrated, the distance, measured from the first character in each string, between the first motif, “A B C”, and the second motif, “D E”, is +5 in the first sequence and +6 in the second sequence. The distance metric that would be assigned to the directed edge in the graph would be the minimum of the two distances , (i.e. +5=min(+5,+6)), in the example shown.

[0057] Directed edge 20 in the graph 7 is exemplary of the dual labeling associated with a directed edge not labeled as a “forbidden pair”. Directed edge 20 is labeled “overlap (2)” indicating that the two motifs overlap each other in the sequence list 8 with a distance metric measure of +2. A distance metric will be computed for each directed edge in the graph not labeled “forbidden pair”. A distance metric is neither required, nor logically assignable to edges labeled “forbidden pair”. Directed edges DE 22, DE 24, and DE 28-38 are each assigned a non-overlap label, with only edge 36 being explicitly illustrated for purposes of clarity of presentation.

[0058] The following pseudo-code generally describes the construction and labeling of directed edges between vertices (motifs) in the directed graph.

[0059] For every motif p_i, do

[0060] For every motif p_j, not equal to p_i, do

[0061] Does there exist a sequence in which both p_iand p_joccur

[0062] No: No directed edge exists between p_iand p_j

[0063] Yes: Are motifs p_iand p_jpairwise compatible?

[0064] No: Assign label forbidden to edge v_iv_j

[0065] Yes: Do motifs p_iand p_joverlap?

[0066] Yes: Assign label overlap to edge v_iv_j

[0067] No: Assign label non-overlap to edge v_iv_j

[0068] Once a directed graph is constructed the method performs the second and third tests for infeasiblity. The second and third tests for infeasibility are distinguishable from the first test in that the first test is directed to identifying offending pairs of motifs. The second and third tests are both directed to identifying groups of three or more offending motifs.

[0069] 2.2 Smallest Cycles Test (Determination of the 2nd Infeasible Set)

[0070] The second test for infeasibility will now be explained. The second test is directed to finding smallest cycles in the graph. A smallest cycle constitutes a closed loop in the graph wherein a closed loop is defined as any path that originates at some starting vertex and follows a path through at least two additional distinct vertices and returns to the starting vertex. Closed loops are distinguishable from closed paths which will be defined with reference to the third test for infeasibility, closed path inconsistencies.

[0071] Referring now to FIG. 2, an example of a cycle (closed loop) is illustrated. In FIG. 2, a directed graph is shown comprising five vertices, V140, V242, V344, V446, and V548 and six directed edges DE 50, DE 52, DE 54, DE 56, DE 58, and DE 60. Each vertex corresponds to a motif identified at stage one, motif discovery. A cycle is illustrated in FIG. 2 by following the path defined by directed edges DE 50, DE 52, and DE 54. The directed edges lead from a starting vertex, V242 to one or more intermediate vertices, defined by V344 and V446 in FIG. 2, and back to the start vertex V242. By contrast, the path defined by directed edges DE 60, DE 58, and DE 56 do not define a cycle or closed loop because starting at any node, the directed edges do not eventually lead back to a start node. If, however, in the directionality of the directed edge DE 56 were reversed then a closed cycle would exist for edges DE 60, DE 58 and DE 56.

[0072] The method identifies any potential closed loops or cycles in the graph by performing a depth first search. Depth first searches are well known in the art and the method steps for performing the search is well documented in the literature, an example of which can be found in chapter 23 of Introduction to Algorithms, herein incorporated by reference in its entirety.

[0073] As the depth first search is being performed, any offending motifs which violate the second test for infeasibility (closed loops) are identified and grouped into an infeasible set of the second type (closed loops).

[0074] 2.3 Closed Path with Inconsistencies Test (Determination of the 3rd Infeasible Set)

[0075] The third test for infeasibility will now be explained. The third test is directed to finding all closed paths with inconsistencies in the graph. A closed path may be generally defined as two or more distinct paths having a common origin and destination vertex. An inconsistent closed path may then be defined as an identified closed path which violates the closed path infeasibility test. The details associated with the infeasibility test for closed paths will be described in detail later.

[0076] Referring to FIG. 3, a directed graph illustrating an inconsistent closed path is shown. The graph is comprised of four vertices; V162, V264, V366, and V468 and four directed edges DE 70, DE 72, DE 74, and DE 76. The graph illustrates two paths P1 and P2 starting from vertex V162 and terminating at vertex V366. Path P1 is defined by directed edges DE 70 and DE 72, and path P2 is defined by edges DE 76 and DE 74.

[0077] The present invention identifies all closed paths in the graph by performing a breadth first search. A breadth first search is performed for each vertex in the graph. That is, for the four vertices illustrated in FIG. 3, four separate breadth first searches will be performed utilizing one of the four vertices as the start vertex in each iteration. Breadth first searches are well known in the art and the method steps for performing the search is well documented in the literature, an example of which can be found in chapter 23 of Introduction to Algorithms, herein incorporated by reference in its entirety. The objective of the algorithm, simply stated, is given a graph and a distinguished source vertex, a breadth first search systematically explores the edges of the graph to discover every vertex that is reachable from the source vertex.

[0078] As the breadth first search is being performed to identify closed paths, the third test for infeasibility (closed path inconsistencies) is performed in parallel. That is, as the breadth first search is being executed, necessary vertex information is provided as input to the third test for infeasibility which is performed in parallel with the breadth first search. In particular, the third test receives information from the breadth first search concerning the identification of any closed paths. A closed path is identified whenever there exists two or more paths, p={1, . . . M}, from any source vertex V1 in the graph to any destination vertex V2. (i.e. p=2 in FIG. 3). A path distance metric, pdm_i, is computed for each path i as the sum of the distance metrics of each directed edge in the path, pdm_i=dme₁+dme₂+dm_i. . .

[0079] Referring again to FIG. 3, a path distance metric would be computed for paths 1 and 2. The path distance metric for path 1, =pdm_i+3, is computed as the sum of the distance metric for directed edge DE 70, (dme₇₀=+1) and the distance metric for directed edge DE 72 (dme₇₂=+2). The path distance metric for path 2, pdm₂=+4. The path distance metric information, along with other data, is used by the infeasibility test for determining whether a closed path is inconsistent. In particular, for a closed path comprising a set S of paths p={1, . . . M}, each path p from the set S will be individually compared with every other path in the set S to determine whether the closed path is inconsistent. For example, assume that three paths, p={1,2,3}, are identified while performing the breadth first search, with a common start and destination vertex. The third test for infeasibility would first compare paths 1 and 2. When path 3 is identified by the search algorithm, the infeasibility test would then compare paths 1 and 3 and 2 and 3. If during any particular path comparison it is determined that the test is violated, no further comparisons are required and the motifs comprising the closed path are determined to be offending motifs of the third type (closed paths with inconsistencies).

[0080] Referring now to FIG. 4, the method steps illustrating the details of performing the third test of infeasiblity for two paths identified by the breadth first search algorithm is shown. At step 80, two paths are identified by the breadth first search algorithm, P1 and P2. At step 82, a determination is made as to whether the first path, P1, is an overlap or a non-overlap path. At step 84, a similar determination is made as to whether the second path, P2, is an overlap path or a non-overlap path. At decision step 86, a determination is made as to whether both paths, P1 and P2, are overlap paths. If both paths P1 and P2 are overlap paths a branch occurs to decision step 92 where a determination is made as to whether the path distance metrics for paths 1, dp1, and 2, dp2, are equal. If not, a branch occurs to step 94 where it is determined that the motifs comprising the closed path constitute an infeasible set. The process then terminates at step 96. If, however, dp1=dp2, the distance metrics for paths 1 and 2 are found to be equal, then a branch occurs to step 93. At step 93 it is determined that the motifs corresponding to the vertices in the closed path do not constitute an infeasible set. The process then terminates at step 96.

[0081] At decision step 86, if either path P1 or P2 is determined not to be an overlap path then a branch is taken to decision step 88. At decision step 88 a determination is made as to whether the distance metric for the path determined to be the overlap path is greater than the distance metric for the path determined to be the non-overlap path. If true, a branch occurs to step 93. Otherwise a branch occurs to step 94. The process terminates at step 96.

[0082] 2.a Pruning Stage

[0083] Upon completing the test for the third type of infeasible set, closed path with inconsistencies, the three infeasible sets of offending motifs are completely determined. Having identified all of the offending motifs, the method of the present invention then prunes or removes at least one offending motif from each identified infeasible set to both allow the sequences to be aligned and to maximize a pre-specified cost function. In the exemplary embodiment, the cost function is defined as the sum of all the matching characters in the columns of the alignment in which at least K characters are identical. Other user defined cost functions include, but are not exclusive to; 1. the sum of the columns of the alignment in which at least K characters are identical, and, 2. retaining a maximum number of motifs for consideration at the alignment stage (i.e. pruning a minimum number of motifs). Other realizable cost functions may be considered by the method.

[0084] To maximize the cost function, as defined by the exemplary embodiment, defined as;

[0085] Cost Function: the sum of all the matching characters in the columns of the alignment in which at least K characters are identical to the method of the present invention utilizes a mathematical algorithm, defined in the prior art as the weighted set-covering algorithm, discussed at length in chapter 37 of Introduction to Algorithms, and herein incorporated by reference in its entirety.

[0086] Application of the weighted set-covering problem to the problem of removing offending motifs will now be explained by way of example. As was described, a number of offending motifs were identified and collected into one or more of the three infeasible sets.

[0087] 1. Pairwise incompatible,

[0088] 2. Closed cycles, and

[0089] 3. Closed path with inconsistencies.

[0090] Assume the first test of infeasibility, forbidden pairs (pairwise compatibility), identified a single pair of offending motifs thereby generating a first infeasible set;

[0091] 1st infeasible set: X={M3, M4}

[0092] Assume the second test of infeasibility, smallest cycles, identified three offending motifs thereby generating a second infeasible set;

[0093] 2nd infeasible set: Y={M1, M2, M7}

[0094] Assume the third test of infeasibility, closed paths with inconsistencies, identified three offending motifs thereby generating a third infeasible set;

[0095] 3rd infeasible set: Z={M3, M5, M7}

[0096] It is notable that the offending motifs which are the respective members of sets X, Y, and Z are not disjoint. That is, motifs M3 and M7 in the example provided, appear in two of the three categories. In general, a motif may appear in one, two, or all three of the sets.

[0097] Application of the weighted set covering problem to the method of the present invention in the context of the present example is as follows. From the three identified infeasible sets X, Y, and Z, derived as a result of performing the infeasibility tests, a union set A is formed where A={X, Y, Z}. In addition to defining union set A, sets B_m1through B_m7are defined, indicating which infeasible set(s) each offending motif, M1-M7, belongs to. For example, set B_m1describes that motif M1 is a member of only infeasible set Y, smallest cycles.

[0098] B_m1={Y}

[0099] B_m2={Y}

[0100] B_m3={X, Z}

[0101] B_m4={Z}

[0102] B_m5={X}

[0103] B_m7={Y, X}

[0104] As a further example, B_m3={X, Z} indicates that offending motif M5 is a member of both the first and third infeasible sets, forbidden pairs and closed paths with inconsistencies.

[0105] In terms of an unweighted set-covering problem solution, the objective of the present method is to find a minimum number of B sets that covers A. That is, what combination B sets will completely characterize set A where A is the set defined as the union of X, Y, and Z. For the present example, one possible solution is the combination of B sets {B_m2and B_m3} where B_m2defines Y and B_m3defines X and Z. Other possible solutions include {B_m7and B_m}}, or {B_m1and B_m3}. These and other solutions are realizable from the example above.

[0106] Any of the above solutions would satisfy the unweighted set-covering problem. To satisfy the weighted set-covering problem, however, removing B sets requires additional considerations. In particular, the weighted set-covering problem associates a weight with each identified B set. Furthermore, the weight is computed differently, depending upon the particular pre-specified cost function that is to be maximized. For the cost function defined by the exemplary embodiment the weight is computed as a multiplication of the number of solid characters of the motif represented by that B set multiplied by the number of sequences in which it appears. For example if motif M2 of set B_m2was “A B . D”, the number of solid characters in the motif is three. If motif M2 appeared in two of the N input sequences then the weight associated with B_m2would be 3*2=6. A weight value is similarly calculated for all identified B sets.

[0107] Having determined weights for each of the B sets, maximizing the pre-specified cost function of the exemplary embodiment then involves removing, or pruning, those B sets from consideration in the alignment stage which minimize the collective weight. For example, three possible solutions were explicitly defined for the present example, additional solutions are identifiable but not cited for the sake of clarity:

[0108] 1. {B_m2and B_m3)with associated weight w1

[0109] 2. (B_m1and B_m3)with associated weight w2

[0110] 3. {B_m4and B_m1)with associated weight w3

[0111] The selected solution, however, would be selected as the one whose associated weight is minimum, min(w1,w2,w3).

[0112] 3.b Alignment Stage

[0113] At the alignment stage, once the offending motifs (B sets) are removed from consideration in accordance with the weighted set-covering method, the N input sequences may then be aligned without conflict. For those offending motifs which are removed from further consideration their associated vertices in the graph will be pruned leaving only feasible vertices (motifs). The resulting directed graph is then used to determine step 1 defined below.

[0114] To render an alignment of the N sequences, each sequence is individually analyzed to determine two things; 1. the relative positions of any feasible motif contained in that sequence, and 2. the filler positions, where a filler position is a position in the sequence not occupied by a feasible motif. Knowing the motif and filler positions, each sequence is then aligned with every other sequence such that the feasible motifs are properly aligned. Fillers or gaps may be necessary to achieve the final alignment.

[0115] The alignment procedure described is further explained by way of example. Assume two sequences, N=2, are input to the method, described as;

[0116] Sequence 1: a b c d e f g h i j k l

[0117] Sequence 2: c d e x y z p q r g h i t u

[0118] Further assume that after the pruning stage, the feasible motifs (vertices) comprise “c d e” and “g h i”. The relative position of Motif “c d e” in sequences 1 and 2 is determined to occupy positions 3-5 in sequence 1 and positions 1-3 in sequence 2. Similarly the relative position of motif “g h i” occupies position 7-9 in sequence 1 and 10-12 in sequence 2. The relative position of the fillers may then be determined. Fillers occupy positions 1-2, 6, 10-12 in sequence 1 and 4-9, and 13-14 in sequence 2. With that information fillers or gaps are appropriately inserted to align the two sequences as shown;

[0119] Sequence 1: a b c d e f . . . g h i j k l

[0120] Sequence 2: c d e x y z p q r g h i t u

[0121] Note that the two sequences align without conflict as a consequence of each of the respective motifs being feasible.

[0122] Referring now to FIG. 5, a sample of some experimental results for the method of the present invention are illustrated. The experimental results pertain to the alignment of a large number of sequences with biological relevance. The input sequences are not shown. In FIG. 5, uppercase letters depict aligned characters and lowercase letters are portions of the input sequence that the method has not taken into account. The dots (‘ . ’) represent don't cares and the dashes represent the gaps enforced by the algorithm. Furthermore, some portions of the input have been removed to avoid clutter in the figure.

[0123] Results are illustrated on three sets of data: the first set is on a highly similar data set (histones) and the second set is on a fairly dissimilar (cytosine-specific DNA methylases list in PROSITE database (Release 13 entry accession #PS00094) data set, and a set of plant hemoglobins (listed in the PROSITE database Release 13 entry accession #PS00208). The first data set has 20 protein sequences of about 200 residues each with an alignment number K=10. The size of the resulting fensible set of motifs was found to be 22 for a maximal motif size of N=648. The second has 46 protein sequences with lengths going up to 500 residues, with an alignment number of K=18. The size of the resulting set of motifs was found to be 132 for a maximal motif size of N=810. A third set has 21 sequences with about 150 residues each, with an alignment number of K=14. A resulting compatible set of 42 motifs were identified from N=873 non-redundant motifs.

[0124] With particular reference to the first set of data, the method successfully identifies COREHISTONE-3 FOLDS. These folds represent sub-structures contained within the proteins. With reference to the second set of data, the alignment identifies the blocks that correspond well to the consensus blocks for the protein as identified in BLOCKS entry number BL0094 of BLOCKS database. With reference to the third set of data, the method identifies the known secondary structures (HELICES) present in the globin family members.

[0125] In the examples, some of the sequences have been truncated in the alignment and the output runs across multiple pages due to the size.

[0126] As will be obvious to those skilled in the art, the present invention can be practiced with a wide variety of modifications and variations, and is not limited except as specified in the accompanying claims.