 Software Article
 Open Access
 Published:
MLDSP: Machine Learning with Digital Signal Processing for ultrafast, accurate, and scalable genome classification at all taxonomic levels
BMC Genomics volume 20, Article number: 267 (2019)
Abstract
Background
Although software tools abound for the comparison, analysis, identification, and classification of genomic sequences, taxonomic classification remains challenging due to the magnitude of the datasets and the intrinsic problems associated with classification. The need exists for an approach and software tool that addresses the limitations of existing alignmentbased methods, as well as the challenges of recently proposed alignmentfree methods.
Results
We propose a novel combination of supervised Machine Learning with Digital Signal Processing, resulting in MLDSP: an alignmentfree software tool for ultrafast, accurate, and scalable genome classification at all taxonomic levels. We test MLDSP by classifying 7396 full mitochondrial genomes at various taxonomic levels, from kingdom to genus, with an average classification accuracy of >97%.
A quantitative comparison with stateoftheart classification software tools is performed, on two small benchmark datasets and one large 4322 vertebrate mtDNA genomes dataset. Our results show that MLDSP overwhelmingly outperforms the alignmentbased software MEGA7 (alignment with MUSCLE or CLUSTALW) in terms of processing time, while having comparable classification accuracies for small datasets and superior accuracies for the large dataset. Compared with the alignmentfree software FFP (Feature Frequency Profile), MLDSP has significantly better classification accuracy, and is overall faster.
We also provide preliminary experiments indicating the potential of MLDSP to be used for other datasets, by classifying 4271 complete dengue virus genomes into subtypes with 100% accuracy, and 4,710 bacterial genomes into phyla with 95.5% accuracy.
Lastly, our analysis shows that the “Purine/Pyrimidine”, “JustA” and “Real” numerical representations of DNA sequences outperform ten other such numerical representations used in the Digital Signal Processing literature for DNA classification purposes.
Conclusions
Due to its superior classification accuracy, speed, and scalability to large datasets, MLDSP is highly relevant in the classification of newly discovered organisms, in distinguishing genomic signatures and identifying their mechanistic determinants, and in evaluating genome integrity.
Background
Of the estimated 8.7 million (±1.3 million) species existing on Earth [1], only around 1.5 million distinct eukaryotes have been catalogued and classified so far [2], leaving 86% of existing species on Earth and 91% of marine species still unclassified. To address the grand challenge of all species identification and classification, a multitude of techniques have been proposed for genomic sequence analysis and comparison. These methods can be broadly classified into alignmentbased and alignmentfree. Alignmentbased methods and software tools are numerous, and include, e.g., MEGA7 [3] with sequence alignment using MUSCLE [4], or CLUSTALW [5, 6]. Though alignmentbased methods have been used with significant success for genome classification, they have limitations [7] such as the heavy time/memory computational cost for multiple alignment in multigenome scale sequence data, the need for continuous homologous sequences, and the dependence on a priori assumptions on, e.g., the gap penalty and threshold values for statistical parameters [8]. In addition, with nextgeneration sequencing (NGS) playing an increasingly important role, it may not always be possible to align many short reads coming from different parts of genomes [9]. To address situations where alignmentbased methods fail or are insufficient, alignmentfree methods have been proposed [10], including approaches based on Chaos Game Representation of DNA sequences [11–13], random walk [14], graph theory [15], iterated maps [16], information theory [17], categorypositionfrequency [18], spacedwords frequencies [19], Markovmodel [20], thermal melting profiles [21], word analysis [22], among others. Software implementations of alignmentfree methods also exist, among them COMET [23], CASTOR [24], SCUEAL [25], REGA [26], KAMERIS [27], and FFP (Feature Frequency Profile) [28]. While alignmentfree methods address some of the issues faced by alignmentbased methods, [7] identified the following challenges they face:

(i)
Lack of software implementation: Most of the existing alignmentfree methods are still exploring technical foundations and lack software implementation, which is necessary for methods to be compared on common datasets.

(ii)
Use of simulated sequences or very small real world datasets: The majority of the existing alignmentfree methods are tested using simulated sequences or very small realworld datasets. This makes it hard for experts to pick one tool over the others.

(iii)
Memory overhead: Scalability to multigenome data can cause memory overhead in wordbased methods, especially when long kmers are used.
To overcome these challenges, we propose MLDSP, a novel combination of supervised Machine Learning with Digital Signal Processing of the input DNA sequences, as a generalpurpose alignmentfree method and software tool for genomic DNA sequence classification at all taxonomic levels.
The main contribution of MLDSP is the feature vector that we propose to be used by the supervised learning algorithms. Given a genomic DNA sequence, its feature vector consists of the pairwise Pearson Correlation Coefficient (PCC) between (a) the magnitude spectrum of the Discrete Fourier Transform (DFT) of the digital signal obtained from the given sequence by some suitable numerical encoding of the letters A, C, G, T into numbers, and (b) the magnitude spectra of the DFT of all the other genomic sequences in the training set. The use of this new feature vector, which has not previously been used in conjunction with machine learning algorithms, allows MLDSP to significantly outperform existing methods in terms of speed, while achieving an average classification accuracy of >97%. This substantial performance improvement allows MLDSP to scale up and successfully classify much larger datasets than existing studies. Indeed, in contrast with previous benchmark datasets, each comprising less than fifty sequences, this study accurately classifies thousands of genomes from a variety of species: eukaryotic (7396 complete mitochondrial genomes), viral (4271 genomes), and bacterial (4710 genomes). In addition, this study provides the first comprehensive analysis and comparison of all thirteen onedimensional numerical representations of DNA sequences used in the Genomic Signal Processing (GSP: digital signal processing applied to genomes) literature for classification purposes. We conclude that the “Purine/Pyrimidine (PP)”, “JustA”, and “Real” numerical representations are the top three performers in terms of classification accuracy of MLDSP for our main dataset. This is surprising given that these three numerical representations do not appear to contain sufficient biological information for the accuracy attained. For example, the numerical representation “JustA” (encoding A as “1”, and G,C,T as “0”) retains the incidence and spacing for A, but not individually for the other three nucleotides.
Numerical representations of DNA sequences
Digital Signal Processing (DSP) can be employed in the context of comparative genomics because genomic sequences can be numerically represented as discrete numerical sequences and hence treated as digital signals. Several numerical representations of DNA sequences, that use numbers assigned to individual nucleotides, have been proposed in the literature [29], e.g., based on a fixed mapping of each nucleotide to a number, without biological significance; using mappings of nucleotides to numerical values deduced from their physiochemical properties; or using numerical values deduced from the doublets or codons that the individual nucleotide was part of [29, 30]. In [31, 32] three physiochemical based representations of DNA sequences (atomic, molecular mass, and ElectronIon Interaction Potential, EIIP) were considered for genomic analysis, and the authors concluded that the choice of numerical representation did not have any effect on the results. A recent study comparing different numerical representation techniques on a small dataset [33] concluded that multidimensional representations (such as Chaos Game Representation) yielded better genomic comparison results than some onedimensional representations. However, in general there is no agreement on whether or not the choice of numerical representation for DNA sequences makes a difference on the genome comparison results, or on which numerical representations are best suited for analyzing genomic data. We address this issue by providing a comprehensive analysis and comparison of thirteen onedimensional numerical representations, for suitability in genome analysis.
Digital signal processing
Following the choice of a suitable numerical representation for DNA sequences, DSP techniques can be applied to the resulting discrete numerical sequences, and the whole process has been termed Genomic Signal Processing (GSP) [30]. DSP techniques have previously been used for DNA sequence comparison, e.g., to distinguish coding regions from noncoding regions [34–36], to align genomic signals for classification of biological sequences [37], for whole genome phylogenetic analysis [38], and to analyze other properties of genomic sequences [39]. In our approach, genomic sequences are represented as discrete numerical sequences, treated as digital signals, transformed via DFT into corresponding magnitude spectra, and compared via Pearson Correlation Coefficient (PCC) to create a pairwise distance matrix.
Supervised machine learning
Machine learning has been used in smallscale genomic analysis studies [40–42], and classification analyses associated with microarray gene expression data [43–45]. In this vein, MLDSP focusses on the use of the primary DNA sequence data for taxonomic classification, and is based on a novel combination of supervised machine learning with feature vectors consisting of the pairwise distances between the magnitude spectrum of the DFT obtained from the digital signal generated from a DNA sequence, and the magnitude spectra of the DFT of the digital signals generated from all other sequences in the training set. The taxonomic labels of sequences are provided for training purposes. Six supervised machine learning classifiers (Linear Discriminant, Linear SVM, Quadratic SVM, Fine KNN, Subspace Discriminant, and Subspace KNN) are trained on these pairwise distance vectors, and then used to classify new sequences. Independently, classical MultiDimensional Scaling (MDS) generates a 3D visualization, called Molecular Distance Map (MoDMap) [46], of the interrelationships among all sequences.
For our computational experiments, we used a large dataset of 7396 complete mtDNA sequences, and six different classifiers, to compare onedimensional numerical representations for DNA sequences used in the literature for classification purposes. For this dataset, we concluded that the “PP”, “JustA”, and “Real” numerical representations were the best numerical representations. We analyzed the performance of MLDSP in classifying the aforementioned genomic mtDNA sequences, from the highest level (domain into kingdoms) to lower level (family into genera) taxonomical ranks. The average classification accuracy of MLDSP was >97% when using the “PP”, “JustA”, and “Real” numerical representations.
To evaluate our method, we compared its performance (accuracy and speed) on three datasets: two previously used small benchmark datasets [47], and a large real world dataset of 4322 complete vertebrate mtDNA sequences. We found that MLDSP had significantly better accuracy scores than the alignmentfree method FFP on all datasets. When compared to the stateoftheart alignmentbased method MEGA7 (with alignment using MUSCLE or CLUSTALW), MLDSP achieved similar accuracy but superior processing times (2250 to 67,600 times faster) for the small benchmark dataset of 41 mammalian genomes. The contrast in running time was even more extreme for the large dataset of 4322 mtDNA genomes, where MLDSP took 28 s, while MEGA7(MUSCLE/CLUSTALW) could not complete the computation after 2 h/6 h and had to be terminated.
Lastly, we provide preliminary computational experiments that indicate the potential of MLDSP to successfully classify viral genomes (4271 complete dengue virus genomes into four subtypes) and bacterial genomes (4710 complete bacterial genomes into three phyla).
Methods and implementation
The main idea behind MLDSP is to combine supervised machine learning techniques with digital signal processing, for the purpose of DNA sequence classification. More precisely, for a given set S={S_{1},S_{2},…,S_{n}} of n DNA sequences, MLDSP uses:

DNA numerical representations to obtain a set N={N_{1},N_{2},…,N_{n}} where N_{i} is a discrete numerical representation of the sequence S_{i}, 1≤i≤n.

Discrete Fourier Transform (DFT) applied to the lengthnormalized digital signals N_{i}, to obtain the frequency distribution; the magnitude spectrum M_{i} of this frequency distribution is then obtained.

Pearson Correlation Coefficient (PCC) to compute the distance matrix of all pairwise distances for each pair of magnitude spectra (M_{i},M_{j}), where 1≤i,j≤n.

Supervised Machine Learning classifiers which take the pairwise distance matrix for a set of sequences, together with their respective taxonomic labels, in a training set, and output the taxonomic classification of a new DNA sequence. To measure the performance of such a classifier, we use the 10fold crossvalidation technique.

Independently, Classical Multidimensional Scaling (MDS) takes the distance matrix as input and returns an (n×q) coordinate matrix, where n is the number of points (each point represents a unique sequence from set S) and q is the number of dimensions. The first three dimensions are used to display a MoDMap, which is the simultaneous visualization of all points in 3Dspace.
DNA numerical representations
To apply digital signal processing techniques to genomic data, genomic sequences are first mapped into discrete numerical representations of genomic sequences, called genomic signals [48]. In our analysis of various numerical representations for DNA sequences (Table 1), we considered only 1D numerical representations, that is, those which produce a single output numerical sequence, called also indicator sequence, for a given input DNA sequence.
We did not consider other numerical representations, such as binary [29], or nearest dissimilar nucleotide [49], because those generate four numerical sequences for each genomic sequence, and would thus not be scalable to classifications of thousands of complete genomes.
Discrete Fourier Transform (DFT)
Our alignmentfree classification method of DNA sequences makes use of the DFT magnitude spectra of the discrete numerical sequences (discrete digital signals) that represent DNA sequences. In some sense, these DFT magnitude spectra reflect the nucleotide distribution of the originating DNA sequences.
To start with, assuming that all input DNA sequences have the same length p, for each DNA sequence S_{i}=(S_{i}(0),S_{i}(1),…,S_{i}(p−1)), in the input dataset, where 1≤i≤n, S_{i}(k)∈{A,C,G,T}, 0≤k≤p−1, we calculate its corresponding discrete numerical representation (discrete digital signal) N_{i} defined as
where, for each 0≤k≤p−1, the quantity f(S_{i}(k)) is the value under the numerical representation f of the nucleotide in the position k of the DNA sequence S_{i}.
Then, the DFT of the signal N_{i} is computed as the vector F_{i} where, for 0≤k≤p−1 we have
The magnitude vector corresponding to the signal N_{i} can now be defined as the vector M_{i} where, for each 0≤k≤p−1, the value M_{i}(k) is the absolute value of F_{i}(k), that is, M_{i}(k)=F_{i}(k). The magnitude vector M_{i} is also called the magnitude spectrum of the digital signal N_{i} and, by extension, of the DNA sequence S_{i}. For example, if the numerical representation f is Integer (row 1 in Table 1), then for the sequence S_{1}=CGAT, the corresponding numerical representation is N_{1}=(1,3,2,0), the result of applying DFT is F_{1}=(6, −1−3i, 0, −1+3i) and its magnitude spectrum is M_{1}=(6, 3.1623, 0, 3.1623).
Figure 1a shows the discrete digital signal (using the “PP” numerical representation, row 6 of Table 1) of the DNA sequence consisting of the first 100 bp of the mtDNA genome of Branta canadensis (Canada goose, NCBI accession number NC_007011.1), and of the DNA sequence consisting of the first 100 bp of the mtDNA genome of Castor fiber (European beaver; NCBI accession number NC_028625.1). Figure 1b shows the DFT magnitude spectra of the same two signals/sequences. As can be seen in Fig. 1b, these mtDNA sequences exhibit different DFT magnitude spectrum patterns, and this can be used to distinguish them computationally by using. e.g., the Pearson Correlation Coefficient, as described in the next subsection. Other techniques have also been used for genome similarity analysis, for example comparing the phase spectra of the DFT of digital signals of full mtDNA genomes, as seen in Fig. 2 and [50, 51].
Note that, with the exception of the example in Fig. 1, all of the computational experiments in this paper use full genomes.
Pearson Correlation Coefficient (PCC)
Consider two variables X and Y (here X and Y are the magnitude spectra M_{i} and M_{j} of two signals), each of length p, that is, X={X_{0},…,X_{p−1}} and Y={Y_{0},…,Y_{p−1}}. The Pearson Correlation Coefficient r_{XY} between X and Y is the ratio of their covariance (measure of how much X and Y vary together) to the product of their standard deviations [52, 53], that is,
where the covariance of X and Y is \(\sigma _{XY}= \sum _{i=0}^{p1} \left (X_{i}\overline {X}\right)\left (Y_{i}\overline {Y}\right)/(p1)\), and the standard deviation is \(\sigma _{X} = \sqrt {\sum _{i=0}^{p1} \left (X_{i}\overline {X}\right)^{2}/(p1)}\), and similarly for σ_{Y}, where the average is defined as \(\overline {X} = \left (\sum _{i=0}^{p1} X_{i}\right)/p\) and similarly for Y. Now the formula for the Pearson Correlation Coefficient becomes
The Pearson Correlation Coefficient between X and Y is a measure of their linear correlation, and has a value between +1 (total positive linear correlation) and −1 (total negative linear correlation); 0 is no linear correlation. We normalized the results, by taking (1−r_{XY})/2, to obtain distance values between 0 and 1 (value 0 for identical signals, and 1 for negatively correlated signals). For our data sets, the PCC values between any two digital signals of DNA sequences ranged between 0 and 0.6.
For each pairwise distance calculation, the Pearson Correlation Coefficient requires the input variables (that is, the magnitude spectra of the two sequences) to have the same length. The length of a magnitude spectrum is equal to the length of corresponding numerical digital signal, which in turn is equal to the length of the originating DNA sequence. Given that genome sequences are typically of different lengths, it follows that their corresponding digital signals need to be lengthnormalized, if we are to be able to use the Pearson Correlation Coefficient. Hoang et al. avoided normalization and considered only the first few mathematical moments constructed from the power spectra for comparison, after applying DFT [54]. The limitation of this method is that one loses information that may be necessary for a meaningful comparison. This is especially important when the genomes compared are very similar to each other.
Different methods for lengthnormalizing digital signals were tested: downsampling [55], upsampling to the maximum length using zero padding [30], even scaling extension [56], periodic extension, symmetric padding, or antisymmetric padding [57]. For example, zeropadding, which adds zeroes to all of the sequences shorter than the maximum length, was used in [30], e.g., for taxonomic classifications of ribosomal S18 subunit genes from twelve organisms. While this method may work for datasets of sequences of similar lengths, it is not suitable for datasets of sequences of very different lengths (our study: fungi mtDNA genomes dataset  1364 bp to 235,849 bp; plant mtDNA genomes dataset  12,998 bp to 1,999,595 bp; protist mtDNA genomes dataset  5882 bp to 77,356 bp). In such cases, zeropadding acts as a tag and may lead to inadvertent classification of sequences based on their length rather than based on their sequence composition. Thus, we employed instead antisymmetric padding, whereby, starting from the last position of the signal, boundary values are replicated in an antisymmetric manner. We also considered two possible ways of employing antisymmetric padding: normalization to the maximum length (where shorter sequences are extended to the maximum sequence length by antisymmetric padding) vs. normalization to the median length (where shorter sequences are extended by antisymmetric padding to the median length, while longer sequences are truncated after the median length).
Supervised machine learning
In this paper we used the Linear discriminant, Linear SVM, Quadratic SVM, Fine KNN, Subspace discriminant and Subspace KNN classifiers from the Classification Learner application of MATLAB (Statistics and Machine Learning Toolbox). The default MATLAB parameters were used.
To assess the performance of the classifiers, we used 10fold cross validation. In this approach, the dataset is randomly partitioned into 10 equalsize subsets. The classifier is trained using 9 of the subsets, and the accuracy of its prediction is tested on the remaining subset. As part of the supervised learning, taxonomic labels are supplied for the DNA sequences in the 9 subsets used for training. The process is repeated 10 times, and the accuracy score of the classifier is then computed as the average of the accuracies obtained in the 10 separate experiments. The standard algorithms were modified so that no information about sequences in the testing set (that is, no distance matrix entries containing distances to/from any sequence in the testing set to any other sequence) was available during the training stage.
Classical multidimensional scaling (MDS)
Classical multidimensional scaling takes a pairwise distance matrix (n×n matrix, for n input items) as input, and produces n points in a qdimensional Euclidean space, where q≤n−1. More specifically, the output is an n×q coordinate matrix, where each row corresponds to one of the n input items, and that row contains the q coordinates of the corresponding itemrepresenting point [11]. The Euclidean distance between each pair of points is meant to approximate the distance between the corresponding two items in the original distance matrix.
These points can then be simultaneously visualized in a 2 or 3dimensional space by taking the first 2, respectively 3, coordinates (out of q) of the coordinate matrix. The result is a Molecular Distance Map [46], and the MoDMap of a genomic dataset represents a visualization of the simultaneous interrelationships among all DNA sequences in the dataset.
Software implementation
The algorithms for MLDSP were implemented using the software package MATLAB R2017A, license no. 964054, as well as the opensource toolbox Fathom Toolbox for MATLAB [58] for distance computation. All software can be downloaded from https://github.com/grandhawa/MLDSP. The user can use this code to reproduce all results in this paper, and also has the option to input their own dataset and use it as training set for the purpose of classifying new genomic DNA sequences.
All experiments were performed on an ASUS ROG G752VS computer with 4 cores (8 threads) of a 2.7GHz Intel Core i7 6820HK processor and 64GB DD4 2400MHz SDRAM.
Datasets
All datasets in this paper can be found at https://github.com/grandhawa/MLDSP in the “DataBase” directory. The mitochondrial dataset comprises all of the 7396 complete reference mtDNA sequences available in the NCBI Reference Sequence Database RefSeq on June 17, 2017. We performed computational experiments on several different subsets of this dataset. The bacteria dataset comprises all 4710 complete bacterial genomes with lengths between 20,000 bp and 500,000 bp, available in the aforementioned NCBI database on the same date. The dengue virus dataset contained all 4721 dengue virus genomes available in the NCBI database on August 10,2017. Note that any letters “N” in these DNA sequences were deleted.
For the performance comparison between MLDSP and other alignmentfree and alignmentbased methods we also used the benchmark datasets of 38 influenza virus sequences, and 41 mammalian complete mtDNA sequences from [47].
Results and discussion
Following the design and implementation of the MLDSP genomic sequence classification tool prototype, we investigated which type of lengthnormalization and which type of distance were most suitable for genome classification using this method. We then conducted a comprehensive analysis of the various numerical representations of DNA sequences used in the literature, and determined the top three performers. Having set the main parameters (lengthnormalization method, distance, and numerical representation), we tested MLDSP’s ability to classify mtDNA genomes at taxonomic levels ranging from the domain level down to the genus level, and obtained average levels of classification accuracy of >97%. Finally, we compared MLDSP with other alignmentbased and alignmentfree genome classification methods, and showed that MLDSP achieved higher accuracy and significantly higher speeds.
Analysis of distances and of length normalization approaches
To decide which distance measure and which length normalization method were most suitable for genome comparisons with MLDSP, we used nine different subsets of full mtDNA sequences from our dataset. These subsets were selected to include most of the available complete mtDNA genomes (Vertebrates dataset of 4322 mtDNA sequences), as well as subsets containing similar sequences, of similar length (Primates dataset of 148 mtDNA sequences), and subsets containing mtDNA genomes showing large differences in length (Plants dataset of 174 mtDNA sequences).
The classification accuracy scores obtained using the two considered distance measures (Euclidean and Pearson Correlation Coefficient) and two different lengthnormalization approaches (normalization to maximum length and normalization to median length) on several datasets are listed in Table 2. The classification accuracy scores are slightly higher for PCC, but sufficiently close to those obtained when using the Euclidean distance to be inconclusive.
In the remainder of this paper we chose the Pearson Correlation Coefficient because it is scale independent (unlike the Euclidean distance, which is, e.g., sensitive to the offset of the signal, whereby signals with the same shape but different starting points are regarded as dissimilar [59]), and the lengthnormalization to median length because it is economic in terms of memory usage.
Analysis of various numerical representations of DNA sequences
We analyzed the effect on the MLDSP classification accuracy of thirteen different onedimensional numeric representations for DNA sequences, grouped as: Fixed mappings DNA numerical representations (Table 1 representations #1, #2, #3, #6, #7, see [29], and representations #10, #11, #12, #13  which are onedimensional variants of the binary representation proposed in [29]), mappings based on some physiochemical properties of nucleotides (Table 1 representation #4, see [29, 32], and representation #5, see [29, 31, 32]), and mappings based on the nearestneighbour values (Table 2 representations #8, #9, see [30]).
The datasets used for this analysis were the same as those in Table 2. The supervised machine learning classifiers used for this analysis were the six classifiers listed in the Methods and Implementation section, with the exception of the datasets with more than 2000 sequences where two of the classifiers (Subspace Discriminant and Subspace KNN) were omitted as being too slow. The results and the average accuracy scores for all these numerical representations, classifiers and datasets are summarized in Table 3.
As can be observed from Table 3, for all numerical representations, the table average accuracy scores (last row: average of averages, first over the six classifiers for each dataset, and then over all datasets), are high. Surprisingly, even using a single nucleotide numerical representation, which treats three of the nucleotides as being the same, and singles out only one of them (“JustA”), results in an average accuracy of 91.9%. The best accuracy, for these datasets, is achieved when using the “PP” representation, which yields an average accuracy of 92.3%.
For subsequent experiments we selected the top three representations in terms of accuracy scores: “PP”, “JustA”, and “Real” numerical representations.
MLDSP for three classes of vertebrates
As an application of MLDSP using the “PP” numerical representation for DNA sequences, we analyzed the set of vertebrate mtDNA genomes (median length 16,606 bp). The MoDMap, i.e., the multidimensional scaling 3D visualization of the genome interrelationships as described by the distances in the distance matrix, is illustrated in Fig. 3. The dataset contains 3740 complete mtDNA genomes: 553 bird genomes, 2313 fish genomes, and 874 mammalian genomes. Quantitatively, the classification accuracy score obtained by the Quadratic SVM classifier was 100%.
Classifying genomes with MLDSP, at all taxonomic levels
We tested the ability of MLDSP to classify complete mtDNA sequences at various taxonomic levels. For every dataset, we tested using the “PP”, “JustA”, and “Real” numerical representations.
The starting point was domain Eukaryota (7396 sequences), which was classified into kingdoms, then kingdom Animalia was classified into phyla, etc. At each level, we picked the cluster with the highest number of sequences and then classified it into the next taxonomic level subclusters. The lowest level classified was family Cyprinidae (81 sequences) into its six genera. For each dataset, we tested all six classifiers, and the maximum of these six classification accuracy scores for each dataset are shown in Table 4.
Note that, at each taxonomic level, the maximum classification accuracy scores (among the six classifiers) for each of the three numerical representations considered are high, ranging from 91.4% to 100%, with only three scores under 95%. As this analysis also did not reveal a clear winner among the top three numerical representations, the question then arose whether the numerical representation we use mattered at all. To answer this question, we performed two additional experiments, that exploit the fact that the Pearson correlation coefficient is scale independent, and only looks for a pattern while comparing signals. For the first experiment we selected the top three numerical representations (“PP”, “JustA”, and “Real”) and, for each sequence in a given dataset, a numerical representation among these three was randomly chosen, with equal probability, to be the digital signal that represents it. The results are shown under the column “Random3” in Table 4: The maximum accuracy score over all the datasets is 96%. This is almost the same as the accuracy obtained when one particular numerical representation was used (1% lower, which is well within experimental error). We then repeated this experiment, this time picking randomly from any of the thirteen numerical representations considered. The results are shown under the column “Random13” in Table 4, with the table average accuracy score being 88.1%.
Overall, our results suggest that all three numerical representations “PP”, “JustA”, and “Real” have very high classifications accuracy scores (average >97%), and even a random choice of one of these representations for each sequence in the dataset does not significantly affect the classification accuracy score of MLDSP (average 96%).
We also note that, in addition to being highly accurate in its classifications, MLDSP is ultrafast. Indeed, even for the largest dataset in Table 2, subphylum Vertebrata (4322 complete mtDNA genomes, average length 16,806 bp), the distance matrix computation (which is the bulk of the classification computation) lasted under 5 s. Classifying a new primate mtDNA genome took 0.06 s when trained on 148 primate mtDNA genomes, and classifying a new vertebrate mtDNA genome took 7 s when trained on the 4322 vertebrate mtDNA genomes. The result was updated with an experiment whereby QSVM was trained on the 4322 complete vertebrate genomes in Table 2, and querried on the 694 new vertebrate mtDNA genomes uploaded on NCBI between June 17, 2017 and January 7, 2019. The accuracy of classification was 99.6%, with only three reptile mtDNA genomes misclassified as amphibian genomes: Bavayia robusta, robust forest bavayia  a species of gecko, NC_034780, Mesoclemmys hogei, Hoge’s toadhead turtle, NC_036346, and Gonatodes albogularis, yellowheaded gecko, NC_035153.
MoDMap visualization vs. MLDSP quantitative classification results
The hypothesis tested by the next experiments was that the quantitative accuracy of the classification of DNA sequences by MLDSP would be significantly higher than suggested by the visual clustering of taxa in the MoDMap produced with the same pairwise distance matrix.
As an example, the MoDMap in Fig. 4a, visualizes the distance matrix of mtDNA genomes from family Cyprinidae (81 genomes) with its genera Acheilognathus (10 genomes), Rhodeus (11 genomes), Schizothorax (19 genomes), Labeo (19 genomes), Acrossocheilus (12 genomes), Onychostoma (10 genomes); only the genera with at least 10 genomes are considered. The MoDMap seems to indicate an overlap between the clusters Acheilognathus and Rhodeus, which is biologically plausible as these genera belong to the same subfamily Acheilognathinae. However, when zooming in by plotting a MoDMap of only these two genera, as shown in Fig. 4b, one can see that the clusters are clearly separated visually. This separation is confirmed by the fact that the accuracy score of the Quadratic SVM classifier for the dataset in Fig. 4b is 100%. The same quantitative accuracy score for the classification of the dataset in Fig. 4a with Quadratic SVM is 91.8%, which intuitively is much better than the corresponding MoDMap would suggest. This is likely due to the fact that the MoDMap is a threedimensional approximation of the positions of the genomerepresenting points in a multidimensional space (the number of dimensions is (n−1), where n is the number of sequences).
This being said, MoDMaps can still serve for exploratory purposes. For example, the MoDMap in Fig. 4a suggests that species of the genus Onychostoma (subfamily listed “unknown” in NCBI) (yellow), may be genetically related to species of the genus Acrossocheilus (subfamily Barbinae) (magenta). Upon further exploration of the distance matrix, one finds that indeed the distance between the centroids of these two clusters is lower than the distance between each of these two clustercentroids to the other clustercentroids. This supports the hypotheses, based on morphological evidence [60], that genus Onychostoma belongs to the subfamily Barbinae, respectively that genus Onychostoma and genus Acrossocheilus are closely related [61]. Note that this exploration, suggested by MoDMap and confirmed by calculations based on the distance matrix, could not have been initiated based on MLDSP alone (or other supervised machine learning algorithms), as MLDSP only predicts the classification of new genomes into one of the taxa that it was trained on, and does not provide any other additional information.
As another comparison point between MoDMaps and supervised machine learning outputs, Fig. 5a shows the MoDMap of the superorder Ostariophysi with its orders Cypriniformes (643 genomes), Characiformes (31 genomes) and Siluriformes (107 genomes). The MoDMap shows the clusters as overlapping, but the Quadratic SVM classifier that quantitatively classifies these genomes has an accuracy of 99%. Indeed, the confusion matrix in Fig. 5b shows that Quadratic SVM misclassifies only 8 sequences out of 781 (recall that, for m clusters, the m×m confusion matrix has its rows labelled by the true classes and columns labelled by the predicted classes; the cell (i,j) shows the number of sequences that belong to the true class i, and have been predicted to be of class j). This indicates that when the visual representation in a MoDMap shows cluster overlaps, this may only be due to the dimensionality reduction to three dimensions, while MLDSP actually provides a much better quantitative classification based on the same distance matrix.
Applications to other genomic datasets
The two experiments in this section indicate that the applicability of our method is not limited to mitochondrial DNA sequences. The first experiment, Fig. 6a, shows the MoDMap of all 4721 complete dengue virus sequences available in NCBI on August 10,2017, classified into the subtypes DENV1 (2008 genomes), DENV2 (1349 genomes), DENV3 (1010 genomes), DENV4 (354 genomes). The average length of these complete viral genomes is 10,595 bp. Despite the dengue viral genomes being very similar, the classification accuracy of this dataset into subtypes, using the Quadratic SVM classifier, was 100%. The second experiment, Fig. 6b, shows the MoDMap of 4710 bacterial genomes, classified into three phyla: Spirochaetes (437 genomes), Firmicutes (1129 genomes), and Proteobacteria (3144 genomes). The average length of these complete bacterial genomes is 104,150 bp, with the maximum length being 499,136 bp and the minimum length being 20,019 bp. The classification accuracy of the Quadratic SVM classifier for this dataset was 95.5%.
Comparison of MLDSP with stateoftheart alignmentbased and alignmentfree tools
The computational experiments in this section compare MLDSP with three stateoftheart alignmentbased and alignmentfree methods: the alignmentbased tool MEGA7 [3] with alignment using MUSCLE [4] and CLUSTALW [5, 6], and the alignmentfree method FFP (Feature Frequency Profiles) [28].
For this performance analysis we selected three datasets. The first two datasets are benchmark datasets used in other genetic sequence comparison studies [47]: The first dataset comprises 38 influenza viral genomes, and the second dataset comprises 41 mammalian complete mtDNA sequences. The third dataset, of our choice, is much larger, consisting of 4,322 vertebrate complete mtDNA sequences, and was selected to compare scalability.
For the alignmentbased methods, we used the distance matrix calculated in MEGA7 from sequences aligned with either MUSCLE or CLUSTALW. For the alignmentfree FFP, we used the default value of k=5 for kmers (a kmer is any DNA sequence of length k; any increase in the value of the parameter k, for the first dataset, resulted in a lower classification accuracy score for FFP). For MLDSP we chose the Integer numerical representation and computed the average classification accuracy over all six classifiers for the first two datasets, and over all classifiers except Subspace Discriminant and Subspace KNN for the third dataset.
Table 5 shows the performance comparison (classification accuracy and processing time) of these four methods. The processing time included all computations, starting from reading the datasets to the completion of the distance matrix  the common element of all four methods. The listed processing times do not include the time needed for the computation of phylogenetic trees, MoDMap visualizations, or classification.
As seen in Table 5 (columns 3, 4, and 6) MLDSP overwhelmingly outperforms the alignmentbased software MEGA7(MUSCLE/CLUSTALW) in terms of processing time. In terms of accuracy, for the smaller virus and mammalian benchmark datasets, the average accuracies of MLDSP and MEGA7(MUSCLE/CLUSTALW) were comparable, probably due to the small size of the training set for MLDSP. The advantage of MLDSP over the alignmentbased tools became more apparent for the larger vertebrate dataset, where the accuracies of MLDSP and the alignmentbased tools could not even be compared, as the alignmentbased tools were so slow that they had to be terminated. In contrast, MLDSP classified the entire set of 4322 vertebrate mtDNA genomes in 28 s, with average classification accuracy 98.3%. This indicates that MLDSP is significantly more scalable than the alignmentbased MEGA7(MUSCLE/CLUSTALW), as it can speedily and accurately classify datasets which alignmentbased tools cannot even process.
As seen in Table 5 (columns 5 and 6), MLDSP significantly outperforms the alignmentfree software FFP in terms of accuracy (average classification accuracy 98.3% for MLDSP vs. 48.3% for FFP, for the large vertebrate dataset), while at the same time being overall faster.
This comparison also indicates that, for these datasets, both alignmentfree methods (MLDSP and FFP) have an overwhelming advantage over the alignmentbased methods (MEGA7 (MUSCLE/CLUSTALW)) in terms of processing time. Furthermore, when comparing the two alignmentfree methods with each other, MLDSP significantly outperforms FFP in terms of classification accuracy.
As another angle of comparison, Fig. 7 displays the MoDMaps of the first benchmark dataset (38 influenza virus genomes) produced from the distance matrices generated by FFP, MEGA7 (MUSCLE), MEGA7 (CLUSTALW), and MLDSP respectively. Figure 7a shows that with FFP it is difficult to observe any visual separation of the dataset into subtype clusters. Figure 7b, MEGA7 (MUSCLE), and Fig. 7c MEGA7 (CLUSTALW) show overlaps of the clusters of points representing subtypes H1N1 and H2N2. In contrast, Fig. 7d, which visualizes the distance matrix produced by MLDSP, shows a clear separation among all subtypes.
Finally Figs. 8 and 9 display the phylogenetic trees generated by each of the four methods considered. Figure 8a, the tree generated by FFP, has many misclassified genomes, which was expected given the MoDMap visualization of its distance matrix in Fig. 7a. Figure 9a displays the phylogenetic tree generated by MEGA7, which was the same for both MUSCLE and CLUSTALW: It has only one incorrectly classified H5N1 genome, placed in middle of H1N1 genomes. Figures 8b and 9b display the phylogenetic tree generated using the distance produced by MLDSP (shown twice, in parallel with the other trees, for ease of comparison). MLDSP classified all genomes correctly.
Discussion
The computational efficiency of MLDSP is due to the fact that it is alignmentfree (hence it does not need multiple sequence alignment), while the combination of 1D numerical representations, Discrete Fourier Transform and Pearson Correlation Coefficient makes it extremely computationally time efficient, and thus scalable.
MLDSP is not without limitations. We anticipate that the need for equal length sequences and use of length normalization could introduce issues with examination of small fragments of larger genome sequences. Usually genomes vary in length and thus length normalization always results in adding (upsampling) or losing (downsampling) some information. Although the Pearson Correlation Coefficient can distinguish the signal patterns even in small sequence fragments, and we did not find any considerable disadvantage while considering complete mitochondrial DNA genomes with their inevitable length variations, length normalization may cause issues when we deal with the fragments of genomes, and the much larger nuclear genome sequences.
Lastly, MLDSP has two drawbacks, inherent in any supervised machine learning algorithm. The first is that MLDSP is a blackbox method which, while producing a highly accurate classification prediction, does not offer a (biological) explanation for its output. The second is that it relies on the existence of a training set from which it draws its “knowledge”, that is, a set consisting of known genomic sequences and their taxonomic labels. MLDSP uses such a training set to “learn” how to classify new sequences into one of the taxonomic classes that it was trained on, but it is not able to assign it to a taxon that it has not been exposed to.
Conclusions
We proposed MLDSP, an ultrafast and accurate alignmentfree supervised machine learning classification method based on digital signal processing of DNA sequences (and its software implementation). MLDSP successfully addresses the limitations of alignmentfree methods identified in [7], as follows:

(i)
Lack of software implementation: MLDSP is supplemented with freely available sourcecode. The MLDSP software can be used with the provided datasets or any other custom dataset and provides the user with any (or all) of: pairwise distances, 3D sequence interrelationship visualization, phylogenetic trees, or classification accuracy scores. A quantitative comparison showed that MLDSP significantly outperforms stateoftheart alignmentbased MEGA7 (MUSCLE/CLUSTALW) and alignmentfree (FFP) software in terms of speed and classification accuracy.

(ii)
Use of simulated sequences or very small realworld datasets: MLDSP was successfully tested on a variety of large realworld datasets, comprising thousands of complete genomes, such as all complete mitochondrial DNA sequences available on NCBI at the time of this study, and similarly large sets of viral genomes and bacterial genomes. MLDSP was tested in different evolutionary scenarios such as different levels of taxonomy (from domain to genus), small dataset (38 sequences), large dataset (4322 sequences), short sequences (1,136 bp), long sequences (1,999,595 bp), benchmark datasets of influenza virus and mammalian mtDNA genomes etc.

(iii)
Memory overhead: MLDSP uses neither kmers nor any compression algorithms. Thus, scalability does not cause an exponential memory overhead, and a high classification accuracy is preserved with large datasets.
In addition, we provided a comprehensive quantitative analysis of all 13 onedimensional numerical representations of DNA sequences used in the Genomic Signal Processing literature and found that, on average, the “PP”, “JustA”, and “Real” representations performed better than others. We also showed that the classification accuracy of MLDSP was significantly higher than the corresponding MoDMap visualizations of the dataset would indicate, likely due to the inherent dimensionality limitations of the latter. Lastly, we showed the potential for MLDSP to be used for classifications of other DNA sequence genomic datasets, such as large datasets of complete viral or bacterial genomes.
Availability and Requirements
Project name: MLDSP
Project home page: https://github.com/grandhawa/MLDSP
Operating system(s): Microsoft Windows
Programming language: MATLAB R2017A, license no. 964054
License: Creative Commons Attribution License
Any restrictions to use by nonacademics: MATLAB license required
Abbreviations
 DFT:

Discrete Fourier transform
 DSP:

Digital signal processing
 FFP:

Feature frequency profile
 GSP:

Genomic signal processing
 MDS:

Classical multidimensional scaling
 MoDMap:

Molecular distance map
 PCC:

Pearson correlation coefficient
 PP:

Purine/Pyrimidine
References
 1
Mora C, Tittensor DP, Adl S, Simpson AGB, Worm B. How many species are there on earth and in the ocean?PLoS Biol. 2011; 9(8):1001127.
 2
May RM. Why worry about how many species and their loss?PLoS Biol. 2011; 9(8):1001130.
 3
Kumar S, Stecher G, Tamura K. MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets. Mol Biol Evol. 2016; 33(7):1870–4.
 4
Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004; 32(5):1792–7.
 5
Thompson JD, Higgins DG, Gibson TJ. CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, positionspecific gap penalties and weight matrix choice. Nucleic Acids Res. 1994; 22(22):4673–80.
 6
Larkin MA, Blackshields G, Brown NP, Chenna R, McGettigan PA, McWilliam H, Valentin F, Wallace IM, Wilm A, Lopez R, Thompson JD, Gibson TJ, Higgins DG. CLUSTAL W and CLUSTAL X version 2.0. Bioinformatics. 2007; 23(21):2947–8.
 7
Zielezinski A, Vinga S, Almeida J, Karlowski WM. Alignmentfree sequence comparison: benefits, applications, and tools. Genome Biol. 2017; 18(1):186.
 8
Vinga S, Almeida J. Alignmentfree sequence comparison—a review. Bioinformatics. 2003; 19(4):513–23.
 9
Schwende I, Pham TD. Pattern recognition and probabilistic measures in alignmentfree sequence analysis. Brief Bioinform. 2014; 15(3):354–68.
 10
Song K, Ren J, Reinert G, Deng M, Waterman MS, Sun F. New developments of alignmentfree sequence comparison: measures, statistics and nextgeneration sequencing. Brief Bioinform. 2014; 15(3):343–53.
 11
Kari L, Hill KA, Sayem AS, Karamichalis R, Bryans N, Davis K, Dattani NS. Mapping the space of genomic signatures. PLoS ONE. 2015; 10(5):0119815.
 12
Hoang T, Yin C, Yau SS. Numerical encoding of DNA sequences by chaos game representation with application in similarity comparison. Genomics. 2016; 108(3):134–42.
 13
Almeida J, Carriço JA, Maretzek A, Noble PA, M F. Analysis of genomic sequences by chaos game representation. Bioinformatics. 2001; 17 5:429–37.
 14
Yao YH, Dai Q, Nan XY, He PA, Nie ZM, Zhou SP, Zhang YZ. Analysis of similarity/dissimilarity of DNA sequences based on a class of 2D graphical representation. J Comput Chem. 2008; 29(10):1632–9.
 15
Qi X, Wu Q, Zhang Y, Fuller E, Zhang CQ. A novel model for DNA sequence similarity analysis based on graph theory. Evol Bioinformatics Online. 2011; 7:149–58.
 16
Almeida JS. Sequence analysis by iterated maps, a review. Brief Bioinform. 2014; 15(3):369–75.
 17
Vinga S. Information theory applications for biological sequence analysis. Brief Bioinform. 2014; 15(3):376–89.
 18
Bao J, Yuan R, Bao Z. An improved alignmentfree model for DNA sequence similarity metric. BMC Bioinformatics. 2014; 15(1):321.
 19
Leimeister CA, Boden M, Horwege S, Lindner S, Morgenstern B. Fast alignmentfree sequence comparison using spacedword frequencies. Bioinformatics. 2014; 30(14):1991–9.
 20
Chang G, Wang H, Zhang T. A novel alignmentfree method for whole genome analysis: Application to HIV1 subtyping and hev genotyping. Inf Sci. 2014; 279:776–84.
 21
Reese E, Krishnan VV. Classification of DNA sequences based on thermal melting profiles. Bioinformation. 2010; 4(10):463–7.
 22
BonhamCarter O, Steele J, Bastola D. Alignmentfree genetic sequence comparisons: a review of recent approaches by word analysis. Brief Bioinform. 2014; 15(6):890–905.
 23
Struck D, Lawyer G, Ternes AM, Schmit JC, Bercoff DP. Comet: adaptive contextbased modeling for ultrafast HIV1 subtype identification. Nucleic Acids Res. 2014; 42(18):144.
 24
Remita MA, Halioui A, Malick Diouara AA, Daigle B, Kiani G, Diallo AB. A machine learning approach for viral genome classification. BMC Bioinformatics. 2017; 18:208.
 25
Kosakovsky Pond SL, Posada D, Stawiski E, Chappey C, Poon AF, Hughes G, Fearnhill E, Gravenor MB, Leigh Brown AJ, Frost SD. An evolutionary modelbased algorithm for accurate phylogenetic breakpoint mapping and subtype prediction in HIV1. PLoS Comput Biol. 2009; 5(11):1000581.
 26
de Oliveira T, Deforche K, Cassol S, Salminen M, Paraskevis D, Seebregts C, Snoeck J, van R EJ, Wensing AMJ, van de Vijver DA, Boucher CA, Camacho R, Vandamme AM. An automated genotyping system for analysis of HIV1 and other microbial sequences. Bioinformatics. 2005; 21(19):3797–800.
 27
SolisReyes S, Avino M, Poon A, Kari L. An opensource kmer based machine learning tool for fast and accurate subtyping of HIV1 genomes. PLoS ONE. 2018; 13(11):0206409.
 28
Sims GE, Jun SR, Wu GA, Kim SH. Alignmentfree genome comparison with Feature Frequency Profiles (FFP) and optimal resolutions. In: Proceedings of the National Academy of Sciences of the USA. USA: National Academy of Sciences: 2009. p. 2677–82. https://doi.org/10.1073/pnas.0813249106.
 29
Kwan HK, Arniker SB. Numerical representation of DNA sequences. In: 2009 IEEE International Conference on Electro/Information Technology. New Jersey: IEEE publishing: 2009. p. 307–10. https://doi.org/10.1109/EIT.2009.5189632.
 30
Borrayo E, MendizabalRuiz EG, VélezPérez H, RomoVázquez R, Mendizabal AP, Morales JA. Genomic signal processing methods for computation of alignmentfree distances from DNA sequences. PLoS ONE. 2014; 9(11):110954.
 31
Adetiba E, Olugbara OO, Taiwo TB. Identification of pathogenic viruses using genomic cepstral coefficients with radial basis function neural network. In: Advances in Nature and Biologically Inspired Computing, Proceedings of the 7th World Congress on Nature and Biologically Inspired Computing: 2016. p. 281–90.
 32
Adetiba E, Olugbara OO. Classification of eukaryotic organisms through cepstral analysis of mitochondrial DNA. In: International Conference on Image and Signal Processing. Berlin: Springer: 2016. p. 243–52. https://doi.org/10.1007/9783319336183_25.
 33
MendizabalRuiz G, RománGodínez I, TorresRamos S, SalidoRuiz RA, Morales JA. On DNA numerical representations for genomic similarity computation. PLoS ONE. 2017; 12(3):0173288.
 34
Chakravarthy N, Spanias A, Iasemidis LD, Tsakalis K. Autoregressive modeling and feature analysis of DNA sequences. EURASIP J Appl Signal Process. 2004; 2004:13–28.
 35
Yu Z, Anh VV, Zhou Y, Zhou LQ. Numerical sequence representation of DNA sequences and methods to distinguish coding and noncoding sequences in a complete genome. In: Proceedings 11th World MultiConference on Systemics, Cybernetics and Informatics. Orlando: International Institute of Informatics and Systemics: 2007. p. 171–6.
 36
AboZahhad M, Ahmed S, AbdElrahman S. Genomic analysis and classification of exon and intron sequences using DNA numerical mapping techniques. Int J Inform Technol Comput Sci. 2012; 4(8):22–36.
 37
Skutkova H, Vitek M, Sedlar K, Provaznik I. Progressive alignment of genomic signals by multiple dynamic time warping. J Theor Biol. 2015; 385:20–30.
 38
Yin C, Yau SST. An improved model for whole genome phylogenetic analysis by Fourier transform. J Theor Biol. 2015; 382:99–110.
 39
LorenzoGinori JV, RodriguezFuentes A, Grau Abalo R, Sanchez Rodriguez R. Digital signal processing in the analysis of genomic sequences. Curr Bioinforma. 2009; 4(1):28–40.
 40
Weitschek E, Cunial F, Felici G. LAF: Logic alignment free and its application to bacterial genomes classification. BioData Mining. 2015; 8:39.
 41
Fiscon G, Weitschek E, Cella E, Lo Presti A, Giovanetti M, BabakirMina M, Ciotti M, Ciccozzi M, Pierangeli A, Bertolazzi P, Felici G. MISSEL: a method to identify a large number of small speciesspecific genomic subsequences and its application to viruses classification. BioData Mining. 2016; 9:38.
 42
Remita MA, Halioui A, Malick Diouara AA, Daigle B, Kiani G, Diallo AB. A machine learning approach for viral genome classification. BMC Bioinformatics. 2017; 18:208.
 43
Lu H, Yang L, Yan K, Xue Y, Gao Z. A costsensitive rotation forest algorithm for gene expression data classification. Neurocomputing. 2017; 228:270–6.
 44
Lu H, Meng Y, Yan K, Gao Z. Kernel principal component analysis combining rotation forest method for linearly inseparable data. Cogn Syst Res. 2018; 53:111–22.
 45
Liu Y, Lu H, Yan K, Xia H, An C. Applying costsensitive extreme learning machine and dissimilarity integration to gene expression data classification. Comput Intell Neurosci. 2016; 2016:1–9.
 46
Karamichalis R, Kari L. MoDMaps3D: an interactive webtool for the quantification and 3D visualization of interrelationships in a dataset of DNA sequences. Bioinformatics. 2017; 33(19):3091–3.
 47
Li Y, He L, Lucy He R, Yau SST. A novel fast vector method for genetic sequence comparison. Sci Rep. 2017; 7(1):1–11.
 48
Cristea PD. Conversion of nucleotide sequences into genomic signals. J Cell Mol Med. 2002; 6(2):279–303.
 49
Afreixo V, Bastos CAC, Pinho AJ, Garcia SP, Ferreira PJSG. Genome analysis with distance to the nearest dissimilar nucleotide. J Theor Biol. 2011; 275(1):52–8.
 50
Cristea PD. Large scale features in DNA genomic signals. Signal Process. 2003; 83(4):871–88.
 51
Skutkova H, Vitek M, Babula P, Kizek R, Provaznik I. Classification of genomic signals using dynamic time warping. BMC Bioinformatics. 2013; 14(10):1.
 52
Asuero AG, Sayago A, González AG. The correlation coefficient: an overview. Crit Rev Anal Chem. 2006; 36(1):41–59.
 53
ElBadawy IM, Aziz AM, Omar Z, Malarvili MB. Correlation between different DNA period3 signals: An analytical study for exons prediction. In: 2017 AsiaPacific Signal and Information Processing Association Annual Summit and Conference. New Jersey: IEEE publishing: 2017. p. 1123–8. https://doi.org/10.1109/APSIPA.2017.8282195.
 54
Hoang T, Yin C, Zheng H, Yu C, He RL, Yau SST. A new method to cluster DNA sequences using Fourier power spectrum. J Theor Biol. 2015; 372:135–45.
 55
Sedlar K, Skutkova H, Vitek M, Provaznik I. Set of rules for genomic signal downsampling. Comput Biol Med. 2016; 69:308–14.
 56
Yin C, Chen Y, Yau SST. A measure of DNA sequence similarity by Fourier transform with applications on hierarchical clustering. J Theor Biol. 2014; 359:18–28.
 57
Strang G, Nguyen T. Wavelets and Filter Banks. Wellesley: WellesleyCambridge Press; 1996.
 58
Jones DL. Fathom Toolbox for MATLAB: software for multivariate ecological and oceanographic data analysis. St. Petersburg: College of Marine Science, University of South Florida; 2017. Available from: https://www.marine.usf.edu/research/matlabresources/.
 59
Lee S, Kwon D, Lee S. Efficient similarity search for time series data based on the minimum distance. In: International Conference on Advanced Information Systems Engineering. Berlin: Springer: 2002. p. 377–91. https://doi.org/10.1007/3540479619_27.
 60
Taki Y. Cyprinid fishes of the genera Onychostoma and Scaphiodonichthys from Upper Laos with remarks on the dispersal of the genera and their allies. Jpn J Ichthyol. 1975; 22(3):143–50.
 61
Zheng L, Yang J, Chen X. Molecular phylogeny and systematics of the Barbinae (Teleostei: Cyprinidae) in China inferred from mitochondrial DNA sequences. Biochem Syst Ecol. 2016; 68:250–9.
Acknowledgements
We thank Michael Pang for an independent Python implementation and reproducing some of the computational results, and Maximillian Soltysiak and Nicholas A. Boehler for comments on the manuscript and for testing the software tool.
Funding
This work was supported by NSERC (Natural Science and Engineering Research Council of Canada) Discovery Grants R2824A01 to L.K., and R3511A12 to K.A.H.
Availability of data and materials
The source code can be downloaded from https://github.com/grandhawa/MLDSP. All datasets can be found at https://github.com/grandhawa/MLDSP/DataBase.
Author information
Affiliations
Contributions
G.S.R. and L.K. conceived the study and wrote the manuscript. G.S.R. designed and tested the software. G.S.R., L.K. and K.A.H. conducted the data analysis and edited the manuscript, with K.A.H. contributing biological expertize. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Randhawa, G., Hill, K. & Kari, L. MLDSP: Machine Learning with Digital Signal Processing for ultrafast, accurate, and scalable genome classification at all taxonomic levels. BMC Genomics 20, 267 (2019). https://doi.org/10.1186/s128640195571y
Received:
Accepted:
Published:
Keywords
 Taxonomic classification
 Whole genome analysis
 Genomic signature
 Alignmentfree sequence analysis
 Machine learning
 Numerical representation of DNA sequences
 Digital signal processing
 Discrete Fourier transform