Proteins can be naturally classified into families of homologous sequences that derive from a common ancestor. The comparison of homologous sequences and the analysis of their phylogenetic relationships provide very useful information regarding the structure, function and evolution of genes. Thanks to the progress of sequencing projects, this comparative approach can now be applied at the whole genome scale in many different taxa, and several databases have been developed to provide a simple access to collections of multiple sequence alignments and phylogenetic trees [1–9]. The building of such phylogenomic databases involves three steps that require important computing resources: 1) compare all proteins to each other to detect sequence similarities, 2) cluster homologous sequences into families (that we will call theclustering step) and 3) compute multiple sequence alignments and phylogenetic trees for each family. With the recent progress of sequencing technologies, there is an urgent need toprepare for the deluge and hence to develop methods able to deal with a huge quantity of sequences. In this paper, we present a new approach for the clustering of homologous sequences, based on single transitive links (single linkage) with alignment coverage constraints and implemented in a software package (calledSiLiX forSIngle LInkage Clustering of Sequences). We model the dataset as asimilarity network where sequences are vertices and similarities are edges [10]. To overcome memory limitations we follow an online framework [11] in which we visit the edges one at a time to update the families dynamically. This approach enables also an incremental procedure where sequences and similarities are added into the dataset so that it would not be necessary to rebuild the families from scratch. Finally, we adopt a divide-and-conquer strategy to deal with the quantity of data [12] and design a parallel algorithm whose theoretical complexity is addressed in this paper.

We evaluated the computational performances and scalability of this method on a very large dataset of more than 3 millions sequences from the HOGENOM phylogenomic database [9]. Our approach presents several advantages over other clustering algorithms: it is extremely fast, it requires only limited memory and can be run on a parallel architecture - which is essential for ensuring its scalability to large datasets.SiLiX outperforms other existing software programs both in terms of speed and memory requirements. Moreover, it allows a satisfying quality of clustering. We discuss the interest ofSiLiX for the clustering of homologous sequences in huge datasets, possibly in combination with other clustering methods.

Modelling

Single linkage and filtering with alignment coverage constraints

The principle of the single-linkage clustering is that if sequence A is considered homologous to sequence B, and B homologous to C, then A, B and C are grouped into the same family, whatever the level of similarity between A and C. The choice of the sequence similarity criteria that is used to infer homology is therefore an essential parameter of the single-linkage clustering approach. Different criteria can be used, separately or in combination (percentage of identity, alignment score or E-value, alignment coveragei.e. percentage of the length of the sequence that is effectively aligned). Then, if a pair of sequences (A, B) does not satisfy the criteria, the pair is not considered for the clustering. The choice of these criteria depends on the goal of the clustering.

The method presented in this paper was motivated by the development of databases of homologous genes (such as HOGENOM or HOVERGEN [9]). The goal of these databases is to allow the study of the evolution of entire proteins considered as a unit, in contrast to databases such as PFAM [13] or PRODOM [14] that aim at studying the domain architecture of proteins. Hence, in HOGENOM, proteins are classified in the same family only if they are homologous over their entire length - or almost. In practice, protein sequences are compared against each other with BLASTP [15]. For each pairwise alignment, the list of High-scoring Segment Pairs (HSPs) is analyzed to exclude HSPs that are not compatible with a global alignment (for details, see [9]). Then, proteins are classified in the same family if the remaining HSPs cover at least a givenpercentage of coverage of the longest protein with apercentage of identity greater or equal to a given threshold (see Figure1). Therefore the first step of the clustering process consists in analyzing pairwise sequence alignments resulting from the all-against-all comparisons (typically a set of alignments obtained with BLAST [15]) in order to obtain a binary information: keeping or excluding pairs whether they meet or not these sequence similarity criteria. This step (that we will refer to as thefiltering step) can be time consuming, but can be easily distributed (see below).

Single linkage clustering with alignment coverage constraints. The four proteins (A, B, C, D) contain some homologous domains (represented by colored boxes). To avoid the clustering in the same family of proteins that do not share any homology (e.g. A and D), pairwise sequence alignments are considered for the clustering only if they cover a minimum threshold of the length of each of the two proteins. This threshold has to be high enough to exclude cases like the alignment (B, C), which would lead to the clustering of A and D.

Full size image

Sequence families are the connected components of the similarity network

Here we consider the second step: given a list of pairs of similar sequences previously positively filtered, group the sequences into families. We define an undirected graphG = (V,E) with the set of verticesV representing sequences and the sets of edgesE representing similarities between these sequences. We definen = |V | andm = |E |. Naturally, finding sequence families consists in computing the connected components ofG. In this paper, we want to address the case of largen andm and we therefore develop a parsimonious approach in terms of memory use. We want to examine the edgesonline [11,16] and avoid storing them into a connectivity matrix. Therefore the classicalDepth-first search algorithm [17] is not adapted. By analogy with external-memory graph algorithms [18], our approach consists in dynamically reducing the connected components into trees. When an edge is examined, we need to execute two operations:find the tree containing each of the two vertices andunion these trees by merging their vertices into a new tree. Consequently, the connected component problem consists in (1) iteratively build a collection of trees representing the connected components of the graphG and (2) transform each resulting tree into a star tree which root is the representative (orleader) of the family. The final formulation of the problem is therefore building a spanning star forestG* = (V,E*).

Using a memory-efficient structure

The connected components ofG actually form a partition ofV into non-overlapping subsets of vertices that we calldisjoint-sets. Initially each vertex is a set by itself. We need to store the information of the partition and be able to update it dynamically. For this purpose, we use thedisjoint-sets data structure [19,20] which is well suited when the graph is discovered edge by edge. This structure allows efficient implementation of the find and union operations by representing each set as a tree. Practically, the forest composed by all the trees is implemented as an arrayparent of sizen. Each elementi of a tree has a parentparent(i) such thatparent(r) =r ifr is the root of the tree. Moreover, it is straightforward and practical to transform each tree into a star tree such that theparent information is a common label for the vertices in a connected component. This will allow to directly retrieve each sequence family by reading theparent information.

Online procedure for a set of similarities

To buildG* from a set of sequence similarities, we develop a two steps procedure. First, we adopt the algorithm calledUnion-Find by Rank with path compression [19,20]. It consists in updating trees of minimal height while discovering the edges of the graphG online. For this purpose, therank of a vertex is basically defined as its height in the tree. Each edge (i,j) is processed as explained in Algorithm 1 (see also Figure2). It is basically based on the FIND function that associates the root of the tree containing a vertex of interest and the PATH COMPRESSION function which connects the vertices in a path to the root of a tree. The time complexity was proved to be in our case almostO(m) [20]. Secondly, we use PATH COMPRESSION for each vertex inO(n) time. This procedure requires the storage ofn parent andn rank values such that the memory requirements areO(n).

An example of the steps involved in the algorithm called Union-Find by rank with path compression [19,20]. Edges (first column, in red) are examined online. The disjoint-sets data structure, represented by trees (third column) and implemented using theparent array (second column), is consequently modified. The two vertices of the current edge of interest are colored in red.

Full size image

Algorithm 1 ADD EDGE(i,j) by UNION- FIND

Function: FIND (i): returns the root of the tree containingi

Function: PATH COMPRESSION(i,r):parent of vertices in the path fromi to the root of the tree containingi are set tor

1:r₁ ← FIND(i);r₂ ← FIND(j)

2:k ← arg max_{l = 1, 2}(rank(r_l ))

3:ifrank(r₁) ==rank(r₂) andr₁ ≠r₂then

4:   rank(r_k)++

5:end if

6: PATH COMPRESSION(i,r_k)

7: PATH COMPRESSION(j,r_k)

Parallelization for multiple sets of similarities

We take advantage of the possibility of exploring series of sets of sequence similarities with a client-server parallel architecture. We assume that it is usually affordable to split a large set intoq sets. For the sake of clarity, we consider here a group ofq processors, which is a reasonable hypothesis in practice. We note that it would also be recommended to have sets of comparable sizes. We adopt a divide-and-conquer strategy where different processors use the previous sequential algorithm to independently obtain a collection of spanning star forests where such that. These subsolutions are successively merged to obtain the final solutionG* [12]. We first design an algorithm to merge two of these forests inO(n) time (see Algorithm 2). It is also based on the disjoint-sets data structure since, for each vertexi, it basically consists in adding in one forest a formal edge betweeni and the root of the tree containingi in the other forest. Then we build a parallel formulation of our approach [21,22] where are obtained with step (1) of the sequential algorithm and iteratively merged (see Algorithm 3). The parallel time complexity can be estimated asO(m/q +nq). We notice that the merge procedure is many orders of magnitude faster than the processing of a single set of similarities. For this reason, we decide not to distribute over the processors the merge procedures that will be consequently performed by the server processor in the order of the availability.

Algorithm 2 MERGE

Function: FIND(i): returns the root of the tree containingi

1:for alli such that FIND(i) ≠i indo

2:   r ← FIND(i) in

3:   ADD EDGE(r,i) in

4:end for

Algorithm 3 ParallelSiLiX

1: each processorr builds with the sequential algorithm

2:ifr > 1 clientthen

3:   MPI_SEND to server processor 1

4:else

5:   fork in 2,...pdo

6:      MPI_RECEIVE among in their order of availability

7:      MERGE

8:   end for

9:   for alli indo

10:      PATH COMPRESSION(i, Find(i))

11:   end for

12:end if

Dealing with partial sequences

Filtering

Because genome sequences are often not 100% complete and hence some genes may overlap with gaps in the genome assembly, it is important to be able to treat somepartial protein sequences (as opposed tocomplete sequences). These partial sequences cannot be classified using the same criteria as the complete ones and are therefore treated separately. In a first step, gene families are built using only complete protein sequences as explained previously. In a second step, partial sequences are added to this classification, using different alignment length thresholds (for details about parameters, see [9]). It is important to note that, if there are several families that meet these alignment coverage criteria, a partial sequence is included in the one with which it shows the strongest similarity score.

Modelling

To allow the treatment of partial sequences, we propose a modified version of our approach. We redefine the previously mentioned graphG = (V_c,E_c) and we define the undirected graphH =G ∪ (V_p,E_p) with two sets of verticesV_c andV_p, the complete and partial sequences respectively, and the set of edgesE_p between complete and partial sequences, each edge inE_p being weighted by the similarity score. We also impose that edges between partial sequences are not allowed. In this case,n_c = |V_c|,n_p = |V_p|,n =n_c +n_p,m_c = |E_c|,m_p = |E_p| andm =m_c +m_p. At this point, we note that sequence families correspond to the connected components of a subgraph ofH obtained by only conserving the edge of maximum weight for each vertex inV_p: this will guarantee that each partial sequence is connected to only one complete sequence and prevent partial sequences to link two connected components. As described previously, the problem consists in building a novel graph that has the following properties:

• H* is a spanning star forest,

• H* is called asemi-bipartite graph, i.e. a graph that can be partitioned into two exclusive and comprehensive parts (V_candV_p) with internal edges (connecting vertices of the same part) only existing within one of the two parts[23]. The particularity is here that edges between the two parts are weighted,

• ∀v ∈V_p,deg(v) = 1.

Online procedure and parallelization

First, it is necessary to insert an additional step between the two steps of the above-mentioned online procedure: build a subset ofE_p by selecting for each vertex the edge of maximal weight, inO(m_p) time. Then we extend the step (2) to all the vertices inV_p for a time complexity inO(n). This procedure runs inO(n) space since it requires the storage ofn parent values. For the parallelized algorithm, we modify the merging of two forests presented in Algorithm 2 to consider vertices ofV_p and once again select edges of maximal weight, such that the overall parallel complexity can be estimated to be inO(m/q+nq).

TheSiLiX software package

All the presented algorithms are implemented into theSiLiX software package which is written in ANSI C++ and uses MPI (Message Passing Interface) and elements of the well-established Boost libraryhttp://www.boost.org.SiLiX can take two kinds of input. First, the user can provide the result file of an all-against-all BLAST search (genomic or protein sequences) in tabular format (option -outfmt 6 in BLAST). In that case,SiLiX performs the filtering step by analyzing BLAST hits to search for pairs of sequences that fulfill similarity criteria (alignment coverage, sequence identity) set by the user to build families. In this mode, partial sequences can be treated separately, as described above. Second, if the user prefers to use other types of criteria for the filtering,SiLiX can simply take as input a list of pairs of sequences IDs and perform the clustering step. Compilation and installation are compliant with the GNU standard procedure. The package is freely available on theSiLiX webpagehttp://lbbe.univ-lyon1.fr/SiLiX. Online documentation and man pages are also available.SiLiX is licensed under the General Public Licensehttp://www.gnu.org/licenses/licenses.html.

SiLiX is faster and more memory efficient than other methods

To testSiLiX and compare it to state-of-the-art programs, we extracted protein sequences from the HOGENOM database (Release 5, [9]). The current release of HOGENOM contains 3,666,568 protein sequences (76% bacteria, 3% archae and 20% eukarya). We selected 3,159,593 non-redundant sequences including about 1% partial sequences. Sequences were compared against each others with BLASTP [15] with an E-value threshold set to 10^-4. The BLAST output file contained 1,905,335,339 pairwise alignments. Then we selected three previously published programs, for which the source code is publicly available:hcluster_sg [24] andMC-UPGMA [25] that are based on hierarchical clustering, andMCL [26] that relies on graph-based heuristics.

The clustering of the protein dataset withSiLiX was very fast (about 2 hours) and required only limited memory capacity (0.4 GB).SiLiX outperformed the 3 other methods, both in terms of speed and RAM usage (see Table1). The programhcluster_sg took 40 times more time thanSiLiX to perform the clustering (about 4 days), and required a very large amount of RAM memory (99 GB). With larger sequence datasets (which are already present in databases), the RAM requirements ofhcluster_sg will certainly exceed computer memory resources presently available.MCL also required a large amount of memory (78 GB) and was very slow (we stopped it after 10 days of calculation, before it finished the clustering).MC-UPGMA is almost as efficient asSiLiX in terms of RAM usage, but requires ample disk space to hold intermediate files (49 GB of HDD). The main problem withMC-UPGMA is that it is too slow on such a large dataset.MC-UPGMA uses an iterative procedure to cluster sequences. On our dataset, the first 20 iterations took 2 days. The authors ofMC-UPGMA tested their method on a smaller dataset and they indicate in their article that 200 iterations were necessary to reach convergence (see [25]). We therefore extrapolated thatMC-UPGMA would take more than 20 days on our dataset, and hence we decided to stop it before it finished the clustering. Note also that to optimize the performance ofMC-UPGMA, the authors recommend using very permissive similarity threshold (E-value = 100, see [25]), which is unaffordable given the number of sequences in our dataset. In conclusion,SiLiX presents the best efficacy to tackle the challenge of huge dataset analysis with CPU and memory requirements equivalent to those of a laptop computer. Note thatSiLiX may also be used in combination with other methods. To test this strategy, we first ranSiLiX with permissive similarity thresholds (sequence identity 25% and alignment coverage 80%), and then we usedhcluster_sg to subdivide families of more than 100 sequences. This combined procedure still runs in a reasonable total CPU time (about 9 hours, i.e. 10 times faster thanhcluster_sg alone, see Table1) and also divides by 4 the RAM usage. Furthermore the second step of this combined procedure can easily be distributed on several computers.

SiLiX is scalable in practice

As the number of available sequences increases dramatically and the number of similarities is quadratic with this number of sequences, the CPU time required for the clustering is expected to increase very rapidly. To ensure the scalability of our method, we designed a parallel implementation ofSiLiX with a low number of inter-processors communications to take advantage of multiple kinds of parallel hardware architectures. This algorithm delocalizes the processing of the sequence similarity dataset, including the filtering step, and merges the results in a last step (see Methods). We designed a divide-and-conquer approach that requires onlyq - 1 communications whereq is the number of processors, with a procedure for merging partial results from two processors that is considerably faster than the independent computations on each processor. For these reasons, we observe practical performances consistent with the theoretical complexity such that the run time decrease is inversely-proportional to the number of processors (see Figure3).

CPU time of the parallelized version ofSiLiX(plain) according to the number of processors on the dataset of similarity pairs extracted from the HOGENOM database [9],compared with theoretical values (dashed). Run on a cluster of 2 octo-bicore Opteron 2.8 Ghz and 2 octo-quadcore Opteron 2.3 GHz.

Full size image

Clustering quality

Although the speed and memory requirements are important parameters for the choice of clustering method, the most important criterion is of course the quality of the results. Single linkage clustering is known to be problematic because spurious similarities can lead to the clustering of non-homologous sequences. Even with stringent sequence similarity criteria, single linkage clustering can lead to erroneous clustering, because of the so-called problem of "domain chaining" [27], as illustrated in Figure1. To avoid this problem,SiLiX performs single linkage clustering with alignment coverage constraints, i.e. pairs of similar sequences are considered for the clustering only if they meet two criteria: i) the alignment should cover at least a given percentage of the longest sequence; ii) sequence similarity within the alignment should exceed a given threshold. To assess the quality ofSiLiX clustering, we used 2 different strategies. First, we compared clustering results to the classification of protein families reported in the InterPro database [28]. Second we assessed the performance ofSiLiX on a set of 13 families of orthologous genes encoded by mitochondrial genomes in 1821 metazoan species.

Evaluation ofSiLiX performances with Interpro

We evaluated the performance ofSiLiX on the HOGENOM dataset, using the procedure proposed by Loewenstein and colleagues [25]: we extracted the most frequent correspondence betweenSiLiX families and protein (not domain) families from InterPro (Release 22, [28]) containing more than 10 sequences, and then we calculated the specificity and sensitivity ofSiLiX classification with respect to the InterPro family. We also computed the Jaccard score (see [25]), which is a standard metric of the trade-off between specificity and sensitivity. To evaluate the impact of sequence similarity criteria, we ranSiLiX with different thresholds for alignment coverage and percentage of identity. The performance ofSiLiX were compared to those obtained withhcluster_sg, used alone or in combination withSiLiX. (NB: the 2 other methods,MC-UPGMA andMCL, could not be evaluated because of their excessive running time). As expected, increasing the alignment coverage and/or the sequence identity thresholds leads to increase specificity (Figure4a), but decreases sensitivity (Figure4b). The best trade-off between specificity and sensitivity was obtained for thresholds of 80% for alignment coverage and 35% for sequence identity (Figure4c). With these parameters, the Jaccard score ofSiLiX is slightly better than that ofhcluster_sg (Table2). Interestingly, the use ofhcluster_sg in combination withSiLiX (with permissive threshold) leads to a better Jaccard score thanhcluster_sg alone. Thus the use ofSiLiX in combination with other methods can strongly decrease computing time without any loss in clustering quality.

Clustering performance evaluation based on InterPro classification. a) specificity, b) sensitivity and c) Jaccard coefficient ofSiLiX, used on similarity pairs extracted from the HOGENOM database, with different values of threshold on the percentage of sequence identity and alignment coverage.

Full size image

Evaluation ofSiLiX performances with metazoan mitochondrial gene families

We used InterPro to evaluate clustering performances because this database is widely recognized for its quality and it has already been used for that purpose [25]. This strategy may however not be optimal because the building of the InterPro database also relies on arbitrary sequence similarity criteria. Hence, some cases that were considered as false positives in the above evaluations might in fact correspond to true homologues (i.e. specificity would be underestimated). Ideally, to evaluate clustering quality, one would need a set of homologous gene families known a priori,i.e. identified without using sequence similarity criteria. The mitochondrial genome of metazoan taxa can provide such an ideal test set. Indeed, in animals the mitochondrial genome contains 13 protein-coding genes. These proteins show different levels of sequence conservation across taxa, but the gene content is extremely conserved: except in very rare cases, all metazoan mitochondrial genomes contain these 13 genes [29,30]. This very strong conservation of synteny allows the identification of orthologous genes, even with very low levels of sequence similarity. We extracted from RefSeq (Release 41, [31]) a set of complete mitochondrial genomes from 1821 different species. The 13 mitochondrial proteins are present in all taxa, except ATP8 that is missing in the genome of 6 species. These mitochondrial proteins were then added to the HOGENOM dataset. Sequence comparisons were performed with BLASTP (using parameters indicated above), and the clustering was performed withSiLiX on the entire dataset, using thresholds of 80% for alignment coverage and 35% for sequence identity. For each of the 13 gene families, the majority of sequences were grouped in a singleSiLiX family. The sensitivity, measured as the number of known proteins that are included in this firstSiLiX family, is generally very high: it is higher than 94% for 11 out of the 13 gene families, and even higher than 99% for 9 of them (Table3). The ND6 family was split into 3 mainSiLiX families corresponding respectively to deuterostomia, protostomia and other metazoa (porifera, placozoa and cnidaria), and containing overall 95% of all known ND6 proteins. The ATP8 protein is very short (about 50 amino-acids) and evolves rapidly. The largestSiLiX family contains only 53% of all known ATP8, and the three largestSiLiX families contain 77% of known sequences, whereas 14% of known ATP8 were not included in anySiLiX family - in most cases because they did not have any BLAST hit with a E-value below 10^-4. Thus, on this test set the sensitivity is generally very good, except for rapidly evolving sequences. To evaluate the specificity, we manually investigated allSiLiX families containing at least one protein of our mitochondrial set. We did not identify a single case where one mitochondrial protein was clustered with non-related proteins. Thus, even though the clustering was performed with the entire HOGENOM dataset (which contained more than 3 million different proteins), we did not detect any false positive clustering. This indicates that single linkage clustering with alignment coverage constraints is robust to spurious similarity matches.

Different methods have been proposed for the clustering of proteins into families of homologous sequences [1,8,9,24–26,32]. These methods differ both in terms of the quality of the clustering, and in terms of the computing resources required to perform the clustering. The single-linkage clustering approach is used in different phylogenomic databases such as EnsemblCompara [8] or HOGENOM [9]. Here we propose a new implementation of the single linkage clustering method with alignment coverage constraints,SiLiX, which is extremely efficient, both in terms of computing time and memory requirements. Moreover, this method can be cost-effectively run on parallel architectures, and hence is easily scalable. Thus, in terms of the computing resource requirements, this method is much more efficient than other available methods for the treatment of huge sequence datasets. In terms of clustering quality,SiLiX performs as well ashcluster_sg, the only other available clustering program that could be run in reasonable time on such a large sequence dataset. Given its speed,SiLiX may also efficiently be used as a first clustering step, before running other algorithms.

• Project name:SiLiX

• Project home page:http://lbbe.univ-lyon1.fr/SiLiX

• Operating system(s): All Unix-like operating systems such as Linux and Mac OS X.

• Programming language: C++

• Other requirements: MPI, the Boost:program options class, and optionally CppUnit and the Boost:unordered_map class.

• License: GNU GPL.