Inken Wohlers
email: inken.wohlers@uni-due.de

Mathilde Le Boudic-Jamin
email: mathilde.leboudic-jamin@irisa.fr

Hristo Djidjev
email: djidjev@lanl.gov

Gunnar W Klau
email: gunnar.klau@cwi.nl

Rumen Andonov
email: rumen.andonov@irisa.fr

Exact protein structure classification using the maximum contact map overlap metric

In this work we propose a new distance measure for comparing two protein structures based on their contact map representations. We show that our novel measure, which we refer to as the maximum contact map overlap (max-CMO) metric, satisfies all properties of a metric on the space of protein representations. Having a metric in that space allows to avoid pairwise comparisons on the entire database and thus to significantly accelerate exploring the protein space compared to nonmetric spaces. We show on a gold-standard classification benchmark set of 6, 759 and 67, 609 proteins, resp., that our exact k-nearest neighbor scheme classifies up to 95% and 99% of queries correctly. Our k-NN classification thus provides a promising approach for the automatic classification of protein structures based on contact map overlap.

Introduction

Understanding the functional role and evolutionary relationships of proteins is key to answering many important biological and biomedical questions. Because the function of a protein is determined by its structure and because structural properties are usually conserved throughout evolution, such problems can be better approached if proteins are compared based on their representations as threedimensional structures rather than as sequences. Databases such as SCOP [START_REF] Murzin | SCOP: a structural classification of proteins database for the investigation of sequences and structures[END_REF] and CATH [START_REF] Orengo | CATH-a hierarchic classification of protein domain structures[END_REF] have been built to organize the space of protein structures. Both SCOP and CATH, however, are constructed partly based on manual curation, and many of the currently over 98, 000 protein structures in the protein data bank (PDB) [START_REF] Bernstein | The protein data bank: A computer-based archival file for macromolecular structures[END_REF] are still unclassified. Moreover, classifying a newly found structure manually is both expensive in terms of human labor and slow. Therefore, computational methods that can accurately and efficiently complete such classifications will be highly beneficial. Basically, given a query protein structure, the problem is to find its place in a classification hierarchy of structures, for example, to predict its family or superfamily in the SCOP database.

One approach to solving that problem is based on having introduced a meaningful distance measure between any two protein structures. Then the family of a query protein q can be determined by comparing the distances between q and members of candidate families and choosing a family whose members are "closer" to q than members of the other families, where the precise criteria for deciding which family is closer depend on the specific implementation. The key condition and a crucial factor for the quality of the classification result is having an appropriate distance measure between proteins.

Several such distances have been proposed, each having its own advantages. Recently, a number of approaches based on a graph-based measure of closeness called contact map overlap (CMO) [START_REF] Godzik | Regularities in interaction patterns of globular proteins[END_REF] have been shown to perform well [START_REF] Andonov | Maximum contact map overlap revisited[END_REF][START_REF] Caprara | 1001 optimal PDB structure alignments: integer programming methods for finding the maximum contact map overlap[END_REF][START_REF] Malod-Dognin | Gr-align: fast and flexible alignment of protein 3d structures using graphlet degree similarity[END_REF][START_REF] Malod-Dognin | Using dominances for solving the protein family identification problem[END_REF][START_REF] Pelta | A simple and fast heuristic for protein structure comparison[END_REF][START_REF] Wohlers | CSA: comprehensive comparison of pairwise protein structure alignments[END_REF][START_REF] Xie | A reduction-based exact algorithm for the contact map overlap problem[END_REF]. Informally, CMO corresponds to the maximum size of a common subgraph of the two contact map graphs, see the next section for the formal definition. Although CMO is a widely used measure, none of the CMO-based distance methods suggested so far satisfies the triangle inequality and, hence, introduces a metric on the space of protein representations. Having a metric in that space establishes a structure that allows much faster exploration of the space compared to non-metric spaces. For instance, all previous CMO-based algorithms require pairwise comparisons of the query with the entire database. With the rapid increase of the protein databases, such a strategy will unavoidably create performance problems even if the individual comparisons are fast.

In this work we propose a new distance measure for comparing two protein structures based on their contact map representations. We show that our novel measure, which we refer to as the maximum contact map overlap (max-CMO) metric, satisfies all properties of a metric. This enables us to describe a given protein database as a metric space where we model each protein family as a ball with a specially chosen protein from the family as center. We exploit this representation to accurately and efficiently classify a query protein according to its k nearest neighbors. We demonstrate that using polynomial-time approximations of max-CMO in terms of lower-bound upper-bound intervals speeds up the classification process significantly, without sacrificing its accuracy. We point out that our approach is not heuristic and guarantees solving the classification problem to provable optimality with respect to our max-CMO metric and that we do so without having to compute all query-target alignments to optimality.

The metric property has been proved for a related graph distance based on the maximum common subgraph of two graphs [START_REF] Bunke | A graph distance metric based on the maximal common subgraph[END_REF]. Here, however, we consider contact map overlap, which is based on an underlying order-preserving alignment. Furthermore, the size of a contact map overlap is the number of common edges, whereas the graph distance is based on the number of common nodes.

We show on a gold-standard classification benchmark set of 6, 759 proteins that our exact k-nearest neighbor scheme classifies up to 224 out of 236 queries correctly, and on a large, extended version of the data set that contains 67, 609 proteins even up to 1361 out of 1369 queries. Our k-NN classification thus provides a promising approach for the automatic classification of protein structures based on flexible contact map overlap alignments. Amongst the other existing (non-CMO) protein structure comparison methods we are aware of only one satisfying the triangle inequality. This so called scaled Gauss metric (SGM) introduced in [START_REF] Rogen | Automatic classification of protein structure by using gauss integrals[END_REF] and further developed in [START_REF] Harder | Fast large-scale clustering of protein structures using Gauss integrals[END_REF] is shown to be very successful in automatic classification. In our work, however, we focus on contact map overlap and a comparison to classification algorithms based on different concepts is outside the scope of this paper.

The maximum contact map overlap metric

We introduce here the notions of contact map overlap (CMO) and the related max-CMO distance between protein structures. A contact map describes the structure of a protein P in terms of a simple, undirected graph G = (V, E) with vertex set V and edge set E. The vertices of V are linearly ordered and correspond to the sequence of residues of P . Edges denote residue contacts, that is, pairs of residues that are close to each other. More precisely, there is an edge (i, j) between residues i and j iff the Euclidean distance in the protein fold is smaller than a given threshold. The size |G| := |E| of a contact map is the number of its contacts. Given two contact maps G 1 (V, E 1 ) and G 2 (U, E 2 ) for two protein structures, let I = (i 1 , i 2 , . . . , i m ) and J = (j 1 , j 2 , . . . , j m ) be subsets of V and U , respectively, respecting the linear order. Vertex sets I and J encode an alignment of G 1 and G 2 in the sense that vertex i 1 is aligned to j 1 , i 2 to j 2 and so on. In other words, the alignment (I, J), is a one-to-one mapping between the sets V and U . Given an alignment (I, J), a shared contact (or common edge) occurs if both (i k , i l ) ∈ E 1 and (j k , j l ) ∈ E 2 exist. We say in this case that the shared contact (i k , i l ) is activated by the alignment (I, J). The maximum contact overlap problem consists in finding an alignment (I * , J * ) that maximizes the number of shared contacts and CMO(G 1 , G 2 ) denotes then this maximum number of shared contacts between the contact maps G 1 and G 2 , see Figure 1.

Computing CMO(G 1 , G 2 ) is NP-hard following from [START_REF] Lathrop | The protein threading problem with sequence amino acid interaction preferences is NP-complete[END_REF]. Nevertheless, maximum contact map overlap has been shown to be a meaningful way for comparing two protein structures [START_REF] Andonov | Maximum contact map overlap revisited[END_REF][START_REF] Caprara | 1001 optimal PDB structure alignments: integer programming methods for finding the maximum contact map overlap[END_REF][START_REF] Malod-Dognin | Gr-align: fast and flexible alignment of protein 3d structures using graphlet degree similarity[END_REF][START_REF] Malod-Dognin | Using dominances for solving the protein family identification problem[END_REF][START_REF] Wohlers | CSA: comprehensive comparison of pairwise protein structure alignments[END_REF][START_REF] Xie | A reduction-based exact algorithm for the contact map overlap problem[END_REF]. Previously, several distances have been proposed based on the maximum contact map overlap, for example, D min [START_REF] Caprara | 1001 optimal PDB structure alignments: integer programming methods for finding the maximum contact map overlap[END_REF][START_REF] Pelta | A simple and fast heuristic for protein structure comparison[END_REF] and D sum [START_REF] Andonov | Maximum contact map overlap revisited[END_REF][START_REF] Malod-Dognin | Gr-align: fast and flexible alignment of protein 3d structures using graphlet degree similarity[END_REF][START_REF] Xie | A reduction-based exact algorithm for the contact map overlap problem[END_REF] with

D min (G 1 , G 2 ) = 1 - CMO(G 1 , G 2 ) min{|E 1 |, |E 2 |} and D sum (G 1 , G 2 ) = 1 - 2CMO(G 1 , G 2 ) |E 1 | + |E 2 | .
These distances have the disadvantage that they are no metrics as the following lemma shows (see the extended version [START_REF] Wohlers | Exact protein structure classification using the maximum contact map overlap metric[END_REF] for a proof).

Lemma 1. Distances D min and D sum do not satisfy the triangle inequality.

Let G 1 (V, E 1 ), G 2 (U, E 2 ) be two contact maps graphs. We propose a new distance D max (G 1 , G 2 ) = 1 - CMO(G 1 , G 2 ) max{|E 1 |, |E 2 |} . (1)
The following claim states that D max is indeed a distance metric on the space of contact maps and we refer to it as the max-CMO metric.

Lemma 2. D max is a metric on the space of contact maps.

Proof. To prove the triangle inequality for the function D max , we consider three contact maps

G 1 (V, E 1 ), G 2 (U, E 2 ), G 3 (W, E 3
), and we want to prove that

D max (G 1 , G 2 ) + D max (G 2 , G 3 ) ≥ D max (G 1 , G 3
). We will use the fact that a similar function d max on sets is a metric, which is defined as

d max (A, B) = 1 - |A ∩ B| max{|A|, |B|} . (2)
The mapping M corresponding to CMO(G 1 , G 2 ) generates an alignment (V ′ , U ′ ), where V ′ ⊆ V and U ′ ⊆ U are ordered sets of vertices preserving the order of V and U , correspondingly. Since M is a one-to-one mapping, we can rename the vertices of U ′ to the names of the corresponding vertices of V ′ and keep the old names of the vertices of U \ U ′ . Denote the resulting ordered vertex set by U and denote by E 2 the corresponding set of edges. Define the graph

G 2 = (U , E 2 ). Note that |E 2 | = |E 2 | and any common edge discovered by CMO(G 1 , G 2 ) has the same endpoints (after renaming) in E 2 as in E 1 ; hence CMO(G 1 , G 2 ) = CMO(G 1 , G 2 ) = |E 1 ∩ E 2 |. Then from (2) D max (G 1 , G 2 ) = 1 - CMO(G 1 , G 2 ) max{|E 1 |, |E 2 |} = 1 - |E 1 ∩ E 2 | max{|E 1 |, |E 2 |} = d max (E 1 , E 2 ) .
Similarly, we compute the mapping corresponding to CMO(G 2 , G 3 ) and generate an optimal alignment (U ′ , W ′ ). As before, we use the mapping to rename the vertices of W ′ to the corresponding vertices of U ′ and denote the resulting sets of vertices and edges by W and E 3 . Similarly to the above case, it follows that

D max (G 2 , G 3 ) = d max (E 2 , E 3 ).
Combining the last two equalities, we get

D max (G 1 , G 2 ) + D max (G 2 , G 3 ) = d max (E 1 , E 2 ) + d max (E 2 , E 3 ) ≥ d max (E 1 , E 3 ). ( 3
)
On the other hand, E 1 ∩ E 3 contains only edges jointly activated by the alignments (V ′ , U ′ ) and (U ′ , W ′ ) and its cardinality is not larger than CMO(G 1 , G 3 ), which corresponds to the optimal alignment between G 1 and G 3 . Hence

|E 1 ∩ E 3 | ≤ CMO(G 1 , G 3 ) and, since |E 3 | = |E 3 |, d max (E 1 , E 3 ) = 1 - |E 1 ∩ E 3 | max{|E 1 |, |E 3 |} ≥ 1 - CMO(G 1 , G 3 ) max{|E 1 |, |E 3 |} = D max (G 1 , G 3 ).
Combining the last inequality with (3) proves the triangle inequality for D max .

The other two properties of a metric, that

D max (G 1 , G 2 ) ≥ 0 with equality if and only if G 1 = G 2 and D max (G 1 , G 2 ) = D max (G 2 , G 1 ), are obviously also true. ⊓ ⊔ If instead of CMO(G 1 , G 2
) one computes lower or upper bounds for its value, replacing those values in (1) produces an upper or lower bound for D max , respectively.

Nearest neighbor classification of protein structures

We suggest to approach the problem of classifying a given query protein structure with respect to a database of target structures based on a majority vote of the k nearest neighbors in the database. Nearest neighbor classification is a simple and popular machine learning strategy with strong consistency results, see for example [START_REF] Altman | An introduction to kernel and nearest-neighbor nonparametric regression[END_REF]. An important feature of our approach is that it is based on a metric and we fully profit from all usual benefits when exploring a metric space [START_REF] Moreno-Seco | A modification of the laesa algorithm for approximated k-nn classification[END_REF].

Finding family representatives

In order to minimize the number of targets with which a query has to be compared directly, i.e., via computing an alignment, we designate a representative central structure for each family. Let d denote any metric. Each family F ∈ C can then be characterized by a representative structure R F and a family radius r F determined by

R F = arg min A∈F max B∈F d(A, B), r F = min A∈F max B∈F d(A, B). (4)
In order to find R F and r F , we compute, during a preprocessing step, all pairwise distances within F. We aim to compute these distances as precise as possible, using a sufficiently long run time for each pairwise comparison. Since proteins from the same family are structurally similar, the alignment algorithm performs favorably and we can usually compute intra-family distances optimally. These distances obtained during preprocessing are later re-used during k-NN classification for computing triangle bounds.

Dominance between target protein structures

In order to find the target structures which are closest to a query q, we have to decide for a pair of targets A and B which one is closer. We call such a relationship between two target structures dominance: Definition 3 (dominance). Protein A dominates protein B with respect to a query q if and only if d(q, A) < d(q, B).

In order to conclude that A is closer to q than B, it may not be necessary to know d(q, A) and d(q, B) exactly. It is sufficient that A directly dominates B according to the following rule.

Lemma 4 (direct dominance). Protein A dominates protein B with respect to a query q if d(q, A) < d(q, B), where d(q, A) and d(q, B) are an upper and lower bound on d(q, A) and d(q, B), respectively.

Proof. Follows from the inequalities d(q, A) ≤ d(q, A) < d(q, B) ≤ d(q, B).

⊓ ⊔

The idea of dominance is crucial for reducing the number of computations in our approach. Based on the relationship of polynomial-time lower and upper bounds a dominated protein is discarded from further consideration. Although the precise distance between proteins and the associated alignment are not computed, which is an NP-hard problem, the accuracy of the classification is not sacrificed. In its simplest form this idea has been first proposed in [START_REF] Malod-Dognin | Using dominances for solving the protein family identification problem[END_REF]. Here we extend it by exploiting the properties of a metric space as shown below.

Given a query q, a representative r and a target A, the triangle inequality provides an upper bound, while the reverse triangle inequality provides respectively a lower bound on the distance from query q to target A d(q, A) ≤ d(q, r) + d(r, A) and d(q, A) ≥ |d(q, r)d(r, A)| .

We define the triangle upper (resp. lower) bound as

d △ (q, A) = min r∈R {d(q, r) + d(r, A)} , d ▽ (q, A) = max r∈R max{d(q, r) -d(r, A), d(r, A) -d(q, r)} . Lemma 5. d ▽ (q, A) ≤ d(q, A) ≤ d △ (q, A) Proof. d ▽ (q, A) = max r∈R max{d(q, r) -d(r, A), d(r, A) -d(q, r)} ≤ max r∈R |d(q, r) - d(r, A)| ≤ d(q, A) ≤ min r∈R d(q, r) + d(r, A) ≤ min r∈R d(q, r) + d(r, A) = d △ (q, A). ⊓ ⊔
Using Lemma 5 we derive supplementary sufficient conditions for dominance, which we call indirect dominances.

Lemma 6 (indirect dominance).

Protein A dominates protein B with respect to query q if d △ (q, A) < d ▽ (q, B).

Proof. d(q, A)

Lemma 5 ≤ d △ (q, A) < d ▽ (q, B) Lemma 5 ≤ d(q, B). ⊓ ⊔

Classification algorithm

K-nearest neighbor classification is a scheme which assigns the query to the class to which most of the k targets belong which are closest to the query. In order to classify, we therefore need to determine the k structures with minimum distance to the query and assign the super-family to which the majority of the neighbors belong. As seen in the previous section, we can use bounds to decide whether a structure is closer to the query than another structure. This can be generalized to deciding whether or not is it possible for a structure to be among the k closest structures in the following way. We construct two priority queues LB and UB whose elements are (t, lb(q, t))) and (t, ub(q, t)), respectively, where q is the query and t the target. The algorithm works with any lower bound lb(q, t) on the distance between q and t, for example d(q, t) or d ▽ (q, t) and with any upper bound ub(q, t) on d(q, t), for example d(q, t) or d △ (q, t). In our current implementation we use D max as a distance while lower and upper bounds are polynomially computed based on Lagrangian relaxation as explained in [START_REF] Andonov | Maximum contact map overlap revisited[END_REF]. The quality of these bounds for the purpose of protein classification has been already demonstrated in [START_REF] Malod-Dognin | Gr-align: fast and flexible alignment of protein 3d structures using graphlet degree similarity[END_REF][START_REF] Malod-Dognin | Using dominances for solving the protein family identification problem[END_REF]. LB and UB are sorted in the order of increasing distance. The k-th element in queue UB is denoted by t UB k . Its distance to the query, d(q, t UB k ), is the distance for which at least k target elements are closer to the query. Therefore we can safely discard all those targets which have a lower bound distance of more than d(q, t UB k ) to query q. That is, t UB k dominates all targets t for which lb(q, t) > ub(q, t UB k ).

Validation setup

We evaluated the classification performance and efficiency of different types of dominance of our algorithm on domains from SCOPCath [START_REF] Csaba | Systematic comparison of SCOP and CATH: a new gold standard for protein structure analysis[END_REF], a benchmark that consists of a consensus of the two major structural classifications SCOP [START_REF] Murzin | SCOP: a structural classification of proteins database for the investigation of sequences and structures[END_REF] (version 1.75) and Cath [START_REF] Orengo | CATH-a hierarchic classification of protein domain structures[END_REF] (version 3.2.0). We use this consensus benchmark in order to obtain a gold-standard classification that very likely reflects structural similarities that are detectable automatically, since two classifications, each using a mix of expert knowledge and automatic methods, agree in their superfamily assignments. SCOPCath has been filtered such that it only contains proteins with less than 50% sequence identity. Since this results in a rather small benchmark with only 6, 759 structures, we added these filtered structures for our evaluation in order to have a benchmark representative of the merged databases SCOP and Cath. There were 264 domains in extended SCOPCath which share more than 50% sequence similarity with a domain in SCOPCath, but do not both belong to the same SCOP family; since their families are perhaps not in SCOPCath and their classification in SCOP and Cath may not agree, we removed them. This way we obtained 60, 850 additional structures. These belong to 1, 348 super-families and 2, 480 families of which 2, 093 families have more than one member. For SCOPCath, there are 1, 156 multi-member families. Structures and families are divided into classes according to Table 1. For super-family assignment, we compared a structure only to structures of the corresponding class since class membership can in most cases be determined automatically, for example by a program that computes secondary structure content. We then computed all-versus-all distances (2) or distance bounds within each family using optimal maximum contact map overlap or the upper Lagrangian bound on it and determined the family representative according to Equation (4). For every pairwise distance computation, we used a maximum time limit of 10 s. Since most comparisons were computed optimally, the average run time is approximately 2 s.

For classification, we randomly selected one query from every family with at least six members. This resulted in 236 queries for SCOPCath and 1, 369 queries for the extended SCOPCath benchmark. For every query, the k = 10 nearest neighbor structures from SCOPCath and extended SCOPCath, respectively, were computed using our k-NN Algorithm. The algorithm is a two-step procedure. First it improves distance bounds by applying several rounds of triangle dominance, in which the maximum contact map overlap bounds from query to representatives are updated, and second it switches to pairwise dominance, for which the distance to any remaining target is computed. In the first step, query representative distances are computed using an initial time limit of τ = 1 s, then triangle dominance is applied to all targets and the algorithm iterates with time limit doubled until a termination criterion is met. This way, bounds on query target distances are improved successively. The computation of triangle dominance terminates if any of the following holds (i) k targets are left (ii) all query-representative distances have been computed optimally or with a time limit of 32 CPU seconds (iii) the number of targets did not reduce from one round to the next. Pairwise dominance terminates if any of the following holds (i) k targets are left or all remaining targets belong to the same super-family (ii) all query-target distances have been computed with a time limit of 32 CPU seconds. The query is then assigned to the super-family to which the majority of the k nearest neighbors belongs. In cases in which the pairwise dominance terminates with more than k targets or more than one super-family remains, the exact k nearest neighbors are not known. In that case we order the targets based on the upper bound distance to the query and assign the super-family using the top ten queries. In the case that there is a tie among the superfamilies to which the top ten targets belong, we report this situation. In order to investigate the impact of k on classification accuracy, we additionally decreased k from 9 to 1, using each time the k + 1 nearest neighbors from the classification result for k + 1. In the case that for a query more than k + 1 queries remained in this classification, we used all of them for searching for the k nearest neighbors, but put an additional termination criterion which prevents extremely long run times for a few queries. Due to the large number of computations, classifications were run on different architectures on clusters with various load and are therefore only used for order of magnitude comparison.

Computational results

Characterizing the distance measure

In a first, preprocessing step we evaluate how well our distance metric captures known similarities and differences between protein structures by computing intra-family and inter-family distances. A good distance for structure comparison should pool similar structures, i.e., from the same family, whereas it should locate dissimilar structures from different families far apart from each other. In order to quantify such characteristics, we compute for each family with at least two members a central, representative structure according to Equation ( 4). Therefore, we compute the distance between any two structures that belong to the same family. Such intra-family distances should ideally be small. We observe that the distribution of intra-family distances differ between classes and are usually smaller than 0.5, except for class c. For the four major protein classes, there is a distance peak close to 0 and another one around 0.2.

We then compute a radius around the representative structure that encompasses all structures of the corresponding family. The number of families with a given radius decreases nearly linearly from 0 to 0.6, with most families having a radius close to zero, and almost no families having a radius greater than 0.6.

Considering that the distance metric is bound to be within 0 and 1, interfamily distances and radii show that our distance overall captures well the similarity between structures. Further, we investigate the distance between protein families by computing their overlap value as defined by d(R F1 , R F2 )r F1r F2 . Most families are not close to each other according to our distance metric. Families of the four most populated classes which belong to different superfamilies overlap in 23-25% of cases for class a, 11-18% for class b, 10-22% for class c and 11-18% for class d. These bounds on the number of overlapping families can be obtained by using the lower and upper bounds on the distances between representatives and the distances between family members appropriately.

Results for the SCOPCath benchmark

When classifying the 236 queries of SCOPCath, we achieve between 89 and 95% correct super-family assignments, see Table 2. Remarkably, the highest accuracy is reached for k=1, so here just classifying the query as belonging to the superfamily of the nearest neighbor is the best choice. Our k-NN classification resulted for any k in a large number of ties, especially for k=2, see Table 2. These currently unresolved ties also decrease assignment accuracy compared to k = 1, for which a tie is not possible. Table 2 further lists the number of queries which have been assigned, where exact denotes that the provable k nearest neighbors have been computed. The percentage of exactly computed nearest neighbors varies between 50 and 99% and increases with decreasing k. A likely reason for this is that the larger k, the weaker is the k-th distance upper bound that is used for domination, especially if the target on rank k is dissimilar to the query. Since SCOPCath domains have low sequence similarity, this is likely to happen. It is also interesting to note that there are for any k quite a few queries which have been assigned exact, but which are nonetheless wrongly assigned, see Table 2. These are cases in which our distance metric fails in ranking the targets correctly with respect to gold standard.

Results for the extended SCOPCath benchmark

Our exact k-NN classification can also be successfully applied to larger benchmarks like extended SCOPCath, which are more representative of databases such as SCOP. Here, the benefit of using a metric distance, triangle inequality and k-NN classification is more pronounced. Remarkably, our classification run time on this benchmark that is about an order of magnitude larger than SCOPCath is for most queries of the same order of magnitude as run times on SCOPCath (except for some queries which need an extremely long run time and finally cannot be assigned exactly). Also here, run time varies extremely between queries, between 0.15 and 85.63 hours for queries of the four major classes which could be assigned exactly. The median run time for all 1120 exactly assigned extended SCOPCath queries is 3.8 hours. The classification results for extended SCOPCath are shown in Table 2. Slightly more queries have been assigned correctly compared to SCOPCath, between 95 and 99%, and significantly more queries have been assigned exactly. Both may reflect that there are now more similar structures within the targets. Further, the number of ties is decreased. Figure 2 displays the progress of the computation. Here, many more target structures are removed by triangle dominance and within the very first iteration of pairwise dominance compared with the SCOPCath benchmark. For example, for most queries, more than 60% of targets are removed by triangle dominance alone. Only very few queries need to explicitly compute the distance to a large percentage of the targets, and most can be assigned after only one round of pairwise dominance.

Conclusion

In this work we introduced a new distance based on the CMO measure and proved that it is a true metric, which we call the max-CMO metric. We analyzed the potential of max-CMO for solving the k-NN problem efficiently and exactly and built on that basis a protein superfamily classification algorithm. Depending on the value of k, our accuracy varies between 89% for k = 10 and 95% for k = 1 for SCOPCath and between 95 and 99% for extended SCOPCath. The fact that the accuracy is highest for k = 1 indicates that using more sophisticated rules than k-NN may produce even better results.

In summary, our approach provides a general solution to k-NN classification based on a computationally intractable metric for which polynomial upper and lower bounds are available that can successfully be applied for exact large-scale protein superfamily classification.

Fig. 1 .

 1 Fig.1. The alignment visualized with dashed lines ((v1 ↔ u1)(v2 ↔ u2)(v3 ↔ u4)(v4 ↔ u5)) maximizes the number of the common edges between the graphs G1 and G2. The alignment activates four common edges that are emphasized in bold (i.e., CMO(G1, G2) = 4).

Fig. 2 .

 2 Fig. 2. Boxplots of the percentage of removed targets at each iteration during triangle and pairwise dominance for the 1369 queries of the extended SCOPCath benchmark.

Table 1 .

 1 For every protein class, the table lists the number of structures in SCOPCath (str) and extended SCOPCath (ext), the corresponding number of families (fam) and superfamilies (sup).

	class	a	b	c	d e	f	g h i j k
	# str 1195 1593 1774 1591 30 103 342 72 11 38 10
	# ext 10796 19215 17497 15679 349 1006 2398 520 43 81 25
	# fam 524 516 548 632 6 59 121 32 5 29 8
	# sup 303 266 191 375 6 52 82 31 5 29 8

Table 2 .

 2 Classification results showing the number of queries out of overall 236 queries for SCOPCath and 1369 queries for extended SCOPCath that have been assigned to a super-family, the number of assignments to the correct superfamily (cor), the number of assignments computed exactly, i.e. queries which terminate with the provable k nearest neighbors (exc), thereof the number of correct classifications (e&c) and the number of ties which do not allow a superfamily assignment based on majority vote.

	SCOPCath ext. SCOPCath
	k cor exc e&c ties cor exc e&c ties
	10 210 117 110 10 1303 1120 1104 35
	9 211 143 134 9 1331 1182 1166 5
	8 213 156 149 11 1334 1228 1215 12
	7 213 165 155 8 1341 1271 1257 6
	6 214 188 178 10 1341 1286 1276 11
	5 217 206 198 10 1346 1339 1329 7
	4 217 204 195 10 1344 1341 1330 9
	3 219 211 205 10 1351 1352 1341 3
	2 213 209 206 20 1348 1347 1343 17
	1 224 234 224 0 1361 1368 1360 0


Guillaume Chapuis

Mathilde Le Boudic -Jamin

Rumen Andonov

Hristo Djidjev

Dominique Lavenier

Parallel seed-based approach to multiple protein structure similarities detection *

Keywords: protein structure comparison, alternative alignments, alignment graph, maximal clique, Streaming SIMD Extensions (SSE), bit-level parallel computations

Finding similarities between protein structures is a crucial task in molecular biology. Most of the existing tools require proteins to be aligned in order-preserving way and only find single alignments even when multiple similar regions exist. We propose a new seed-based approach that discovers multiple pairs of similar regions. Its computational complexity is polynomial and it comes with a quality guaranteethe returned alignments have both Root Mean Squared Deviations (coordinate-based as well as internal-distances based) lower than a given threshold, if such exist. We do not require the alignments to be order preserving (i.e. we consider non-sequential alignments), which makes our algorithm suitable for detecting similar domains when comparing multi-domain proteins as well as to detect structural repetitions within a single protein. Because the search space for non-sequential alignments is much larger than for sequential ones, the computational burden is addressed by extensive use of parallel computing techniques: a coarse-grain level parallelism making use of available CPU cores for computation and a fine-grain level parallelism exploiting bit-level concurrency as well as vector instructions.

Introduction

A protein's three dimensional structure tends to be evolutionarily better preserved than its sequence. Therefore, finding structural similarities between two proteins can give insights into whether these proteins share a common function or whether they are evolutionarily related. Structural similarity between two proteins is usually defined by two functions -a one-to-one mapping (also called alignment or correspondence [START_REF] Kolodny | Approximate protein structural alignment in polynomial time[END_REF]) between two subchains of their three dimensional representations and a specific scoring function that assesses the alignment quality. The structural alignment problem is to find the mapping that is optimal with respect to the scoring function. Hence, the complexity of the protein structural alignment problem and the quality of the found solution strongly depend on the way that scoring function is defined.

The most commonly used among the various measures of alignment similarity are the internal-distances root mean squared deviation (RM SD d ) and the coordinate root mean squared deviation (RM SD c ) (see [START_REF] Chapuis | Parallel seed-based approach to protein structure similarity detection[END_REF] and (2) respectively for the exact definitions). Tightly related to these measures are the two main approaches for solving the structural alignment problem. The similarity score in the first approach is based on the internal distances matrix, where a set of distances between elements in the first protein is matched with a set of distances in the second protein. The second approach uses the actual Euclidean distances between corresponding atoms in two proteins and aims to determine the rigid transformation that superimposes the two structures.

A huge majority of the algorithms representing these approaches are heuristics [START_REF] Holm | Protein structure comparison by alignment of distance matrices[END_REF][START_REF] Malod-Dognin | Gr-align: fast and flexible alignment of protein 3d structures using graphlet degree similarity[END_REF][START_REF] Shindyalov | Protein structure alignment by incremental combinatorial extension (ce) of the optimal path[END_REF][START_REF] Zemla | LGA: a method for finding 3D similarities in protein structures[END_REF][START_REF] Zhang | TM-align: a protein structure alignment algorithm based on the TM-score[END_REF] (excellent reviews can be found in refs. [START_REF] Koehl | Protein structure similarities[END_REF][START_REF] Liu | A mathematical framework for protein structure comparison[END_REF]) and as such, do not guarantee finding an optimal alignment with respect to any scoring function. The fact that finding exact solutions in this field is computationally hard is related to the fact that computing the longest alignment of protein structures is typically modeled as an NP-hard problem, e.g., the protein threading problem [START_REF] Lathrop | The protein threading problem with sequence amino acid interaction preferences is NP-complete[END_REF], the problem of enumerating all maximal cliques [START_REF] Gibrat | Surprising similarities in structure comparison[END_REF][START_REF] Schmitt | A new method to detect related function among proteins independent of sequence and fold homology[END_REF], or finding a maximum clique [START_REF] Konc | ProBiS algorithm for detection of structurally similar protein binding sites by local structural alignment[END_REF][START_REF] Malod-Dognin | Maximum cliques in protein structure comparison[END_REF][START_REF] Strickland | Optimal protein structure alignment using maximum cliques[END_REF].

These results have been generalized by Kolodny and Linial [START_REF] Kolodny | Approximate protein structural alignment in polynomial time[END_REF], who showed that protein structural alignment is NP-hard if the similarity score is distance based. They also point out that a correct and efficient solution of the structural-alignment problem must exploit the fact that the proteins lie in three-dimensional Euclidean space.

In this paper we present an algorithm that avoids the fundamental in-tractabilities pointed out in [START_REF] Kolodny | Approximate protein structural alignment in polynomial time[END_REF]. Our algorithm is both internal-distances based and Euclidean-coordinates based (i.e., it uses a rigid transformation to optimally superimpose the two structures). Its computational complexity is polynomial and it comes with a quality guarantee -for a given threshold τ , it guarantees to return alignments that have RM SD c as well as RM SD d less than 2τ , if such exist.

Our algorithms is motivated by a class of exact structural-alignments algorithms that look for the largest clique in the so-called product (or alignment) graphs [START_REF] Konc | ProBiS algorithm for detection of structurally similar protein binding sites by local structural alignment[END_REF][START_REF] Malod-Dognin | Maximum cliques in protein structure comparison[END_REF][START_REF] Strickland | Optimal protein structure alignment using maximum cliques[END_REF]. The edges in such graphs encode information about pairs of residues in the two proteins that match based on internal distances between them, namely, if the difference between corresponding distances does not exceed some fixed parameter τ . Then a clique of size k would correspond to subsets of k residues in both proteins that match.

Here, we relax this condition and accept cliques such that edges correspond to matching of similar internal distances up to 2τ . For this relaxed problem, we propose a polynomial algorithm that takes advantage of internal-distance similarities among both proteins to search for an optimal transformation to superimpose their structures. We also replace the goal of finding the largest clique by the one of returning several very dense "nearclique" subgraphs. This choice is strongly justified by the observation that distinct solutions to the structural-alignment problem that are close to the optimum are all equally viable from the biological perspective, and hence are all equally interesting from the computation standpoint [START_REF] Berbalk | Accuracy analysis of multiple structure alignments[END_REF][START_REF] Kolodny | Approximate protein structural alignment in polynomial time[END_REF].

To the best of our knowledge, our tool is unique in its capacity to generate multiple alignments with "guaranteed good" both RM SD c and RM SD d values. We do not require the alignments to be order preserving which makes our algorithm suitable for detecting similar domains when comparing multidomain proteins. Thanks to this property, the tool is able to find both sequential and non-sequential alignments, as well to detect structural repetitions within a single protein and between related proteins.

However, to enumerate exhaustively multiple similar regions requires a more systematic approach than those developed in other existing heuristicbased tools. The computational burden is addressed by extensive use of parallel computing techniques: a coarse grain level parallelism making use of available CPU cores for computation and a fine grain level parallelism exploiting bit-level optimization as well as vector instructions.

Other non-sequential structure alignment methods have been recently proposed (excellent review on this topic can be found in the very recent reference [START_REF] Minami | Mican : a protein structure alignment algorithm that can handle multiple-chains, inverse alignments, c alpha only models, alternative alignments, and nonsequential alignments[END_REF]). None of them is close to the approach proposed here. As they are all heuristic and do not guarantee finding an optimal alignment, a detailed comparison with algorithms based on different concepts requires extensive numerical experiments and is outside the scope of this study.

Here we present a significantly improved and expanded version of a paper originally presented at the PPAM 2013 conference [START_REF] Chapuis | Parallel seed-based approach to protein structure similarity detection[END_REF]. In comparison to [START_REF] Chapuis | Parallel seed-based approach to protein structure similarity detection[END_REF], the current version contains detailed description and explanation of all steps of the algorithms, all pseudo codes, supplementary figures illustrating the algorithms and the experimental results, and extended reference section. Additional sections are added including a comparison between the straightforward and the bit-vector implementations based on complexity analysis as well as detailed analysis of the work from the point of view of future performance improvements and additional possible applications.

Preliminaries

Measures for protein alignments

Consider a protein P of n atoms, P = (p 1 , . . . , a n ), with p i ∈ R 3 . Many measures have been proposed to assess the quality of a protein alignment. These measures include additive scores based on the distance between aligned residues such as the TM-score [START_REF] Zhang | Scoring function for automated assessment of protein structure template quality[END_REF], the DALI score [START_REF] Holm | Searching protein structure databases has come of age[END_REF][START_REF] Wohlers | DALIX: Optimal DALI protein structure alignment[END_REF], the PAUL score [START_REF] Wohlers | PAUL: Protein structural alignment using integer linear programming and lagrangian relaxation[END_REF] and the STRUCTAL score [START_REF] Subbiah | Structural similarity of dnabinding domains of bacteriophage repressors and the globin core[END_REF] and Root Mean Square Deviation (RMSD) based scores, such as RMSD100, SAS and GSAS [START_REF] Kolodny | Comprehensive evaluation of protein structure alignment methods: scoring by geometric measures[END_REF]. Given a set of n deviations S = s 1 , s 2 , ..., s n , its Root Mean Square Deviation is

RM SD(S) = 1 n n i=1 s 2 i . ( 1
)
Two different RMSD measures are used for protein structure comparison. The first one, RM SD c , takes into account deviations consisting of the euclidean distances between matched residues after optimal superposition of the two structures and is defined as

RM SD c (P ) = min P 1 n n i=1 (p i -pi ) 2 . ( 2
)
where P is the image of protein P under a rigid transformation.

The second one, denoted here by RM SD d , takes into account deviations consisting of absolute differences of internal distances within the matched structures. The measured deviations are |d(i, j)-d(k, l)|, for all couples of matching pairs "i ↔ k, j ↔ l." Let M be the latter set and N m , its cardinality. We have that

RM SD d (M ) = 1 N m (ij,kl)∈M (|d(i, j)-d(k, l)| 2 ). (3)

Alignment graphs

An undirected graph G = (V, E) is represented by a set V of vertices and a set E of edges between these vertices. In this paper, we focus on a subset consisting of grid-like graphs, referred to as alignment graphs.

An m × n alignment graph G = (V, E) is a graph in which the vertex set V is depicted by an m × n array T , where each cell T [i][k] contains at most one vertex (i, k) from V . An example of such an alignment graph for protein comparison is given in Figure 1.

A good matching of two proteins P 1 and P 2 can be found by analyzing their alignment graph G = (V, E), where

V = {(v 1 , v 2 ) | v 1 ∈ V 1 , v 2 ∈ V 2 } and V 1 (resp. V 2 )
is the set of atoms of interest in protein P 1 (resp. protein P 2 ). A vertex (I, I ′ ) is present in V only if atoms I ∈ V 1 and I ′ ∈ V 2 are compatible. An example of incompatibility could be different electrostatic properties of the two atoms. An edge ((I, I ′ ), (J, J ′ )) is in E if and only if the distance between atoms I and J in protein P 1 , d(I, J), is similar to the distance between atoms I ′ and J ′ in protein P 2 , d(I ′ , J ′ ). In our case, these distances are considered similar if |d(I, J) -d(I ′ , J ′ )| < τ , where τ is a given threshold.

We arbitrarily order the vertices of the alignment graph and assign a corresponding index to each of them. In subsequent sections, the notion of successors of a vertex v refers to all the vertices that are adjacent to v in the alignment graph and have a higher index than v with respect to those fixed indices. The notion of neighbors of a vertex v refers to the set of all the vertices that are adjacent to v in the alignment graph regardless of their respective indices. Formally, we define Figure 1: Example of an alignment graph used here to compare the structures of two proteins. The presence of an edge between vertex (1, 1) and vertex [START_REF] Chapuis | Parallel seed-based approach to protein structure similarity detection[END_REF][START_REF] Berbalk | Accuracy analysis of multiple structure alignments[END_REF] means that the distance between atoms 1 and 3 of protein 1 is similar to the distance between atoms 1 and 2 of protein 2. The clique (2, 1) (3, 2) (4, 3) indicates that RM SD d of structures [START_REF] Berbalk | Accuracy analysis of multiple structure alignments[END_REF][START_REF] Chapuis | Parallel seed-based approach to protein structure similarity detection[END_REF][START_REF] Dagum | OpenMP: an industry standard API for shared-memory programming[END_REF] and (1, 2, 3) is less than τ [START_REF] Malod-Dognin | Maximum cliques in protein structure comparison[END_REF].

successors(v i ∈ V ) ={v j ∈ V |(v i , v j ) ∈ E & i < j} neighbors(v i ∈ V ) ={v i } ∪ {v j ∈ V |(v i , v j ) ∈
In an alignment graph of two proteins P 1 and P 2 , a subgraph with high density of edges denotes similar regions in both proteins. One can, therefore, find similarities between two proteins by searching in the corresponding alignment graph for subgraphs with high edge density. The highest possible edge density is found in a clique, a subset of vertices that are all connected to each other. Then a clique of size k would correspond to subsets of k residues in both proteins that match since all k(k- 1) of similar internal distances up to 2τ . We also replace the goal of finding the largest clique by the one of returning several very dense "near-clique" subgraphs.

Our method uses geometric properties of the 3-d space. Instead of testing for the presence of all edges among a subset S of vertices, in order to determne if S defines a clique, we only test for the presence of edges between every vertex of the subset and an initial 3-clique (triangle), referred to as seed. This reduction of the performed comparisons is crucial and leads to a polynomial algorithm that takes advantage of internal-distance similarities among both proteins to search for an optimal transformation to superimpose their structures.

Overview of the algorithm

Algorithm 1 gives an overview of our approach. The algorithm consists of the following three steps:

• Seeds in the alignment graph are enumerated. In our case, a seed is a set of three points in the alignment graph that correspond to two triangles (one in each protein) with similar internal distances. This step is detailed in Section 3.3.

• Each seed is then extended. Extending a seed consists in adding all pairs of atoms, for which distances to the seed are similar in both proteins, to the set of three pairs of atoms that make up the seed. Seed extension is detailed in Section 3.4.

• Each seed extension is filtered -cf. lines 10 through 15. Extension filtering is detailed in Section 3.5 and consists in removing pairs of atoms that do not match correctly.

Filtered extensions are then ranked according to their size -number of aligned pairs of atoms -and a user-defined number of best matches are returned. This process is explained in Section 3.7. For very large alignment graphs, the graph can be partitioned into a user-defined number of parts to speed up computations. The graph is partitioned using a min-cut type of heuristic to preserve the quality of the results. Each subgraph is then processed independently. This process is explained in Section 5. The overall worst-case complexity of the whole algorithm without partitioning is

O(|V | * |E| 3/2 ).

Seed enumeration

A seed consists of three pairs of atoms that form similar triangles in both proteins. A triangle IJK in protein P 1 is considered similar to a triangle

I ′ J ′ K ′ in protein P 2 if the following conditions are met: |d(I, J) -d(I ′ , J ′ )| < τ , |d(I, K) -d(I ′ , K ′ )| < τ and |d(J, K) -d(J ′ , K ′ )| < τ
, where d denotes the Euclidean distance and τ is a user-defined threshold parameter. The default value for τ is 2.0 Ångstroms. In the alignment graph terminology, these conditions for a seed

(v i = (I, I ′ ), v j = (J, J ′ ), v k = (K, K ′ )) in graph G(V, E) translate to the following: (v i , v j ) ∈ E, (v i , v k ) ∈ E and (v j , v k ) ∈ E.
A seed thus corresponds to a 3-clique in the alignment graph; i.e., three vertices that are connected to each other. Enumerating all the seeds is therefore equivalent to enumerating every 3-clique in the input alignment graph.

Not all 3-cliques, however, are relevant. Suitable 3-cliques are composed of triangles for which a unique transformation can be found to optimally superimpose them. Namely, 3-cliques composed of triangles that appear to be too "flat" will not yield a useful transformation. We thus ensure that the triangles in both proteins defined by a potential seed are not composed of colinear points (or points which are close to being colinear). The method is detailed in Section 3. The worst-case complexity of this step is O(|E| 3/2 ) using, e.g., the algorithms from [START_REF] Schank | Finding, counting and listing all triangles in large graphs, an experimental study[END_REF].

Algorithm

Seed extension

Extending a seed consists in finding the set of vertices that correspond to pairs of atoms that potentially match well (see Section 3.5 for details) when the two triangles defined by the seed are optimally superimposed. Finding a superset of pairs of atoms that match well is performed by triangulation with the three pairs of atoms composing the seed. Let (v i = (I, I ′ ), v j = (J, J ′ ), v k = (K, K ′ )) be a seed of an alignment graph G(V, E) as defined in Section 3.3. Then the extension of the seed

(v i , v j , v k ) is defined as the set extension(v i , v j , v k ) ={v l = (L, L ′ )| |d(L, I) -d(L ′ , I ′ )| < τ ∧ |d(L, J) -d(L ′ , J ′ )| < τ ∧ |d(L, K) -d(L ′ , K ′ )| < τ }.
In the alignment graph terminology, the previous definition translates to:

extension(v i , v j , v k ) = neighbors(v i ) ∩ neighbors(v j ) ∩ neighbors(v k ).
The detail of seed extension is given in Algorithm 3. The sequential computational complexity associated to this step is O(|V |). For the complexity of our parallel implementation refer to subsection 4.

Extension filtering

The triangulation performed when extending a seed is not sufficient to find alignments with good RM SD measures. Indeed, in most cases, knowing the distance of a point in a three dimensional space to three other nonaligned points yields two possible locations. These locations are symmetrical with respect to the plane defined by the three reference points. A vertex in a seed extension represents a pair of atoms, one in each studied proteins. By construction, these atoms have similar distances to the three points of their respective triangles. It may happen that one of the two points, say L, is located, in protein P 1 , on one side of the plane defined by its reference triangle, while the second point, say L ′ , in protein P 2 , lies on the other side of the plane defined by the two optimally superimposed reference trianglessee Figure 2.

I ′ J ′ K ′ .
Using quadruplets of vertices as seeds does improve the quality of seed extensions but greatly increases the computational cost of seed enumeration and degeneration check on the corresponding tetrahedra. Moreover, larger seeds do not completely ensure the quality of extensions. Namely, in cases where, for a vertex v l = (L, L ′ ), atom L (resp. L ′ ) is very distant from atoms I, J and K (resp. atoms I ′ , J ′ and K ′ ) of a seed (v i = (I, I ′ ), v j = (J, J ′ ), v k = (K, K ′ )), distance similarities to the atoms of the seed do not ensure similar positions of atoms L and L ′ in the two proteins.

In order to remove issues with symmetry (where the atoms in the extending pair are roughly symmetrical with respect to the plane determined by the seed atoms) and distance from the seed, we implement a method to filter seed extensions. This method consists in computing the optimal transformation T to superimpose the triangle from the seed corresponding to the first protein onto the triangle corresponding to the second. The optimal transformation is obtained in constant time with respect to the size of the alignment by using the fast, quaternion-based method of [START_REF] Liu | Fast determination of the optimal rotational matrix for macromolecular superpositions[END_REF]. For each pair of atoms (L,L') composing the extension of a seed (v i = (I,

I ′ ), v j = (J, J ′ ), v k = (K, K ′ )),

Guarantees on resulting alignments' RMSD scores

By construction, the filtering method ensures that the RMSD for a resulting alignment is less than τ : the distance between two aligned residues after superimposition of the two structures is guaranteed to be less than τ . Internal distances between any additional pair of atoms and the seed is also guaranteed, by construction to be less than τ . Concerning internal distances between two additional pairs of atoms, we show that in the worst possible case the difference does not exceed 2τ , see Figure 3. The worst possible case happens when two additional pairs of atoms v l = (L, L ′ ) and v m = (M, M ′ ), added to the extension of a seed (v i , v j , v k ), have atoms L, L ′ , M and M ′ aligned, after superimposition, and atoms from one protein lie within the segment defined by the two other atoms. In such a case, the filtering step ensures that d(L, L ′ ) < τ and d(M, M ′ ) < τ ; it follows that

|d(L, M ) -d(L ′ , M ′ )| < 2τ . v i =(I , I ' ) v j =(J , J ') v k =(K , K ' ) v l = L, L'  v m = M , M ' d  L, L'  d  M , M '

Result ranking

When comparing two proteins, we face a double objective: finding alignments that have both long overlaps and low RM SD scores. The methodology described in Section 3.5 ensures that any returned alignment will have a RM SD d lower or equal to twice the user-defined parameter τ . We can therefore leave the responsibility to the user to define a threshold for RM SD scores of interest. However, ranking alignments that conform to this RM SD threshold simply based on their lengths is not an acceptable solution. In a given alignment graph, several seeds may lead to very similar transformations and thus very similar alignments. The purpose of returning multiple alignments is to find distinct similar regions in both proteins. Therefore, when two alignments are considered similar, we discard the shorter of the two.

Two alignments are considered similar, when they share a defined number of pairs of atoms. This number can be adjusted depending on the expected length of the alignments or even set to a percentage of the smaller of the two compared alignments. This methodology of ranking results ensures that no two returned alignments match the same region in the first protein to the same region in the second protein.

Graph splitting

Large protein alignment graphs can contain millions of edges. In order to reduce the computations induced by such large graphs, a graph splitting scheme is implemented.

Graph splitting is performed using a min-cut like heuristic, also known as multi-level graph partitioning, provided by the METIS library [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF]. This heuristic partitions the graph in k components of similar number of vertices and aims at minimizing the number of inter-component edges -edges between vertices that belong to distinct components. In order to further minimize the number of inter-component edges, we allow the sizes in terms of numbers of vertices of the components to vary up to an order of magnitude. The assumption is that such partitions will keep each area of interest in the graph within a single component.

Once a partition is obtained, subgraphs corresponding to the k components are sorted according to their respective numbers of vertices. Each subgraph is then processed in decreasing order of sizes starting with the largest one. The list of best results is transmitted from one subgraph to another, in order to be able to discard seeds whose extensions are smaller than the best results found so far.

In practice, partitioning the graph tends to group vertices of each of the best results within a single component. However, several of these vertices may be placed in different components. To address this issue, seeds yielding the best results in a subgraph are extended and filtered once more using atoms from the initial global graph.

This second extension and filtering phase significantly improves the length of resulting alignments but does not guarantee to provide the same results as without partitioning. However, experimental results show that a given graph could be partitioned in up to 10 components with only a 2% loss in terms of alignment length and a four fold overall speedup.

The graph splitting scheme is described in Algorithm 5.

Overview of the implemented parallelism

The overall computational complexity of our algorithm being O(|V | * |E| 3/2 ), handling large protein comparison with a decent level of precision -i.e., using alignment graphs with a large number of edges -can prove time-consuming.

Our approach is however parallelizable at multiple levels.

Figure 4 shows an overview of our parallel implementation. Multiple seeds are treated simultaneously to form a coarse-grain level of parallelism, while a finer grain parallelism is used when extending a single seed.

neighbors i 1 ∩ neighbors j 1 ∩ neighbors k 1 seed 1 i 1 , j 1 ,k 1 : Filtering neighbors i 2 ∩ neighbors j 2 ∩ neighbors k 2 seed 2 i 2 , j 2 ,k 2 : Filtering neighbors i n ∩ neighbors j n ∩ neighbors k n seed n i n , j n ,k n : Filtering . . .

Coarse-grain parallelism

Fine-grain parallelism

Coarse-grain parallelism

This level of parallelism is implemented using the open MP standard [START_REF] Dagum | OpenMP: an industry standard API for shared-memory programming[END_REF]. A user-defined number of threads is spawned to handle, in parallel, computations for the seeds generated in the seed enumeration procedure.

In order to output only significant (sufficiently large) alignments we perform here a kind of pruning strategy. It is based on the size of a seed extension which is in fact an easily computed upper bound of the alignment to be produced by the given seed. Whenever this size is smaller than the length of any of the alignments present in the result list (those are the best alignments, so called records), the corresponding seed is discarded, thus avoiding the cost of a filtering step.

In this regard, our pruning strategy is similar to the one applied in a branch-and-bound (B&B) kind of algorithms. A sequential implementation of this strategy is relatively easy since there is only one list of records. The order in which seeds (tasks to be performed) are treated, is crucial for the efficiency of the approach. As sooner as the best alignments are computed, more tasks can be discarded and many computations can be avoided. However, parallelizing the approach induces the same challenges that parallel B&B implementations face [START_REF] Lai | Anomalies in parallel branch-andbound algorithms[END_REF] which can be resumed as follows.

Each thread has now its own (local) list of records which are only lower bounds of the global records found so far by all threads. The pruning strategy is less efficient and could even result in an increase of the amount of computations. One option to avoid this is to make threads share the list of records and to keep it updated. However, inserting frequently new entries in this global-result list, or often updating its content, would prove rather inefficient, because thread safety would need to be ensured by using locks around accesses to this result list. With such locks, threads would often stall whenever inserting a new alignment and the time lost on these accesses would only increase with the number of threads in use. In order to avoid any bottleneck when inserting a new alignment in the result list, in our implementation each thread has its own private list. These lists are merged at the end of the computations to form a final result list. This implementation choice prevent the need for a synchronization mechanism and allow threads to be completely independent.

Fine-grain parallelism

Seed extension makes extensive use of set intersection operations. We define for each vertex v a set set(v) containing all neighbors of v. Then, for each edge (v, w), the set of seeds containing (v, w) can be computed by finding the intersection of set(v) with set(w). In this section, we compare two alternative implementations for computing set intersection.

Straight-forward implementation

The straight-forward implementation of a set of vertices is to store in a sorted array the indices of the vertices that are present in the set. Computing the intersection of two such sets, S 1 and S 2 , therefore consists in traversing both arrays and adding common vertices to the resulting intersection. The complexity of such an approach is thus O(M + N ), where M is the length of S 1 and N is the length of S 2 . Once the intersection has been found, the length of the resulting array is obtained in constant time.

Bit-vector implementation

In order to speed up set intersection operations, we implemented a bit-vector representation of the neighbors set of each vertex of the alignment graph. To any vertex v i , we associate a vector neighbors(v i ) where the bit at position j is 1 if and only if vertexes v i and v j are connected in the alignment graph (see Fig. 5 for illustration).

Word size

(i , j)∈E |V| (i ,k )∉ E j k neighboors ( v i ) : Figure 5: Bit-vector representation of the neighbors of vertex v i in an align- ment graph G(V, E). In this example, v j unlike v k is a neighbor of v i .
This representation of the neighbors sets allows bit parallel computations of set intersection. A simple logic and operation over every word element of the two sets yields the intersection. For faster traversal of the neighbors set a traditional list representation is also kept. This list representation allows easy access to the first and the last elements of the neighbors set. Knowing the first and last elements of the sets allows to restrict the area of interest for intersection operations as described below. Let f i (resp. f j ) be the first nonzero bit in neighbors(v i ) (resp. neighbors(v j )), while l i (resp. l j ) denote the last non-zero bit in neighbors(v i ) (resp. neighbors(v j )). We then apply the logic and operator only over the interval [b, e] where b = max{f i , f j } and e = min{l i , l j }.

Intersection operations also benefit from SSE1 instructions. A number of atomic operations equal to the size of the SSE registers available on the machine (typically 128 or 256) can be computed simultaneously.

However, this bit-vector approach to computing set intersections increases the number of atomic operations to perform. Namely, vertices, which are not neighbors of any of the two vertices for which the intersection is computed, will induce atomic operations; provided such vertices reside in the area of interest. Such vertices would not be considered in the first approach to set intersection.

In order to efficiently compute the size of the intersection in case of a bit-vector implementation, we use a built-in population count instruction (POPCNT) available in SSE4. This operation returns, in constant time, the number of bits set in a single machine word. For architectures without a built-in population count instruction, a slower yet optimized alternative is provided.

Straight-forward versus bit-vector implementations

As mentioned in section 4.3.1, the complexity of an intersection operation with the straight-forward implementation is in O(M + N ), where M and N denote the lengths of the sets. The size of the resulting set is known after performing the intersection operation, i.e. its complexity is also O(M + N ).

With the bit-vector implementation, the complexity of a set intersection operation is in O(length(A)/SSE_SIZE), where A is the area of interest for a given set intersection, length(A) is the number of vertices in the area of interest and SSE_SIZE is the size of the SSE registers available on the machine2 . It is to be noted that the area of interest contains vertices that are present as well as vertices that are absent from a given set. The complexity of computing the size of the resulting set with the optimized implementation is also in O(length(A)/SSE_SIZE) using POPCNT instructions.

Comparing both implementations requires comparing O(M + N ) versus O(length(A)/SSE_SIZE). The later is not straightforward, but it is obvious that for enough dense graphs (such as more of the alignment graphs and especially when they model atoms on proteins surfaces) the values of M and N tend to increase, and bit-vector implementation incline to be faster.

Additionally, the bit-vector implementation induces regular data access patterns, making it a more cache friendly implementation. This property becomes crucial, when considering a future GPU implementation of the algorithm. Given these observations, we directly implemented the bit-vector alternative.

Note that the size of the intersections of the neighbors of two vertices is an upper bound of the cliques than contains these vertices. If this upper bound is less than the size of a previously found clique, any seed containing these vertices will be discarded. These tests are evaluated when considering a new seed for extension and filtering.

Experimental Evaluation

We evaluated our algorithm with respect to accuracy and speed. In order to test the effectiveness of our approach to detect multiple regions of interest, we considered two proteins (PDB IDs 4clna and 2bbma). These proteins are each composed of two similar domains -named A and B (resp. C and D) for the first (resp. second) protein, separated by a flexible bridge (see Fig. 6).

The corresponding alignment graph contains 21904 vertices and 40762390 edges (17% of the total number of possible edges). Existing approaches for global protein structure comparison, such as PAUL [START_REF] Wohlers | PAUL: Protein structural alignment using integer linear programming and lagrangian relaxation[END_REF] and ones based on contact map overlap (CMO) [START_REF] Andonov | Maximum contact map overlap revisited[END_REF] tend to match both proteins integrally, yielding larger alignments but poorer RMSD scores. TM_align [START_REF] Zhang | TM-align: a protein structure alignment algorithm based on the TM-score[END_REF], the reference tool for protein comparison, only matches domain A onto domain C. The four top results of our tool correspond to all four possible combinations of domain matching, (see Fig. 7). Our tool was run using 12 cores of an Intel(R) Xeon(R) CPU E5645 @ 2.40GHz and the distance threshold τ was fixed to 3 Ångströms in the alignment graph. Scores corresponding to these alignments are displayed in Table 1. In order to test the efficiency of our coarse-grain parallel implementation, we compare run times obtained with various numbers of threads on a single instance. The input alignment graph for this instance contains 4378 vertices and 525547 edges. 3 Computations were run using a varying number of cores of an Intel(R) Xeon(R) CPU E5-2670 0 @ 2.60GHz. Table 2 shows run times and speedups with respect to the number of CPU cores. Run times displayed in this table are averages out of 100 runs. The table also indicates the standard deviation of each set of 100 runs.

CMO

Using more threads for computations provides substantial speedups, but also induces different and unpredictable explorations paths of the search space of the alignment graph. Since finding the optimal set of solutions allows us to prune the search space, the order in which seeds of the graph are considered has an impact on the overall performance of the algorithm. With a single thread, the exploration path of the graph is fixed and run times are homogeneous: the standard deviation with a single thread is 0.9% of the average run time. With more than one thread, the exploration path potentially changes at each run and impacts the total run time. This is particularly true for 2 and 4 threads with standard deviations of respectively 7.6% and 8.4% of the average run time. Further increasing the number of threads reduces the unpredictability by increasing the odds of finding optimal solutions early. This behavior is similar to that exhibited in parallel branch and bound algorithms [START_REF] Lai | Anomalies in parallel branch-andbound algorithms[END_REF].

Conclusion and perspectives

In this paper, we introduce a novel approach to find similarities between protein structures. Resulting alignments are guaranteed to score well for both RM SD d and RM SD c , while remaining polynomial. This approach takes advantage of internal distance similarities, described in an alignment graph, to narrow down the search for an optimal transformation to superimpose two substructures of the proteins. We consider two main possible directions for extending the results reported here: i) improving the performance (accelerating the solver); ii) diversifying the application area.

Performance improvement

Though not implemented, an even higher level of parallelism could be considered when graph splitting is performed. Computations for each subgraph are also independent and could therefore be run in parallel. Since a multicore parallelism implementation is already provided, a cluster level parallelism could be implemented. Each subgraph would be sent to a single cluster node using for example using an MPI approach (for Message Passing Interface [START_REF] Gropp | Using MPI-: Portable Parallel Programming with the Message Passing Interface[END_REF]). However, load balancing would be a challenging task due to the limited number of subgraphs that can be generated without a prohibitive loss of accuracy and the difference in terms of numbers of vertices of these subgraphs. Moreover, the total amount of computations would increase if subgraphs were treated in parallel, since the optimal lower bound found in one subgraph could not be used to solve other subgraphs. This issue would also be similar to that observed in parallel branch and bound algorithms and first described in [START_REF] Lai | Anomalies in parallel branch-andbound algorithms[END_REF].

The structure of this algorithm as well as the required operations make it suitable for a GPU implementation, which could speed up the computations. A bit-vector implementation for set intersection operations allows regular data access patterns. These access patterns make set intersection operations more cache-friendly and could thus be efficiently mapped to the GPU paradigm. Moreover, GPUs provide all the necessary bit-level instructions such as population cound and bit scanning. Seed listing and result ranking operations are however too irregular to be efficiently computed on a GPU; therefore, seeds could be listed by the CPU and sent to the GPU for extension and filtering operations. Results could then be transferred back to the CPU for ranking operations.

Diversifying the applications

Detection of structural repeats. A big advantage of the approach is its capacity to find multiple/alternative alignments with good RM SD c and RM SD d property. This allows a deeper analysis of protein structures. One promising perspective is the investigation of repeats inside a protein structure. Structural repeats are common in protein structures [START_REF] Murray | Toward the detection and validation of repeats in protein structure[END_REF][START_REF] Parra | Detecting repetitions and periodicities in proteins by tiling the structural space[END_REF]. However, they are currently unsatisfactorily studied because of the lack of suitable algorithms. Preliminary results show that our approach returns the repetitions in several reliable alignments, so further investigations are in progress.

Combining local alignments into global ones.

The idea here to further analyze the returned multiple local non-overlapping alignments and to combine them in a new global alignment. Such an approach allows to introduce flexibility and non-sequentiality in protein structure alignments. Similar methods already exist such as LGA [START_REF] Zemla | LGA: a method for finding 3D similarities in protein structures[END_REF] or FlexSnap [START_REF] Berbalk | Accuracy analysis of multiple structure alignments[END_REF][START_REF] Lackner | ProSup: a refined tool for protein structure alignment[END_REF][START_REF] Salem | FlexSnap: flexible non-sequential protein structure alignment[END_REF] and we are currently testing our approach on the corresponding datasets.

Algorithm 1

 1 Overview of the algorithm f u n c t i o n f i n d _ a l i g n m e n t s ( graph ) INPUT : graph , an a l i g n m e n t graph between atoms from two p r o t e i n s OUTPUT: r e s L i s t , a l i s t o f th e l a r g e s t d i s t i n c t a l i g n m e n t s found R e s u l t L i s t r e s L i s t = e m p t y _ r e s u l t _ l i s t ( ) S e e d L i s t s e e d s = enumerate_seeds ( graph ) For each s e e d i n s e e d s V e r t e x S e t s e t = extend_seed ( s e e d ) V e r t e x S e t r e s u l t = empty_set ( ) For each v e r t e x i n s e t I f ( i s _ v a l i d ( v e r t e x ) ) r e s u l t . add ( v e r t e x ) End I f r e s L i s t . i n s e r t _ i f _ b e t t e r ( r e s u l t ) End For End For

For

  each v e r t e x 2 i n g e t _ s u c c e s s o r s ( v e r t e x 1 ) For each v e r t e x 3 i n g e t _ s u c c e s s o r s ( v e r t e x 2 ) I f i s _ e d g e ( v e r t e x 1 , v e r t e x 3 ) I f c o l l i n e a r i t y _ c h e c k ( v e r t e x 1 , v e r t e x 2 , v e r t e x 3 ) s e e d L i s t . add ( v e r t e x 1 , v e r t e x 2 , v e r t e x 3 )

3 . 3

 33 AlgorithmSeed extension f u n c t i o n extend_seed ( vertexA , vertexB , vertexC ) INPUT : a s e e d r e p r e s e n t e d by t h r e e v e r t i c e s ( o r p a i r s o f atoms ) from t h e a l i g n m e n t graph OUTPUT: r e s , a s e t o f p a i r s o f atoms t h a t p o t e n t i a l l y match w e l l when atoms from th e s e e d a r e o p t i m a l l y superimposed ; s i z e , t h e s i z e o f th e r e t u r n e d s e t B i n a r y V e r t e x S e t setA = g e t _ n e i g h b o r s ( vertexA ) B i n a r y V e r t e x S e t setB = g e t _ n e i g h b o r s ( vertexB ) B i n a r y V e r t e x S e t setC = g e t _ n e i g h b o r s ( vertexC ) B i n a r y V e r t e x S e t tmp = i n t e r s e c t i o n ( setA , setB ) B i n a r y V e r t e x S e t r e s = i n t e r s e c t i o n ( tmp , vertexC ) int s i z e = pop_count ( r e s )

Figure 2 :

 2 Figure 2: Example of symmetry issues. Even though, vertex v l = (L, L ′ ) belongs to the extension of seed(v i , v j , v k ), points L and L ′ lie on different sides of the plane defined by optimally superimposed triangles IJK and I ′ J ′ K ′ .

Algorithm 4

 4 we compute the Euclidean distance between T (L) and L ′ . If the distance is greater than a given threshold τ , the pair is removed from the extension. The filtering method is detailed in Algorithm 4. The complexity of this step is O(|V |) per seed. Extension filtering algorithm f u n c t i o n f i l t e r _ e x t e n s i o n ( e x t e n s i o n ) INPUT : e x t e n s i o n , a s e t o f p a i r s o f atoms OUTPUT: r e s u l t , a s u b s e t o f th e e x t e n s i o n c o n t a i n i n g o n l y p a i r s o f atoms t h a t match w e l l V e r t e x S e t r e s u l t = empty_set ( ) Matrix t r a n s f o r m a t i o n = g e t _ o p t i m a l _ t r a n s f o r m a t i o n ( s e e d ) For each v e r t e x i n e x t e n s i o n Point L = g e t _ c o o r d i n a t e s _ i n _ f i r s t _ p r o t e i n ( v e r t e x ) Point L_prime = g e t _ c o o r d i n a t e s _ i n _ s e c o n d _ p r o t e i n ( v e r t e x ) Point L_transformed = a p p l y _ t r a n s f o r m a t i o n (L , t r a n s f o r m a t i o n ) F l o a t d i s t a n c e = d i s t ( L_transformed , L_prime ) I f ( d i s t a n c e < t h r e s h o l d ) r e s u l t . i n s e r t ( v e r t e x )

Figure 3 :

 3 Figure3: Illustration of the guarantee on the similarity of internal distances between two pairs of atoms v l = (L, L ′ ) and v m = (M, M ′ ), here represented in yellow, added to a seed (v i , v j , v k ) represented in blue. Dashed lines represent internal distances, the similarity of which is tested in the alignment graph.

Algorithm 5

 5 Graph splitting algorithm f u n c t i o n s p l i t _ a n d _ s o l v e ( globalGraph ) INPUT : globalGraph , an a l i g n m e n t graph between atoms from two p r o t e i n s OUTPUT: g l o b a l R e s , a l i s t o f th e l o n g e s t d i s t i n c t a l i g n m e n t s found i n t h e graph R e s u l t L i s t g l o b a l R e s = e m p t y _ r e s u l t _ l i s t ( ) Graph [ ] subGraphs = s p l i t ( globalGraph ) s o r t ( subgraphs ) For each subGraph i n subGraphs S e e d L i s t b e s t _ s e e d s = e m p t y _ l i s t ( ) S e e d L i s t s e e d s = enumerate_seeds ( subGraph ) For each s e e d i n s e e d s V e r t e x S e t c u r r e n t _ r e s = e x t e n d _ a n d _ f i l t e r ( subGraph , s e e d ) b e s t _ s e e d s . i n s e r t _ i f _ b e t t e r ( s e e d ) End For For each s e e d i n b e s t _ s e e d s V e r t e x S e t c u r r e n t _ r e s = e x t e n d _ a n d _ f i l t e r ( globalGraph , s e e d ) g l o b a l R e s . i n s e r t _ i f _ b e t t e r ( c u r r e n t _ r e s )

Figure 4 :

 4 Figure 4: Overview of the implemented parallelism.

Figure 6 :

 6 Figure 6: These two proteins are both composed of two similar domainsnamed A and B for 4clna (left), and C and D for 2bbma (right). These domains are separated by a a flexible bridge.

Figure 7 :

 7 Figure 7: Visualizations of the results for the comparison of proteins 4clna and 2bbma returned by CMO, PAUL and the four top alignments of our approach.

Fig. 8

 8 Fig. 8 shows run times for graphs with a varying number of edges and the same number of vertices -21904. Computations were run using 12 cores of an Intel(R) Xeon(R) CPU E5645 @ 2.40GHz. Input alignment graphs were all generated from the same two proteins and different parameters to allow a varying number of edges. The diagram shows dependence of the run time on the number of edges consistent with the theoretical O(|E| 2/3 ), where the dependence on |V | is absorbed into the big O factor.

Figure 8 :

 8 Figure 8: Evolution of run times with respect to # of edges in an alignment graph of 21904 vertices.

  E}

		4	4,1		4,3	4,4
	Atoms of	3		3,2	3,3	3,4
	interest in
	protein P 1	2	2,1			2,4	2,5
		1	1,1	1,2		1,4
			1	2	3	4	5
				Atoms of interest in
				protein P	2

2

  Seed enumeration

	f u n c t i o n enumerate_seeds ( graph )
	INPUT : graph , an a l i g n m e n t graph between atoms from two
	p r o t e i n s
	OUTPUT: s e e d L i s t , a l i s t o f s u i t a b l e 3-c l i q u e s ( i . e .

t r i p l e t s o f v e r t i c e s t h a t a r e c o n n e c t e d t o each o t h e r and c o r r e s p o n d t o non-d e g e n e r a t e d t r i a n g l e s i n both p r o t e i n s ) S e e d L i s t s e e d L i s t = empty_seed_list ( ) For each v e r t e x 1 i n graph

Table 1 :

 1 Details of the alignments returned by other tools -columns 2 through 4 -and our method -columns 5 through 8. Best scores are in italics.

			PAUL	TMAlign	AC	BD	AD	BC
	# of aligned residues	148	148	79	72	70	66	64
	% of aligned residues	100	100	53.4	48.7	47.3	44.6	43.2
	RM SDc	14.781	14.781	2.935	2.048	1.731	1 .592	2.210
	RM SD d	10.838	10.838	2.627	1.797	1.475	1 .414	1.770
	TM_score	0.161	0.161	0 .422	0.411	0 .422	0.405	0.358

Table 2 :

 2 Average run times, standard deviation and speedups for varying # of cores.

	# of cores		1	2	4	6	8	10	12	14	16
	Run time (ms)	49744	25873	14269	10315	8555	7425	6668	6189	5679
	Stand. deviation (ms)	458	1969	1192	604	401	284	257	197	193
	Stand. deviation (%)	0.9	7.6	8.4	5.9	4.7	3.8	3.9	3.2	3.4
	Speedup		1	1.9	3.5	4.8	5.8	6.7	7.5	8.0	8.8
				Evolution of run time w.r.t. # of edges
		800
		700
		600
	Run time (s)	200 300 400 500
		100
		0
		0		5000000	1 0000000	15000000		2 0000000
					# of edges

internal distances are similar. For these reason approaches like ProBis[START_REF] Konc | ProBiS algorithm for detection of structurally similar protein binding sites by local structural alignment[END_REF] and DAST[START_REF] Malod-Dognin | Maximum cliques in protein structure comparison[END_REF] aim at finding the maximum clique in an alignment graph.

Methods3.1 Our approachSince finding a maximum clique in a graph is NP-hard and any exact solver faces prohibitively long run times for some instances, here we relax the above definition of clique and accept cliques such that edges correspond to matching

Streaming SIMD Extensions

SSE_SIZE is 128 on the xeon machine

In order to be able to evaluate more accurately the impact of the number of the threads on the computation time, we have chosen a larger instance than the one used for the same experiment in the conference version of this paper[START_REF] Chapuis | Parallel seed-based approach to protein structure similarity detection[END_REF].

Acknowledgment

We are grateful to Noël Malod-Dognin and Frédéric Cazals for motivating discussions in the initial stage of this study.