Parallel Seed-Based Approach to Multiple Protein Structure Similarities Detection

,


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 [12]) 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 (3) 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 [7,20,27,32,33] (excellent reviews can be found in refs.[10,18]) 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 [16], the problem of enumerating all maximal cliques [5,26], or finding a maximum clique [13,19,28].
These results have been generalized by Kolodny and Linial [12], 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 [12].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 [13,19,28].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 [2,12].
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 [21]).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 [3].In comparison to [3], 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.

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 [34], the DALI score [8,30], the PAUL score [31] and the STRUCTAL score [29] and Root Mean Square Deviation (RMSD) based scores, such as RMSD100, SAS and GSAS [11].Given a set of n deviations S = s 1 , s 2 , ..., s n , its Root Mean Square Deviation is 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 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

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 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 (3,2) 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 (2,3,4) and (1, 2, 3) is less than τ [19].
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) 2 internal distances are similar.For these reason approaches like ProBis [13] and DAST [19] aim at finding the maximum clique in an alignment graph.

Our approach
Since 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 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 Algorithm 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 End For End For

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: , 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 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 [25].

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 In the alignment graph terminology, the previous definition translates to: 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.3.

Algorithm 3 Seed 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 he 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 he 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 )

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. 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 [17].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 )), 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 End I f End For

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

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 [9].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.

Coarse-grain parallelism
This level of parallelism is implemented using the open MP standard [4].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 [15] 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
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 [31] and ones based on contact map overlap (CMO) [1] tend to match both proteins integrally, yielding larger alignments but poorer RMSD scores.TM_align [33], 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. 3Computations 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 [15].

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 [6]).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 [15].
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 [22,23].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 [32] or FlexSnap [2,14,24] and we are currently testing our approach on the corresponding datasets.

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 )

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 )

Figure 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 .

4 .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 :
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
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 he 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 :
Figure 4: Overview of the implemented parallelism.

Figure 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 :
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
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 :
Figure 8: Evolution of run times with respect to # of edges in an alignment graph of 21904 vertices. E}

2
Seed enumeration 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 :
Details of the alignments returned by other tools -columns 2 through 4 -and our method -columns 5 through 8. Best scores are in italics.

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