On the maximal cliques in c-max-tolerance graphs and their application in clustering molecular sequences

Viewing the subject sequences aligned to a query sequence that result from a BLAST-based [1] comparison, in many cases one can identify groups of sequences clustering around different subintervals of the query sequence. Often, the decision by eye to which cluster a certain sequence belongs, is strongly depending on the order in which the sequences are presented. Fig. 1a) shows a schematic sketch of aligned sequences in random order. The sequences seem to form two, or maybe three groups. The same sequences in Fig. 1b) are ordered according to how many positions they have in common and colors indicate those sequences that share a large part of their sequence. The algorithm finds three different clusters of sequences. A cluster in this sense can be defined as a subset of aligned sequences that have approximately the same length and that are aligned to approximately the same subinterval of the query sequence. As we have argued above, the human eye may be fooled by the ordering of the presented sequences and humans are limited in the number of sequences they can group into clusters, and thus the automatic and objective computation of a clustering is an important task.

One kind of clustering can be obtained by using BlastClust from the NCBI-BLAST package [2]. It is a clustering tool designed especially to cluster protein or DNA sequences based on pairwise matches returned by the BLAST algorithm. It uses BLAST scores to assign statistical significance, matches pairs that reach that level of significance, then constructs clusters using a simple, greedy, single-linkage clustering method. However, empirically it is known that BlastClust is very limited with respect to the size of the input set. Furthermore, in many instances, the application of this heuristic is also limited by the fact that one unique clustering does not exist and that in fact many clusterings are possible and reasonable for different applications: for example, when annotating putative domains of a protein sequence, it is important to find common regions that are shared by large groups of sequences, while for the computation of phylogenetic trees it is important to find groups that share the largest possible part of their sequences.

In any case, under a given constraint, the groups should be as large as possible. Both constraints can be formulated as follows: Given two sequences, they are said to tolerate each other if both overlap at least an absolute amount t, the tolerance, or a relative tolerance c% of the length of the other sequence. This leads naturally to the family of max-tolerance graphs [3], where vertices represent intervals, and two vertices are connected by an edge if the corresponding sequences tolerate each other. The question for maximal clusters now reduces to finding maximal cliques in the max-tolerance graph.

Interestingly, the history of a special max-tolerance graph with t = 1 is deeply intertwined with the study of the DNA: Benzer was the first to raise the question, whether the sub-elements of the DNA are arranged linearly, or not [4]. To answer this question, sequences of DNA are represented as vertices, where two vertices are connected by an edge if the corresponding sequences are overlapping on at least one position. The resulting graph is an interval graph [5, 6]. It could be shown that a graph generated in this way can only be derived from fragments of a linear sequence if it has a certain property. Since the overlapping graph derived from DNA fragments showed this linearity property, this was yet another confirmation that the genome of organisms consists of long, linear DNA molecules. Searching for maximal cliques in interval graphs can be done efficiently in time linear to the size of the graph. However, for computing meaningful biological clusters, an absolute tolerance of t = 1 is not appropriate. For this purpose, it is necessary to increase the tolerance to a reasonable value and in addition to use the relative tolerance c of the length of the sequences. Thus, for c = 0.5 two sequences will tolerate each other if both have approximately the same length and share half of their sequence, or if one is not longer than twice the length of the other and overlaps the second completely, and anything in between. We call this kind of max-tolerance graph with a fixed relative tolerance c for all sequences a c-max-tolerance graph. The question of finding maximal groups of sequences that pairwise tolerate each other now reduces to the question of finding all maximal cliques, in the respective c-max-tolerance graph.

For general graphs, the computation of all maximal cliques is an NP-hard problem, since it can be reduced to the maximum clique problem which is again a classical NP-complete graph problem [7]. A well-known and popular branch-and-bound based software to compute maximum and maximal cliques in general graphs is Cliquer [8]. Cliquer's runtime is exponential. It will turn out that the maximal clique problem for the graphs considered here is not NP-hard.

Our aim here is now to show that, based on results of general max-tolerance graphs, an efficient algorithm can be applied to find all maximal cliques in a c-max-tolerance graph. In addition, we show how this set of all maximally cliques can then be used to find reasonable clusterings for biological sequences aligned to each other.

Focusing on the biological application of the idea of c-max-tolerance graphs, our approach is two-fold: First, we will give a practical algorithm to report all maximal cliques in polynomial time. This algorithm is conceptually simple, thus ensuring that it is easy to implement it correctly and numerically stable. We will also show in the experimental section that its implementation is very fast for typical data sizes in biological applications. Second, based on the general principle of the simple algorithm, we will present a more involved, output-sensitive algorithm in O(n² log n + out), where out denotes the size of the output. This is a considerable improvement compared with the corresponding bound for general max-tolerance graphs obtained in [9].

With this two-fold approach we follow the ideas of the field of algorithm design [10], that stresses the point that an easy algorithm should always be preferred for an implementation unless the data is so big that the more efficient algorithm has to be used.

The result of the algorithm can then be used to cluster the data in an automated and objective fashion. We are facing two problems here: First, the number of maximal cliques is strongly depending on the chosen c-value. Second, for many instances most of the sequences are elements of more than one maximal clique. Thus, we will give a heuristic based on the set of all maximal cliques that is able to find a reasonable clustering. We illustrate the algorithm using artificial as well as three different biological examples. The first of the biological data sets is an example of a repeat analysis, where a transposon from the nematode Pristionchus pacificus is compared with reads of BAC clone data covering about half of the genome of that organism. The second and third example are proteins that are compared with the UniProt database [11]. The article is structured as follows: in Section 2 we give a short review on results that are essential for the understanding of the two algorithms together with required definitions. Section 3 presents the first, simpler algorithm, followed by the description of the algorithm with a new, lower runtime complexity for c-max-tolerance graphs in Section 4. Section 5 gives some results on experiments and Section 6 summarizes and discusses our results.

In this section we will define the formal basis and repeat some earlier derived results that are essential for the understanding of the rest of the article.

Let S denote a set of n intervals I _i = [x _i , y _i ], 1 ≤ i ≤ n, x _i , y _i ∈ ℕ, x _i <y _i , where x _i denotes the coordinate of the left endpoint and y _i denotes the coordinate of the right endpoint. The length |I _i | of an interval is defined as y _i - x _i . Let I _i = [x _i , y _i ] and I _j = [x _j , y _j ] be two intervals. Their overlap length |I _i ⋂ I _j | is given by:

|I _i ⋂ I _j | = min{y _i , y _j } - max{x _i , x _j }.

A c-max-tolerance graph G = (V, E) is associated with a tolerance parameter 0 ≤ c ≤ 1 and it is then given by V = {v₁, ...,v _n } representing the set of intervals and

E = {{v _i , v _j }||I _i ⋂ I _j | ≥ c · max{|I _i |, |I _j |}}       (1)

for 1 ≤ i,j ≤ n. Pairs of intervals satisfying this condition are said to tolerate each other. More generally, interval I _i tolerates interval I _j if |I _i ⋂ I _j | > c · |I _i |.

In a general max-tolerance graph each interval I _i is associated with a tolerance 0 ≤ t _i ≤ |I _i | and here, two intervals are said to tolerate each other, if |I _i ⋂ I _j | ≥ max{t _i , t _j }.

Given a graph G = (V, E), a clique is a subset C ⊆ V of pairwise-connected vertices. A clique C is called a maximal clique if there is no other clique C' ⊃ C. A clique is a maximum clique if its cardinality is of largest possible size.

The following lemma has been derived for general max-tolerance graphs and applies also to c-max-tolerance graphs [9]:

Lemma 1 The number of maximal cliques in (c-)max-tolerance graphs is O(n³).

This result has been derived by a fundamental equivalence, for which we will need the following definitions: Let S = {[x₁, y₁], [x₂, y₂],..., [x _n , y _n ]} be a set of intervals and let 0 <c ≤ 1.0 denote the same relative tolerance for all intervals. Let T(S) denote a set of triangles Δ _i = {A _i , B _i , C _i }, A _i , B _i , C _i ∈ (ℕ × ℕ) with the following coordinates:

A _i = (x _i , |I _i | - ⌊|I _i | * (1 - c)⌋)       (2)

B _i = (x _i + ⌊|I _i | * (1 - c)⌋, |I _i | - ⌊|I _i | * (1 - c)⌋)       (3)

C _i = (x _i , |I _i |)       (4)

Thus, these triangles look like squares of length ⌊|I _i | * (1 - c)⌋, located at (x _i , ⌊|I _i | * (1 - c)⌋), cut in half by the diagonal from their upper left to their lower right corner. Therefore, we will call these triangles semi-squares. The next lemma states the following equivalence:

Lemma 2 The c-max-tolerance graph of (S,c) is equivalent to the intersection graph of T(S).

The lemma states that an edge between two intervals is drawn if and only if their corresponding semi-squares intersect. The formal proof is given in [9] but we want to motivate it here: an edge between two intervals A and B is drawn if they tolerate each other. It follows, that A can tolerate B if B does not start right of ⌊|A| (1 - c)⌋; otherwise |A ⋂ B <c · |A|. Note that the y-coordinate of a semi-square corresponding to an interval is equal to ⌊|A| (1 - c)⌋ and that the semi-square ends at x _A + ⌊|A| (1 - c)⌋. On the other hand, any overlap with A cannot be longer than |A|, and this only in the case where B fully overlaps A. If B starts at x _A + 1, their overlap can be of length at most |A| - 1, etc. The semi-square illustrates exactly this: The y-coordinate of the hypotenuse of any semi-square at position x gives a bound on the length of an overlap the corresponding interval can have with any other interval starting at x. Thus, two semi-squares will intersect if and only if the corresponding intervals tolerate each other: Let A be the longer interval, i.e., it determines the minimal length of the intersection, denoted by the y-coordinate of the basis of its corresponding semi-square. This semi-square can only be intersected by those semi-squares whose basis are below that of A and have a hypotenuse that is not lower than the basis of A at least at some point within the interval of A. Fig. 3 shows a set of intervals and their corresponding set of semi-squares.

The semi-squares can be ordered according to three different criteria: either by their left sides x _i , or their y-coordinate of their basis at |I _i | - ⌊|I _i | * (1 - c)⌋ or by their hypotenuse. Hypotenuses can be easily ordered according to their height at a given x-position by noting that the slope is -1 for all of them. Thus, for every x-coordinate with x between x _i and x _i + ⌊|I _i | * (1 - c)⌋, the following equation is valid:

const = y _i - (x - x _i )       (5)

where y _i is the endpoint of the corresponding interval I _i . Thus, for a given x-position the height of the hypotenuse of each active semi-square is determined by y _i + x _i - x (s. Fig. 2).

Since x is the same for all active semi-squares we can order the semi-squares by y _i + x _i . An important observation is that the order by the y-coordinate of the basis of the semi-squares or by their hypotenuse is not the same in general.

A very general tool in geometric algorithms is that of a sweep-line [12]: Here, so called events are ordered and then processed in this order, thus, we process the events in a sweep, where the invariant is that all events up until a certain point in the ordered list are already processed and that the events after this point still have to be processed. In geometrical problems such as finding the intersection points of a set of lines, an event is typically something as the begin or birth of a geometric structure and the end or death of it.

The sweep-line algorithm presented in this section is based on the semi-square intersection graph corresponding to a given c-max-tolerance graph. We will show that the theoretical runtime is a bit higher than the one given in the following section, but this algorithm is very easy to implement and it will lay the conceptual basis for the more involved algorithm presented in the next section. We will show in Section 5 that the implementation based on this simpler algorithm is nonetheless very fast, it will compute all maximal cliques in real data sets within a few seconds. The virtual sweep is done along the x-axis, handling two kinds of events: The birth of semi-squares and their death, i.e., start and end of a semi-square, respectively. As long as the sweep handles events with an x-position within the basis of a semi-square X, X is said to be active. The so-called y-structure Q holds all maximal cliques of active members at any time and every maximal clique is contained only once (y-structure invariant), where maximality refers only to the set of active semi-squares and not necessarily to the set of all semi-squares. In the first step all x-events are sorted non-decreasingly according to their x-position. If there are death and birth events at the same x-position, the birth events are first inserted into the sorted list. Then, the events are extracted from the sorted list which enables the virtual x-axis sweep.

If the event signals the death of a semi-square, it has to be deleted from every clique C it is an element of. If any of these cliques has been increased by the addition of a semi-square since the last deletion of a semi-square, it is a maximal clique and has to be printed out before the semi-square is deleted. Else, the semi-square is just deleted from C, without printing out the clique. A deletion can yield a reduced clique that is not unique anymore, or that is only a subset of some other clique in Q. Thus, for every clique in which some semi-square has been deleted, one has to verify its uniqueness and extensibility, and it will be deleted from Q if it is not unique or if it is extendable. Note that after checking these two properties, the y-structure invariance holds, because the y-structure just holds all maximal cliques of active semi-squares. Thus, not every maximal clique in the y-structure is a maximal clique with respect to the whole set of semi-squares.

The y-structure suffices to guarantee the correctness of this algorithm, because it guarantees that all maximal cliques of active semi-squares are kept in Q. Since a semi-square X is only non-active if it has not yet begun or if the rest of its corresponding interval is not long enough for any new semi-square to overlap it enough, it cannot be an element of a clique when it is non-active. For its active time, the y-structure guarantees that all its maximal cliques are contained in Q. This algorithm is thus correct because the y-structure invariance is maintained in all steps. The pseudocode is sketched in Fig. 4.

Runtime analysis

There are exactly 2n x-events that are sorted in O(n log n). For every birth of a semi-square X we have to check all cliques in Q whether they contain at least one semi-square that is intersecting with X. In the worst case, the y-structure can hold up to O(n³) cliques (see Lemma 1) with O(n) elements each. For the further analysis we will denote the maximal number of cliques in the y-structure by |Q|. Checking whether X can be added to a clique or building a candidate clique can thus be done in O(n * |Q|) for all cliques in Q. Since every newborn semi-square can either be added to a clique, split it, or not intersect any of its elements, there can be at most O(|Q|) candidate cliques with size O(n). For each of them we have to check whether it is extendable by one of the active semi-squares and whether there is a second candidate clique that is equal to that one. To check for equality we have to match each candidate clique with each other candidate clique which takes O(|Q|² * n²) in the naive approach. To check for extensibility we have to test each active semi-square whether it can extend the candidate clique. This can be done in O(|Q| * n), naively. This results in a theoretical runtime of O(|Q|² * n²) per birth. For death events, the clique has to be printed if it is a maximal clique. Altogether the size of out is equal to the number of symbols to print all maximal cliques. Since there are at most n elements per clique, the size of out is O(n⁴). Furthermore, the deletion of a semi-square in a clique can result in a clique that is no longer maximal with respect to all active semi-squares or that is equal to another clique in the y-structure. Thus, we have to check for equality and extensibility again. The clique will be deleted from the y-structure if its size is 0, if it is equal to some other clique in the y-structure, or if it is extendable by one of the other active semi-squares. Also this step requires a runtime of O(|Q|² * n²). The theoretical runtime for the whole algorithm is thus given by O(|Q|² * n³ + out), which is in the worst case O(n⁹).

We will now present some small improvements in the data structure that help to reduce the work for most real data sets. First, as a fast look-up we will compute the adjacency matrix adj for the intersection of all pairs of semi-squares in O(n²).

As we noted before, each semi-square can be characterized by y _i + x _i , the diagonal on which the hypotenuse lies. A clique C is now represented by two ordered lists: H(C), a list of lists of semi-squares for each diagonal on which at least one semi-square has its hypotenuse, and B(C) a list of lists of semi-squares for each y-coordinate on which at least one basis of any semi-square in the clique is positioned. A clique at a given x-position is now an assembly of ordered hypotenuses and bases as depicted in Fig. 5. A new semi-square X can intersect the set of bases and hypotenuses in the following ways:

1. 1.

   X intersects none of the bases and none of the hypotenuses. Since the bases are ordered, this only occurs if the upper point of X is lower than the lowest basis, or if the basis of X is higher than the highest hypotenuse. This can be calculated in O(1) because B(C) and H(C) are sorted. This results in a runtime of O(|Q|) for all cliques per birth event.

    
2. 2.

   X intersects either all bases and/or all hypotenuses. This is the case if X' s basis is lower than the lowest basis in C and X' s upper point is higher than the highest basis, and/or X' s basis is lower than the lowest hypotenuse and its upper point is higher than C' s highest hypotenuse. Both can be checked in O(1). If this is the case, X has just to be added to C. To insert X' s basis and hypotenuse at the correct place, a binary search with X' s height of the basis and the height of its hypotenuse is accomplished in O(log n) for every clique it is added to, i.e. with a runtime of at most O(|Q| log n) per birth event.

    
3. 3.

   X neither intersects all bases nor all hypotenuses, but at least one basis or hypotenuse in C. In this case, X will split the clique. It is easiest to just run through both lists, B(C) and H(C), and look in adj which semi-squares of C and X are intersecting. By running through the ordered lists, we can just add all intersecting semi-squares in the same order to the new candidate clique, maintaining the correct sorting. Each clique has at most n elements, look-up of adjacency is done in O(1) and adding to the list is also done in O(1) per intersecting semi-square. Last, X has to be inserted into the candidate clique, again in O(log n). Thus, generating a candidate clique can be done in O(n), resulting in |Q| * n for all cliques per birth event. To find out whether this candidate clique C' is maximal we use the following trick which is quite fast: Given the semi-square X _s representing the shortest interval, calculate the set of all active semi-squares that are not in C' but do intersect X _s . This set Z can be determined in O(n) by checking X' _s row in adj. Then we have to examine the extensibility of C' for only those semi-squares in Z. This takes O(n²), resulting in O(|Q| * n²) for every birth event. The check for equality is also quite cheap, because Q can again be ordered by for example the height of the lowest basis in a clique. This equality search for any other clique in the y-structure with the same height of the lowest basis can be done by a binary search in O(|Q| log |Q|). Since all cliques are ordered, a check on equality can be done in O(n), resulting in O(|Q| * n). Inserting the clique in the correct place in Q can again be done in O(|Q| log |Q|).

    

Again, the deletion of a semi-square will result in an addend of out to the runtime, and we have to check for equality and extensibility. Thus, our theoretical runtime is O(|Q| * n³+ out) which in the worst case is given by O(n⁶). Note that for real data sets this theoretical runtime will not be reached: First of all, intervals have very different lengths and therefore not all of them will be active at the same time. Second, most of the time the set Z of candidates for the extensibility test is very small. The same is true for the set of cliques with the same height of their lowest basis. Since the test for extensibility and uniqueness are the most expensive steps in the algorithm, this shows why the implementation can be so fast in real data sets, as will be demonstrated in section 5. But first we will give a more elaborated algorithm that shows that the runtime complexity of finding all maximal cliques in c-max-tolerance graphs is considerably lower than the corresponding runtime in max-tolerance graphs.

This section shows some results on using our implementation for finding all maximal cliques in artificial and biological data sets.

Runtime on artificial and biological datasets

In this section we want to demonstrate how fast the implementation of our simple algorithm is. We compare the runtime to that of Cliquer [8]. For the test we use both artificial as well as biological data sets. The artificial datasets were generated in the following way: For a given interval of some maximal length ℓ _max we computed all possible sub-intervals within this interval. From this set of possible sub-intervals, a given percentage number of sub-intervals was chosen uniformly at random for the further evaluations. Note that the number of possible sub-intervals in an interval with length ℓ is roughly ℓ², and thus the number of maximal cliques is O(ℓ⁶). We chose two values for ℓ _max , namely 20 and 100. Then, the overall runtime to calculate all maximal cliques was determined for all constraints within 0.05 ≤ c ≤ 0.95 in steps of 0.05. For Cliquer we could only use an interval of length ℓ _max = 20 and a randomly chosen set of 30% of all possible sub-intervals. For other test sets Cliquer had runtimes longer than many hours or even days. For ℓ _max = 20 and c = 0.05 Cliquer took nearly three hours to compute all maximal cliques, for c = 0.1 it took nearly one hour to compute all maximal cliques, for c = 0.15, c = 0.2, c = 0.25 it took approximately 8, 4, and 1 minutes, respectively. For all c ≥ 0.3 it took under a minute to compute all maximal cliques and it took below 0.01 seconds for all c ≥ 0.5. Our implementation needed no more than 50 ms for any of the c-values.

A second experiment was made on a much larger data set where the interval had a length of ℓ _max = 100. We used 10%, 20%, 30%, 40% and 50% of the total set of possible sub-intervals chosen uniformly at random. Again, Cliquer could not compute the maximal cliques for any of the data sets in any reasonable time. The results for the runtimes of our algorithm for these data sets are visualized in Fig. 9. Note that the data sets contained 1000, 2000, 3000, 4000 and 5000 intervals, respectively. We see that in all cases and for all chosen constraint values c, the runtimes were below 30 seconds. Even for the largest data set with 5000 sub-intervals the computation of the full statistics needed less than 4 minutes.

It is an interesting observation that the runtime seems to follow the number of cliques in the set, in contrast to the runtime of the Cliquer that is mainly determined by how interwoven the cliques are. The runtime of our algorithm thus indicates that for these random data sets the runtime is bound by out. Biological data sets are very different in structure from random data sets, which is caused by the clusteredness of their sequences. Since genes are often composites of different functional domains, aligned sequences have a high probability to center around those domains and build groups of 'similar' sequences with respect to the position and length of their alignment to the query.

For the test of our algorithm on biological sequences three BLAST-derived data sets were chosen as follows: From a bacterial artificial chromosome (BAC) and fosmid library of the genome of Pristionchus pacificus [13] that contains 78690 sequences, we reconstructed a genetic element that we characterized as a transposon from the maT family (S. Steigele, unpublished). Experimental hybridisation assays proved that this transposon is highly repetitive in the Pristionchus pacificus genome (A. Breit, unpublished). A BLASTN search (E-value < 10^-6) of this transposon against all sequences of the BAC and fosmid library was conducted, resulting in 126 hits below the E-value threshold. We refer to this data set by PpmaT. Furthermore, two multi-domain proteins, both containing repeated domains, were subject to BLASTP searches (E-value < 10^-6) versus the UniProt database [11]. The protein SH3 and multiple ankyrin repeat domains protein 1, short Shank1, contains 7 ANK repeats, 1 PDZ domain, 1 SAM domain and 1 SH3 domain. The BLAST search resulted in 86 hits below the E-value threshold. We refer to this data set by Shank1.

The protein Cadherin EGF LAG seven-pass G-type receptor 2, short Celsr2, is a receptor that may have an important role in cell-cell signaling during nervous system formation. It belongs to the G-protein coupled receptor 2 family, LN-TM7 subfamily, and contains 4 cadherin domains, 7 EGF-like domains, 1 GPS domain, 1 laminin EGF-like domain, 2 laminin G-like domains and 1 transmembrane receptor domain. The BLAST search resulted in 249 hits below the E-value threshold. We refer to this data set by Celsr2. Fig. 10 shows the times needed to compute all maximal cliques for the given data sets for all c between 0.05 and 0.95 in steps of 0.05. Generally we observe in all three cases that the runtime within one data set is almost independent of the constraint value c, differing only by a factor of about 2 across each curve. For all cases runtime was well below one second, while Cliquer needed for some of the c-values many hours and even longer. We have therefore again not presented Cliquer's runtime results for these data sets. The computation of the full curves needed less than one second.

For the same data sets, Fig. 11a) shows the number of cliques found at these constraint values and Fig 11b) shows the average number of cliques a sequence is member of. In comparison to the artificially derived data sets, runtime in these cases seems to be more dependent on the input size than on the number of maximal cliques. Furthermore, there is also no obvious dependency between the chosen c-value and the number of maximal cliques. We also observe for quite a large range of constraint values c on average a sequence is contained in up to 30% of all maximal cliques. Even for large c-values (c > 0.6) on average a sequence is contained in up to 10% of the maximal cliques.

Computing a partition of S from maximal cliques

Since the membership of a sequence is almost never unambiguous, a heuristic is needed in order to partition the set of sequences into meaningful clusters. The set of all maximal cliques for a given constraint c guarantees us that every two members of a maximal clique will at least share a part of their sequence of length (|I _i | * c). Let I _c (C _i ) denote that part of the query sequence that is overlapped by all sequences in clique C _i , i.e.,

I _c (C _i ) = [max {x _j |I _j ∈ C _i }, min {y _j |I _j ∈ C _i }]       (6)

It is an important observation that this common part of the sequence does not have to overlap each sequence in C _i by c* 100%, as is shown in Fig. 12. If now a sequence X is a member of many maximal cliques, it should be assigned to that clique C _i whose shared sequence part I _c (C _i ) overlaps X best, i.e. is maximal.

The other problem is that the number and structure of maximal cliques is of course depending on the chosen constraint value c. There are two opposing situations that make it hard to find the best c-value: The first is sketched in Fig 13a) where three different clusters of sequences are visible, but they will only be found if c is larger than 0.6. In Fig. 13b) there is only one cluster, but if c is larger than 0.6 then two maximal cliques will emerge.

Thus, we have two contradicting goals: To get a higher specificity of the maximal cliques it is good to use a large c-value. But if it is too large it will split 'good' cliques into too many parts. Our aim was to find a heuristic that is based on a c-value that is small enough as not to split 'good' clusters but can find subsets of the maximal cliques that are more specific than the whole clique.

For this, the following heuristic is used:

1. 1.

   Add all maximal cliques to a list L.

    
2. 2.

   For all cliques C in L compute the shared overlap interval I _C , i.e., the maximal interval that is overlapped by all of the intervals in the clique.

    
3. 3.

   For every interval I that is in more than one clique, compute the overlap between I _C and I and multiply this value with the number of intervals in this clique.

    
4. 4.

   Assign every interval I to that clique where the product of overlap with I _C and number of intervals in that clique is maximal.

    

This basic variant yields good results but it can be improved by the following procedure: For every pair A, B of cliques compute the intersection A ⋂ B. If there are many intervals in this intersection set, e.g., |A ⋂ B| > 3, and only a few intervals that are not, the intersection of the two cliques represents intervals that are very similar to each other. These intervals would be in one clique if the c-value would have been more restrictive, i.e., larger. But due to the low c-value, there are other intervals to the left and right of the intervals in the intersection set that interact with the intervals in A ⋂ B but not with each other. As discussed above, a low c-value has the advantage of keeping good clusters together, and thus c should not be too high in the beginning. However, the effect of a larger c-value can be mimicked by adding those intersection sets to the list L in the first step of the above described mechanism and to allow the intervals to assign them to any set in L be it a clique or an intersection set.

The given heuristic guarantees that every two intervals in a set tolerate each other, because it is built on the maximal cliques of tolerating intervals. The decision of membership has been built on a composite of the overlap with the shared sequence and the size of the set to reduce the number of sets globally. The size of a set is of course volatile because not all of its designated members will finally be assigned to it. But in each case the overlap with the shared sequence of any set will never decrease if another member is deleted from the set, so the heuristic is stable in this respect.

We have already indicated at the beginning of this article (s. Fig. 1), that many partitions of a set S of sequences seem to be reasonable by regarding only their begin and ends. In order to facilitate a manual revision of the calculated partition of each set of sequences these sets are displayed in the order of the diagonal of their smallest member, i.e., the sum of its x-position and its length.

The results of this heuristic applied to the three biological data sets are given in Fig. 14a–c. The colorization helps to identify the members of one clique. When comparing the cliques with the annotation of the input query sequence individual domains as well as combination of several domains were retrieved and corresponded almost perfectly to the annotated domains of the input protein sequence.

BLAST is a very common algorithm in bioinformatics used for annotation of an unknown sequence or for sampling sequences from a database sharing homologous regions with the input sequence. One goal of BLAST is to identify regions of conserved DNA or protein sequence. Another goal is to detect groups of sequences that are highly similar with respect to the position and the length of conservation with the query sequence. For a given query sequence, BLAST returns the hits ordered statistical significance of each hit. When viewed, this ordering does not necessarily represent a clustering with respect to length of the hits and position in the query sequence. Thus, extraction of possible clusterings of the sequences is not always straightforward. One goal of this article was therefore to provide a method that allows an automatic clustering of sequences returned from a BLAST run, such that the user can decide whether to maximize the lengths of the common region of the sequences within a cluster or whether to maximize the size of the clusters.

In this article we have shown that finding groups from BLAST reports can be reduced to computing maximal cliques in so-called c-max-tolerance graphs. We have presented two algorithms that allow the enumeration of all maximal cliques in c-max-tolerance graphs in polynomial time. The first algorithm has an easy implementation to compute all maximal cliques in c-max-tolerance graphs, but it has a theoretical runtime of O(n⁶). The second algorithm is technically very elaborate, and shows that enumeration of all maximal cliques in c-max-tolerance graphs can be computed in a runtime of O(n² log n + out), where out denotes the size of the output. The algorithm improves a result found for general max-tolerance graphs, for which the computation of all maximal cliques has a runtime of O(n³) [9].

Though our simple algorithm in theory has a higher runtime than the more technical one, when applied to real data, runtimes were nonetheless impressive. Comparison with the popular application Cliquer [8], another exact algorithm to compute all maximal cliques in general graphs, showed that for many instances of our test cases, Cliquer could not report results within hours, while our algorithm finished most computations within milliseconds.

The maximal cliques represent a clustering of the sequences detected with BLAST. Thus with the help of this clique-finding algorithm, clusters can now be quickly and automatically identified even from very large BLAST reports. The result of the clustering can then be further processed for example for phylogenetic reconstructions of the sequences in the clusters. At this point one might discuss the biological usefulness to compute all maximal cliques. As we have seen from the biological data sets used for our experiments, it is very common to have largely interlinked sequences, i.e., where one unique clustering and/or maximal clique does not exist. This is for example a common feature in the case of multidomain proteins, where the query protein contains say k domains, and BLAST reported subject proteins will show a variety of combinations of these k domains. A goal of applying the maximal cliques search algorithm is to cluster individual domains as well as all possible combinations of domains. The examples in figure 14 illustrate the result of these clusterings. Here one also clearly sees that other reasonable partitions are possible.

Another application for this algorithm is the prediction of gene boundaries within a genomic sequence. For this, mapping of cDNA/ESTs to a genome sequence is used. The goal of this strategy is to predict the gene boundaries and/or its exon/intron structure using cDNA/ESTs. Applying our simple maximal clique algorithm to a large test case (with more than 2800 intervals), for which a human genomic locus was compared against the EST database from the NCBI, yielded results within a few seconds with a perfect clustering of EST loci (data not shown).

The use of the c-value is an important parameter for biological applications. Again we think of the situation of reconstructing phylogenetic trees from sequences reported from a BLAST run. Here different c parameters can be used to either optimize the length of the overlaps or the size of the maximal cliques. A model where the factor c is not a tolerance relative to the length but a fixed constant independent of the length might be of some biological relevance. An extension towards such a model or even a mixture seems plausible.

Examples from biological data sets indicated that most of the sequences are contained in several cliques. Again, it is not straightforward to deduce a disjoint clustering of the sequences from the maximal cliques. We have described a heuristic that computes a partitioning of sequences from the set of all maximal cliques in a c-max-tolerance graph. This heuristic when applied to protein sequences containing several domains partitions the set of all maximal cliques into clusters, where the sequences within the clusters reflect the domain structure of the input sequence.

After cliques and clusters have been computed postprocessing such as visualization of the groups or integration into pipelines that need clustered sequence data is now very easy and automizable. And therefore we hope that the algorithms described here have further impact for the bioinformatics community.

The implementation of the first algorithm is available on request from the first author.