BACKGROUND
the amount of biological information being gathered is growing faster than the rate at which it can be analyzed.
data clustering, which compresses the problem space by reducing redundancy, is one viable tool for managing the explosive growth of data.
in general, clustering algorithms are designed to operate on a large set of related values, eventually generating a smaller set of elements that represent groups of similar data points.
a central data element may then be used as the sole representative of a group.
significant clustering work relating to bioinformatics may be traced to the late  <dig> s when methods for quick generation of nonredundant  protein databases were developed.
these combined identical or nearly identical protein sequences into single entries  <cit> .
the primary benefits of these methods include faster searches of the nr protein databases and reduced statistical bias in the query results  <cit> .
similarly, computer programs such as those in icatools  <cit>  were developed for compressing dna databases by removing redundant sequences found via clustering resulting in faster database queries.
note that the use of the term "clustering" in these applications differs from another use often found in the literature where clustering refers to generating a phylogenetic distance matrix, such as in  <cit> .
the operation of clustering used in this work identifies groups of sequences related by phylogeny; and it additionally applies to redundancy removal by identifying a sequence that suitably represents similar sequences.
recently, dna/rna clustering has attracted attention for a variety of reasons.
the drive to lower the expense of genome sequencing has led to the development of high-throughput sequencing technologies capable of generating millions of sequence fragments simultaneously.
a clustering preprocessing step can be used to remove a great amount of fragment redundancy which, in turn, allows for quicker fragment reassembly.
one of the more popular dna/rna clustering algorithms is cd-hit-est  <cit>  which was based on the protein clustering methods of  <cit>  and was developed for clustering dna/rna database data such as non-intron-containing expressed sequence tags .
a major application of cd-hit has been for clustering large data sets from microbiota analysis , often as a preprocessing step to create sets of highly related sequences representing operational taxonomic units .
these otus are subsequently used as a basis for estimating species diversity between treatment groups or quantitative relationships of taxa between treatment groups.
alternatively, representative sequences from the otus are used for phylogeny-based analyses.
a recent effort in  <cit>  to develop software tools which reduce the time required by blast  <cit>  to search large biological databases has resulted in a set of programs, including ublast and usearch, that reduce the search time by orders of magnitude.
as part of the work, an additional clustering program called uclust was created which utilizes the heuristic algorithm provided by usearch.
uclust generates results that dramatically improve upon the time required by cd-hit.
this work presents gramcluster, a fast and accurate algorithm for clustering large data sets of 16s rdna sequences based on the inherent grammar of dna and rna sequences.
lempel-ziv parsing  <cit>  is used to estimate the grammar of each sequence to provide a distance metric among sequences.
the implementation of this algorithm allows for fast and accurate clustering of biological information.
the following sections describe the algorithm and present results, including comparisons with the cd-hit-est algorithm and the recently developed uclust algorithm.
RESULTS
grammar
necessary concepts for understanding how a grammar model is specified are briefly reviewed in this section.
an alphabet, Σ, is a finite, nonempty set of symbols from which finite-length sequences, or strings, are formed.
strings are constructed via the binary operation of concatenation which begins with a copy of the left string and appends a copy of the right string.
a language, l, is then defined as a subset of strings selected from the set of all strings over an alphabet, and a problem is defined as the question of deciding whether a given string is a member of some particular language.
that is, given a string, w ∈ Σ*, and l, a language over Σ, decide if w ∈ l.
as l may be infinite, it is useful to have a compact description of the strings in l. such an abstract model is called a grammar, g. typically, a grammar is specified by the 4-tuple g = , where v is the set of variables and t is the set of terminals which are symbols that form the strings of l. p is the set of productions, each of which represent the recursive definition of l; and s ∈ v is the start symbol, which is the variable that defines l. each production consists of a head variable followed by the production operator → and a body string of zero or more terminals and variables.
each production represents one way to form strings in l from the head variable.
given g = , the language, l, is defined by
l={ w|w∈t*, and s ⇒*w}.
that is, l is the set of all strings derived from s.
it was observed in  <cit> , that a grammar, g, used to model a string can be converted to an lz <dig> representation in a simple way.
the term lz <dig> refers to lempel-ziv dictionary-based lossless compression detailed in  <cit>  and  <cit> .
subsequently, an algorithm was presented in  <cit>  to use an inverted process to map an lz77-compressed sequence into a grammar.
while the inverted process is more involved, it demonstrates the fact that lempel-ziv compression can be thought of as inferring a grammar from the sequence it compresses.
the original concept behind abstract grammars is that a grammar, g, is meant to completely describe the underlying structure of a corpus of sequences.
because most naturally occurring sequences contain repetition and redundancy, grammars are often able to describe sequences efficiently.
algorithm
a general overview of the gramcluster algorithm is shown in figure  <dig>  the set of sequences, s, is regarded as input to the algorithm with s = {s <dig> ...,sn}, where si is the ith sequence and i ∈ { <dig> ..., n}.
the goal of the algorithm is to partition s where each sequence is grouped with similar sequences from s such that all sequences within each resulting cluster are more similar to each other than sequences from other clusters.
the final partition is represented by the set of clusters, c = {c <dig> ..., cm}, where cj is the jth cluster and j ∈ { <dig> ..., m}.
the algorithm initially generates a suffix tree, ti, and grammar dictionary, di, associated with each sequence, si.
for each sequence, si, these data structures are used to determine if an existing cluster contains sufficiently similar sequences to si or if a new cluster needs to be created.
if a cluster, cj ∈ c, already exists with similar sequences, the sequence si is added to cj.
however, if no cluster contains similar sequences, a new cluster containing only si is added to c. this clustering continues for all sequences in s. the algorithm is described in more detail below with reference to the various blocks in figure  <dig>
dictionary creation
one of the core processes of the clustering algorithm is the formation of a distance estimate between an unprocessed sequence, si, and each cluster, cj, already in the partition, c. to this end, one sequence, called the representative sequence, is used to represent all other sequences within each cluster.
the distance between si and srj∈cj, where srj represents cj, is used to determine if si should be added to cj.
each sequence, si, is compared with, at most, the set of representative sequences, {srj|srj represents cj∈c}, to discover the correct cluster for si.
the distance metric relies on the structural rules necessarily present in all information-containing sequences.
gramcluster uses the grammar estimation method based on lempel-ziv  parsing  <cit>  as used in  <cit>  for language-phylogeny inference, in  <cit>  for phylogeny reconstruction, and in  <cit>  to construct a guide tree for multiple sequence alignment.
a similar grammar-based distance is also the focus of  <cit>  which analyzes the quality of the distance metric as a function of the length of the sequences.
the primary aspects of lz dictionary creation are shown in figure  <dig> where a set of grammar rules for each sequence is calculated.
initially, the dictionary, di1=∅, is empty, a fragment, f <dig> = si, is set to the first residue of the corresponding sequence, and only the first element, si, is visible to the algorithm.
at the kth iteration of the procedure, the kth residue is appended to the fragment resulting from the th step; and the visible sequence is checked.
if fk ∉ si, then fk is considered a new rule and so added to the dictionary, dik=dik−1∫{fk}; and the fragment is reset, fk=∅.
however, if fk ∈ si, then the current dictionary contains enough rules to reproduce the current fragment, i.e., dik=dik− <dig>  in either case, the iteration completes by appending the kth residue to the visible sequence.
this procedure continues until the visible sequence is equal to the entire sequence, at which time the size of the dictionary, |di| is determined for use in the metric calculation.
the correlation of the lz-based distance with phylogenetic distance was exploited in  <cit>  to obtain phylogenies for a set of mammalian species using complete mitochondrial dna and for the superfamily cavioidea using exon# <dig> of the growth hormone receptor  gene, the transthyretin  gene, and the  <dig> s rrna gene.
in  <cit> , the same distance metric was used to obtain phylogenies for fungal species using the cytochrome b gene and internal transcribed spacer regions of the rdna gene complex.
suffix tree construction
as shown in figure  <dig> the algorithm also constructs a suffix tree for the sequence.
suffix trees are data structures designed to contain all l suffix substrings of a length-l sequence  <cit> .
for example, a suffix tree for the sequence "gagacat" is schematically shown in figure  <dig>  all seven suffixes {gagacat, agacat, gacat, acat, cat, at, t} are found by tracing a unique path from the root node to one of the seven leaf nodes along solid lines.
one valuable use of suffix trees is searching for substrings which can be thought of as the preffix of a suffix.
by using a suffix tree, a length-l sequence can be completely scanned for a length-f fragment in o time as opposed to o  for a brute force search.
also depicted in figure  <dig> are the dashed-line suffix links which are a fundamental feature for linear-time construction of the suffix tree  <cit> .
a sequence, si, can be converted into a suffix tree, ti, in linear time and then searched for substrings in linear time based on the fragment length.
as will be shown, suffix tree sequence representation is important for reducing the time required for gramcluster to complete all necessary grammar-based comparisons.
clustering
the final component of the algorithm depicted in figure  <dig> is represented by the block labeled, "add to cluster."
the procedure for adding a sequence to a cluster is shown in greater detail in figure  <dig>  the algorithm checks each cluster, cj ∈ c, until a cluster is found where the distance between the representative sequence, srj and si,dj=dist is less than a user-defined threshold, t. once this condition is met, the cluster is updated, cj = cj ∪ {si}; and processing in this block terminates.
if no clusters meet the condition of d < t, a new cluster is created with si as its first member.
the following sections describe the cluster data structure, the representative sequence selection method, and the grammar-based distance calculation.
cluster data structure
in order to follow the cluster classification process, it is helpful to understand the data structure used to represent each cluster.
in particular, every cluster uses a list of suffix trees, ti, and dictionary sizes, |di|, to identify its set of sequences.
the remaining components contained in the data structure are used to determine and specify the representative sequence, srj, of the cluster, cj.
a good selection for srj is a sequence that appears grammatically similar to all other sequences within the cluster.
this implies the need to estimate the grammar-based distance between all sequences of the cluster, a computationally expensive task.
to avoid this cost, gramcluster selects only a few specific sequences in the cluster, that we will call "basis sequences," to which all others are compared.
the representative sequence, srj, can be determined by considering the sets of relative distances between all sequences and each basis sequence.
the centroid of the cluster is then defined as the vector containing the mean values of each set of relative distances.
the sequence with relative distances nearest to the centroid is selected as srj.
to see why this method is effective, consider that clustering is often performed in vector spaces where each element being classified is specified by a vector.
the points spatially near each other are placed into the same cluster, and the representative is typically selected as the point that is closest to the center of the cluster.
this idea is adapted in gramcluster, with an example depicted in figure  <dig>  the example in the figure contains forty sequences plotted in a two-dimensional space.
each dimension represents the grammar-based distance between the plotted sequence point and a basis sequence.
the data set used in this example contained forty 16s rdna sequences each from four genera .
of the two initially selected basis sequences, one came from acetobacter and the second from flavobacterium.
then, the pair of distances between each sequence and the basis sequences was computed and plotted.
as can be seen from the plot, the sequences group into clusters which correspond to their genus.
note that the basis sequences are not orthogonal; however, use is made of the fact that the grammar-based distances tend to obey the transitive property such that if
db=distdc=dist
and if db is close to dc, then sb and sc tend to be grammatically similar to each other.
the example in figure  <dig> demonstrates this by the use of basis sequences from acetobacter  and flavobacterium .
one would expect that comparing all sequences to one sequence would provide separation between the sequences from the same genus as the basis sequence and the rest.
however, sequences from the other genera also form into clusters as a result of sequences being compared to a single basis sequence.
in our example, all forty sequences are compared to just two sequences; and four clear clusters appear.
the method presented here for building vectors of distances relative to basis sequences is similar to the concept of embedding presented in  <cit> .
the work of  <cit>  details an algorithm called mbed that operates on a set of sequences to generate a distance matrix representing a phylogenetic guide tree, a process that is closely related to the data clustering problem presented here.
the mbed algorithm selects a subset of t seed reference sequences that are not close together relative to a distance metric.
then each sequence has a t-dimensional vector associated with it where each coordinate value is the distance between the sequence and the respective reference sequence.
the distance used in  <cit>  was selected to be the k-tuple distance measure of  <cit>  and implemented in clustalw  <cit> .
the basis sequence concept used in this work is similar, with the grammar-based distance metric replacing the k-tuple distance measure being the primary difference.
additionally, a single reference subset is used in  <cit>  to build all vectors.
the algorithm presented here creates vectors for each sequence contained in a cluster relative to basis sequences also sampled from the same cluster.
representative sequence selection
as shown in figure  <dig> the clustering process begins by comparing sequence si to the representative sequence of cluster cj ∈ c. for clusters containing many sequences, a representative sequence is determined using the basis sequence method described above.
in this case, only the representative sequence, srj, is compared to si
d=dist.
however, the progressive addition of sequences to clusters means there are clusters containing only a few sequences.
these clusters do not contain a large enough sample set to yield a reliable representative.
thus, until a cluster is large enough, all sequences are considered representative and compared to si
dk=dist ∀sk∈cj.
the minimum distance, min⁡k{dk}, is used as the classification metric.
grammar-based distance calculation
the distance metric used in gramcluster is a modified form of the grammar-based distance metric introduced in  <cit>  and used in  <cit> .
the original distance metric is computed by concatenating the two sequences being compared into a single sequence and then performing the operations detailed in figure  <dig>  formally, consider the process of comparing sequences sm and sn.
initially, the dictionary, dm,n1=dm, is set to that of sequence sm, a fragment, f <dig> = sn, is set to the first residue of the nth sequence, and the visible sequence is all of sm.
the algorithm operates as described previously, resulting in a new dictionary size, |dm,n|.
when complete, more grammatically similar sequences will have a new dictionary size with fewer entries as compared to sequences that are less grammatically similar.
therefore, the size of the new dictionary, |dm,n|, will be close to the size of the original dictionary, |dm|.
the distance between the sequences is estimated using the dictionary sizes, in particular
d={|dm,n|−|dm|dm×|sm||sn|if|sm| > |sn|,|dn,m|−|dn||dn|×|sn||sm|if|sm| ≤ |sn|.
this particular metric accounts for differences in sequence lengths and normalizes accordingly.
smaller values of d indicate a stronger similarity.
intuitively, sequences with a similar grammar should be clustered with each other.
while this grammar-based distance metric works well, it requires that the extended sequence be rescanned for every residue in the second sequence.
this means that sm will be rescanned completely for every character in sn.
this process is repeated as many times as the number of sequences compared to sm.
as a result, approximately 75% of the computation is devoted to string searching and concatenation.
to improve the execution time, we introduce two significant modifications described below.
fragment markers
the original distance calculation would simply repeat the process depicted in figure  <dig> on the concatenation of two sequences being compared.
thus, for the kth character in the second sequence, the first sequence is completely scanned along with the initial k -  <dig> portion of the second sequence.
however, this is quite unnecessary since many fragments formed from the second sequence were already found in the second sequence during the initial scan.
formally, consider sequences sm and sn which have already had their own dictionaries created in a previous step.
now suppose the concatenated sequence sm·n is being processed for the kth character in sn, at which point there is a nonempty fragment, fk.
the process begins with the fragment completely composed of consecutive letters from sn, which means that this fragment has already been created once before when sn was processed by itself.
as long as fk was previously found within sn, there will be no new information gained by scanning sm·n, because it is certain to be there since sn ⊂ sm·n.
so, there is no need to scan for fragments that have been previously found during any distance calculation.
the inverse statement is also true: fragments not previously found do need to be scanned for during a distance calculation.
this is implemented as shown in figure  <dig> in which fragment fk ∉ = si, so k is added to a list of marked fragment indices.
the same distance metric given by  is used, but there is no longer a need to perform string concatenation; and only the first string is scanned for the marked fragments from the second string.
formally, consider the process of comparing sequences sm and sn.
initially, the dictionary, dm,n1=dm, is set to that of sequence sm, a fragment, f marked, is set to the first marked substring of the nth sequence, and the visible sequence is always just sm.
the algorithm simply scans sm for an occurrence of the fragment and adds one to the dictionary if the fragment is not found.
either way, the fragment is updated to the next marked substring of sn; and sm is scanned again.
this continues for all marked fragments from sn resulting in a new dictionary size, | dm,n|.
this fragment marking process significantly reduces the total number of substring searches performed, as well as the character concatenations that would be otherwise required.
the second optimization involves a time-efficient method of searching a string for a substring of characters, a very relevant problem for suffix trees.
suffix tree searches
as stated previously, a length-l sequence stored in a suffix tree data structure can be completely scanned for a length-f fragment in o time.
to see why this is true, consider the simple example depicted in figure  <dig>  every suffix is represented in the data structure as a unique path beginning at the root node and traversing along a solid line to a leaf node.
any substring occurring in this string has to be the start of a suffix, so searching for a substring amounts to finding a suffix that begins with the substring.
consider searching "gagacat" for the substring fragment "gac" which is present in the string.
the first step is to find a branch beginning with "g" leaving the root, which is found as the third entry in the data structure.
following the branch to the internal node indicates that all suffixes in this tree that begin with "g" are always followed by an "a," which is also true of the fragment.
at the internal node, the next step is to search for any branch that begins with "c," which is found as the second entry in the data structure, concluding the search.
next, consider searching for the substring fragment "gact," which follows the previous search to the internal node and includes identifying the branch beginning with "c." the final step is looking at the subsequent character along the branch, which is "a," and does not match.
this search finishes having determined that "gact" is not a substring of "gagacat."
the use of the suffix tree in this context means that the time necessary for identifying whether previously marked fragments from sequence sn are present in sequence sm is o.
algorithm complexity
the algorithm complexity of gramcluster may be broken into three pieces, beginning with the generation of each sequence grammar dictionary, di for i ∈ { <dig>  ..., n}, where n is the number of sequences.
suppose the average sequence length is l, then each di results in complexity o , so all dictionaries are generated with complexity o .
next, each suffix tree, ti, has a complexity o , so all sequences are converted into trees with complexity o. finally, suppose the average number of clusters is m. as an upper bound, all clusters are scanned until each sequence is classified and each scanning process has complexity o. the result is a total scanning complexity of o. thus, the entire time complexity for gramcluster is o, which simplifies to o.
regarding the memory complexity of gramcluster and continuing with n as the number of sequences, suppose the average sequence header length in the fasta file is h. because every header line is stored for subsequent file output, this memory complexity is o. as before, if the average sequence length is l, then each sequence is stored in o. the worst-case memory usage for the clusters themselves occurs if every cluster created has an incomplete set of basis sequences.
in this case, each cluster has a memory complexity of o where c is the number of sequences held within the cluster and b is the number of basis sequences per cluster.
because there are n sequences stored in memory during this worst-case scenario, a final upper bound on the memory complexity is on) in which the most significant component has a memory complexity of o.
testing
we performed several clustering experiments to validate the proposed algorithm, gramcluster version  <dig>  .
the training procedures for obtaining the default parameters are described in the methods section.
in particular, we used gramcluster to cluster sets of 16s rdna sequences.
as detailed in the methods section, the resulting clusters were analyzed for correctness whereby the genus of each sequence was compared to that of all other sequences in the data set.
correct classification is considered when sequences belonging to the same genus fall into the same cluster.
likewise, incorrect classification occurs when sequences belonging to different genera are placed into the same cluster.
each output set was analyzed using several statistical quality metrics described in the methods section.
for comparison, cd-hit-est   <cit>  and uclust version  <dig> .
<dig>  <cit>  were also used to cluster the same 16s rdna sequences and analyzed using the same quality metrics.
experiments with moderate-sized data set
the proposed algorithm was evaluated using the folkes and mallows index, the jaccard coefficient, and rand statistic measures  <cit> , along with in-cluster classification and sequence differentiation percentages, all defined in the methods section.
the results for gramcluster, cd-hit-est, and uclust are presented in figure  <dig>
results indicate that cd-hit-est achieved  <dig> % in-cluster classification and  <dig> % sequence differentiation out of the  <dig>  total clusters determined.
that is, for sequences that were supposed to be in the same cluster, cd-hit-est placed them together  <dig> % of the time; and for sequences that were not supposed to be in the same cluster, it correctly kept them in different clusters  <dig> % of the time.
improved results for uclust show  <dig> % and  <dig> % in-cluster classification and sequence differentiation out of the  <dig>  total clusters determined.
by comparison, gramcluster achieved  <dig> % in-cluster classification and  <dig> % sequence differentiation out of the  <dig>  total clusters identified.
clearly, gramcluster provides a significant improvement in clustering sequences correctly.
this improvement can be further observed using common statistical measures for evaluating the performance of clustering algorithms  <cit>  described in the methods section.
these measures are shown for gramcluster, cd-hit-est, and uclust operating on a set of  <dig>  16s rdna genes obtained from  <dig>  different genera.
the jaccard coefficient and folkes and mallows index exceed those of cd-hit-est four-fold and over two-fold, respectively.
the cpu execution time of gramcluster  is on the same order as that of cd-hit-est , which is considered ultra-fast  <cit> .
the uclust cpu execution time  is much faster than gramcluster, however its quality metrics fall significantly short of those provided by gramcluster.
experiments with large data set
in order to simulate the application of clustering a large set of unknown fragments that typically result from  <dig> pyrosequencing, the previous fasta file was modified such that every sequence was reduced to only the first  <dig> bases and then repeated  <dig> times for a total of  <dig> , <dig> sequences from  <dig>  genera.
figure  <dig> contains data covering the same categories as in the previous experiment.
cd-hit-est achieved only  <dig> % in-cluster classification and  <dig> % sequence differentiation of the  <dig>  clusters found.
so, for sequences that were supposed to be in the same cluster, cd-hit-est placed them together  <dig> % of the time; and for sequences that were not supposed to be in the same cluster, it correctly kept them in different clusters  <dig> % of the time.
as in the previous experiment, results for uclust show  <dig> % and  <dig> % in-cluster classification and sequence differentiation out of the  <dig>  total clusters determined.
by comparison, gramcluster achieved  <dig> % and  <dig> % out of the  <dig>  clusters identified.
gramcluster continues to show a significant improvement in terms of clustering sequences correctly with each other.
this improvement can be seen further with the higher statistical measures, especially in the jaccard coefficient and folkes and mallows index which are over six and two times those of cd-hit-est.
perhaps most interestingly, gramcluster identified a more accurate number of clusters at  <dig> , even though the length of the sequences was significantly reduced, while both cd-hit-est and uclust reported identifying over  <dig>  clusters.
we also tested blastclust  <cit>  on 16s sequences.
the program was too slow for classifying the original set of  <dig>  sequences so we tested it using only 10% of the sequences.
the results are shown in figure  <dig>  as can be seen, the results of cd-hit-est, uclust, and gramcluster all tend to match those of figure  <dig>  as can be seen in figure  <dig> blastclust resulted in lower statistical metric scores in all categories, a high number of clusters compared to the number of genera.
it is clear that the exclusion of blastclust from the other experiments due to its inability to operate on the size of the input data set has not diminished the results.
varying command line options
next, we consider the effect of varying the command line options primarily responsible for affecting the resulting data set partition.
we ran two additional clustering experiments on the original set of sequences with gramcluster and uclust.
the gramcluster experiments had both grammar-based distance thresholds altered from the default setting of  <dig>  to  <dig>  and  <dig> .
similarly, the uclust experiments had the identity threshold altered from the default setting of 90% to 85% and 95%.
experiments clustering on species
the final experiment operated on the original set of sequences, but the partitioning was based on the sequence species instead of their genus.
CONCLUSIONS
the primary goal of this work was to introduce a computationally efficient clustering algorithm which can be used for clustering large datasets with high accuracy.
the algorithm introduced was validated against a specific class of datasets containing 16s rdna sequences but was designed to cluster any set of rna, dna, or protein sequences.
the grammar-based distance work introduced in  <cit>  and previously used in  <cit>  was modified to generate an estimation of the proper classification in which sequences are to be grouped.
results from clusters generated were presented in an attempt to study the overall quality of the resultant classifications as well as the computation time necessary to achieve the outputs.
accurate clustering of large numbers of biological sequences in an efficient amount of time is an important and challenging problem with a wide spectrum of applications.
in this work, we adapted existing ideas in a novel way and introduced significant improvements.
the proposed algorithm achieved higher-quality clusters compared to existing methods while operating at similar, high-speed execution times.

 the amount of biological information being gathered is growing faster than the rate at which it can be analyzed. data clustering, which compresses the problem space by reducing redundancy, is one viable tool for managing the explosive growth of data. in general, clustering algorithms are designed to operate on a large set of related values, eventually generating a smaller set of elements that represent groups of similar data points. a central data element may then be used as the sole representative of a group. significant clustering work relating to bioinformatics may be traced to the late s when methods for quick generation of nonredundant protein databases were developed. these combined identical or nearly identical protein sequences into single entries. the primary benefits of these methods include faster searches of the nr protein databases and reduced statistical bias in the query results. similarly, computer programs such as those in icatools were developed for compressing dna databases by removing redundant sequences found via clustering resulting in faster database queries. note that the use of the term "clustering" in these applications differs from another use often found in the literature where clustering refers to generating a phylogenetic distance matrix, such as in. the operation of clustering used in this work identifies groups of sequences related by phylogeny; and it additionally applies to redundancy removal by identifying a sequence that suitably represents similar sequences. recently, dna/rna clustering has attracted attention for a variety of reasons. the drive to lower the expense of genome sequencing has led to the development of high-throughput sequencing technologies capable of generating millions of sequence fragments simultaneously. a clustering preprocessing step can be used to remove a great amount of fragment redundancy which, in turn, allows for quicker fragment reassembly. one of the more popular dna/rna clustering algorithms is cd-hit-est which was based on the protein clustering methods of and was developed for clustering dna/rna database data such as non-intron-containing expressed sequence tags. a major application of cd-hit has been for clustering large data sets from microbiota analysis, often as a preprocessing step to create sets of highly related sequences representing operational taxonomic units. these otus are subsequently used as a basis for estimating species diversity between treatment groups or quantitative relationships of taxa between treatment groups. alternatively, representative sequences from the otus are used for phylogeny-based analyses. a recent effort in to develop software tools which reduce the time required by blast to search large biological databases has resulted in a set of programs, including ublast and usearch, that reduce the search time by orders of magnitude. as part of the work, an additional clustering program called uclust was created which utilizes the heuristic algorithm provided by usearch. uclust generates results that dramatically improve upon the time required by cd-hit. this work presents gramcluster, a fast and accurate algorithm for clustering large data sets of 16s rdna sequences based on the inherent grammar of dna and rna sequences. lempel-ziv parsing is used to estimate the grammar of each sequence to provide a distance metric among sequences. the implementation of this algorithm allows for fast and accurate clustering of biological information. the following sections describe the algorithm and present results, including comparisons with the cd-hit-est algorithm and the recently developed uclust algorithm.  grammar necessary concepts for understanding how a grammar model is specified are briefly reviewed in this section. an alphabet, Σ, is a finite, nonempty set of symbols from which finite-length sequences, or strings, are formed. strings are constructed via the binary operation of concatenation which begins with a copy of the left string and appends a copy of the right string. a language, l, is then defined as a subset of strings selected from the set of all strings over an alphabet, and a problem is defined as the question of deciding whether a given string is a member of some particular language. that is, given a string, w ∈ Σ*, and l, a language over Σ, decide if w ∈ l. as l may be infinite, it is useful to have a compact description of the strings in l. such an abstract model is called a grammar, g. typically, a grammar is specified by the 4-tuple g =, where v is the set of variables and t is the set of terminals which are symbols that form the strings of l. p is the set of productions, each of which represent the recursive definition of l; and s ∈ v is the start symbol, which is the variable that defines l. each production consists of a head variable followed by the production operator → and a body string of zero or more terminals and variables. each production represents one way to form strings in l from the head variable. given g =, the language, l, is defined by l={ w|w∈t*, and s ⇒*w}. that is, l is the set of all strings derived from s. it was observed in, that a grammar, g, used to model a string can be converted to an lz representation in a simple way. the term lz refers to lempel-ziv dictionary-based lossless compression detailed in and. subsequently, an algorithm was presented in to use an inverted process to map an lz77-compressed sequence into a grammar. while the inverted process is more involved, it demonstrates the fact that lempel-ziv compression can be thought of as inferring a grammar from the sequence it compresses. the original concept behind abstract grammars is that a grammar, g, is meant to completely describe the underlying structure of a corpus of sequences. because most naturally occurring sequences contain repetition and redundancy, grammars are often able to describe sequences efficiently. algorithm a general overview of the gramcluster algorithm is shown in figure the set of sequences, s, is regarded as input to the algorithm with s = {s...,sn}, where si is the ith sequence and i ∈ {..., n}. the goal of the algorithm is to partition s where each sequence is grouped with similar sequences from s such that all sequences within each resulting cluster are more similar to each other than sequences from other clusters. the final partition is represented by the set of clusters, c = {c..., cm}, where cj is the jth cluster and j ∈ {..., m}. the algorithm initially generates a suffix tree, ti, and grammar dictionary, di, associated with each sequence, si. for each sequence, si, these data structures are used to determine if an existing cluster contains sufficiently similar sequences to si or if a new cluster needs to be created. if a cluster, cj ∈ c, already exists with similar sequences, the sequence si is added to cj. however, if no cluster contains similar sequences, a new cluster containing only si is added to c. this clustering continues for all sequences in s. the algorithm is described in more detail below with reference to the various blocks in figure dictionary creation one of the core processes of the clustering algorithm is the formation of a distance estimate between an unprocessed sequence, si, and each cluster, cj, already in the partition, c. to this end, one sequence, called the representative sequence, is used to represent all other sequences within each cluster. the distance between si and srj∈cj, where srj represents cj, is used to determine if si should be added to cj. each sequence, si, is compared with, at most, the set of representative sequences, {srj|srj represents cj∈c}, to discover the correct cluster for si. the distance metric relies on the structural rules necessarily present in all information-containing sequences. gramcluster uses the grammar estimation method based on lempel-ziv parsing as used in for language-phylogeny inference, in for phylogeny reconstruction, and in to construct a guide tree for multiple sequence alignment. a similar grammar-based distance is also the focus of which analyzes the quality of the distance metric as a function of the length of the sequences. the primary aspects of lz dictionary creation are shown in figure where a set of grammar rules for each sequence is calculated. initially, the dictionary, di1=∅, is empty, a fragment, f = si, is set to the first residue of the corresponding sequence, and only the first element, si, is visible to the algorithm. at the kth iteration of the procedure, the kth residue is appended to the fragment resulting from the th step; and the visible sequence is checked. if fk ∉ si, then fk is considered a new rule and so added to the dictionary, dik=dik−1∫{fk}; and the fragment is reset, fk=∅. however, if fk ∈ si, then the current dictionary contains enough rules to reproduce the current fragment, i.e., dik=dik− in either case, the iteration completes by appending the kth residue to the visible sequence. this procedure continues until the visible sequence is equal to the entire sequence, at which time the size of the dictionary, |di| is determined for use in the metric calculation. the correlation of the lz-based distance with phylogenetic distance was exploited in to obtain phylogenies for a set of mammalian species using complete mitochondrial dna and for the superfamily cavioidea using exon# of the growth hormone receptor gene, the transthyretin gene, and the s rrna gene. in, the same distance metric was used to obtain phylogenies for fungal species using the cytochrome b gene and internal transcribed spacer regions of the rdna gene complex. suffix tree construction as shown in figure the algorithm also constructs a suffix tree for the sequence. suffix trees are data structures designed to contain all l suffix substrings of a length-l sequence. for example, a suffix tree for the sequence "gagacat" is schematically shown in figure all seven suffixes {gagacat, agacat, gacat, acat, cat, at, t} are found by tracing a unique path from the root node to one of the seven leaf nodes along solid lines. one valuable use of suffix trees is searching for substrings which can be thought of as the preffix of a suffix. by using a suffix tree, a length-l sequence can be completely scanned for a length-f fragment in o time as opposed to o for a brute force search. also depicted in figure are the dashed-line suffix links which are a fundamental feature for linear-time construction of the suffix tree. a sequence, si, can be converted into a suffix tree, ti, in linear time and then searched for substrings in linear time based on the fragment length. as will be shown, suffix tree sequence representation is important for reducing the time required for gramcluster to complete all necessary grammar-based comparisons. clustering the final component of the algorithm depicted in figure is represented by the block labeled, "add to cluster." the procedure for adding a sequence to a cluster is shown in greater detail in figure the algorithm checks each cluster, cj ∈ c, until a cluster is found where the distance between the representative sequence, srj and si,dj=dist is less than a user-defined threshold, t. once this condition is met, the cluster is updated, cj = cj ∪ {si}; and processing in this block terminates. if no clusters meet the condition of d < t, a new cluster is created with si as its first member. the following sections describe the cluster data structure, the representative sequence selection method, and the grammar-based distance calculation. cluster data structure in order to follow the cluster classification process, it is helpful to understand the data structure used to represent each cluster. in particular, every cluster uses a list of suffix trees, ti, and dictionary sizes, |di|, to identify its set of sequences. the remaining components contained in the data structure are used to determine and specify the representative sequence, srj, of the cluster, cj. a good selection for srj is a sequence that appears grammatically similar to all other sequences within the cluster. this implies the need to estimate the grammar-based distance between all sequences of the cluster, a computationally expensive task. to avoid this cost, gramcluster selects only a few specific sequences in the cluster, that we will call "basis sequences," to which all others are compared. the representative sequence, srj, can be determined by considering the sets of relative distances between all sequences and each basis sequence. the centroid of the cluster is then defined as the vector containing the mean values of each set of relative distances. the sequence with relative distances nearest to the centroid is selected as srj. to see why this method is effective, consider that clustering is often performed in vector spaces where each element being classified is specified by a vector. the points spatially near each other are placed into the same cluster, and the representative is typically selected as the point that is closest to the center of the cluster. this idea is adapted in gramcluster, with an example depicted in figure the example in the figure contains forty sequences plotted in a two-dimensional space. each dimension represents the grammar-based distance between the plotted sequence point and a basis sequence. the data set used in this example contained forty 16s rdna sequences each from four genera. of the two initially selected basis sequences, one came from acetobacter and the second from flavobacterium. then, the pair of distances between each sequence and the basis sequences was computed and plotted. as can be seen from the plot, the sequences group into clusters which correspond to their genus. note that the basis sequences are not orthogonal; however, use is made of the fact that the grammar-based distances tend to obey the transitive property such that if db=distdc=dist and if db is close to dc, then sb and sc tend to be grammatically similar to each other. the example in figure demonstrates this by the use of basis sequences from acetobacter and flavobacterium. one would expect that comparing all sequences to one sequence would provide separation between the sequences from the same genus as the basis sequence and the rest. however, sequences from the other genera also form into clusters as a result of sequences being compared to a single basis sequence. in our example, all forty sequences are compared to just two sequences; and four clear clusters appear. the method presented here for building vectors of distances relative to basis sequences is similar to the concept of embedding presented in. the work of details an algorithm called mbed that operates on a set of sequences to generate a distance matrix representing a phylogenetic guide tree, a process that is closely related to the data clustering problem presented here. the mbed algorithm selects a subset of t seed reference sequences that are not close together relative to a distance metric. then each sequence has a t-dimensional vector associated with it where each coordinate value is the distance between the sequence and the respective reference sequence. the distance used in was selected to be the k-tuple distance measure of and implemented in clustalw. the basis sequence concept used in this work is similar, with the grammar-based distance metric replacing the k-tuple distance measure being the primary difference. additionally, a single reference subset is used in to build all vectors. the algorithm presented here creates vectors for each sequence contained in a cluster relative to basis sequences also sampled from the same cluster. representative sequence selection as shown in figure the clustering process begins by comparing sequence si to the representative sequence of cluster cj ∈ c. for clusters containing many sequences, a representative sequence is determined using the basis sequence method described above. in this case, only the representative sequence, srj, is compared to si d=dist. however, the progressive addition of sequences to clusters means there are clusters containing only a few sequences. these clusters do not contain a large enough sample set to yield a reliable representative. thus, until a cluster is large enough, all sequences are considered representative and compared to si dk=dist ∀sk∈cj. the minimum distance, min⁡k{dk}, is used as the classification metric. grammar-based distance calculation the distance metric used in gramcluster is a modified form of the grammar-based distance metric introduced in and used in. the original distance metric is computed by concatenating the two sequences being compared into a single sequence and then performing the operations detailed in figure formally, consider the process of comparing sequences sm and sn. initially, the dictionary, dm,n1=dm, is set to that of sequence sm, a fragment, f = sn, is set to the first residue of the nth sequence, and the visible sequence is all of sm. the algorithm operates as described previously, resulting in a new dictionary size, |dm,n|. when complete, more grammatically similar sequences will have a new dictionary size with fewer entries as compared to sequences that are less grammatically similar. therefore, the size of the new dictionary, |dm,n|, will be close to the size of the original dictionary, |dm|. the distance between the sequences is estimated using the dictionary sizes, in particular d={|dm,n|−|dm|dm×|sm||sn|if|sm| > |sn|,|dn,m|−|dn||dn|×|sn||sm|if|sm| ≤ |sn|. this particular metric accounts for differences in sequence lengths and normalizes accordingly. smaller values of d indicate a stronger similarity. intuitively, sequences with a similar grammar should be clustered with each other. while this grammar-based distance metric works well, it requires that the extended sequence be rescanned for every residue in the second sequence. this means that sm will be rescanned completely for every character in sn. this process is repeated as many times as the number of sequences compared to sm. as a result, approximately 75% of the computation is devoted to string searching and concatenation. to improve the execution time, we introduce two significant modifications described below. fragment markers the original distance calculation would simply repeat the process depicted in figure on the concatenation of two sequences being compared. thus, for the kth character in the second sequence, the first sequence is completely scanned along with the initial k - portion of the second sequence. however, this is quite unnecessary since many fragments formed from the second sequence were already found in the second sequence during the initial scan. formally, consider sequences sm and sn which have already had their own dictionaries created in a previous step. now suppose the concatenated sequence sm·n is being processed for the kth character in sn, at which point there is a nonempty fragment, fk. the process begins with the fragment completely composed of consecutive letters from sn, which means that this fragment has already been created once before when sn was processed by itself. as long as fk was previously found within sn, there will be no new information gained by scanning sm·n, because it is certain to be there since sn ⊂ sm·n. so, there is no need to scan for fragments that have been previously found during any distance calculation. the inverse statement is also true: fragments not previously found do need to be scanned for during a distance calculation. this is implemented as shown in figure in which fragment fk ∉ = si, so k is added to a list of marked fragment indices. the same distance metric given by is used, but there is no longer a need to perform string concatenation; and only the first string is scanned for the marked fragments from the second string. formally, consider the process of comparing sequences sm and sn. initially, the dictionary, dm,n1=dm, is set to that of sequence sm, a fragment, f marked, is set to the first marked substring of the nth sequence, and the visible sequence is always just sm. the algorithm simply scans sm for an occurrence of the fragment and adds one to the dictionary if the fragment is not found. either way, the fragment is updated to the next marked substring of sn; and sm is scanned again. this continues for all marked fragments from sn resulting in a new dictionary size, | dm,n|. this fragment marking process significantly reduces the total number of substring searches performed, as well as the character concatenations that would be otherwise required. the second optimization involves a time-efficient method of searching a string for a substring of characters, a very relevant problem for suffix trees. suffix tree searches as stated previously, a length-l sequence stored in a suffix tree data structure can be completely scanned for a length-f fragment in o time. to see why this is true, consider the simple example depicted in figure every suffix is represented in the data structure as a unique path beginning at the root node and traversing along a solid line to a leaf node. any substring occurring in this string has to be the start of a suffix, so searching for a substring amounts to finding a suffix that begins with the substring. consider searching "gagacat" for the substring fragment "gac" which is present in the string. the first step is to find a branch beginning with "g" leaving the root, which is found as the third entry in the data structure. following the branch to the internal node indicates that all suffixes in this tree that begin with "g" are always followed by an "a," which is also true of the fragment. at the internal node, the next step is to search for any branch that begins with "c," which is found as the second entry in the data structure, concluding the search. next, consider searching for the substring fragment "gact," which follows the previous search to the internal node and includes identifying the branch beginning with "c." the final step is looking at the subsequent character along the branch, which is "a," and does not match. this search finishes having determined that "gact" is not a substring of "gagacat." the use of the suffix tree in this context means that the time necessary for identifying whether previously marked fragments from sequence sn are present in sequence sm is o. algorithm complexity the algorithm complexity of gramcluster may be broken into three pieces, beginning with the generation of each sequence grammar dictionary, di for i ∈ {..., n}, where n is the number of sequences. suppose the average sequence length is l, then each di results in complexity o, so all dictionaries are generated with complexity o. next, each suffix tree, ti, has a complexity o, so all sequences are converted into trees with complexity o. finally, suppose the average number of clusters is m. as an upper bound, all clusters are scanned until each sequence is classified and each scanning process has complexity o. the result is a total scanning complexity of o. thus, the entire time complexity for gramcluster is o, which simplifies to o. regarding the memory complexity of gramcluster and continuing with n as the number of sequences, suppose the average sequence header length in the fasta file is h. because every header line is stored for subsequent file output, this memory complexity is o. as before, if the average sequence length is l, then each sequence is stored in o. the worst-case memory usage for the clusters themselves occurs if every cluster created has an incomplete set of basis sequences. in this case, each cluster has a memory complexity of o where c is the number of sequences held within the cluster and b is the number of basis sequences per cluster. because there are n sequences stored in memory during this worst-case scenario, a final upper bound on the memory complexity is on) in which the most significant component has a memory complexity of o. testing we performed several clustering experiments to validate the proposed algorithm, gramcluster version. the training procedures for obtaining the default parameters are described in the methods section. in particular, we used gramcluster to cluster sets of 16s rdna sequences. as detailed in the methods section, the resulting clusters were analyzed for correctness whereby the genus of each sequence was compared to that of all other sequences in the data set. correct classification is considered when sequences belonging to the same genus fall into the same cluster. likewise, incorrect classification occurs when sequences belonging to different genera are placed into the same cluster. each output set was analyzed using several statistical quality metrics described in the methods section. for comparison, cd-hit-est and uclust version. were also used to cluster the same 16s rdna sequences and analyzed using the same quality metrics. experiments with moderate-sized data set the proposed algorithm was evaluated using the folkes and mallows index, the jaccard coefficient, and rand statistic measures, along with in-cluster classification and sequence differentiation percentages, all defined in the methods section. the results for gramcluster, cd-hit-est, and uclust are presented in figure results indicate that cd-hit-est achieved % in-cluster classification and % sequence differentiation out of the total clusters determined. that is, for sequences that were supposed to be in the same cluster, cd-hit-est placed them together % of the time; and for sequences that were not supposed to be in the same cluster, it correctly kept them in different clusters % of the time. improved results for uclust show % and % in-cluster classification and sequence differentiation out of the total clusters determined. by comparison, gramcluster achieved % in-cluster classification and % sequence differentiation out of the total clusters identified. clearly, gramcluster provides a significant improvement in clustering sequences correctly. this improvement can be further observed using common statistical measures for evaluating the performance of clustering algorithms described in the methods section. these measures are shown for gramcluster, cd-hit-est, and uclust operating on a set of 16s rdna genes obtained from different genera. the jaccard coefficient and folkes and mallows index exceed those of cd-hit-est four-fold and over two-fold, respectively. the cpu execution time of gramcluster is on the same order as that of cd-hit-est, which is considered ultra-fast. the uclust cpu execution time is much faster than gramcluster, however its quality metrics fall significantly short of those provided by gramcluster. experiments with large data set in order to simulate the application of clustering a large set of unknown fragments that typically result from pyrosequencing, the previous fasta file was modified such that every sequence was reduced to only the first bases and then repeated times for a total of, sequences from genera. figure contains data covering the same categories as in the previous experiment. cd-hit-est achieved only % in-cluster classification and % sequence differentiation of the clusters found. so, for sequences that were supposed to be in the same cluster, cd-hit-est placed them together % of the time; and for sequences that were not supposed to be in the same cluster, it correctly kept them in different clusters % of the time. as in the previous experiment, results for uclust show % and % in-cluster classification and sequence differentiation out of the total clusters determined. by comparison, gramcluster achieved % and % out of the clusters identified. gramcluster continues to show a significant improvement in terms of clustering sequences correctly with each other. this improvement can be seen further with the higher statistical measures, especially in the jaccard coefficient and folkes and mallows index which are over six and two times those of cd-hit-est. perhaps most interestingly, gramcluster identified a more accurate number of clusters at, even though the length of the sequences was significantly reduced, while both cd-hit-est and uclust reported identifying over clusters. we also tested blastclust on 16s sequences. the program was too slow for classifying the original set of sequences so we tested it using only 10% of the sequences. the results are shown in figure as can be seen, the results of cd-hit-est, uclust, and gramcluster all tend to match those of figure as can be seen in figure blastclust resulted in lower statistical metric scores in all categories, a high number of clusters compared to the number of genera. it is clear that the exclusion of blastclust from the other experiments due to its inability to operate on the size of the input data set has not diminished the results. varying command line options next, we consider the effect of varying the command line options primarily responsible for affecting the resulting data set partition. we ran two additional clustering experiments on the original set of sequences with gramcluster and uclust. the gramcluster experiments had both grammar-based distance thresholds altered from the default setting of to and. similarly, the uclust experiments had the identity threshold altered from the default setting of 90% to 85% and 95%. experiments clustering on species the final experiment operated on the original set of sequences, but the partitioning was based on the sequence species instead of their genus.  the primary goal of this work was to introduce a computationally efficient clustering algorithm which can be used for clustering large datasets with high accuracy. the algorithm introduced was validated against a specific class of datasets containing 16s rdna sequences but was designed to cluster any set of rna, dna, or protein sequences. the grammar-based distance work introduced in and previously used in was modified to generate an estimation of the proper classification in which sequences are to be grouped. results from clusters generated were presented in an attempt to study the overall quality of the resultant classifications as well as the computation time necessary to achieve the outputs. accurate clustering of large numbers of biological sequences in an efficient amount of time is an important and challenging problem with a wide spectrum of applications. in this work, we adapted existing ideas in a novel way and introduced significant improvements. the proposed algorithm achieved higher-quality clusters compared to existing methods while operating at similar, high-speed execution times. 
the performance of the proposed algorithm is validated via comparison to the popular dna/rna sequence clustering approach, cd-hit-est, and to the recently developed algorithm, uclust, using two different sets of 16s rdna sequences from genera. its statistical clustering quality is validated by clustering large datasets containing 16s rdna sequences. the validation results are especially striking for large datasets. we propose a sequence clustering algorithm and compare the partition quality and execution time of the proposed algorithm with those of a popular existing algorithm. the proposed algorithm maintains a comparable cpu execution time with that of cd-hit-est which is much slower than uclust, and has successfully generated clusters with higher statistical accuracy than both cd-hit-est and uclust. the proposed clustering algorithm uses a grammar-based distance metric to determine partitioning for a set of biological sequences. if comparison fails to identify a suitable cluster, a new cluster is created.
 we propose a sequence clustering algorithm and compare the partition quality and execution time of the proposed algorithm with those of a popular existing algorithm. the proposed clustering algorithm uses a grammar-based distance metric to determine partitioning for a set of biological sequences. the algorithm performs clustering in which new sequences are compared with cluster-representative sequences to determine membership. if comparison fails to identify a suitable cluster, a new cluster is created.  the performance of the proposed algorithm is validated via comparison to the popular dna/rna sequence clustering approach, cd-hit-est, and to the recently developed algorithm, uclust, using two different sets of 16s rdna sequences from genera. the proposed algorithm maintains a comparable cpu execution time with that of cd-hit-est which is much slower than uclust, and has successfully generated clusters with higher statistical accuracy than both cd-hit-est and uclust. the validation results are especially striking for large datasets.  we introduce a fast and accurate clustering algorithm that relies on a grammar-based sequence distance. its statistical clustering quality is validated by clustering large datasets containing 16s rdna sequences. 