Genomic signal processing for DNA sequence clustering

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 ${\hat{X}}^{\alpha }$ by applying Eq. (1) (1)${\displaystyle \begin{array}{c}{\hat{\mathbf{X}}}_1\left(i\right)=\left\{\begin{array}{cc}1\hfill & \text{if}X\left(i\right)=A\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}\right.\hfill \\ {\hat{\mathbf{X}}}_2\left(i\right)=\left\{\begin{array}{cc}1\hfill & \text{if}X\left(i\right)=G\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}\right.\hfill \\ {\hat{\mathbf{X}}}_3\left(i\right)=\left\{\begin{array}{cc}1\hfill & \text{if}X\left(i\right)=C\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}\right.\hfill \\ {\hat{\mathbf{X}}}_4\left(i\right)=\left\{\begin{array}{cc}1\hfill & \text{if}X\left(i\right)=T\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}\right.\hfill \end{array}}$

By applying the Discrete Fourier transform to the DNA signal ${\hat{X}}^{\alpha }$, it is possible to compute the power spectral density (PSD) ${\hat{S}}^{\alpha }$ which describes how power of a signal or time series is distributed over frequency.¹ 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 $d\left({\hat{S}}^{\alpha },{\hat{S}}^{\beta }\right)$ 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 {O₁, O₂, …, O_m} represented as a n-dimensional vector, into K ≤ m clusters. Each cluster is represented by a centroid C_j with $j\in \left[1,2,\dots ,k\right]$, 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 C_j 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 Ω = [ω₁, ω₂, …, ω_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:

1. Compute a main centroid point M in the n- dimensional space corresponding to the geometrical center of the K centroids location computed as: (2)${\displaystyle M\left[i\right]=\frac{1}{k}\sum _{j=1}^k{C}_j\left[i\right]}$ where i ∈ [1, 2, …, n].

2. For each cluster j, compute the Euclidean distance d_j of its centroid C_j with respect to the main centroid M: (3)${\displaystyle d_j=\sqrt[]{\sum _{i=1}^n{\left(C_j\left[i\right]-M\left[i\right]\right)}^2}.}$

3. Each centroid of the k clusters is sorted according to its distance to the main centroid and an angle is assigned to them, according to its index $\iota \in \left[0,1,\dots ,k\right]$ in the sorted array: (4)${\displaystyle {\theta }_{\iota }=\iota \frac{2\pi }{k}.}$

4. The main centroid M and the clusters centroids C_ι are mapped into a two dimensional space ϕ, where the main centroid corresponds to the origin (i.e., the point with coordinates (x = 0, y = 0)).

5. Each centroid C_ι is plotted as a point around the main centroid M point, according to its distance and its angle as computed by: (5)${\displaystyle x_{\iota }=d_{\iota }cos\left({\theta }_{\iota }\right)}$ (6)${\displaystyle y_{\iota }=d_{\iota }sin\left({\theta }_{\iota }\right).}$

6. We sort each set of DNA sequences in Ω assigned to a specific centroid C_ι, according to the distance δ_z of each sequence z, with respect to its assigned centroid. The angle θ_z is also computed similarly to step 3.

7. Finally, each sequence z is then plotted into ϕ by computing their correspondent coordinates as: (7)${\displaystyle x_z={\delta }_{z}cos\left({\theta }_z\right)+x_{\iota }}$ (8)${\displaystyle y_z={\delta }_{z}sin\left({\theta }_z\right)+y_{\iota }.}$

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.² 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.