Ultrafast sequence clustering from similarity networks with SiLiX
 Vincent Miele^{1}Email author,
 Simon Penel^{1} and
 Laurent Duret^{1}
DOI: 10.1186/1471210512116
© Miele et al; licensee BioMed Central Ltd. 2011
Received: 28 October 2010
Accepted: 22 April 2011
Published: 22 April 2011
Abstract
Background
The number of gene sequences that are available for comparative genomics approaches is increasing extremely quickly. A current challenge is to be able to handle this huge amount of sequences in order to build families of homologous sequences in a reasonable time.
Results
We present the software package SiLiX that implements a novel method which reconsiders single linkage clustering with a graph theoretical approach. A parallel version of the algorithms is also presented. As a demonstration of the ability of our software, we clustered more than 3 millions sequences from about 2 billion BLAST hits in 7 minutes, with a high clustering quality, both in terms of sensitivity and specificity.
Conclusions
Comparing stateoftheart software, SiLiX presents the best uptodate capabilities to face the problem of clustering large collections of sequences. SiLiX is freely available at http://lbbe.univlyon1.fr/SiLiX.
Background
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 the clustering 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 to prepare 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 (called SiLiX for SIngle LInkage Clustering of Sequences). We model the dataset as a similarity 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 divideandconquer 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 of SiLiX for the clustering of homologous sequences in huge datasets, possibly in combination with other clustering methods.
Implementation
Modelling
Single linkage and filtering with alignment coverage constraints
The principle of the singlelinkage 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 singlelinkage clustering approach. Different criteria can be used, separately or in combination (percentage of identity, alignment score or Evalue, alignment coverage i.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.
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 graph G = (V, E) with the set of vertices V representing sequences and the sets of edges E representing similarities between these sequences. We define n = V  and m = E . Naturally, finding sequence families consists in computing the connected components of G. In this paper, we want to address the case of large n and m and we therefore develop a parsimonious approach in terms of memory use. We want to examine the edges online [11, 16] and avoid storing them into a connectivity matrix. Therefore the classical Depthfirst search algorithm [17] is not adapted. By analogy with externalmemory 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 and union 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 graph G and (2) transform each resulting tree into a star tree which root is the representative (or leader) of the family. The final formulation of the problem is therefore building a spanning star forest G* = (V, E*).
Using a memoryefficient structure
The connected components of G actually form a partition of V into nonoverlapping subsets of vertices that we call disjointsets. 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 the disjointsets 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 array parent of size n. Each element i of a tree has a parent parent(i) such that parent(r) = r if r is the root of the tree. Moreover, it is straightforward and practical to transform each tree into a star tree such that the parent information is a common label for the vertices in a connected component. This will allow to directly retrieve each sequence family by reading the parent information.
Online procedure for a set of similarities
Algorithm 1 ADD EDGE(i, j) by UNION FIND
Function: FIND (i): returns the root of the tree containing i
Function: PATH COMPRESSION(i, r): parent of vertices in the path from i to the root of the tree containing i are set to r
1: r_{1} ← FIND(i); r_{2} ← FIND(j)
2: k ← arg max_{l = 1, 2}(rank(r_{ l } ))
3: if rank(r_{1}) == rank(r_{2}) and r_{1} ≠ r_{2} 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 clientserver parallel architecture. We assume that it is usually affordable to split a large set into q sets. For the sake of clarity, we consider here a group of q 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 divideandconquer 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 solution G* [12]. We first design an algorithm to merge two of these forests in O(n) time (see Algorithm 2). It is also based on the disjointsets data structure since, for each vertex i, it basically consists in adding in one forest a formal edge between i and the root of the tree containing i 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 as O(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 containing i
1: for all i such that FIND(i) ≠ i in do
2: r ← FIND(i) in
3: ADD EDGE(r, i) in
4: end for
Algorithm 3 Parallel SiLiX
1: each processor r builds with the sequential algorithm
2: if r > 1 client then
3: MPI_SEND to server processor 1
4: else
5: for k in 2,...p do
6: MPI_RECEIVE among in their order of availability
7: MERGE
8: end for
9: for all i in do
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 some partial protein sequences (as opposed to complete 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 graph G = (V_{ c }, E_{ c }) and we define the undirected graph H = G ∪ (V_{ p }, E_{ p }) with two sets of vertices V_{ c } and V_{ p }, the complete and partial sequences respectively, and the set of edges E_{ p } between complete and partial sequences, each edge in E_{ 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 } and m = m_{ c } + m_{ p }. At this point, we note that sequence families correspond to the connected components of a subgraph of H obtained by only conserving the edge of maximum weight for each vertex in V_{ 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 a semibipartite graph, i.e. a graph that can be partitioned into two exclusive and comprehensive parts (V_{ c }and V_{ 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 abovementioned online procedure: build a subset of E_{ p } by selecting for each vertex the edge of maximal weight, in O(m_{ p }) time. Then we extend the step (2) to all the vertices in V_{ p } for a time complexity in O(n). This procedure runs in O(n) space since it requires the storage of n parent values. For the parallelized algorithm, we modify the merging of two forests presented in Algorithm 2 to consider vertices of V_{ p } and once again select edges of maximal weight, such that the overall parallel complexity can be estimated to be in O(m/q+ nq).
The SiLiX software package
All the presented algorithms are implemented into the SiLiX software package which is written in ANSI C++ and uses MPI (Message Passing Interface) and elements of the wellestablished Boost library http://www.boost.org. SiLiX can take two kinds of input. First, the user can provide the result file of an allagainstall 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 the SiLiX webpage http://lbbe.univlyon1.fr/SiLiX. Online documentation and man pages are also available. SiLiX is licensed under the General Public License http://www.gnu.org/licenses/licenses.html.
Results and Discussion
SiLiX is faster and more memory efficient than other methods
To test SiLiX and compare it to stateoftheart 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 nonredundant sequences including about 1% partial sequences. Sequences were compared against each others with BLASTP [15] with an Evalue 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] and MCUPGMA [25] that are based on hierarchical clustering, and MCL [26] that relies on graphbased heuristics.
CPU time and memory requirements for SiLiX and three stateoftheart programs on the dataset of similarity pairs extracted from the HOGENOM database [9].
SiLiX is scalable in practice
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 nonhomologous sequences. Even with stringent sequence similarity criteria, single linkage clustering can lead to erroneous clustering, because of the socalled problem of "domain chaining" [27], as illustrated in Figure 1. 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 of SiLiX 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 of SiLiX on a set of 13 families of orthologous genes encoded by mitochondrial genomes in 1821 metazoan species.
Evaluation of SiLiX performances with Interpro
Comparison of clustering performances of SiLiX and hcluster_sg (used alone or in combination with SiLiX).
method (%identity)  Jac.  Spec.  Sens. 

SiLiX (0.25)  0.69  0.73  0.95 
SiLiX (0.35)  0.78  0.88  0.89 
SiLiX (0.45)  0.74  0.94  0.78 
SiLiX (0.25) + hcluster_sg^{(100)}  0.77  0.85  0.92 
hcluster_sg  0.76  0.84  0.91 
Evaluation of SiLiX performances with metazoan mitochondrial gene families
Evaluation of SiLiX performances on mitochondrial genes of metazoan taxa.
Gene  Nb. Seq.  Nb. SiLiX families  Nb. Seq. 1^{ st } fam. (%)  Nb. Seq. 2^{ nd } fam. (%)  Nb. Seq. 3^{ rd } fam. (%)  Nb. Singletons (%) 

ATP8  1815  26  959 (52.8)  294 (16.2)  144 (7.9)  258 (14.2) 
ATP6  1821  2  1814 (99.6)  2 (0.1)    5 (0.3) 
COX1  1821  1  1820 (99.9)      1 (0.1) 
COX2  1821  1  1818 (99.8)      3 (0.2) 
COX3  1821  1  1821 (100)       
CYTB  1821  1  1820 (99.9)      1 (0.1) 
ND1  1821  1  1821 (100)       
ND2  1821  11  1714 (94.1)  53 (2.9)  3 (0.2)  34 (1.9) 
ND3  1821  2  1813 (99.6)  2 (0.1)    6 (0.3) 
ND4  1821  2  1812 (99.5)  2 (0.1)    7 (0.4) 
ND4L  1821  7  1758 (96.5)  4 (0.2)  3 (0.2)  47 (2.6) 
ND5  1821  2  1815 (99.7)  2 (0.1)    4 (0.2) 
ND6  1821  16  1366 (75.0)  313 (17.2)  55 (3.0)  45 (2.5) 
Conclusion
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 singlelinkage 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 costeffectively 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 as hcluster_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.
Availability and requirements

Project name: SiLiX

Project home page: http://lbbe.univlyon1.fr/SiLiX

Operating system(s): All Unixlike 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.
Declarations
Acknowledgements and Funding
The authors would like to thank Bastien Boussau, Daniel Kahn, Vincent Lacroix, MarieFrance Sagot, Franck Picard and Eric Tannier for helpful discussions and comments, Bruno Spataro for the computing facilities and Yanniv Loewenstein, Jan Baumbach and Antje Krause for their answers about the availability and use of their programs. This work has been supported by the French Agence Nationale de la Recherche under grant NeMo ANR08BLAN030401.
Authors’ Affiliations
References
 Petryszak R, Kretschmann E, Wieser D, Apweiler R: The predictive power of the CluSTr database. Bioinformatics 2005, 21: 3604–3609. 10.1093/bioinformatics/bti542View ArticlePubMed
 Meinel T, Krause A, Luz H, Vingron M, Staub E: The SYSTERS Protein Family Database in 2005. Nucleic Acids Res 2005, 33: 226–229. 10.1093/nar/gki471View Article
 Dehal PS, Boore JL: A phylogenomic gene cluster resource: the Phylogenetically Inferred Groups (PhIGs) database. BMC Bioinformatics 2006., 7:
 Li H, Coghlan A, Ruan J, Coin LJ, Heriche JK, Osmotherly L, Li R, Liu T, Zhang Z, Bolund L, Wong GK, Zheng W, Dehal P, Wang J, Durbin R: TreeFam: a curated database of phylogenetic trees of animal gene families. Nucleic Acids Res 2006, 34: 572–580. 10.1093/nar/gkj118View Article
 Hartmann S, Lu D, Phillips J, Vision TJ: Phytome: a platform for plant comparative genomics. Nucleic Acids Res 2006, 34: 724–730. 10.1093/nar/gkj045View Article
 Tian Y, Dickerman AW: GeneTrees: a phylogenomics resource for prokaryotes. Nucleic Acids Res 2007, 35: 328–331. 10.1093/nar/gkl905View Article
 Wall PK, LeebensMack J, Muller KF, Field D, Altman NS, dePamphilis CW: PlantTribes: a gene and gene family resource for comparative genomics in plants. Nucleic Acids Res 2008, 36: 970–976.View Article
 Vilella AJ, Severin J, UretaVidal A, Heng L, Durbin R, Birney E: EnsemblCompara GeneTrees: Complete, duplicationaware phylogenetic trees in vertebrates. Genome Res 2009, 19: 327–335.PubMed CentralView ArticlePubMed
 Penel S, Arigon AM, Dufayard JF, Sertier AS, Daubin V, Duret L, Gouy M, Perriere G: Databases of homologous gene families for comparative genomics. BMC Bioinformatics 2009, 10(Suppl 6):S3. 10.1186/1471210510S6S3PubMed CentralView ArticlePubMed
 Atkinson HJ, Morris JH, Ferrin TE, Babbitt PC: Using sequence similarity networks for visualization of relationships across diverse protein superfamilies. PLoS ONE 2009, 4: e4345. 10.1371/journal.pone.0004345PubMed CentralView ArticlePubMed
 Fiat A, Woeginger JG, (Eds): Online Algorithms, The State of the Art. Volume 1442. Lecture Notes in Computer Science, Springer; 1998.
 Das SK, Narsingh D: Divide and conquerbased optimal parallel algorithms for some graph problems on EREW PRAM model. IEEE transactions on circuits and systems 1988, 35(3):312–322. 10.1109/31.1744View Article
 Finn RD, Mistry J, Tate J, Coggill P, Heger A, Pollington JE, Gavin OL, Gunasekaran P, Ceric G, Forslund K, Holm L, Sonnhammer EL, Eddy SR, Bateman A: The Pfam protein families database. Nucleic Acids Res 2010, 38: D211–222. 10.1093/nar/gkp985PubMed CentralView ArticlePubMed
 Servant F, Bru C, Carrere S, Courcelle E, Gouzy J, Peyruc D, Kahn D: ProDom: automated clustering of homologous domains. Brief Bioinformatics 2002, 3: 246–251. 10.1093/bib/3.3.246View ArticlePubMed
 Altschul SF, Madden TL, Schaffer AA, Zhang J, Zhang Z, Miller W, Lipman DJ: Gapped BLAST and PSIBLAST: a new generation of protein database search programs. Nucleic Acids Res 1997, 25: 3389–3402. 10.1093/nar/25.17.3389PubMed CentralView ArticlePubMed
 Vishwanathan S: Randomized online graph coloring. Journal of algorithms 1992, 13: 657–669. 10.1016/01966774(92)90061GView Article
 Tarjan R: Depthfirst search and linear graph algorithms. SIAM Journal on Computing 1972, 1(2):146–160. 10.1137/0201010View Article
 Katriel I, Meyer U: Elementary Graph Algorithms in External Memory. Algorithms for Memory Hierarchies 2002, 62–84.
 Tarjan RE: Efficiency of a Good But Not Linear Set Union Algorithm. Journal of the ACM 1975, 22(2):215–225. 10.1145/321879.321884View Article
 Alsuwaiyel MH: Algorithms: Design Techniques and Analysis. World Scientific Publishing Company; 1998.
 Krishnamurthy A, Lumetta SS, Culler DE, Yelick K: Connected Components on Distributed Memory Machines. Parallel Algorithms, DIMACS Series in Discrete Mathematics and Theoretical Computer Science 1997.
 Han Y, Wagner RA: An efficient and fast parallelconnected component algorithm. Journal of the ACM 1990, 37(3):626–642. 10.1145/79147.214077View Article
 Bramoulle Y, LopezPintado D, Goyal S, VegaRedondo F: Network formation and anticoordination games. International Journal of Game Theory 2004, 33(1):1–19. 10.1007/s001820400178View Article
 Ruan J, Li H, Chen Z, Coghlan A, Coin LJ, Guo Y, Heriche JK, Hu Y, Kristiansen K, Li R, Liu T, Moses A, Qin J, Vang S, Vilella AJ, UretaVidal A, Bolund L, Wang J, Durbin R: TreeFam: 2008 Update. Nucleic Acids Res 2008, 36: D735–740.PubMed CentralView ArticlePubMed
 Loewenstein Y, Portugaly E, Fromer M, Linial M: Efficient algorithms for accurate hierarchical clustering of huge datasets: tackling the entire protein space. Bioinformatics 2008, 24: i41–49. 10.1093/bioinformatics/btn174PubMed CentralView ArticlePubMed
 Enright AJ, Van Dongen S, Ouzounis CA: An efficient algorithm for largescale detection of protein families. Nucleic Acids Res 2002, 30: 1575–1584. 10.1093/nar/30.7.1575PubMed CentralView ArticlePubMed
 Joseph JM, Durand D: Family classification without domain chaining. Bioinformatics 2009, 25: 45–53. 10.1093/bioinformatics/btp207View Article
 Hunter S, Apweiler R, Attwood TK, Bairoch A, Bateman A, Binns D, Bork P, Das U, Daugherty L, Duquenne L, Finn RD, Gough J, Haft D, Hulo N, Kahn D, Kelly E, Laugraud A, Letunic I, Lonsdale D, Lopez R, Madera M, Maslen J, McAnulla C, McDowall J, Mistry J, Mitchell A, Mulder N, Natale D, Orengo C, Quinn AF, Selengut JD, Sigrist CJ, Thimma M, Thomas PD, Valentin F, Wilson D, Wu CH, Yeats C: InterPro: the integrative protein signature database. Nucleic Acids Res 2009, 37: D211–215. 10.1093/nar/gkn785PubMed CentralView ArticlePubMed
 Boore JL: Animal mitochondrial genomes. Nucleic Acids Res 1999, 27: 1767–1780. 10.1093/nar/27.8.1767PubMed CentralView ArticlePubMed
 Signorovitch AY, Buss LW, Dellaporta SL: Comparative genomics of large mitochondria in placozoans. PLoS Genet 2007, 3: e13. 10.1371/journal.pgen.0030013PubMed CentralView ArticlePubMed
 Pruitt KD, Tatusova T, Klimke W, Maglott DR: NCBI Reference Sequences: current status, policy and new initiatives. Nucleic Acids Res 2009, 37: D32–36. 10.1093/nar/gkn721PubMed CentralView ArticlePubMed
 Wittkop T, Baumbach J, Lobo FP, Rahmann S: Large scale clustering of protein sequences with FORCE A layout based heuristic for weighted cluster editing. BMC Bioinformatics 2007, 8: 396. 10.1186/147121058396PubMed CentralView ArticlePubMed
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.