Genomic Signal Processing Methods for Computation of Alignment-Free Distances from DNA Sequences

Genomic signal processing (GSP) refers to the use of digital signal processing (DSP) tools for analyzing genomic data such as DNA sequences. A possible application of GSP that has not been fully explored is the computation of the distance between a pair of sequences. In this work we present GAFD, a novel GSP alignment-free distance computation method. We introduce a DNA sequence-to-signal mapping function based on the employment of doublet values, which increases the number of possible amplitude values for the generated signal. Additionally, we explore the use of three DSP distance metrics as descriptors for categorizing DNA signal fragments. Our results indicate the feasibility of employing GAFD for computing sequence distances and the use of descriptors for characterizing DNA fragments.

The second group includes methods for which the numerical values are defined according to certain biophysical or biochemical properties of the DNA molecules. Examples of this type of mapping include the use of electron-ion interaction potentials (EIIP) (i.e., A = 0.1260, C = 0.1340, T = 0.1335, G = 0.0806) [12] and the use of single atomic numbers (i.e., A = 70, C = 58, T = 66, G = 78) [13]. Other examples include paired nucleotide representations that consider nucleotide complementarity (i.e., A = T = 0, C = G = 1) [14] and graphical approaches such as the DNA-walk model, in which a step is taken upwards (+1) if S(t) is a pyrimidine (C or T) or downwards (21) if it is a purine (A or G). Finally, this category also includes Z-curve representation [15], which maps a DNA sequence into a 3-dimensional sequence where Dx t distinguishes between purines and pyrimidines; Dy t distinguishes between amino-type and keto-type molecules, and Dz t distinguishes between weak and strong hydrogen bonds.
Most GSP methods reported in the literature are focused on the detection of coding regions (e.g., [11,[16][17][18][19][20][21][22][23]). In general, these methods consist of performing DNA-to-signal mapping and obtaining the power spectrum of sections of the signal employing the short time discrete Fourier transform (STFT) using a sliding window of fixed length. When a period-3 frequency peak is detected in the power spectra, the section of sequence corresponding to that window is labeled as a coding region. Other applications of GSP include searching for genomic repeats using STFT [24] and determining the structural, thermodynamic, and bending properties of DNA by Fourier analysis [25].
Determining the distance between different genomic sequences is one of the most common types of analysis. In this scope, phylogenetic trees are one of the most essential tools in DNA analysis because they provide structured classification of DNA sequences and enable organization of our growing knowledge of biological diversity. Moreover, this method provides insight into events that occur during evolution. A phylogenetic tree may be constructed from a distance matrix M containing a set of values that represents the pairwise distance d(S i ,S j ) of a set of sequences S~½S 1 ,S 2 ,:::,S m . Examples of distance matrix-based methods for phylogenetic tree construction include neighbor-joining [26] and the Fitch-Margoliash method [27], among others.
A distance-matrix corresponding to a set of DNA, RNA, or protein sequences is commonly determined by assessing the distance based on alignment of sequence pairs. Alignment methods have also been used to identify domains, assemble genome contigs, and study sequence variations. Techniques for determining the alignment of a pair of sequences include dotmatrix, dynamic programming, and k-tuple methods. In dot matrix-based methods, a recurrence plot is generated by comparing all elements of both sequences to form a twodimensional matrix in which a dot is placed at the intersection where characters match. Dynamic programming methods compute the optimal alignment between two sequences by considering possible differences due to mutations, insertions, and deletions. This method can also be used for global alignments via the Needleman-Wunsch (NW) algorithm [28] or local alignments via the Smith-Waterman [29] algorithm. The k-tuple (word) method attempts to identify sub-sequences of length k in the query sequence. Although this method does not guarantee an optimal solution, it is significantly more efficient than dynamic programming, making it suitable for the analysis of large-scale databases. Two of the most popular local alignment methods are FASTA [30] and BLAST [31]. A GSP method has also been proposed for aligning multiple sequences (i.e., MAFFT [32]). In this method, amino acid sequences of different proteins are converted into two numerical vectors consisting of values that correspond to the volume and polarity of the components. The correlation between the two amino acid sequences is computed by the fast Fourier transform (FFT) using a sliding window of fixed length. By assessing the correlation score of both sliding windows, it is possible to detect regions of matching sequences.
Several alignment-free methods for DNA distance computation have been proposed. In general, these methods are based on statistics of word frequencies (i.e., k-tuples) using metrics such as weighted Euclidean distance, correlation, co-variance, information theory-based measurements, and angle metrics. [33]. However, other methods based on graphical DNA representations apply dinucleotide (doublet) histograms [34], graph theory [35], trinucleotide (triplet) curves [36], or the average bandwidth of distance/distance (D/D) matrices [37]. A widely used tool for computing phylogenetic trees is the phylogeny inference package (Phylip) [38], which applies different methods such as parsimony, jackknife, bootstrapping, and consensus trees using molecular sequences, gene frequencies, restriction sites and fragments, distance matrices, and discrete characters.
In this paper, we present a novel GSP alignment-free method (GAFD) for determining the distance between two DNA sequences. We introduce a new DNA-to-signal mapping tool that is based on using doublets with a mapping function inspired by Kstrings [39], which increases the number of possible amplitude values for the generated signal. Additionally, we explore three GSP distance metrics that may be used as descriptors for categorizing dissimilarities between pairs of DNA signal fragments, and that set the basis for developing methods for domain search and characterizing sections of DNA. Our results demonstrate that GAFD performs similarly to the NW and Phylip methods for computing distances among a set of DNA sequences. Moreover, the results obtained using the proposed descriptors show the feasibility of this method in characterizing the types of differences present between sections of sequences. All the methods and algorithms were implemented in MATLAB R2010b. We employed the ARfit module for the autoregressive model computation and signal processing toolbox of MATLAB for FFT computation. NW analysis and phylogenetic tree construction was computed using the bioinformatics toolbox of MATLAB. Source code is available for download at: http://hypatia.cucei. udg.mx/invteorica/DNASignals/.

DNA sequence-to-signal mapping
Our proposed DNA sequence-to-signal mapping method was inspired by the alignment-free distance method, which employs the nearest-neighbor method (NN) [40,41]. NN was originally developed for determining the double strand melting temperature [42,43] and was based on the rationale that the interaction between bases on different strands depends to some extent on the neighboring bases. The model assumes that, under specific environmental conditions, the stability of hydrogen bridges between strands of a nucleic acid duplex for a given doublet and its complementary pairs depends on the identity of its neighboring bases.
Our mapping tool requires that numerical values for all possible combinations of two consecutive bases (doublets) are defined. Let S(t) indicate the nucleotide at position t with respect to the beginning of a sequence S of length m. For each doublet, we define a numeric value H(S(t), S(tz1)). Then, a one-dimensional discrete signalŜ S(t) is generated in a manner similar to the K-strings approach [39] by combining the values generated by doublets inside a window of magnitude 2az1: After mapping the DNA sequence to a discrete signal, a noise reduction method is applied. A typical solution for noise reduction from non-stationary signals is the wavelet denoising method. Wavelet-based noise decomposition of a signal using orthogonal discrete wavelet transform (DWT) can ''concentrate'' the informative signal into a few wavelet coefficients with large absolute values without modifying the random distribution of noise. Then, DWT-based denoising can be achieved by limiting the number of wavelet coefficients that represent the signal. Consider the model whereŜ S(t) represents the original discrete signal, c(t) represents the noiseless unknown version ofŜ S(t), and n(t) represents the noise. Since DWT is a linear transformation, the wavelet coefficient vectors for each term in Eq. (2) (i.e., wŜ S , w c , and w n ) are related by: Denoising is performed by computing the wavelet transform of a signal and then removing the coefficients that correspond to high frequencies by applying a threshold T. The wavelet coefficients corresponding to low frequencies remain unchanged. The main challenge of this denoising method is determining the T value between small and large wavelet coefficients. To determine this threshold, several algorithms have been proposed recently [44,45]. In this work, we employed the SureShrink algorithm (i.e., Stein's unbiased risk estimator [46]) because this wavelet denoising method has been successfully applied to biological signals without the loss of important information [47].

Alignment-free distance computation
Each de-noised signalŜ S i corresponding to a DNA sequence in a set was converted into its frequency representation by applying discrete Fourier transform (DFT) followed by computation of its power spectrumF F i . To perform a direct comparison of sequences of different length, we employed zero padding to compute the DFT using the maximum length of all DNA signals in the set. Subsequently, for a given pair of DNA signals, we calculated the alignment-free distance d(i,j) by computing the mean squared error (MSE) of their corresponding power spectrum: Finally, a distance matrix M was created by performing pairwise comparison of all sequences in the set.

Cluster overlapping score
Because the distance matrices generated with GAFD and NW cannot be directly compared, we evaluated these methods by employing a modified version of the cluster overlapping score [48], which is based on the splitting method by Robinson and Foulds. An NJ tree is a dichotomic hierarchic classifier that can be defined by a set of leaves T~l 1 ,l 2 , . . . ,l t f gand a set of branches B~B 1 ,B 2 , . . . ,B t{1 f g . Every branch of the tree is represented in the subset B~l 1 ,l 2 , . . . ,l b f g , where t~b for the root, and twb §2 for all other branches of the tree. In this report, a cluster G is defined as a branch that is not the root and lies beyond the immediate vicinity of a leaf twg §3. Note that the set of clusters G~G 1 ,G 2 , . . . ,G g È É is a subset of B. Given two clusters from two different trees (T 1~T2^G1 =G 2 ) where G T1 \G T2 =0 = , cluster overlapping can be measured using the Dice coefficient as: The cluster overlapping similarity score k is calculated from the arithmetic mean of all k occurrences. Comparison of a tree against itself yields a value of 1.

GSP distance descriptors
This proposed method for distance computation can determine the global similarity of two DNA sequences. However, this metric may not be capable of determining local differences between the sequences. Therefore, we explored three DSP-based metrics used to generate descriptors that may be useful for characterizing and classifying differences between sections of the sequences.
Correlation coefficient. Correlation can be used to measure the dependency of a signal on itself or another signal. For two signals s i and s j with the same length n, the correlation coefficient is defined as: Depending on the data to be evaluated, the correlation coefficient will be r*1 for signals that are highly correlated, r*0 for non-correlated signals, and r*{1 for signals that are inversely correlated. Because we were only concerned with the degree of correlation and not its type, we discarded the sign and defined the descriptor as R~DrD.
Coherence. Coherence, a relationship measurement used to estimate the degree of linear association between two signals is defined as: where s ij is the cross-spectral density that describes the common power distribution between the two signals, while s ii and s jj denote the auto-spectral density of s i and s j respectively, at a frequency t. We defined the descriptor C as the mean of the coherences c(t) Vt). A C value close to 0 indicates that signals at this frequency are linearly independent, whereas a value close to 1 represents a very high linear correlation. In this work, the spectral densities used for determining the coherence between two signals were computed using an autoregressive (AR) model, which is one of the most widely used tools in DSP [49]. For a given interval, the multidimensional AR model is given by: where A(k)~a 1 (k),a 2 (k), . . . ,a n (k) ½ T is the n|n AR coefficients matrix, n the number of channels, x(t{k) the time-delayed values vector, p the model order, and e(t) the error vector. To solve Equation 8, it is necessary to fit p (smaller than the sequence length) and then estimate the AR coefficient matrix [50].
Derivative comparison. Given that our proposed DNA sequence-to-signal mapping accounts for neighboring nucleotides, small differences between s i and s j due to indels and mutations produce a shift in the intensity of the resulting signal (Fig. 1). To account for these changes, we compared the derivatives of the two signals by using finite differences and computed the mean slope D as a descriptor representing the degree of similarity between the two signals. A value D*0 indicates strong similarity.
Similarity space. For a given pair of sequences, the three descriptors were used to generate a point with coordinates (R,C,D) in a three-dimensional space. We hypothesize that it may be possible to characterize sub-sequences by employing clustering or classification methods in this space.

Results and Discussion
DNA sequence-to-signal mapping Our DNA sequence-to-signal mapping tool requires that different values be set for every possible doublet (i.e., 16 different values). For all the experiments presented in this section, we employed the values listed in Table 1. The proposed DNA sequence-to-signal mapping was designed to consider the nucleotides within a window defined by a. An example of the effect of a on the proposed mapping is depicted in Figure 1. As the value of a increases, the resulting DNA signal becomes smoother as the values corresponding to nucleotides within the window are combined. This indicates that the value of a determines how far a change is propagated through the signal. Note that a single nucleotide substitution produces a vertical shifting effect depending on the value of a with respect to a signal corresponding to a similar sequence. As a increases, a substitution has less impact on the signal shift. Indels in the DNA sequence are reflected as a horizontal shift with respect to another similar sequence proportional to the number of deleted or inserted bases. Figure 2 depicts the distance as computed by GAFD with respect to different numbers of changes in a given sequence employing different values of a. Note that, compared to methods that perform DNA sequence-to-signal mapping using individual nucleotides, a determines the robustness of our method with respect to subtle differences between the sequences being evaluated. In this work, we chose to employ a~3 since this value allows us to distinguish between different numbers of signal changes.
In a sense, substitution matrices may be considered equivalent to simple-DNA-mapping functions. Where the former assigns a value to the difference between two nucleotides on different DNA sequences, the latter replaces the nucleotide with a number. Therefore, when comparing two mapped sequences, the difference between two nucleotides would be represented by their corresponding numbers. It becomes evident then that simple-mapping functions reduce the degrees of freedom when element-wise comparisons are made. The main advantage when mapping DNA is the ability to treat it as a series, which allows for the use of DSP and other concepts such as context information-dependent entities.
In this paper, we gathered contextual information using the NN algorithm, where each nucleotide value is considered along with its neighboring nucleotides. This approach, albeit still unidimensional, increases the degrees of freedom with respect to simple-DNAmappings while also including contextual information. Although this is only a ''proof of concept'' study, we hypothesize that the referring context may not be only local, but also distant. Moreover, it may contain the sequence information itself as well as data from different annotation levels according to related known ontology or metadata. Adding the contextual elements will improve analysis by encompassing the DNA grammar into the structural information.

GSP distance computations using GAFD
To evaluate the performance of GAFD with respect to existing methods for computing sequence distances, we assessed the similarity of unrooted phylogenetic trees generated by the NJ method [51,52] (equal variance and independence of evolutionary distance estimates) using distance matrices computed with GAFD and NW (nuc44 scoring matrix, gap penalty of 8, and use of Jukes-Cantor for the maximum likelihood estimate of the number of substitutions) of various DNA sequences belonging to different organisms. In addition, we computed phylogenetic trees employing the Phylip method using ordinary parsimony and without randomization, with a search for the best 100 trees. The Phylip method was fed with sequences aligned using ClustalW with gap open penalty = 10, gap extension penalty = 0.05, and no weight transition. The resulting tree typologies were compared using the previously described cluster overlapping score k.
Examination of the ribosomal S18 subunit gene. Two experiments were performed by analyzing two sets of DNA sequences corresponding to the ribosomal S18 subunit (KEGG orthology K02964). This gene was selected because it is the broadest evolutionary marker discernable between all eukaryotes. In the first experiment, three basic clusters were built, namely mammals, insects, and plants, according to general taxonomy. The resulting phylogenetic trees generated from the distance matrices computed by the three methods are depicted in Figure 3. Note that the eutheryan (a mammal subgroup) were grouped in GAFD, NW, and Phylip. However, the insects were grouped differently by the three methods (e.g., Nasonia vitripennis was located far outside the other insects according to GAFD and Phylip). These results are consistent with the known complexity of insect genetics due to horizontal transference, spurious recombination, and high variability rate. Note that NW represented the outside eukaryote Saccharomyces cerevisiae appropriately, while GAFD placed it incorrectly among the plant group. Phylip placed this sequence in an outer group next to Monodelphis domestica and N. vitripennis. Although M. domestica was expected to be placed  in an external group within mammals, it was placed in the outer branch of all trees. Lastly, with the exception of S. cerevisiae in GAFD, all plants were properly clustered. In the second experiment, all entries of the aforementioned orthology were compared. A total of 149 organisms and 231 entries were analyzed, resulting in mean similarity scores of k~93:12 between GAFD and NW, k~94:19 for GAFD and Phylip, and k~95:75 for NW and Phylip.
Assessment of other evolutionary markers. In this experiment, we selected evolutionary markers corresponding to coding (i.e., 21 tRNA synthetases and 2 ribosomal proteins) and non-coding (i.e., 20 tRNAs and 2 rRNAs) genes. We included species present in all KEGG orthologies and then selected all entries belonging to these organisms. We constructed and compared the phylogenetic trees generated using GAFD, NW, and Phylip. Figure 4 depicts two examples of trees generated by NW and GAFD for two selected orthologies (tRNA-Asp and tRNA-GLU). Note the similarity in gene clustering by GAFD and NW. Tables 2, 3, and 4 list the similarity scores k for the noncoding tRNAs, coding tRNA synthetases, and coding/non-coding ribosomal genes, respectively. The mean scores for the non-coding genes were k~0:85+0:06, while k~0:82+0:06 was exhibited for the coding genes. In general, the cluster overlapping scores   between the methods were relatively high, indicating that GAFD can group similar sequences effectively. We assessed statistical significance by applying the nonparametric Wilcoxon Signed-Rank test, which does not require any assumptions regarding the normality of the data distribution. The null hypothesis for our tests is that the median differencek k between the similarity of pairs of evaluations is not significant (H 0 :k k 1~k k 2 ), while the alternative hypothesis is the statistically significant difference between both medians (H 1 :k k 1 =k k 2 ). The resulting p-values for each test at a significance level of 0:05 are listed in Table 5. No significant differences were observed between the three methods when examining coding sequences. However, for non-coding genes, GAFD performed differently from the other two methods. This may be related to the fact that coding genes appear to have a certain periodic structure [11,[16][17][18][19][20][21][22][23] which will affect GAFD since it also considers the frequency content of the mapped sequence. Figure 5 depicts the times required to compute the distance matrices using NW and GAFD on a desktop PC (i-Core 7, 2GHz, 6 GB RAM) for different numbers of sequences. GAFD performed faster than NW despite the implementation of high level MATLAB code. We believe that this performance could be improved by employing low level coding (e.g., C++) and tools such as GPU and parallel computing. A comparison of computer times for Phylip was not necessary because this method does not compute a similarity matrix.
Several elements may impact GAFD results to different extents. Even if they are related, these elements can be grouped according to the source of two critical phenomena, namely those related to the method, such as sequence-to-signal mapping and zero padding, and those related to the nature of the genomic grammar. Regarding the former, GAFD sequence-to-signal mapping is based on uniformly-euclidean unidimensional mappings. Even though comparisons will sometimes yield different results under particular sequence conditions, higher-order NN mappings are very difficult to implement and analyze. Additionally, since biologically meaningful sequences may not be of the same length, power spectra comparisons are influenced by zero padding. Regarding the latter, genomic information is typically full of challenging sequences, i.e. palindromes, inversions, translocations, repeats, duplications, and indels. All of these will exhibit distinct characteristics in power spectra that will in turn lead to inconsistencies when several sequence comparisons are performed. For example, differences among inversions, translocations, and palindromes may not be observed, while repetitions and duplications will display specific frequency peaks. How the interaction between all of these elements affects GAFD analysis is outside the scope of this paper. Moreover, NN mapping using contextsensitive information and DNA distance determination through power spectra comparison should be explored in the future. GAFD is not intended for sequence alignment, but rather for comparing them in another domain and rendering a similarity value. We believe that, after refinement, this approach will enable us to discover relationships between sequences that are not bound to the sequence itself, but to specific underlying patterns in the genomic grammar that is yet to be fully understood.

GSP distance descriptors
To explore the three-dimensional space generated by the proposed descriptors (R,C,D), we performed an experiment in which we perturbed a randomly generated DNA sequence S r that generates a DNA signalŜ S r of length b~20. Using S r as the ''mother sequence'', we generated all the DNA sequences and signals corresponding to all possible combinations of one, two, and three changes, considering all possible types of changes (i.e., substitutions, deletions, and insertions). Every pair of signals generated a point (R,C,D) in the this space ( Figure 6). Our results from the comparisons corresponding to one change were located near the origin, while those corresponding to two or three changes were positioned at increasing distance from the origin according to the number of changes. Additionally, the points corresponding to substitutions were well-separated from those corresponding to insertions and deletions ( Figures 6 D and 6 E). These results demonstrate that GAFD can characterize the type of change present using a classification technique that combines several descriptors. However, coherence exhibited the poorest results since a lack of specificity for detecting insertions and substitutions was observed. This result is supported by Sims, et al. [53], where it was reported that optimal resolutions (length of b) proved critical for genomic comparisons. Moreover, studies have shown that coherence AR models depend highly on the parameters employed [54].
As a preliminary domain search assay, we conducted another experiment using real data (i.e., ribosomal S18 subunit sequences from the previous 26 selected species). The signal corresponding to the Homo sapiens sequence was segmented into non-overlapping fragments of length b~20 to generate a ''signal dictionary''. From the dictionary, seven entries were selected at random and compared against the complete signal set employing a sliding window of length b. For each position within the sliding window, we computed the proposed similarity descriptors. We considered the segment of signal contained within the sliding window as similar to that from the dictionary if the correlation and coherence descriptors were larger than 0.9 and the comparison derivative was less than 0.8. The resulting alignment schematic is depicted in Figure 7. Even when the fragments were selected randomly, our results provide evidence that most mammals share similar fragments. Note that the number of shared fragments decreases as the sequences become less related to the original sequence (i.e., H. sapiens). Interestingly, insects shared the least number of fragments. These data suggest that it may be possible to determine biologically significant elements among compared sequences. Figure 6. Depiction of the similarity space created with the three GSP distance descriptors. A random 20 nt sequence was created. Using this sequence as a template, all possible combinations for up to three substitutions were created and measured against the template using the three distance descriptors. The dots in A, B, and C correspond to the distances for one (red), two (blue), and three (green) substitutions, respectively. As expected, the more substitutions present, the farther they scattered along the frequency peak. Subsequently, starting with the same template, all possible combinations of insertions, deletions, and substitutions were created and measured similarly as aforementioned. The dots in D, E, and F correspond to the distances for insertions (yellow), deletions (brown), and substitutions (green). The distance scatters shift between substitutions and indels, which is especially evident in the Correlation and Derivative descriptors. The blue scatter on A through C is equal to the green scatter on D through F. doi:10.1371/journal.pone.0110954.g006