Introduction

Cluster analysis is one of the most common and useful tools in pattern recognition, statistical data analysis, and exploratory data mining. It has many applications such as image segmentation, recognition of objects, document retrieval and others (Jain, Murty & Flynn, 1999). The main advantage of employing clustering techniques is the possibility of finding a hidden structure in the data without the requirement of prior information or knowledge about it. A clustering task consists of dividing a dataset into groups (i.e., clusters) that share common properties or that are related in some way, according to given criteria and similarity metrics (Baikey, 1994).

The most popular method used to perform cluster analysis is the K-means algorithm (Jain, 2010). K-Means clustering is an iterative partition technique which finds mutual exclusive spherical groups (Joshi & Kaur, 2013). The main advantage of the K-means algorithm is its ease of implementation and its linear time complexity (Jain, Murty & Flynn, 1999). However, the K-means algorithm rely on the frequent computation of similarity metrics between all of the elements to be clustered and the proposed centroids of each of the k-clusters. Therefore, its application in practice is limited to the type of data for which those similarity scores can be computed in a efficient way.

In bioinformatics, traditional methods for computing the similarity scores between sequences consist of applying DNA and amino acid sequence alignment methods, whose main objective is to identify portions of successive nucleotide or amino acids that are common in two or more sequences. They are then rearranged to easily visualize those similar portions (White et al., 2010). The comparison of two sequences is known as pairwise sequence alignment (PSA). When more than two sequences are compared, the process is known as multiple sequence alignment (MSA) (Sharma, 2008).

One of the most popular applications of PSA is phylogenetic analysis. It consists of establishing an evolutionary relationship among nucleic acid or protein families sequences. It is generally depicted by the use of dichotomous trees, for which the branches represent organism separations. Branches that are close to each other, suggest a similar organism. By contrast, the farthest branches indicates large differences (Mount, 2004). Some of the most popular algorithms for MSA are ClustalW (Thompson, Higgins & Gibson, 1994), Muscle (Edgar, 2004), T-COFFEE (Notredame, Higgins & Heringa, 2000), MAFFT (Katoh et al., 2005), and K-Align (Lassmann & Sonnhammer, 2005).

However, since these methods require large computational times for determining similarity among sequences, the use of K-means is not feasible for this application. Therefore, other approaches for DNA clustering have been proposed based on the use of these similarity computation methods. Two of the most popular algorithms for clustering biological sequences are the CD-HIT (Li & Godzik, 2006) and the UCLUST (Edgar, 2010). Both algorithms use a greedy approach for identifying representative sequences that can be used as a “seed” to group all of the sequences that have a similarity score above a certain threshold. However, the computational resources necessary to perform the multiple sequence alignments remain the main challenge which limits the number of sequences that can be clustered.

More recently, an approach for the analysis of genomic data that has captured the attention of researchers in recent years, is the use of genomic signal processing (GSP) which is based on the use of digital signal processing (DSP) theory and algorithms to analyze DNA or protein sequences. GSP methods require the transformation or mapping of the biological sequences, usually represented as a string of characters (i.e., A, T, G and C) to a numeric representation (i.e., a signal) that can be processed using mathematical functions (Kwan & Arniker, 2009). Examples of the use of GSP methods include the identification of protein-coding regions in DNA sequences (Das & Turkoglu, 2017; Mabrouk, 2017; Das & Turkoglu, 2015; Inbamalar & Sivakumar, 2012; Marhon & Kremer, 2011; Akhtar, Epps & Ambikairajah, 2008; Akhtar, Epps & Ambikairajah, 2007; Rushdi & Tuqan, 2006; Yin & Yau, 2005; Kotlar, 2003; Anastassiou, 2000), finding for genomic repeats (Sharma et al., 2004), determining the structural, thermodynamic, and bending properties of DNA (Gabrielian & Pongor, 1996), biological sequence querying (Ravichandran et al., 2010), estimating of DNA sequence similarity (Mendizabal-Ruiz et al., 2017; Hoang, Yin & Yau, 2016; Yin, Yin & Wang, 2014; Borrayo et al., 2014; Cheever et al., 1989), and sequence alignment (Skutkova et al., 2015).

One of the main advantages of GSP methods is that the analysis of the genomic data can be performed very quickly because of the optimal coding of the algorithms and the processors that have been designed specifically for those tasks.

Cluster analysis of DNA signals through the use of GSP methods have been previously proposed by Zhao, Duan & Yau (2011) and Hoang et al. (2015). However, these methods are based on the computation of a number of features from the Fourier spectrum which may reduce the dimensionality of the data and perhaps its discriminative power as compared with the use of the whole raw spectrum. Moreover, those works employed a hierarchical clustering algorithm instead of the K-means. Comparatively, K-means properties allow us to generate plots that are different from the traditional dendrograms and that facilitate the exploration of the results.

In this paper, we propose an approach for performing cluster analysis of DNA sequences that is based on the use of GSP methods and the K-means algorithm. We also present a visualization method that allows us to easily inspect and analyze the results. Our results indicate the feasibility of employing the proposed method to find and easily visualize interesting features of sets of DNA data.

DNA sequence to signal

In order to be able to employ the DSP methods in genomic data, it is necessary to first perform a transformation or mapping of the DNA sequences to be analyzed into numerical values representing the information contained by them. There are several proposed DNA numerical representations. However, one of the most popular of this DNA to signal mapping is the Voss representation, which employs four binary indicator vectors, each meant to denote the presence of a nucleotide of each type at a specific location within the DNA sequence (Voss, 1992).

Given a DNA sequence α (e.g., α = ATTCGCAT...) we can employ the Voss representation to compute its corresponding fourth-dimensional DNA signal \documentclass[12pt]{minimal}

\usepackage{amsmath}

\usepackage{wasysym}

\usepackage{amsfonts}

\usepackage{amssymb}

\usepackage{amsbsy}

\usepackage{upgreek}

\usepackage{mathrsfs}

\setlength{\oddsidemargin}{-69pt}

\begin{document}

}{}${\hat {X}}^{\alpha }$\end{document}X ˆα by applying Eq. (1)

(1)\documentclass[12pt]{minimal}

\usepackage{amsmath}

\usepackage{wasysym}

\usepackage{amsfonts}

\usepackage{amssymb}

\usepackage{amsbsy}

\usepackage{upgreek}

\usepackage{mathrsfs}

\setlength{\oddsidemargin}{-69pt}

\begin{document}

}{}\begin{eqnarray*}\begin{array}{@{}l@{}} \displaystyle {\hat {\mathbf{X}}}_{1}(i)= \left\{ \begin{array}{@{}ll@{}} \displaystyle 1&\displaystyle \text{if} X(i)=A\\ \displaystyle 0&\displaystyle \text{otherwise} \end{array} \right. \vskip{.5pc}\\ \displaystyle {\hat {\mathbf{X}}}_{2}(i)= \left\{ \begin{array}{@{}ll@{}} \displaystyle 1&\displaystyle \text{if} X(i)=G\\ \displaystyle 0&\displaystyle \text{otherwise} \end{array} \right. \vskip{.5pc}\\ \displaystyle {\hat {\mathbf{X}}}_{3}(i)= \left\{ \begin{array}{@{}ll@{}} \displaystyle 1&\displaystyle \text{if} X(i)=C\\ \displaystyle 0&\displaystyle \text{otherwise} \end{array} \right. \vskip{.5pc}\\ \displaystyle {\hat {\mathbf{X}}}_{4}(i)= \left\{ \begin{array}{@{}ll@{}} \displaystyle 1&\displaystyle \text{if} X(i)=T\\ \displaystyle 0&\displaystyle \text{otherwise} \end{array} \right. \end{array}\end{eqnarray*}\end{document}X ˆ1i=1ifXi=A0otherwiseX ˆ2i=1ifXi=G0otherwiseX ˆ3i=1ifXi=C0otherwiseX ˆ4i=1ifXi=T0otherwise

By applying the Discrete Fourier transform to the DNA signal \documentclass[12pt]{minimal}

\usepackage{amsmath}

\usepackage{wasysym}

\usepackage{amsfonts}

\usepackage{amssymb}

\usepackage{amsbsy}

\usepackage{upgreek}

\usepackage{mathrsfs}

\setlength{\oddsidemargin}{-69pt}

\begin{document}

}{}${\hat {X}}^{\alpha }$\end{document}X ˆα, it is possible to compute the power spectral density (PSD) \documentclass[12pt]{minimal}

\usepackage{amsmath}

\usepackage{wasysym}

\usepackage{amsfonts}

\usepackage{amssymb}

\usepackage{amsbsy}

\usepackage{upgreek}

\usepackage{mathrsfs}

\setlength{\oddsidemargin}{-69pt}

\begin{document}

}{}${\hat {S}}^{\alpha }$\end{document}S ˆα which describes how power of a signal or time series is distributed over frequency.1

1For further details regarding the formal definition of a PSD refer to Stoica & Moses (2005).In our case, the PSD is a descriptor of the nucleotide patterns that may be present within the DNA sequence (Borrayo et al., 2014).

The relatedness or similarity score of any two given DNA sequences α and β, can then be estimated by comparing the components of their PSDs \documentclass[12pt]{minimal}

\usepackage{amsmath}

\usepackage{wasysym}

\usepackage{amsfonts}

\usepackage{amssymb}

\usepackage{amsbsy}

\usepackage{upgreek}

\usepackage{mathrsfs}

\setlength{\oddsidemargin}{-69pt}

\begin{document}

}{}$d({\hat {S}}^{\alpha },{\hat {S}}^{\beta })$\end{document}dS ˆα,S ˆβ using a similarity metric (Mendizabal-Ruiz et al., 2017).

DNA signal clustering

K-means is a two step algorithm which performs the partitioning of a given set of observations {O1, O2, …, Om} represented as a n-dimensional vector, into K ≤ m clusters. Each cluster is represented by a centroid Cj with \documentclass[12pt]{minimal}

\usepackage{amsmath}

\usepackage{wasysym}

\usepackage{amsfonts}

\usepackage{amssymb}

\usepackage{amsbsy}

\usepackage{upgreek}

\usepackage{mathrsfs}

\setlength{\oddsidemargin}{-69pt}

\begin{document}

}{}$j\in \left[ 1,2,\ldots ,k \right] $\end{document}j∈1,2,…,k, which is defined as a point in a n-dimensional space generated by computing the average of each element of the vectors of the observations that belong to that cluster. In the first step, an observation is assigned to the cluster Cj that scores the highest similarity to the point represented by the observation’s vector, according to a specific metric. In the second step, the centroids of the k clusters are updated, according to the observations assigned to them in the previous step. The best groups and their centroids are obtained by the minimization of the total sum of the distances between the observations and their corresponding centroids.

Consider a set of PSD Ω = [ω1, ω2, …, ωm] corresponding to a number m of different DNA sequences. The K-means algorithm is applied to the data in Ω by considering the power spectra as the vector that describes the DNA sequence in a n-dimensional space. In this work, we chose the Euclidean distance between these vectors as the similarity metric to be employed by the K-means algorithm. Since the K-means results depends on the initial labels assigned to each entry, which are assigned randomly, we repeat the computation 50 times and keep the best convergence score. As a result, we obtain a label for each element of Ω which defines the assigned cluster.

DNA clusters visualization

The raw results of the clustering procedure may be difficult to analyze and interpret. Therefore, we propose to produce graphical representations of the results that can easily provide an insight into the DNA sequence clustering results. The generation of the proposed graphical representation (Fig. 1) from the K-means clustering result, consists of the following steps:

Experimental data

To assess our DNA signal clustering method and the proposed visualization technique, we employed a set of 141 DNA sequences corresponding to the Cytochrome c oxidase I gene (COXI) marker belonging to 131 different species obtained from the Kyoto Encyclopedia of Genes and Genomes (KEGG) K02256 (Kanehisa et al., 2017; Kanehisa et al., 2016; Kanehisa & Goto, 2000). We selected the COXI marker because it performs a fundamental role in the terminal oxidative step for energy metabolism (Adkins & Honeycutt, 1994) and is a very well known marker commonly used for the identification of species (Patwardhan, Ray & Roy, 2014). In the selected set, a total of 112 organisms have only one copy. However, there are some species represented by more than one sequence. This is the case of the Yak, Bos grunniens (bom:102267288, bom:102278784, bom:22161768), a Bat, Myotis davidii (myd:102771221, myd:22203924), the Spotted green pufferfish, Tetraodon nigroviridis (tng:BAE79219, tng:GSTEN00036010G001), the Pacific giant oyster, Crassostrea gigas (crg:109618508, crg:109618509, crg:808829), Yarrowia lipolytica (yli:YalifMp03, yli:YalifMp05, yli:YalifMp06), Loa loa, the parasite responsible for filariasis disease (loa:COX1, loa:LOAG_19059), the Castor oil tree, Ricinus communis (rcu:10221395, rcu:8272741), and the Picoplanktons, Ostreococus tauri (ota:OstapMp24, ota:OstapMp40), Bathycoccus prasinos (bpg:BathyMg00110, bpg:BathyMg00240), and Micromonas commoda (mis:MicpuN_mit45, mis:MicpuN_mit7). It is important to note that all gene copies were considered during the experiments and that the selected organisms belong to the total spectrum of the Eukaryote domain.

The selected organisms were manually organized according to their respective taxon, based on the Catalogue of Life (Roskov et al., 2017), and were divided into seven kingdoms, 17 phyla, and 35 classes. To easily identify the different categories, we employed different colors and symbols as described in Fig. 2

Comparison with other cluster methods

We evaluated the performance of the MATLAB implementation of proposed algorithm “Signal Tool for the Analysis of the Relationship between Sequences” (i.e., STARS) in terms of computational time with respect to ClustalW and UCLUST. While ClustalW is not strictly a clustering method, we used it for comparison because it is one of the most commonly used tools to evaluate the similarity of multiple sequences. We employed a CPU Intel XEON E5-1650 at 3.50 GHz with 16 GB RAM.2

2When analyzing the processing times of the compared methods, it is important to consider that the STARS was implemented in MATLAB without any parallelization, in comparison with the highly optimized implementations of ClustalW and parallelized UCLUST.Table 1 list the processing time in seconds for the three methods for sets of 8, 17, 35, 70, and 141 sequences of COXI. The time required to transform the 141 sequences from strings of characters to their corresponding PSDs was 0.921 s and it is not considered in Table 1 since this is performed only one time. Note that the time required by STARS is significantly smaller with respect to ClustalW. UCLUST is time-constant at 1 s for every experiment, however, note that the number of clusters generated by this method was practically the same number of sequences (i.e., the method assigns a cluster to each sequence). This is because UCLUST requires a sequences identity range of at least 40% for amino acids and 65% for nucleotides (Edgar, 2010).

Table 2 list the processing times of five datasets of different number of sequences with different length where UCLUST generated a number of clusters different from one cluster for each sequence. Dataset A consisted of 31 sequences of Mammals with average length of 16,695 nucleotides labeled into to seven groups; Dataset B consisted of 38 sequences of Influenza A viruses with average length of 1,407 nucleotides labeled into to six groups; Dataset C consisted of 116 sequences of Human Rhinovirus with average length of 7,154 nucleotides labeled into four groups; Dataset D consisted of 34 Coronavirus sequences with average length of 27,567 nucleotides labeled into six groups; and Dataset E consisted of 30 sequences of Bacteria with average length of 3,361,393 labeled into eight groups. Note that the computational time required for performing the clustering of the sequences’ PSD data is smaller when compared with UCLUST for the same number of clusters for datasets A, B, and D. In the case of dataset E, we could not achieve a result using UCLUST (i.e., the program throws a fatal error) indicating that the data was too big.

Discussion

Numerous reports have discussed the use of different molecular markers to determine the appropriate phylogenetic divergence at many levels of the tree of life. For our experiments, we considered molecular markers previously employed in phylogenetic analysis for the evaluation of the differentiation capability of DNA sequences (Hoang et al., 2015) (e.g., COXI, mtDNA, influenza A virus, human rhinovirus, coronavirus and bacterial genomes). In this work, we focused on the Mitochondrial Cytochrome C Oxidase Subunit I (COXI) coding gene as a marker to evaluate our approach for group clustering of relevant similar sequences. The COXI gene has been proposed as one of the most relevant marker genes for molecular taxonomy (Patwardhan, Ray & Roy, 2014). While no single gene is even close to establishing the systematic classification of organisms, the COXI gene is one of the most closely related to the consensus evolutionary divergence. Therefore, it is important to wrap our results not as how the Catalogue of Life (Roskov et al., 2017) classification should become, but as how accurate our marker could be. The selection (COXI) was based on three criteria: (i) the marker must code for proteins since it has been already proven that these type of markers have steadier mutation rates, (ii) the marker should have already been employed in a wide range of the tree of life, at least for eukaryotes, and should be able to discriminate for the intended groups, and (iii) the marker should have a homogeneous length and have a minimum number of reported copies in the selected database, since both duplication events and large indels may bias new cluster formation. To rule out any bias in the clustering of organisms with respect to their downloaded sequences, we incorporated all of their stored copies.

A possible expected result was that the clusters generated for each selected value of k would correspond to the organization depicted in Fig. 2. However, since the K-means method promote the generation of centroids in highly populated regions of the feature space, it is more likely to obtain clusters of organisms that are highly related among them, instead of organisms related by possible common ancestors or groups with a small number of less homogeneous organisms (e.g., primates formed a cluster early in small k values and kept together at larger numbers of k).

The COXI gene is one of the most accepted general markers to establish divergence (Patwardhan, Ray & Roy, 2014). It spans from Phylum to Class, and when using introns in selected species, it has been shown to properly classify Genera and Species (Zardoya & Meyer, 1996). Our results were remarkably good in clustering up to the Family level by using only the coding region and without the need to pretreat or manually curate the sequences.

Hebert, Ratnasingham & De Waard (2003) established a divergence rate from 0.01% to 64% with a median of 8% across a number of species on 11 Metazoan phyla. In that study, Arthropod (with the exception of Lepidoptera and Diptera class) and Plathelminth phyla displayed the greatest divergences, while Chordata showed the second lowest divergence. Our results showed high cohesiveness, particularly for Chordates, where they quickly established stable, compact clusters, predominantly with their own classes. Since the downloaded data for each phylum or class was not balanced, we had the opportunity to evaluate how the sequences are clustered in a real-life condition. For example, when sampling whole ecosystems (i.e., microbiomes), bacterial populations will not be balanced across their species, but will show predominant phylogenetic diversity toward certain groups. Our results show that our presented method is very sensitive to both, the relative abundance of tight clusters, and the K-number. Far from being a disadvantage, we found that changing the number of clusters in an experiment may provide new insights about the relationship between the various sequences.

Sequence mutations of COXI coding regions have not been shown to distribute bias towards any segment or region, something like what happens to other markers such as the ribosomal 16S gene, where changes on highly conserved regions are very few and slow, while changes on hypervariable regions show rapid changes that can determine divergence along several phylogenetic groups, according to which hypervariable region is being evaluated. If this would be the case, K-Means clustering may be adapted to steps of low mutation rates before high mutation rate regions. COXI gene mutations spanning all of the sequence may increase the amount of spurious clustering due to converging hotspots.

The presence of a spurious cluster that is gathered together by their size, is an indication of the need to filter out sequences with large indels. Despite such mishaps, the proposed method is capable of performing an analysis of relationships between multiple DNA sequences with minimum handling and without the need of sequence alignment, which results in less human and computational time compared to traditional methods. We tested this method with a number of markers (i.e., mammal mtDNA, influenza A virus, human rhinovirus, coronavirus, and bacterial genomes) previously employed for Fourier DNA spectra phylogenetic analysis (Hoang et al., 2015), the results are shown in Supplemental Information 1 of this article. Briefly, most sequences evaluated under our method cluster properly and consistently with previous reports (Hoang et al., 2015). Also consistent with COXI results, the most evident aspect is the tendency to prioritize division of heavily populated groups.

The proposed method may be used to evaluate the capability of a marker or gene to differentiate between organisms at different levels, to identify subgroups within a set of organisms, and perform classification of organisms with respect to known sequences or classification of sections of a DNA sequence. Furthermore, this method can also be used to perform similar analysis with amino acid sequences.

We have demonstrated that it is possible to group DNA sequences based on their frequency components. It is the subject of future work to identify whether distinct frequency bands amount to greater weight in the clustering of sequences.

The proposed method has been coded and executed in MATLAB. The source code and the datasets employed for the results presented in this paper are available at Github

Conclusion

We have presented a method for performing cluster analysis of DNA sequences that is based on the use of GSP methods and the K-means algorithm. We also proposed a visualization method that allows us to easily inspect and analyze the results and possible nontrivial relationships. Our results indicate the feasibility of employing the proposed method to find and easily visualize interesting features of sets of DNA data.