Abstract
Background
Clustering DNA sequences into functional groups is an important problem in bioinformatics. We propose a new alignmentfree algorithm, mBKM, based on a new distance measure, DMk, for clustering gene sequences. This method transforms DNA sequences into the feature vectors which contain the occurrence, location and order relation of ktuples in DNA sequence. Afterwards, a hierarchical procedure is applied to clustering DNA sequences based on the feature vectors.
Results
The proposed distance measure and clustering method are evaluated by clustering functionally related genes and by phylogenetic analysis. This method is also compared with BlastClust, CDHITEST and some others. The experimental results show our method is effective in classifying DNA sequences with similar biological characteristics and in discovering the underlying relationship among the sequences.
Conclusions
We introduced a novel clustering algorithm which is based on a new sequence similarity measure. It is effective in classifying DNA sequences with similar biological characteristics and in discovering the relationship among the sequences.
Background
With the development of advanced biotechnology, more and more biological sequence information has been generated. The amount of genetic data is growing faster than the rate at which it can be analyzed. Clustering techniques provide a viable solution for handling and analyzing such rapidly growing genetic data. Clustering algorithms partition sequences into different biologically meaningful groups, facilitating therefore the prediction of functions of genes [1]. When a new gene is assigned to a cluster, the biological function of this cluster can be attributed to this gene with high confidence. On the other hand, clustering gene sequences into groups may also help with analyzing evolutionary relationships among the sequences in a cluster [2].
Clustering of gene sequences requires calculation of similarity between sequences. There are two clustering approaches according to the similarity measure used in a clustering method. One is based on sequence alignment. The similarity between two gene sequences is measured by the scores obtained from an alignment algorithm such as BLAST [3] or FASTA [4]. Although sequence alignment gives good solutions, it is relatively difficult to cluster a large number of sequences because of its computational complexity. Moreover, if the sequences in the set vary in length, a satisfactory alignment is hard to achieve, resulting in a low accuracy of clustering.
The other approach for similarity measure is to use alignmentfree methods [510]. In recent years, several alignmentfree measures have been proposed. The wordbased measure is one of the most widely used methods [1114]. This method chooses a short word length k, maps each sequence onto an ndimensional vector according to its klength tuple (also called ktuple or kword) properties, and then assesses the similarity of any two vectors by measures such as Euclidean distance [15], Mahalanobis distance [16], Kullback–Leibler discrepancy [17], cosine distance [18] or Pearson’s correlation coefficient [19]. In recent years, several novel alignmentfree measures [20,21] have been designed for DNA sequences analysis. Yang et al. [22] extended the ktuple distance, which is based on the difference in tuple frequencies, to clustering gene sequences. Their tuplebased method determines the similarity of sequences by considering only tuple frequencies and ignoring the positional information within a sequence.
Major algorithms used in gene sequence clustering can be divided into two categories according to the result format: hierarchical clustering algorithms and partitional clustering algorithms [23]. Hierarchical clustering is widely used for detecting clusters in genomic data. It generates a set of partitions forming a cluster hierarchy. According to linkage criteria, there are three hierarchical clustering methods including singlelinkage clustering (SL), completelinkage clustering (CL) and averagelinkage clustering (AL) [24]. With SL, clusters may be merged together due to single sequences being close to each other, even though many of the sequences in each cluster may be very distant to each other [25]. CL tends to find compact clusters of approximately equal diameters [25]. With CL, all objects in a cluster are similar to each other. AL can be seen as an intermediate between single and complete linkage clustering, resulting in more homogeneous clusters than those obtained by the singlelinkage method [26]. For instance, BlastClust [27] and GeneRage [28] employ single linkage clustering approach; SWORDS [29] is based on word frequencies as profiles to merge clusters hierarchically; and Uchiyama [30] use average linkage clustering algorithm to classify genes. Hierarchical approaches may yield fairly good results, but they require the similarity of all pairs of sequences and quickly arrive at a bottleneck in terms of computational time and memory usage for largescale data sets [31].
Partitioning algorithms have also been used. Partitional clustering obtains a partition of data objects by optimizing some clustering criterion. Partitional clustering algorithms are simple and wellsuited for clustering large datasets [32]. Kmeans (KM) [33,34] is a commonly used method of partitional clustering methods. KM has a lower order of computational complexity and demands less physical memory than the hierarchical method. It is suitable for clustering large gene data. Some KMbased algorithms, such as those introduced by Wan et al. [33], Kelarev et al. [34], Tseng et al. [35] and Ashlock et al. [36], have been developed to group DNA sequences. The major drawback of KM compared to hierarchical clustering algorithms is the lack of hierarchical relationships in its results. To remedy the problem, bisecting Kmeans (BKM), a hierarchical variation of KM, was proposed to build a tree of clusters in a topdown fashion by splitting the least homogeneous cluster into two more homogeneous ones. BKM can produce either a flat clustering or a hierarchical clustering by recursively applying KM. It has a linear complexity and is relatively efficient and scalable. Recent study [37] concluded that BKM outperforms KM and performs equally well or better than hierarchical methods when it partitions the dataset based on a homogeneity criteria. The bisecting approach is very attractive for genomic studies [38].
Hierarchical clustering produces a nested series of partitions, where the results are usually depicted as a dendrogram while partitional clustering produces a flat partition. BlastClust [27] is a hierarchical clustering method based on BLAST scores as the measure of sequence similarity. BlastClust computes pairwise similarity of all sequences by BLAST alignment and then clusters sequences by the single linkage clustering method which produces clusters of linear topology. The performance of BlastClust is limited by the size of the input data. CDHITEST [39], a partitional approach, is also widely used to cluster DNA sequences. CDHITEST uses an incremental clustering process and avoids the unnecessary alignments by a short word filtering mechanism, which detects similar sequences by counting the number of identical short words between them. The purpose of filters is to decide whether the identity between two sequences is above or below a threshold without aligning them, therefore speeding up the clustering process. Though CDHITEST is based on alignment, it can avoid too many pairwise alignments by using a filter, thus it is faster than BlastClust, and can handle larger datasets.
Recent studies reveal also that BlastClust is less effective for clustering divergent sequences [40], and its performance strongly depends on the choice of optimal BLAST parameters including similarity threshold, percent identity, and alignment length [41]. CDHITEST, on the other hand, does not provide hierarchical relationships between clusters of sequences. In many situations both CDHITEST and BlastClust yield clusters with only one sequence [41]. All the traditional clustering methods based on sequence alignment encounter computational difficulties in dealing with large biological databases.
The approach presented in this paper involves a new alignmentfree distance measure based on ktuples, DMk (Distance Measure based on ktuples) [42], and a modified bisecting Kmeans clustering algorithm, mBKM (modified Bisecting KMeans algorithm). mBKM aims to speed up the clustering process by using the alignmentfree similarity measure, and is able to produce either a hierarchical clustering or a partition clustering result. We have applied mBKM with DMk in clustering gene sequences and performing phylogenetic analysis. DMk shows better performance than the ktuple distance in our experiments, and mBKM outperforms SL, CL, AL, BKM and KM when tested on public gene sequence datasets. Furthermore, the proposed method also outperforms alignmentbased methods such as BlastClust and CDHITEST.
Methods
A gene is a stretch of DNA that codes for a single polypeptide chain [43]. A gene sequence is a succession of four symbols {A, C, G, T}. Because the similarity between the genes of two species indicates their evolutionary relationship, it is used in many clustering algorithms. The goal of sequence clustering is to partition biological sequences into meaningful/functional groups according to the similarity information, which is calculated using either an alignmentbased method or an alignmentfree method.
The traditional approach for clustering DNA sequences requires allbyall comparisons from alignment [4446]. Given two sequences: S_{1} = AGCACACA and S_{2} = ACACAGTA, S_{1}^{P} and S_{2}^{P} are used to represent the p^{th} characters in S_{1} and S_{2}, respectively. The alignment score [45] for (S_{1}S_{2}) is given by
where E is the cost of an alignment operation: deletion, substitution, or insertion. However this distance measure relies on sequence alignment. Since sequence alignment suffers in computational aspect with regard to large biological databases, clustering methods relying on sequence alignment have difficulties in dealing with the large gene data. An alignmentfree similarity measure helps avoid the computational complexity of multiple sequence alignment for similarity computation. In this paper we propose a new alignmentfree similarity measure, DMk, based on which we developed mBKM to cluster gene sequences.
In the follows, we will present DMk first, and then describe mBKM algorithms.
A new similarity measure: DMk
In this section, we introduce a new similarity measure which takes into account the occurrence, location and order relation of ktuple in a DNA sequence.
Sequences are numerically transformed to feature vectors that can be processed by data mining algorithms. Let Σ be the alphabet set of nucleotides (Σ = {A, C, G, T}). A sequence of length s, S, is defined as a linear succession of s symbols from Σ. A segment of k consecutive symbols in sequence S (k ≤ s) is designated as a ktuple. There is a set of 4^{k} possible ktuples, W_{k}. The number of occurrences of a ktuple w, N_{w}, is counted by moving a sliding window of length k over the sequence with k  1 bp overlapping step size.
To explore the correlation properties of DNA, Nair et al. [47] provided a presentation of genomic data using the internucleotide distance sequence. Based on a similar idea, we utilize the gaps between the locations where ktuple occur in the sequence to explore the sequence structure. For a DNA sequence Sp_{r} is the location of the r^{th} occurrence of ktuple w, where p_{0} = 0. And α_{r} is given as,
in which m stands for the number of occurrences of w. α_{r} reflects the density of w and is closely related to the location where w occurs in the sequence. Each w begins at the 1/α_{1} position, and {α_{1,}α_{2},…,α_{m}} for repetition of w forms an array whose r^{th} element indicates the relative position of two neighboring w in the sequence. This array allows us to find all subsequent repeats of w.
To characterize the order of α_{r}, we define β_{j} as a partial sum of {α_{r}}. β_{j} is calculated by the following formula:
{α_{r}} is a list of nonnegative real numbers, and β_{j} is totally ordered by ≤, so β_{1}, β_{2}, …, β_{m} is also an ordered set. {α_{1}, α_{2},…, α_{m}} and {β_{1}, β_{2},…, β_{m}} determine each other uniquely. β_{j} is only dependent of the number and positions of w and independent on other ktuples. Given the set of {β_{1}, β_{2},…, β_{m}}, one can obtain where w occurs and how many times w occurs in the sequence.
Shannon’s entropy[48], which illuminates the total information measure of source on the average, is a measure of order/disorder. According to [49], when using the totally ordered set {β_{1}β_{2},…, β_{m}} to calculate the probabilities, the Shannon entropy reflects the degree of importance of position in a sequence. We construct a discrete probability distribution , where , and . The Shannon entropy of the discrete probability distribution is calculated by
For each ktuple w in the sequence, not only the information of tuple numbers but also the information of tuple positions is involved in the definition of H. We take H as the feature of w in the sequence, and then construct a vector consisted of H of all possible ktuples in the given sequence.
For a fixed k, there are 4^{k} distinct ktuples to be considered. These ktuples in a fixed 4^{k}dimension feature vector are denoted by , where H_{i} means the feature representation of the ith ktuple. This feature vector based on H can be regarded as an index for its corresponding sequence.
Cluster analysis algorithms partition objects into groups based on the distances between objects. Euclidean distance is the square root of the summation of the squares of the differences between all pairs of corresponding objects. The ktuple distance is the sum of the differences in frequency over all possible ktuples; on the other hand, we use Euclidean distance between Shannon entropy of ktuples in sequences to measure the similarity. This distance measure method is referred as DMk. For any two sequences X and Y, DMk can be calculated as:
where and represent the Shannon entropy values of the i^{th}ktuple in sequences X and Y, respectively. DMk can be calculated from following algorithm:
Algorithm Name: DMk for similarity measure
Input: sequences {S_{1}, S_{2},…, S_{N}}.
Output: similarity matrix, (d(X,Y))_{N*N}.
Steps:
1. For each sequence, search and locate each ktuple;
1.1 For each ktuple, use Equation (1) to calculate
1.2 For each ktuple, use Equation (2) to calculate ;
1.3 For each ktuple, use Equation (3) to calculate H;
2. For each sequence, construct 4^{k} component vector by H of all ktuples.
3. For any two sequences, use Equation (4) to calculate the distance between the two sequences.
4. Return {d}.
A new clustering algorithm: mBKM
KM can be used to obtain a hierarchical clustering solution using a repeated bisecting approach [50,51]. BKM is such an algorithm and it can produce either a partitional or a hierarchical clustering.
BKM has a linear time complexity in each bisecting step. Recent study [51] concludes BKM outperforms KM as well as the agglomerative approach in terms of accuracy and efficiency. Consequently, the bisecting approach is very attractive in many applications for clustering and genomic data analysis.
BKM initially regards the whole data set as a cluster, and splits one cluster into two subclusters at each bisecting step using KM until singleton clusters are obtained at the leafs or until K clusters are obtained. The outcome is structured as a binary tree. There are two key steps in a typical BKM. The first one is the selection of initial centroids. Generally the initial centroids are chosen randomly in BKM. The second key step is the rule, ζ, for selection of a existing cluster to be split in each bisecting step. ζ is typically given by the following three approaches [50]:
1) Choosing the cluster with largest size;
2) Selecting the cluster with the overall similarity
The overall similarity is either minimized or maximize, depending on the definition of d(s, s^{’}). C is a cluster;
3) Using a criterion based on both size and overall similarity.
Because the differences between these methods are small in terms of the final clustering result, the way of splitting the largest remaining cluster is recommended [50].
There are two problems in BKM algorithm:
1. Randomly choosing the initial centroids in BKM may result in too adjacent elements selected. If the initial centroids are too close, the algorithm will reach a local optimization. Moreover, different sets of initial cluster centroids can lead to different final clustering results.
2. The algorithm for choosing one existing cluster to split in each bisecting step usually selects the cluster with the largest size. Although this leads to reasonably good and balanced clustering solution, it cannot gracefully work for datasets where the natural clusters are of different sizes, as it will tend to partition larger clusters first. In real biological data, the number of elements in every cluster may not always be similar.
To address the above two problems and obtain more natural hierarchical solutions, we develop a modified bisecting Kmeans, mBKM, which choose the initial centroids by the maximum and minimum principle and select the cluster to split based on the compactness of clusters.
1) Selecting Initial Cluster Centroids
In order to achieve stable and reliable clustering results, we use the maximum distance, which can avoid obtaining adjacent elements, to select the initial centroids. For a set of sequences, {s_{1}, s_{2}, …, s_{N}}, let d(s_{i}, s_{j})(i, j = 1, 2, …N) be the distance between any two sequences in the dataset. We choose the sequence and as the cluster centroid according the following rule:
2) Selecting the Cluster to Split
BKM algorithm usually partitions the largest size cluster into two smaller ones and yields clusters with similar size. However, a cluster with large number is not always the loose one. If one existing cluster is a loose one, in which its members are not closely related to each other, the cluster will be selected to be split.
Variance is a measure of how far a set of numbers are spread out from each other, and it can measure the compactness of the clusters. So we select the cluster to split on the basis of the compactness of clusters measured by variance. The variance of cluster C_{j} is defined as following:
where μ_{j} is the centroid of sequences in C_{j}, d (s_{i}, μ_{j}) is the distance between s_{i} and μ_{j}, and n_{j} is the number of sequences in the cluster.
A small variance of a cluster indicates that the members in the cluster tend to be closely related to the mean. In other words, the smaller the variance is, the more compact the cluster is, and vice versa.
Based on the above idea, we outline mBKM algorithm as follows.
Algorithm Name: mBKM for clustering sequences
Input: sequences {s_{1}, s_{2}, …, s_{N}}, a distance function d between sequences, the number of clusters K.
Output: Set of K clusters.
Steps:
1. Initialization: Regard the whole dataset {s_{1}, s_{2}, …, s_{N}} as a single cluster.
2. Pick a cluster to split.
3. Find two subclusters:
3.1 Select two initial centroids using Equation (6);
3.2 Assign the sequences to the closest centroid;
3.3 Recalculate two centroids based on the sequences assigned to the cluster;
3.4 Repeat steps 3.2 and 3.3 until no change in cluster centroid calculation.
4. Calculate the variance of each cluster according Equation (7) and take the split that produces the clustering result with the highest variance.
5. Repeat steps 2, 3 and 4 until the desired number K is reached.
This algorithm outputs a binary tree of sequences, where each leaf represents a sequences and each node represents a sequence collection.
Results and discussion
The proposed method is evaluated by clustering functionally related gene sequences and by phylogenetic analysis. We present our evaluation results in two parts. The first one aims at testing the efficiency of our similarity measure, DMk. The second one is to illustrate the efficiency of the proposed clustering method, mBKM.
To measure the quality of the clustering results, our experiments adopt Fmeasure [52] to evaluate the clustering performance. For cluster j and class iF (ij) is defined as:
where i =1, 2, …, ej = 1, 2, …, fprecision(i, j) = n_{ij}/n_{j}recall(i, j) = n_{ij}/n_{i}e is the number of classes, and f is the number of clusters. n_{ij} is the number of the sequences of class i in cluster jn_{i} is the number of the sequences of class i, and n_{j} is the number of the sequences of cluster j.
The Fmeasure of the whole clustering result is defined as:
where N is the total number of sequences in the data set. Clearly, an Fmeasure has a value between 0 and 1. The larger the Fmeasure is, the better the clustering result is.
Evaluation of similarity measure
To evaluate the proposed similarity measure, we test DMk on gene sequence data sets and compare it with the ktuple distance. We also verify the effectiveness of DMk by assessing how well it performs on phylogenetic analysis.
Gene sequences clustering
Genes of the same family usually share similar sequences, functional domains, and even interacting partners. When a new gene is assigned to a cluster, the biological function of this cluster can be attributed to this gene with high confidence.
Four data sets are extracted from different gene repositories as shown in Table 1. The sequences of DS1 are downloaded from NCBI (http://www.ncbi.nlm.nih.gov webcite). The other three datasets, DS2, DS3 and DS4, are taken from PBIL (http://pbil.univlyon1.fr/ webcite). DS2 is taken from HOVERGEN of PBIL, a database of homologous vertebrate genes. DS3 is taken from HOGENOM, which contains homologous gene families from microbial organisms. DS4 is randomly selected from HOMOLENS, a database of homologous genes from Ensembl organisms and Ensembl families.
Table 1. Description for the Data Sets
Four widely used clustering algorithms, including KM, singlelinkage clustering (SL), completelinkage clustering (CL) and averagelinkage clustering (AL), have been chosen in the experiments. For comparison, we perform the clustering tests on all data sets using the ktuple distance and DMk distance. In this paper, we set k value to 3. For protein coding genes, a tuple size of 3 is a good choice according to reference [22]. We also tested the clustering performance on different k values, and the result confirms that a small k value is preferred, see Additional file 1: Table S1. For larger k values, there are more tuples with zero frequencies and less information is captured by the algorithm.
Additional file 1. Supplementary Data.
Format: DOC Size: 340KB Download file
This file can be viewed with: Microsoft Word Viewer
KM algorithm would yield different results during multiple executions due to its stochastic feature for initialization. We examine KM in ten runs and report the average performance. The AL, CL and SL hierarchical algorithms generate one solution for each of them. We obtain the result of hierarchical clustering algorithms by analyzing the hierarchical tree using the expected number of cluster as input parameters.
According to Table 2, the Fmeasure values for each of the data sets using DMk are clearly higher than those obtained with the ktuple distance. In our experiments, on average, the value of the Fmeasure given by DMk is 18% better than by the ktuple distance (p = 0.0165, onesided paired ttest) in KM, 49.7% better in SL (p = 0.0028), 24.9% better in CL (p = 0.016), and 35.8% better in AL (p = 0.01885). Clearly, DMk provides a significant improvement in clustering sequences. On the four data sets, the Fmeasure of DMk is improved more than 20% compared with that of the ktuple distance during the same clustering process in most cases. DMk outperforms the ktuple distance in the experiments. This is because DMk considers the occurrence, location and order relation of tuples in sequence and can capture more information in the sequence, while the ktuple distance considers frequency alone and ignore the position of tuples in a sequence. In addition, we have tested DMk and ktuple measures on protein sequences with a k value of 2, and the results indicate that DMk performs better than ktuple distance (data not shown). Thus in practical DMk measure can also be applied in clustering protein sequences after tuning current algorithm.
Table 2. The Fmeasures of the Data Sets
Phylogenetic analysis
In this experiment, the proposed similarity measure DMk is further tested by phylogenetic analysis. In order to evaluate the similarity measures, we use UPGMA in the PHYLIP package, a widely used clustering algorithm in phylogenetic analysis. The tree is drawn by TREEVIEW program [53].
The selected data set includes the full βglobin gene sequences of 10 species reported by Feng et al. [54], which are downloaded from NCBI (http://www.ncbi.nlm.nih.gov webcite). Their names, accession numbers, locations and lengths are listed in the Additional file 1: Table S2. The similarity/dissimilarity matrices for the full sequences of βglobin gene of the 10 species using DMk are shown in Table 3, respectively. The smaller the distance is, the more similar the two sequences are.
Table 3. The similarity/dissimilarity matrix for the 10 full βglobin gene sequences based on DMk
In Table 3, the most similar species pairs are humangorilla, humanchimpanzee and gorillachimpanzee, which are expected from their evolutionary relationship. A slightly less similar species pair is goatbovine. On the other hand, gallus is separated from the rest, this coincides with the fact that gallus is the only nonmammalian species among these 10 species. We can also find that opossum is far away from the remaining mammals. These results are consistent with biological morphology.
The quality of the constructed tree shows the quality of the distance matrix and the method of abstracting information from DNA sequences. In Figure 1v, we show the phylogenetic tree of 10 βglobin gene sequences based on DMk, generated by UPGMA. For comparison, the phylogenetic tree of the ktuple distance is shown in Figure 1(a).
Figure 1. The phylogenetic trees for 10 species using the full DNA sequences of βglobin.
The tree in Figure 1 (a) has some consistencies with biological morphology. Although it supports the separation of gallus relative to other species, its obvious drawback is that it fails to separate (mouse, rat) and (goat, bovine) from opossum. From Figure 1 (b), gallus is separated from the rest and opossum is far away from the other species. This topology is in good agreement with that presented by Feng et al. [54] and Cao et al. [55] except for the relative position of rodents.
DMk measures the similarity between DNA sequences more effective than the ktuple distance. This is because DMk measures the distance between DNA sequences based on sequence structure and composition. Through evaluation on gene families and constructing phylogenetic trees of full gene sequences of 10 species, we find that DMk gives more competitive results compared to the ktuple distance.
Evaluation of clustering methods
To evaluate the effectiveness of the proposed clustering algorithm, mBKM, we apply mBKM in clustering gene sequences and compare it with several clustering algorithms. Moreover, we use our method, mBKM with similarity measure DMk, in phylogenetic analysis to show how well the genes are grouped together and how well the resulting trees agree with existing phylogenies.
Performance comparison of clustering methods
In order to illustrate the efficiency of mBKM in gene sequence clustering, we ran mBKM with the ktuple distance and DMk on real data sets listed in Table 1. The clustering results are compared with those of KM, SL, CL, AL and BKM algorithms. For BKM, the number of iterations for each bisecting step is set to 5. We ran BKM 10 times to obtain the average Fmeasure. By combing the six clustering algorithms with two similarity measures, we have 12 combinations of clustering algorithm for performance assessment. The combinations are KM with ktuple, SL with ktuple, CL with ktuple, AL with ktuple, BKM with ktuple, mBKM with ktuple, KM with DMk, SL with DMk, CL with DMk, AL with DMk, BKM with DMk and mBKM with DMk.
The clustering performance of different clustering methods is the result of a combination of factors, including the types of sequence distances used for clustering and the choice of clustering algorithms. Table 2 shows the clustering performance on the data sets for all 12 clustering methods. For each data set, we set the number of cluster as the real number of class during the clustering run. For example, the real number of cluster is 8 in DS1 and 6 in DS2.
From Table 2, we observe that mBKM using DMk achieves best result and clearly outperforms other methods for the four data sets. The average Fmeasure of mBKM with ktuple is about 2.2% higher than KM with ktuple (p = 0.036), 45% higher than SL with ktuple (p = 0.00195), 11.4% higher than CL with ktuple (p = 0.0424), 19% higher than AL with ktuple (p = 0.08615) and 2.3% higher than BKM (p = 0.0141). For mBKM with DMk, Fmeasures for DS1, DS2, DS3, and DS4 are 0.808, 0.9645, 0.9143, and 0.9587 respectively. On average, the value of Fmeasure given by mBKM is 14.2% better than KM (p = 0.00025), 21.3% better than SL (p = 0.0105), 15.4% better than CL (p = 0.02835), 10.1% better in AL (p = 0.0686), and 2.3% higher than BKM (p = 0.0015) respectively. These results show that our method, combining mBKM with DMk, is able to achieve high quality results on all the data sets.
Because the clustering methods listed in Table 2 use the numbers of cluster as input parameters, we analyze the effects of varying the number of clusters on the clustering performance. This analysis is applied to DS1, DS2, DS3 and DS4 datasets and all 12 combinations. Figures 2 and 3 show the results of these runs based on the ktuple distance and DMk, respectively. The data used for generating these figures are included in Additional file 1: Tables S3S10.
Figure 2. The distribution of Fmeasure as a function of the number of clusters based on the ktuple distance (The real numbers of DS1, DS2, DS3 and DS4 are 8, 6, 6, and 6, respectively).
Figure 3. The distribution of Fmeasure as a function of the number of clusters based on DMk (The real numbers of DS1, DS2, DS3 and DS4 are 8, 6, 6, and 6, respectively).
Figure 2 illustrates the results of the six clustering algorithms with the ktuple distance. From Figure 2 and Additional file 1: Tables S3S6, mBKM achieves better Fmeasures than other five clustering algorithms for the real number of clusters on all the data sets. Although the other clustering algorithms give slightly better results in terms of Fmeasure in some cases, mBKM performs better than the other clustering algorithms in terms of the average of the Fmeasures values (average values are shown in Additional file 1: Tables S3S6). This result shows that on average, mBKM performs better than other clustering algorithms for a range of cluster numbers, in the vicinity of real number of clusters. It also implies that varying the number of clusters as input for these clustering algorithms does not affect the performance.
Figure 3 shows the results of clustering algorithms with DMk. mBKM obtains the highest Fmeasure values among the six clustering algorithms at the real number of clusters. On average, mBKM achieves better results than the other clustering algorithms for DS2, DS3, and DS4. For DS1, the average value of mBKM is very close to that of AL and higher than those of the other clustering algorithms. Overall mBKM produces consistently high quality clusters in the neighborhood of the real number of cluster (data shown in Additional file 1: Tables S7S10). The Fmeasures given by mBKM are higher than those of other clustering methods at the corresponding number of clusters in most cases.
From Figures 2 and 3, we can see that DMk achieves better cluster quantity than the ktuple distance in terms of Fmeasure. Using same clustering algorithm on the same data set, DMk achieves higher average of the Fmeasure values than the ktuple distance, and DMk also obtains higher Fmeasures at corresponding number of clusters (data shown in Additional file 1: Tables S3S10). From both Figures, we find that Fmeasure changes as the number of cluster changes. As it is known, Fmeasure is a balanced measure of precision and recall. It is an ideal condition when the number of cluster is equal to the real number. When the number of cluster is greater than or less than the real number, the Fmeasure will be affected.
With regard to clustering algorithms, SL performs poorly in many cases, and this may be because that SL uses the nearest pair of sequences and may lead to bad splits of one cluster if two or more clusters show different pattern densities. For KM and BKM, the results of many runs are lower than those of mBKM. On the whole, mBKM achieves better results than other clustering algorithms, and mBKM combining with DMk achieves best results among these clustering methods in our experiments.
The task of sequence clustering is to group given sequences into clusters. The similarity measure, DMk, measures the similarity between DNA sequences based solely on the ktuple. It is more effective than the ktuple distance, which is one of the most widely used methods. The clustering algorithm, mBKM, can obtain better clustering results and can reveal the relationships among clusters in hierarchical manner. In the next experiments, we combine mBKM with DMk to clustering DNA sequences.
In order to further illustrate the efficiency of our method, combining mBKM and DMk, we compare mBKM with DMk to two other clustering programs: BlastClust [27] and CDHITEST [39]. BlastClust is an alignmentdependent clustering algorithm. BlastClust is from NCBI Blast package. BlastClust accepts a number of parameters that can be used to control the clustering stringency including thresholds for score density (−S parameter), and alignment length (−L parameter). CDHITEST is a popular DNA clustering program based on greedy incremental clustering method. CDHITEST groups DNA sequences into clusters that meet a userdefined similarity threshold (−c parameter) and uses shortword filters to rapidly determine that if two sequences are similar, which reduces the number of full alignments necessary.
We perform tests using BlastClust and CDHITEST on the data sets listed in Table 1. In order to obtain the best possible performance of BlastClust, we set p as F (input type is nucleotide sequence) and vary the input parameters, S and –L, to evaluate the results. The score density, –S parameter, varies between 10 and 90 with step size 10, and the alignment length, –L parameter, varies between 0.1 and 0.9 with step size 0.1. Other parameters are kept default. For CDHITEST, because the sequence identity threshold, c parameter, should be greater than or equal to 0.8 in the program, we vary c parameter between 0.8 and 1 with step size 0.02, and set the word length as default value. The best results from different parameter combination are recorded. For mBKM with DMk, we set the size of ktuple as 3 and use the real number of clusters as input. As BlastClust and CDHITEST do not use the number of clusters as input, we choose the resulting class i, which has the max F(i,j) for cluster j, to calculate the Fmeasures. The results, which contain the corresponding Fmeasures and the execution time, are summarized in Table 4.
Table 4. Clustering results on the data sets listed in Table 1
Table 4 demonstrates that mBKM with DMk produces good results relative to each original cluster set in terms of Fmeasure. Every Fmeasure of mBKM with DMk is higher than 0.8 and the highest is 0.9645. It is also seen in the table that mBKM with DMk outperforms BlastClust and CDHITEST on all the data sets. BlastClust and CDHITEST tend to give more clusters than the real numbers of classes, therefore, BlastClust and CDHITEST give high precision and low recall value. But neither of these two performs well in terms of Fmeasure. The execution times reported in Table 4 for algorithm comparison show mBKM with DMk is faster than BlastClust and CDHITEST.
For the cases that the real number of clusters is unknown, the performance of our algorithm will be affected. In order to compare with BlastClust and CDHITEST on a relatively fair ground, we can vary the number of clusters and take the average of the Fmeasure values over the different numbers of clusters. For instance, we run mBKM with DMk with the range of 3–20 numbers and the average values of Fmeasure are 0.7065, 0.8533, 0.8205 and 0.8429 for DS1, DS2, DS3 and DS4, respectively. As shown in Additional file 1: Tables S7S10, these values are also higher than the corresponding Fmeasure of BlastClust and CDHITEST.
Phylogenetic analysis
In this experiment, we used mBKM with DMk to construct phylogenetic trees.
1) The clustering result of 10 species
We apply mBKM with DMk to the 10 DNA sequences of βglobin gene in Table 4. The clustering result is shown in Figure 4(a). Using the same data set, we also build the phylogenetic tree using CLUSTALW [56] and MUSCLE [57] for alignment, and UPGMA and Maximum Likelihood (ML) method (in the PHYLIP package) for presenting the tree. Figure 4(b) and 4(c) shows the tree built by CLUSTALW with UPGMA and MUSCLE with ML respectively. The trees built by MUSCLE with UPGMA and CLUSTALW with ML are provided in Figure 1 of Additional file 1.
Figure 4. The phylogenetic trees for 10 species using the full DNA sequences of βglobin.
In Figure 4(a), human, gorilla, chimpanzee and lemur are closer to bovine and goat than to mouse and rat, this topology is in complete agreement with Feng et al. [54] and Cao et al. [55] confirming the outgroup status of rodents relative to ferungulates and primates. Moreover, the tree in Figure 4(a) is identical to the tree in Figure 4(b)4(c) and the tree built MUSCLE with UPGMA. In experiment, the branch (bovine, goat) is not classified well by CLUSTALW with ML. Furthermore, it took about 0.1 second for our method. However, UPGMA with CLUSTALW and MUSCLE for the same data set took 5.1 and 1.2 seconds to build the tree, respectively, and ML with CLUSTALW and MUSCLE took 8 and 4.1 seconds to build the tree, respectively.
2) The Clustering result of 60 H1N1 viruses
H1N1 is subtype of the influenza A virus which can cause illness in humans and many other animal species. Analysis of H1N1 is critical for preparing a strategy to prevent and to control influenza epidemics and pandemics. The H1N1 avian influenza is characterized by its continuous antigen variation, which is mainly caused by the HA and NA proteins in which HA protein has highest rate of mutation. HA protein plays a critical role in identifying and adsorbing the host cell receptor in the infection process, and it is the decisive factor of host specific. We use our method to verify the phylogenetic relationships of H1N1, and the result is included in Additional file 1. The clustering result using mBKM with DMk is shown in Figure 5(a). As a comparison, we also use CLUSTALW with UPGMA and MUSCLE with ML to construct the phylogenetic tree and they are presented in Figure 5(b) and 5(c).
Figure 5. The phylogenetic trees for 60 H1N1 viruses.
As is seen from Figure 5(a), 60 H1N1 viruses are distinctly divided into four main groups using our method. The four groups, include European swine older than 2009 (G1), the avian older than 2009 (G2), American swine older than 2009 (G3) and the new 2009 viruses from human, swine and avian (G4). The result shows that the new 2009 human H1N1 viruses have closer relationship with old American swine than old avian and European swine. This grouping result is generally consistent with the topology given by CLUSTALW with UPGMA, which is shown in Figure 5(b), and the one presented by MUSCLE with UPGMA, which is provided in the Additional file 1, as well as the result suggested by zhao et al. [58]. Figure 5(c), built by MUSCLE using ML method, also shows the new 2009 human H1N1 viruses have close relationship with old American swine except the position of the group (old avian swine, European swine) is different from the positions in Figure 5(a) and 5(b). CLUSTALW with ML (in Additional file 1) also classifies the 60 H1N1 viruses into four groups except that swine/Wisconsin/1961 and swine/Wisconsin/1961 are not classified well.
Our method analyzed the 60 H1N1 viruses within 1 second, while UPGMA with CLUSTALW and MUSCLE of the same data set took 460 and 60.1 seconds to build the tree, and ML with CLUSTALW and MUSCLE took 571 and 188.1 seconds to build the tree, respectively.
Our method, mBKM with DMk, performs well when clustering 10 species and 60 H1N1 viruses. It obtains similar results to the alignmentbased method. Furthermore, our method is much faster than the alignmentbased methods.
In order to compare the speed of our method with the multiple sequence alignment based methods, CLUSTALW and MUSCLE, we performed the test on two sets of sequences. The first set consists of six datasets. All the six datasets include 100 sequences. The lengths of all sequences in the six datasets are around 1000, 2000, 3000, 4000, 5000 and 6000 respectively. Another set also consists of six datasets. The number of sequences in each dataset is 20, 40, 60, 80, 100, 120 respectively; the lengths of all the sequences are around 3000. Because ML method is slower than UPGMA, we use UPGMA to build the phylogenetic tree of the results from CLUSTALW and MUSCLE and record the time used for each method. The results in Figure 6 show that our method is much faster than the other two methods. The actual time differences are much higher than the visual differences in the figure since we are using the log(time) as the label of yaxis.
Figure 6. The time comparison of three methods.
Scalability test
For DMk, the time complexity of transforming the gene sequence s_{1}⋯s_{l} to a vector is O (l4^{K}), thus the time complexity of generating the vectors for the whole sequence database is , where is the average length of the sequences and N is the number of sequences. The value of k set to 3 yields good results in our experiments, and we fix k to 3 as the size of ktuple. DMk have linear time complexity with respect to both and N.
The time consumed for mBKM calculation is primarily determined by choosing the initial cluster centroids. For N sequences, this step has a time complexity of O (N^{2}). The time complexity of clustering step in mBKM is O (N log K). The following scalability test on our method, mBKM with DMk, confirms that our method has linear time complexity with respect to the average length of the sequences. The scalability test uses theoretical model sequences composed of the four symbols ‘A’, ‘C’, ’G’ and ‘T’. The method is implemented in Java and on a computer with 3.00 GHz CPU and 2 GB RAM.
Figure 7(a) illustrates the relationships between the runtime and the number of sequences (implemented on a computer with 8 GB RAM). To test the scalability with respect to the number of sequences, we use five data sets which consist of 5000, 10000, 15000, 20000, 25000, 30000, 35000 and 40000 sequences. Each data set contains 10 clusters and all the sequences have the same length, 100. The curve in Figure 7(a) is primarily consistent with the time complexity of mBKM with O (N^{2}). The scalability with respect to the length of sequences was tested on five datasets with five different sequence lengths: 10000, 20000, 30000, 40000, 50000 and each set consists of 4 clusters and 100 sequences. The sensitivity with respect to the length of the sequence is illustrated in Figure 7(b), from which we can see that the time of our method increases linearly when the length of sequences increases.
Figure 7. The relationship between the runtime and different numbers of sequences and length of sequences.
Conclusions
In this paper, we presented a novel approach for DNA sequence clustering, mBKM, based on a new sequence similarity measure, DMk, which is extracted from DNA sequences based on the position and composition of oligonucleotide pattern. The experimental results show the method of combining mBKM with DMk is effective in classifying DNA sequences with similar biological characteristics and in discovering the underlying relationship among the sequences. In addition, DMk can achieve comparable or better accuracy than the frequencybased distance measure. Our proposed method can be applied to study gene families and it can also help with the prediction of novel genes. Furthermore, mBKM with DMk can generate cluster trees that are useful to understand the processes governing the gene evolution. In addition, our method may be extended for protein sequence analysis and metagenomics of identifying source organisms of metagenmic data. Our method has limitations too. For example, the method did not consider edge length, and has not address problems with long repeated sequences or long insertions. In future we will try to address these problems.
Competing interests
The authors declare that there are no competing interests.
Authors’ contributions
DW designed the algorithm, conducted the experiments, and wrote the manuscript. QJ supervised the project and proposed data mining algorithm. YW guided the experiments, wrote the manuscript and analyzed the results. SW guided the experiment analysis, and proposed ideas for sequence clustering algorithm. All authors read and approved the final manuscript.
Acknowledgements
This work is supported by the National Natural Science Foundation of China under Grant No.61175123 and No.10771176, and the Shenzhen New Industry Development Fund under grant No.CXB201005250021A.
References

Demuth JP, De Bie T, Stajich JE, Cristianini N, Hahn MW: The evolution of mammalian gene families.
PLoS One 2006, 1:110. Publisher Full Text

Zhao B, Duan V, Yau SS: A novel clustering method via nucleotidebased Fourier power spectrum analysis.
JTheor Biol 2011, 279:8389. Publisher Full Text

Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ: A basic local alignment search tool.

Pearson WR, Lipman DJ: Improved tools for biological sequence comparison.
ProcNatlAcad Sci USA 1988, 85(8):24442488. Publisher Full Text

Vinga S, Almeida J: Alignmentfree sequence comparisona review.
Bioinformatics 2003, 19(4):513523. PubMed Abstract  Publisher Full Text

Haubold B, Reed FA, Pfaffelhuber P: Alignmentfree estimation of nucleotide diversity.
Bioinformatics 2011, 27(4):449455. PubMed Abstract  Publisher Full Text

Liu Z, Meng J, Sun X: A novel featurebased method for whole genome phylogenetic analysis without alignment: application to HEV genotyping and subtyping.
Biochem Biophys Res Commun 2008, 368(2):223230. PubMed Abstract  Publisher Full Text

DomazetLoso M, Haubold B: Efficient estimation of pairwise distances between genomes.
Bioinformatics 2009, 25(24):32213227. PubMed Abstract  Publisher Full Text

DomazetLoso M, Haubold B: Alignmentfree detection of local similarity among viral and bacterial genomes.
Bioinformatics 2011, 27(11):14661472. PubMed Abstract  Publisher Full Text

Kelil A, Wang S, Brzezinski R, Fleury A: CLUSS: Clustering of protein sequences based on a new similarity measure.
BMC Bioinformatics 2007, 8:286. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Reinert G, Chew D, Sun FZ, Waterman MS: Alignmentfree sequence comparison (I): statistics and power.

Dai Q, Liu X, Yao Y, Zhao F: Numerical characteristics of word frequencies and their application to dissimilarity measure for sequence comparison.
JTheor Biol 2011, 276(1):174180. Publisher Full Text

Lu G, Zhang S, Fang X: An improved string composition method for sequence comparison.
BMC Bioinformatics 2008, 9(Suppl 6):S15. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Aita T, Husimi Y, Nishigaki K: A mathematical consideration of the wordcomposition vector method in comparison of biological sequences.
BioSystems 2011, 106:6775. PubMed Abstract  Publisher Full Text

Blaisdell BE: A measure of the similarity of sets of sequences not requiring sequence alignment.
ProcNatlAcad Sci USA 1986, 83:51555159. Publisher Full Text

Wu TJ, Burke JP, Davison DB: A measure of DNA sequence dissimilarity based on Mahalanobis distance between frequencies of words.
Biometrics 1997, 53(4):14311439. PubMed Abstract  Publisher Full Text

Wu TJ, Hsieh YC, Li LA: Statistical measures of DNA dissimilarity under Markov chain models of base composition.
Biometrics 2001, 57(2):441448. PubMed Abstract  Publisher Full Text

Stuart GW, Moffett K, Baker S: Integrated gene and species phylogenies from unaligned whole genome protein sequences.
Bioinformatics 2002, 18(1):100108. PubMed Abstract  Publisher Full Text

Fichant G, Gautier C: Statistical method for predicting protein coding regions in nucleic acid sequences.

Wang J, Zheng X: WSE, a new sequence distance measure based on word frequencies.
Math Biosci 2008, 215(1):7883. PubMed Abstract  Publisher Full Text

Zheng X, Qin Y, Wang J: A Poisson model of sequence comparison and its application to coronavirus phylogeny.
Math Biosci 2009, 217(2):159166. PubMed Abstract  Publisher Full Text

Yang K, Zhang L: Performance comparison of gene family clustering methods with expect curated gene family data set in Arabidposis thaliana.
Planta 2008, 228:439447. PubMed Abstract  Publisher Full Text

Dong G, Pei J: Classification, clustering, features and distances of sequence Data.
Sequence Data Mining 2007, 33:4765. Publisher Full Text

Sokal RR, Rohlf FJ: Biometry: The Principles and Practice of Statistics in Biological Research. 3rd edition. W. H. Freeman and Company, New York; 1995.

Everitt BS, Landau S, Leese M: Cluster Analysis. Oxford University Press, London; 2001.

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(13):i41i49. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

National Center for Biotechnology Information (NCBI):
Documentation of the BLASTCLUSTalgorithm.
ftp://ftp.ncbi.nih.gov/blast/documents/blastclust.html webcite

Enright AJ, Ouzounis CA: GeneRAGE: a robust algorithm for sequence clustering and domain detection.
Bioinformatics 2000, 16(5):451457. PubMed Abstract  Publisher Full Text

Chaudhuri P, Das S: SWORDS: A statistical tool for analyzing large DNA sequences.
J Biosci 2002, 27(1):16. PubMed Abstract  Publisher Full Text

Uchiyama I: Hierarchical clustering algorithm for comprehensive orthologousdomain classification in multiple genomes.
Nucleic Acids Res 2006, 34(2):647658. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hao X, Jiang R, Chen T: Clustering 16 S rRNA for OTU prediction: a method of unsupervised Bayesian clustering.
Bioinformatics 2011, 27(5):611618. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Sperisen P, Pagni M: JACOP: a simple and robust method for the automated classification of protein sequences with modular architecture.
BMC Bioinformatics 2005, 6:216. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Wan X, Bridges SM, Boyle JA, Boyle AP: Interactive Clustering for Exploration of Genomic Data.

Kelarev A, Kang B, Steane D: Clustering Algorithms for ITS Sequence Data with Alignment Metrics.
Lect Notes ComputSci 2006, 4304:10271031. Publisher Full Text

Tseng GC: Penalized and weighted Kmeans for clustering with scattered objects and prior information in highthroughput biological data.
Bioinformatics 2007, 23(17):22472255. PubMed Abstract  Publisher Full Text

Ashlock D, Warner E: Classifying Synthetic and Biological DNA Sequences with Side Effect Machines.
IEEE Symposium on Computational Intelligence in Bioinformatics and Computational Biology: 15–17 Sept. 2008; Sun Valley, ID 2008, 2229.

Criterion functions for document clustering: experiments and analysis. 2001.
Technical Report

Kashef R, Kamel MS: Enhanced bisecting kmeans clustering using intermediate cooperation.
Pattern Recognit 2009, 42(11):25572569. Publisher Full Text

Li W, Godzik A: Cdhit: a Fast Program for Clustering and Comparing Large Sets of Protein or Nucleotide Sequences.
Bioinformatics 2006, 22(13):16581659. PubMed Abstract  Publisher Full Text

Alam I, Cornell M, Soanes DM, Hedeler C, Wong HM, Rattray M, Hubbard SJ, Talbot NJ, Oliver SG, Paton NW: A Methodology for Comparative Functional Genomics.

Picardi E, Mignone F, Pesole G: EasyCluster: a fast and efficient geneoriented clustering tool for largescale transcriptome data.
BMC Bioinformatics 2009, 10(Suppl 6):S10. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Wei D, Jiang Q: A DNA Sequence Distance Measure Approach for Phylogenetic Tree Construction.
In 5th IEEE International Conference on BioInspired Computing: Theories and Applications: 23–26 Sept. 2010 Changsha Edited by Li K, Tang Z, Li R, Nagar AK, Thamburaj R. 2010, 204212.

NeumannHeld EM: The gene is deadLong live the gene: Conceptualizing genes the constructionist way. In Sociobiology and Bioeconomics: the Theory of Evolution in Biological and Economic Theory. Edited by Koslowski P. Springer, Berlin; 1999:105137.

White JR, Navlakha S, Nagarajan N, Ghodsi M, Kingsford C, Pop M: Alignment and clustering of phylogenetic markers implications for microbial diversity studies.
BMC bioinformatics 2010, 11:152. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Waterman MS: Introduction to Computational Biology: Maps, Sequences, and Genomes. Chapman and Hall, Lodon; 1995.

Durbin R, Eddy SR, Krogh A, Mitchison G: Biological Sequence Analysis: probabilistic models of proteins and nucleic acids. Cambridge University Press, Cambridge; 1998.

Nair ASS, Mahalakshmi T: Visualization of genomic data using internucleotide distance signals.
Proceedings of IEEE Genomic Signal Processing: 11–13 July 2005; Bucharest 2005.

Schmitt AO, Herzel H: Estimating the Entropy of DNA sequences.
JTheor Biol 1997, 188(3):369377. Publisher Full Text

Steinbach M, Karypis G, Kumar V: A comparison of document clustering techniques.

Zhao Y, Karypis G: Hierarchical Clustering Algorithms for Document Datasets.
Data Mining Knowl Discov 2005, 10:141168. Publisher Full Text

Larsen B, Aone C: Fast and effective text mining using lineartime document clustering.
In Proceedings of the fifth ACM SIGKDD international conference on Knowledge discovery and data mining: 15–18 August 1999; San Diego Edited by Fayyad U, Chaudhuri S, Madigan D. 1999, 1622.

Page RD: TreeView: an application to display phylogenetic trees on personal computers.
Bioinformatics 1996, 12:357358. Publisher Full Text

Feng J, Hu Y, Wan P, Zhang A, Zhao W: New method for comparing DNA primary sequences based on a discrimination measure.
JTheor Biol 2010, 266(4):703707. Publisher Full Text

Cao Y, Janke A, Waddell PJ, Westerman M, Takenaka O, Murata S, Okada N, Paabo S, Hasegawa M: Conflict among individual mitochondrial proteins in resolving the phylogeny of eutherian orders.

Chenna R, Sugawara H, Koike T, Lopez R, Gibson TJ, Higgins DG, Thompson JD: Multiple sequence alignment with the Clustal series of programs.
Nucleic Acids Res 2003, 31(13):34973500. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Edgar RC: MUSCLE: multiple sequence alignment with high accuracy and high throughput.
Nucleic Acids Res 2004, 32(5):17921797. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Zhao B, He RL, Yau SS: A new distribution vector and its application in genome clustering.