Alignment-Free Sequence Analysis and Applications
Jie Ren Xin Bai Yang Young Lu Kujin Tang Ying Wang Gesine Reinert Fengzhu Sun











Domain transfer convolutional attribute embedding
    
Fang Sua^∗
^∗Corresponding author.,
Jing-Yan Wangb

a
School of Economics and Management, Shaanxi University of Science & Technology, Xi'an, ShaanXi Province, P.R.C, 710021

b
New York University Abu Dhabi, Abu Dhabi, United Arab Emirates
v5.0 released July 2015

    .
============================================================================================================================================================================================================================================================================










Genome and metagenome comparisons based on large amounts of next generation sequencing (NGS) data pose significant challenges for alignment-based approaches due to the huge data size and the relatively short length of the reads. Alignment-free approaches based on the counts of word patterns in NGS data do not depend on the complete genome and are generally computationally efficient.  Thus, they contribute significantly to genome and metagenome comparison. Recently, novel statistical approaches have been developed for the comparison of both long and shotgun sequences. These approaches have been applied to many problems including the comparison of gene regulatory regions, genome sequences,  metagenomes, binning contigs in metagenomic data, identification of virus-host interactions, and detection of horizontal gene transfers.  We provide an updated review of these applications and other related developments of word-count based approaches for alignment-free sequence analysis.

keywords: sequence comparison, alignment, alignment-free, phylogeny, metagenomics, virus-host interaction, Markov chain, horizontal gene transfer



















§ INTRODUCTION




Molecular sequence comparison is one of the most basic and fundamental problems in computational biology, and has been widely used to study the evolution of whole genome sequences and gene regulatory regions, gene function prediction, sequence assembly, and finding the relationships among microbial communities. The most widely used methods for molecular sequence comparison are alignment-based including the Smith-Waterman algorithm <cit.>, BLAST <cit.>, BLAT<cit.>, etc. Although alignment-based approaches are most accurate and powerful for sequence comparison when they are feasible, their applications are limited in some situations. First, for whole genome comparison, there are many duplications, translocations, large insertions/delections, and horizontal gene transfers in the genomes.
This situation makes it difficult to use alignment-based methods to investigate the relationship among whole genome sequences.
Second, in the current next generation sequencing (NGS) era, investigators can sequence the genomes using NGS efficiently and economically. However, some parts of the genomes may not be sequenced due to the stochastic distribution of the reads along the genomes and the difficulties of sequencing some parts of the genomes, especially when the coverage is relatively low. Even if we can assemble the reads into long contigs, these contigs may not share long homologous regions making it challenging to study the relationships among the genomes using alignment in such situations.
Third, noncoding regions such as gene regulatory regions are not highly conserved except for some functional regions such as transcription binding sites, and cannot be reliably aligned.
Therefore, alignment-based approaches are not well suited to study the evolution of gene regulatory regions.
Fourth, alignment is not suitable to compare sequences of large divergence.
When we investigate the relationship between viruses and their hosts, infecting virus-host pairs may only share a tiny fraction of their genomes such as CRISPR regions, and thus alignment-based approaches can potentially identify the hosts of only a small fraction of viruses.
Fifth, many large genome and metagenome data sets from shotgun NGS sequencing are available and alignment-based methods are too time consuming. For all these scenarios, alignment-free methods for genome and metagenome comparison provide promising alternative approaches.




Alignment-free approaches for sequence comparison can be divided into several different groups: a) word-counts <cit.>, b) average longest common substrings <cit.>, shortest unique substrings <cit.>, or a combination of both <cit.>, c) sequence representation based  on chaos theory <cit.>, d) the moments of the positions of the nucleotides <cit.>, e) Fourier transformation <cit.>, f) information theory  <cit.>, and g) iterated maps <cit.>. Several excellent reviews on various alignment-free sequence comparison methods have been published <cit.>


In this review, we concentrate on methods that can be applied to the comparison of sequences based on NGS data. Since the word-count-based approaches are the most adaptable to NGS reads data, we deal with word-count-based approaches as in <cit.>. These methods first count the number of occurrences of word patterns (k-mers, k-grams, k-tuples) along a sequence or in a NGS sample using different algorithms such as Jellyfish <cit.>, DSK <cit.>, and KMC 2 <cit.>. Secondly, a similarity/dissimilarity measure is defined between any pair of sequences based on the word-count frequencies. Finally, various clustering algorithms such as hierarchical clustering and neighbor-joining are used to group the sequences. In the rest of the review, we use “word" and “k-mer" interchangeably.

The use of k-mer frequencies to compare molecular sequences traces back to the early work of Carl Woese and colleagues from the early  1970s to the mid 1980s when they generated oligonucleotide catalogs of 16S rRNA sequences from about 400 organisms <cit.>.
They used a similarity measure, S_AB, for two sequences A and B  using k-mers similar to the Bray-Curtis dissimilarity <cit.>.  When the whole 16S rRNA sequences for many organisms were available, they showed a positive correlation between the dissimilarity of two sequences using k-mers with the distance calculated by alignment although the correlation is not very high (0.40) <cit.>. Ragan et al. <cit.> gave an excellent review of these early efforts to study the relationships among 16S sequences using oligonucleotide patterns and compared the dendrograms derived using multiple sequence alignment, the similarity measure S_AB, and the newly developed d_2^S statistic <cit.>. It was shown that the tree constructed based on d_2^S for k from 6 to 16 yielded the dendrogram that was most consistent with the maximum likelihood tree using multiple sequence alignment.

Many word-count-based methods for sequence comparison have been developed including the un-centered  correlation of word count vectors between two sequences <cit.>, χ^2-statistics <cit.>, composition vectors <cit.>, nucleotide relative abundances <cit.>, and the recently developed d_2^* and d_2^S statistics <cit.>.
It was shown that alignment-free methods are more robust than alignment-based methods especially against genetic rearrangements and horizontal gene transfers <cit.>. Since word frequencies are generally stable across different genomic regions, alignment-free methods work well even with sequences coming from different regions of the genomes.
Song et al. <cit.> presented an review of the  development and applications of these methods before 2013. In the current review, we provide further developments of d_2^* and d_2^S and their applications in recent years including a) how to determine the background Markov chain model of the sequences, b) genome, metagenome, and transcriptome comparison using Markov chains, c) inference of virus-bacterial host infectious associations, d) identification of horizontal gene transfers, and e) integrated software for alignment-free genome and metagenome comparison. We will also present an review of other developments related to d_2^* and d_2^S in recent years. For a recent review of other alignment-free sequence comparison methods and their applications, see <cit.>.













§ DETERMINATION OF THE BACKGROUND MARKOV CHAIN MODELS OF THE GENOMES



Alignment-free sequence comparison methods using k-mers generally involve counting the number of occurrences of words of length k in genomic sequences and comparing sequences using dissimilarity measures defined based on k-mer frequencies. Different dissimilarity measures have been developed using a number of principles.
The measures can be broadly classified into two groups: measures that require background word frequencies and those that do not. Lu et al. <cit.> developed a one-stop platform for computing a suite of 28 different alignment-free measures and provided various forms of visualization  tools including dendrograms, heatmaps, principal coordinate analysis and network display. The definitions of the 28 measures can be found in the supplementary material for <cit.>.


For measures that do not require background word frequencies, the observed word frequency or word presence (or absence) are directly used to compute the dissimilarity measures. The measures include but are not limited to,
Euclidian distance (Eu),
Manhattan distance (Ma),
d_2 <cit.>,
Feature Frequency Profiles (FFP) <cit.>,
Jensen-Shannon divergence (JS) <cit.>,
 Hamming distance, and Jaccard index. 
For measures that take background word frequency into account, dissimilarity between sequences is computed using the normalized word frequencies, where the expected word frequencies estimated using a background model are subtracted from the observed word frequencies to eliminate the background noise and enhance the signal.
This group of measures includes
d_2^*, d_2^S <cit.> and their variants <cit.>, CVTree <cit.>, Teeling <cit.>, EuF <cit.>
and Willner <cit.>, where different forms of sequence background models are incorporated.

The second group of measures requires the knowledge about the approximate distribution of word counts in the background sequences.
Markov chains (MC) are widely used to model genomic sequences <cit.> with many applications including
the study of dependencies between bases <cit.>,
the enrichment and depletion of certain word patterns <cit.>,
prediction of occurrences of long word patterns from short patterns
<cit.>,
and the detection of signals in introns <cit.>.
The defining feature of a MC model is the “memorylessness” property, that implies that the future state of the sequence can be well predicted solely based on its latest history without knowing the full history.
In particular, an r-th order MC assumes that the distribution of the future state only depends on the states of the past r positions regardless of the earlier history, i.e. P(X_t | X_1… X_t-1) = P(X_t | X_t-r… X_t-1),
where X_1, X_2,  … X_t, are the states in the sequence X, and X_i takes its values from a finite alphabet of size L.
For DNA sequences, the alphabet set is 𝒜={A, C, G, T}.
The MC can be represented in the form of a L^r × L matrix, where the element in the matrix corresponds to the transition probability P(w | w_1 w_2 … w_r), w ∈𝒜.
A 0-th order MC is the simplest case; in this case the positions in the sequence are independent and identically distributed (i.i.d.).
 




 §.§ INFERENCE OF MC PROPERTIES FOR A LONG GENOMIC SEQUENCE


For a long genomic sequence, efficient statistics are available to determine the order of the MC <cit.>.
For reviews on the application of MCs to molecular sequence analysis, see <cit.>.
In particular,
under the hypothesis that the long sequence follows a (k-2)-th order MC,
it holds that twice the log-likelihood ratio of the likelihood of the sequence
under a (k-1)-th order MC versus that under the (k-2)-th order
MC model follows approximately a χ^2-distribution with
df_k=(L-1)^2 L^k-2 degrees of freedom.
The log-likelihood ratio can be approximated by the Pearson-type statistic
 
    S_k = ∑_𝐰∈𝒜^k (  N_𝐰 - E_𝐰)^2 /  E_𝐰,

where 𝐰 = w_1 w_2⋯w_k
denotes a k-mer consisting of letters w_i ∈𝒜,
^-𝐰 = w_2⋯w_k, 𝐰^- = w_1
w_2⋯w_k-1, and ^-𝐰^- = w_2 ⋯w_k-1,
N_𝐰 denotes the count of the word
𝐰 in the sequence,
and E_𝐰 = N_^-𝐰 N_𝐰^-/N_^-𝐰^-
is the estimated
expected count of 𝐰 if the sequence is generated by a
MC of order (k-2), for k ≥ 3.
For k=2, N_^-𝐰^- is replaced by the total number of bases in the sequence.












Several estimators for the order of MC have been proposed based on the above results of the hypothesis testing.
Menndez et al. <cit.> proposed a procedure for estimating the order by performing a sequence of tests for increasing orders until the null hypothesis is accepted.
Papapetrou and Kugiumatzis <cit.> similarly used sequential hypothesis tests to find the optimal order of MC based on the significance of the conditional mutual information (CMI) of different orders.
Moray and Weiss <cit.>, Peres and Shields <cit.> and Dalevi et al. <cit.> developed methods to estimate the order of a MC based on the observation of a maximal sharp transition of |N_𝐰 - E_𝐰| at the true order. Baigorri et al. <cit.> estimated the order of MC by considering the change of χ^2-divergence involving S_k.
For the cases where a χ^2-test fails due to inefficient data, Besag and Mondal <cit.> provided exact goodness-of-fit tests for Markov chains.


Model selection approaches have also been widely used in the determination of the order of MC.
The Akaike information criterion (AIC) <cit.>, AICc <cit.>, the Bayesian information criterion (BIC) <cit.>, and the Efficient Determination Criterion (EDC) <cit.> were proposed to estimate the order of MC, and their consistency were studied in Katz <cit.> and Peres and Shields <cit.>.
 

All of these model selection methods were formulated using the logarithm of the maximum likelihood of the sequence and a penalty term related to the number of parameters in the model.
Let X be a sequence under the r-th order Markov model ℳ_r.
Then the log-maximum likelihood of the data under the model ℳ_r is

    l (X; ℳ_r) = ∑_w_1 w_2 … w_r ∈𝒜^r, w ∈𝒜 N_w_1 w_2 … w_r wlog( P̂(w | w_1 w_2 … w_r)),

where P̂(w | w_1 w_2 … w_r)=N_w_1 w_2 … w_r w/N_w_1 w_2 … w_r is the estimated transition probability.
Then the optimal order r^* of the MC is found by minimizing various criteria as follows.

    AIC(r)    = -2 l (X; ℳ_r) + 2 |ℳ_r|, 
    
    AICc(r)    = AIC(r) + 2|ℳ_r| ( |ℳ_r| + 1)/( |X_r| - |ℳ_r| - 1), 
    
    BIC(r)    = -2 l (X; ℳ_r) + |ℳ_r| log |X_r|, 
    
    EDC(r)    = -2 l (X; ℳ_r) + |ℳ_r| c(|X_r|),

where |X_r| is the data size of X_r, i.e. the total number of (r+1)-words in the sequence, |ℳ_r| is the number of parameters in the model (L^r × L in this case), and c(·) is a general increasing function.
Narlikar et al. <cit.> evaluated the AIC, AICc and BIC methods for estimating the order of a MC of a genomic sequence.
The results showed that the order of a MC had marked effects on the performance of sequence clustering and classifications.
The MC order obtained based on the BIC optimality criterion yielded the best performance among all the model selection criteria.




 §.§ INFERENCE OF MC PROPERTIES BASED ON NGS DATA


One successful application of alignment-free methods is comparing different genomes using NGS reads data for which each sample contains millions of short reads randomly sampled from different parts of the genomes.
For NGS reads data, it is challenging to assemble short reads to recover the original genomic sequences.
Ren et al. <cit.> developed an assembly-free method to estimate background MCs based solely on short reads.
The NGS reads data are modeled as generated by a two-layer stochastic process: first, a (un-observed) long MC sequence is generated, and second, short reads are randomly sampled from the long MC sequence.

The classic statistic S_k defined in equation <ref> for the long sequence was extended to S_k^R for the NGS data by replacing the word frequencies in a long sequence with that in NGS short reads.
Let N_𝐰^R be the count of the k-word 𝐰 in the NGS short reads (the superscript R refers to the “read” data). Define

    S_k^R = ∑_𝐰∈𝒜^k (  N_𝐰^R - E_𝐰^R )^2 /  E_𝐰^R .

Due to the additional randomness introduced in the process of sampling short reads from genomic sequences, the new statistic S_k^R no longer follows the classic χ^2- distribution.
Instead, it was shown that S_k^R follows a gamma distribution when the reads are sampled based on the Lander-Waterman model <cit.>.
In particular,
let f_i be the fraction of the genome that is covered by exactly i reads, i=1,2,⋯. Define the effective coverage

    d = ∑_i i^2 f_i/∑_i i f_i.

The statistic S_k^R/d has an approximate χ^2-distribution
    with df_k = (L-1)^2 L^k-2 degrees-of-freedom; equivalently, the
    statistic S_k^R has an approximate gamma distribution with shape
    parameter df_k/2 and scale parameter 2d.
Several estimators for the order of a MC based on NGS data using various criteria, such as the sharp transition of S_k^R, AIC and BIC, were proposed and compared in <cit.>, by extending the classical order estimators for long genomic sequences to those for NGS data.





§ APPLICATIONS OF THE ALIGNMENT-FREE METHODS TO COMPARATIVE GENOMICS


Among the various alignment-free sequence comparison methods, the measures using normalized k-mer counts, d_2^* and d_2^S <cit.>, have been shown to have
superior performance for comparing genomic sequences.
Wan et al. <cit.> and Burden et al. <cit.> studied the theoretical statistical properties of the d_2^* and d_2^S measures.
Song et al. <cit.> extended the definition of d_2^* and d_2^S from two long genomic sequences to the comparison of two samples based on NGS reads data, and investigated theoretically the properties of the measures.
As an application, the relationship of 13 tropical tree species in <cit.> were revealed without assembly using d_2^* and d_2^S.
Ren et al. <cit.> clustered genomic sequences of 28 vertebrate species based on NGS reads using d_2^* and d_2^S under different MC models.
Using the appropriate order of MC, the pairwise dissimilarity scores using d_2^* and d_2^S are highly correlated (Spearman's rank correlation coefficient 0.92) with the true pairwise evolutionary distances inferred based on the multiple sequence alignment of homologous genes in <cit.>.
Compared to d_2^*, d_2^S is less affected by the order of the MC model. For example, the Spearman's rank correlation coefficient using d_2^S  is 0.86 even under the i.i.d model.

Bernard et al. <cit.> and  Chan et al. <cit.>  systematically assessed the performance of various alignment-free measures under different evolutionary scenarios using simulations and empirical data.
The results showed that the alignment-free methods are sensitive to sequence divergence, less sensitive to lateral genetic transfer, and robust against genome rearrangement, among-site rate heterogeneity and compositional biases.
Chan et al. <cit.> performed phylogenetic inference using alignment-free measures for 4,156 nucleotide sequences.
The topology obtained using d_2^S was most congruent with the phylogeny inferred using multiple sequence alignment.
Similarly, the relationship among 143 bacteria and archaea genomes <cit.>, 63 Enterobacteriaceae genomes <cit.>, 27 Escheriachia coli and Shigella genomes <cit.>, 21 primate genomes <cit.> and 27 primate mitochondrial genomes <cit.>, 14 plants <cit.>, and 8 Yersinia genomes <cit.> were inferred using d_2^S and compared with the evolutionary tree built based on multiple sequence alignment. Despite some incongruence, the clustering results in general had highly similar structures with the classical evolutionary trees.

To evaluate the robustness of the clustering, different resampling methods, including jackknife <cit.> and bootstrap <cit.>, were applied for resampling sequences to provide a measure of robustness for the branches in the inferred clustering tree.
The studies showed that alignment-free methods can accurately recover phylogenetic relationship even with low sequencing coverage.
The time complexity for alignment-free methods was significantly lower compared to the traditional maximum likelihood and Bayesian methods based on multiple sequence alignment <cit.>.
It was estimated that alignment-free methods are approximately 140-fold faster than the traditional methods <cit.>.
Normalization of the background and including inexact matches increases the time complexity.
Alignment-free methods based on k-mers lend themselves to parallel algorithms, and parallel computational methods have been applied to achieve speedup and scalability for alignment-free methods <cit.>.
When k is large, memory is a main limitation for storing k-mer counts and computing alignment-free measures <cit.>.






§ PREDICTION OF VIRUS-PROKARYOTIC HOST INTERACTIONS USING ALIGNMENT-FREE METHODS


It is widely recognized that bacteria and archaea (prokaryotes) play important roles in many ecosystems and significantly impact the health of humans, animals, and plants <cit.>.
However, much less is known about the viruses that infect prokaryotes.

Since viral infections can lead to lysis of host cells, viruses consequently can indirectly impact ecological processes by regulating and controlling the abundance of prokaryotes.
Metagenomic sequencing, that uses NGS to recover genetic material of microbial organisms from environment samples, can be used for high-throughput identification of bacteria, archaea, and viruses regardless of culturability.
Increasing numbers of new viruses have been discovered by assembling short reads from various environments including human gut <cit.>, ocean <cit.>, and soil <cit.>.
Yet, their biological functions and prokaryotic hosts cannot be directly inferred from the metagenomic data.

A few computational approaches have been developed recently for predicting the host given a viral sequence.
The most straight forward method is alignment-based gene homology search and CRISPR search between virus and host genomes <cit.>.
However, not many viruses share regions with hosts and not many hosts have CRISPR spacers.
In contrast, alignment-free methods can be powerful for revealing virus-host interaction relationships, because it is observed that viruses share highly similar k-mer usage with their hosts, possibly due to the fact that virus replication is dependent on translational machinery of its host <cit.>.

Edwards et al. <cit.> and Roux et al. <cit.> used Euclidean and Manhattan distances based on tetramers (k=4) to measure the distance between viruses and hosts, and predicted the host as the one with the smallest distance to the query virus.

Ahlgren et al. <cit.> conducted a comprehensive evaluation of alignment-free dissimilarity measures over various k-mer lengths for host prediction.
The study evaluated a suite of 11 measures including those based on the observed word frequencies such as Euclidean and Manhattan distances and those based on normalized word frequencies such as d_2^* and d_2^S.
The prediction accuracy of the measures were assessed based on the largest benchmark dataset containing 1,427 virus isolate genomes whose true hosts are known and ∼32,000 prokaryotic genomes as host candidates.

The measures based on normalized frequencies in general have better discriminatory power of separating true interacting virus-host pairs from random pairs than those based on observed word frequencies.

Increasing k-mer length from 4 to 6 also improves the discriminatory power.

Among the 11 measures, d_2^* at k= 6 and a second order MC yielded the highest host prediction accuracy (Figure <ref>).
Requiring a minimum dissimilarity score for making predictions (thresholding) and taking the consensus
of the 30 most similar hosts further improved accuracy.
While prediction accuracy decreases for shorter contigs,
the method is able to make decent predictions on contigs as short as 5 kbp.
A software called VirHostMatcher was developed for predicting hosts of viruses and visualizing the predicted results using alignment-free methods.




Following the same principle that the virus-host genomes tend to have high similarity, Galiez et al. <cit.> developed a program, WIsH, that computes the likelihood of the query viral sequence under each of the Markov models for candidate bacteria genomes, and predicts the host as the one whose model yields the highest likelihood.

Since the program only relies on the Markov models for bacteria complete genomes, the method achieves decent accuracy even for viral contigs as short as 3 kbp, and it is generally faster than VirHostMatcher.
WIsH uses a fixed 8th order MC to model the bacteria genomes,
so the method may not be readily applicable for metagenomic contigs where the host contigs are so short that no sufficient data is available for estimating a high order MC.

Another group of host prediction methods is based on the observation that similar viruses often share the same host range.
Different virus-virus similarity measures have been investigated using various principles <cit.>, and the clusters in the gene-based virus-virus similarity network show high association with the host classes <cit.>.
Villarroel et al. <cit.> developed a host prediction tool, HostPhinder, that predicts the host of a query virus as the host of the most similar reference virus. The similarity was defined based on the proportion of the shared k-mers between the query and the reference virus genomes.
Zhang et al. <cit.> developed machine learning based classifiers to predict if a query virus can infect a particular host genus, based on the common k-mer features learned from the existing infectious viruses.
However, the method is only applicable to hosts that have a relatively large number of known infecting viruses.




§ GENOME AND TRANSCRIPTOME COMPARISON USING ALIGNMENT-FREE APPROACHES WITH VARIABLE LENGTH MARKOV CHAINS


Using Fixed Order Markov Chains (FOMC) to model the background sequence has several potential limitations.  First, the MC order needs to be set manually. However, for most sequences of interest, there is no prior knowledge available for setting the correct MC order. Second, FOMC is not structurally rich. The number of parameters in an r-th order MC is (L-1) L^r where L is the alphabet size, and there are no MC models with number of parameters between (L-1) L^r and (L-1) L^r+1.
Third, the number of parameters grows exponentially with the  MC order r. When the length of the sequence is short or sequencing depth is relatively low, the parameters cannot be accurately estimated.

Therefore, Liao et al. <cit.> investigated the use of  the data-driven Variable Length Markov Chain (VLMC) <cit.> model as an alternative to FOMC to model background sequences. VLMC was originally designed for modeling one long sequence and was represented as a context tree structure <cit.>. Liao et al. <cit.> designed a three-step approach for prunning a tree based on NGS short reads data. First, a full prefix tree  based on 1, 2, ⋯, 10-mer frequency vectors was built. However, the tree usually overfits the data.
Second, the full prefix tree was pruned to remove the redundant branches based on the Kullback-Leibler divergence <cit.>. The pruned tree is called a context tree <cit.>. The threshold value K for the Kullback-Leibler divergence determines the complexity of the pruned tree. The value of K was chosen by optimizing the Akaike Information Criterion (AIC) <cit.> designed for the high-throughput sequencing data. AIC measures the relative quality of statistical models for a given set of data. Third, transition probabilities were estimated with respect to the VLMC from the context tree, and the probabilities of words were then computed accordingly.


Using VLMC for background modeling, beta diversity measures d^S_2and d^*_2 were applied to one simulated dataset and four real datasets to compare bacterial transcriptomes or metatranscriptomes. According to the previous experiments, generally the optimal k is 6-9. For comparison, d^S_2and d^*_2with 0-4th order FOMC, three L_p-norm measures andd_2were also applied to the testing datasets. Triplet based distances was selected to evaluate the topological consistency between the reference and the clustering trees due to its robust and fine-grained measure.

The simulated metatranscriptomic dataset was composed of 90 samples belonging to 3 different groups with 5000 genes from 5 microbes. d^S_2 on VLMC depicted clear groups of three abundance distributions among samples. Four real datasets were transcriptomic data from marine microbial eukaryotes, global ocean metatranscriptomic samples, different depths of ocean metatranscriptomic and metagenomic samples, and metatranscriptomic samples from different iron-rich microbial mats. The first two datasets were used to evaluate the effect of VLMC-based measures in identifying group relationships. The last two datasets were to study the performance of VLMC-based measures in revealing environmental gradient relationships. (1)The transcriptomic dataset is composed of 18 RNA-Seq data from marine microbial eukaryotes of Phylum Chlorophyta. The molecular phylogeny reference tree <cit.> was reconstructed based on the 18S rRNA genes with maximum likelihood (ML) method. The clustered tree on VLMC represented very similar topological structure of the reference phylogenetic tree, which demonstrated that VLMC is able to depict the genetic relationships among the species under same Phylum. (2) The 88 metatranscriptomic samples were collected from 7 geological locations in global ocean by 12 study projects. Except two samples, all other samples were consistently grouped according to the marine locations. VLMC reveals clear location relationships among the samples. (3) The four depths of ocean samples were composed of 8 metatranscriptomic and 8 metagenomic samples from underwater 25m, 75m, 125m and 500m. Samples from the same depth were clustered first, then the samples belonging to the photic zone (25m, 75m and 125m) were merged, and finally, samples belonging to the mesopelagic zone (500m). (4) For 14 metatranscriptomic samples from two depths of 0.03m and 0.08m within a typical iron-rich microbial mat, samples from surface water and from deeper regions were clearly separated. The PCA plot based on VLMC reflected the gradient information for collection depths and sites as the first and second principal component. In contrast, FOMC mixed surface samples with deeper samples successively.

Liao et al. <cit.> evaluated the performance of  d^S_2 and d^*_2 using both simulations and real data. It was shown that VLMC outperformed FOMC to model the background sequences in transcriptomic and metatranscriptomic samples.
Moreover, d^S_2 based on VLMC background model can identify underlying relationships among metatranscriptomic samples from different microbial communities, and can reveal a gradient relationship among the metatranscriptomic samples. VLMC is easier to apply than FOMC because of being free from MC order selections. The flexible number of parameters in VLMC avoids estimating the vast number of parameters of high-order MC under limited sequencing depth. 
In contrast, the VLMC model does not work as well as FOMC for investigating the relationship among whole genome or metagenome data. It was hypothesized that whole genomes and metagenomes contain mixtures of coding and noncoding regions and are too complex to be modeled by relatively concise VLMC models.
Yet, the coding regions are more homogeneous than the whole genome.
The clustering performance can be improved for metatranscriptomic data using the VLMC to model the background sequence, but not for whole genome or metagenomic data.

For the comparison of metagenomes, Jiang et al. <cit.> showed that  d_2^S with the i.i.d background model and k-mer length between 6 to 9 bps generally performs well compared to other measures.

It is time-consuming to model VLMC due to the generation and the pruning of the prefix tree. Behnam and Smith <cit.> measured the dissimilarity between metagenomic samples with dot product distance based on the i.i.d. model, and they integrated a randomized hashing strategy based on locality-sensitive hashing and the regular nearest neighbor graph to reach logarithmic query time for identifying similar metagenomes even as the database size reaches into the millions. Meanwhile, also focusing on fast comparisons among large-scale multiple metagenomic samples, Benoit et al. <cit.> developed the program,  Simka, to compute 16 standard ecological distances by a parallel k-mer counting strategy on multiple data sets. Simka was able to compute in a few hours both qualitative and quantitative ecological distances based on hundreds of metagenomic samples.





§ IMPROVING METAGENOMIC CONTIG BINNING USING D_2^S 


 Wang et al. <cit.> used d_2^S to improve contig binning. Assigning assembled contigs into discrete clusters, known as bins, is a key step toward investigating the taxonomic structure of microbial communities <cit.>. Contig binning using k-mer composition is based on the observation that relative sequence compositions are similar across different regions of the same genome, but differ between distinct genomes <cit.>. Contigs in the same bin are expected to come from the same taxonomic group.  Three different types of strategies have been used to bin contigs: sequence composition, abundance and a hybrid between the two.  Sequence composition based methods use k-mer frequencies with k=2-6 as genomic signatures of contigs <cit.>. Abundance based methods use the relative abundance levels of species and the distribution of the number of reads containing certain k-mers to bin contigs  <cit.>. The hybrid approaches use both composition and abundance of k-mers to bin contigs <cit.>. Most of the currently available binning methods used the frequency of k-mers directly, but this represented absolute, not relative, sequence composition. Here, “absolute" frequency refers to the number of occurrences of a k-mer over the total number of occurrences of all k-mers. On the other hand, “relative" frequency refers to the difference between the observed frequency of a k-mer and the corresponding expected frequency under a given background model.
The dissimilarity measures d^S_2 based on relative frequencies of k-mers have been successfully used for sequence comparison as reviewed above.
Therefore, we expected that calculating the dissimilarity between contigs using d_2^S would improve contig binning compared to other contig binning methods that are based on the difference of absolute k-mer frequencies. However, directly using d_2^S for contig binning is too time consuming and is impractical for most metagenomic data.

Instead of binning contigs directly using d_2^S, Wang et al. <cit.> developed d^S_2Bin that uses d_2^S to improve reasonable contig binning results using other fast and efficient programs such as MetaCluster3.0 <cit.>, MetaWatt <cit.>, SCIMM <cit.>, MaxBin1.0 <cit.>, and MyCC <cit.>. Each contig was modeled with a MC based on its k-mer frequency vector. The center of the bin was represented by the average k-mer frequency vectors of all contigs in this bin and was also modeled with a MC. Then, d^S_2 was used to measure the dissimilarity between a contig and the center of a bin based on relative k-mer composition. Finally, a K-means clustering algorithm was applied to cluster the contigs based on the d^S_2 dissimilarities.  Recall, precision and Adjusted Rand Index (ARI) were used to evaluate the binning performance. Wang et al. <cit.> showed that d^S_2Bin consistently achieved the best performance with k=6-mers under the i.i.d. background model.
d^S_2Bin improves the binning performance in 28 out of 30 testing experiments.
Experiments showed that d_2^S accurately measures the dissimilarity between contigs of metagenomic reads and that measures defined based on relative sequence composition are more suitable for contig binning.
Also, d^S_2Bin can be applied to any existing contig-binning tools for single metagenomic samples to improve binning results.

To evaluate the performance of d^S_2Bin, five representative binning tools on sequence composition (MetaCluster3.0 <cit.>, MetaWatt <cit.> and SCIMM <cit.>) and the hybrid of sequence composition and abundance (MaxBin1.0 <cit.>, MyCC <cit.>) were selected to bin six testing datasets, after which d^S_2Bin was applied to adjust the binning results. The six testing datasets included five synthetic and one real datasets, which were generated, or used by previous binning tools. These datasets covered various species diversity, species dissimilarity, sequencing depth, and community complexity to yield a comprehensive evaluation of d^S_2Bin. The five synthetic metagenomic datasets had been used by MaxBin1.0 <cit.>, with different setting of number of genomes, coverages and complexities in abundance. The real dataset had been applied to test the binning tools COCACOLA <cit.> and CONCOCT <cit.>, composed of a time-series of 11 fecal microbiome samples from a premature infant <cit.>. All metagenomic sequencing reads were merged and assembled, and contigs unambiguously aligned were set as ground truth <cit.>. d^S_2Bin improves the binning performance in 28 out of 30 testing experiments (6 datasets with 5 binning tools). In most cases, the three criteria were improved by 1%-22%. The three indexes increased significantly on the first iteration and reached steady state quickly. Irrespective of the different strategies employed by the contig-binning tools, d^S_2Bin was able to achieve better performance for all tools tested. d^S_2Bin can be applied to any existing contig-binning tools for individual metagenomic samples to obtain better binning results. d^S_2Bin is available at https://github.com/kunWangkun/d2SBin.


d^S_2Bin gives credence to the relative sequence composition model over the direct application of absolute sequence composition. Current "Taxonomy-independent" binning tools extract features from contigs to infer bins based on the three strategies. (1) Sequence composition. It is usually denoted as frequencies of k-mers with k=2-6 as genomic signatures of contigs <cit.>. (2) Abundance. It estimated the relative abundance levels of species on Poisson distributions of k-mers <cit.>. (3) Hybrid of composition and abundance. It combined k-mer frequencies, marker genes and scaffold coverage levels <cit.>. Most of them used the frequency of k-mers directly, but this represented absolute, not relative, sequence composition. Here, "absolute" frequency refers to the number of occurrences of a k-mer over the total number of occurrences of all k-mers. On the other hand, "relative" frequency refers to the difference between the observed frequency of a k-mer and the corresponding expected frequency under a given background model. Therefore, relative sequence composition is more reasonable to bin the contigs.




§ IMPROVING THE IDENTIFICATION OF HORIZONTAL GENE TRANSFER USING D_2^* OR CVTREE



Horizontal gene transfer (HGT) or lateral gene transfer (LGT) describe the transmission of genetic material between organisms that are not in a parent-offspring relationship.
HGT plays an important role in the evolution of microbes and is responsible for metabolic adaption <cit.> and the spread of antibiotic resistance <cit.>. Existing computational methods for HGT inference can be broadly separated into two groups: alignment-based and alignment-free methods.


Alignment-based, or phylogenetic methods for detecting HGT rely on phylogenetic conflicts; that is, finding genes whose phylogenetic relationships among multiple organisms differ significantly from that of other genes <cit.>. Although alignment-based methods are considered to be the gold standard <cit.> for HGT detection because of their explicit models, finding topological incongruences is computationally demanding, requires large memory, and requires that genomes of interest are annotated and their phylogenetic relationships are known.
In addition, alignment-based methods can only be applied to coding sequences and thus have limited ability to detect horizontal transfer in non-coding regions.


Instead, alignment-free methods, also called compositional parametric methods, can be used to avoid these limitations. Alignment-free methods infer horizontal gene transfer by detection of regions in a genome with atypical word pattern composition based on the observation that sequences transferred from donor genomes have different composition signatures from that of the host genome <cit.>. Recently, Cong et al. <cit.> introduced TF-IDF as a scalable alignment-free approach for HGT detection in large molecular-sequence data sets by combining multiple genomes and k-mer frequencies. However, these methods require the phylogenetic relationship among a group of genomes and they can only detect HGT within this group of genomes. More widely used alignment-free methods apply a sliding window approach to scan a single genome and calculate the dissimilarity between each window and the whole genome. Consecutive windows with dissimilarity higher than a threshold are inferred as HGT. The performances of k-mer-based alignment-free methods depend largely on the choice of dissimilarity measures between a genomic region and the whole genome, on the k-mer length, on the sliding window size, and on the evolutionary distance between host and donor genomes. Manhattan and Euclidean distances between the k-mer frequency vector of a genomic region and that of the whole genome are the most frequently used measures for detecting HGTs because of their simplicity. For example, Dufraigne et al. <cit.> analyzed HGT regions of 22 genomes by using Euclidean distance with k-mer length of 4 bps.
Rajan et al. <cit.> used Manhattan distance with k-mer length of 5 bps to detect HGT in 50 diverse bacterial genomes.


Several papers compared the performances of different dissimilarity measures for HGT detection. Because the true HGT history is unknown, the evaluation and benchmarking of HGT detection methods typically relies on simulated artificial genomes, for which the true simulated history is known.
Tsirigos and Rigoutsos <cit.>  investigated several dissimilarity measures between the relative frequencies of a genomic region and the whole genome under the i.i.d. model including correlation, covariance, Manhattan distance, Mahalanobis distance, and Kullback–Leibler (KL) distance for HGT detection. They showed that k-mers of length 6-8 bps with covariance dissimilarity perform the best under their simulated situations.
Becq et al. <cit.> reviewed alignment-free methods on horizontal gene transfer detection and showed that k-mer-based methods with a 5 kbps sliding window outperformed other alignment-free methods based on features such as GC content <cit.>, codon usage <cit.> and dinucleotide content <cit.>. However, they only tested Euclidean distance with k-mer length 4 bps as genomic signature <cit.> for k-mer-based methods.  


Recently, we evaluated the performance of different dissimilarity measures including Manhattan, Euclidean,  CVtree, d_2, d_2^*, d_2^s with different choices of k-mer length and Markov order. We also studied the influence of window size and evolutionary distance between host and donor genomes on HGT detection by both simulation and real data in terms of F_1 score and precision-recall curve (PRC).
We showed that none of these dissimilarity measures work well when the donor and host genomes are within the same order level since the donor and host genomes are too similar and it is challenging to distinguish the transferred regions.
All dissimilarity measures perform well when the donor and host genomes are in different class levels since the host and donor genomes are highly different and most of these methods can identify their differences.
For HGT between genomes  from different order levels but in the same class level, background adjusted dissimilarity measures that consider Markov order of sequences, such as CVtree with k = 4 and d_2^* with k = 3 and Markov order 1 can achieve significantly better performance than the other methods. The PRC results for different scenarios are shown in Figure <ref>.


Therefore, k-mer-based alignment-free methods for HGT detection are suitable when host and donor genomes are in different order levels and HGT length is greater than 5 kbps.
Therefore, alignment-free methods should not replace alignment-based methods in all cases. Instead, they are complimentary as each has unique advantages in different scenarios and they also tend to find complimentary sets of HGT regions <cit.>. Alignment-free methods are preferred when no evolutionary trees are available or genomes are not well annotated. Our study suggests that CVTree with word length of 4, d_2^* with word length 3, Markov order 1, and d_2^* with word length 4, Markov order 1 all perform well in most situations.
















§ OTHER WORD-COUNT BASED APPROACHES FOR SEQUENCE COMPARISON


Many other sequence dissimilarity measures based on k-mer frequencies have been developed in recent years. 
Liu et al. <cit.> proposed local alignment-free measures by summing up the maximal pairwise scores between any sub-fragments of a fixed length in the sequence.
Ren et al. <cit.> developed a suite of alignment-free multiple sequence comparison methods to enable measuring similarity among a set of more than two sequences. Several alignment-free  methods incorporating potential mismatches, sequencing errors,  or spaced word patterns have been developed for sequence comparison
<cit.>.
Fan et al. <cit.> developed a method called Assembly and Alignment-Free (AFF) that defines the distance based on the proportion of shared k-mers as an indication of the amount of divergence between the species.


In most of the dissimilarity measures reviewed above,  the k-mers are treated equally. Differential weighting of the k-mers may help study the relationship among the sequences.
Patil and McHardy <cit.> generalized Euclidean distance to a weighted Euclidean distance where the weights are learned from the training data, and evaluated on the independent test data. The learned weighted Euclidean distances specified for a group of species increase the accuracy for inferring taxonomic relationships of a new species from the same group.

Qian and Luan <cit.> developed an alternative approach for weighting the different k-mers by maximizing the weighted L_1 norm between the frequency vectors among all the sequences with c_𝐰 being the weight for the word 𝐰.
Qian and Luan <cit.> proposed to maximize

    ∑_𝐰∈𝒜^k∑_i,j = 1^n  c_𝐰 | f_i 𝐰 - f_j 𝐰|,

with the constraint of  ∑_𝐰∈𝒜^k c_𝐰 = 1,
where n is the number of sequences to be compared. Once the values of c_𝐰 were determined,  they modified the definitions of d_2, d_2^* and d_2^S by putting the weight c_𝐰 in front of the corresponding terms. Applications to the identification of homologous genes and cis-regulatory modules (CRM) showed that the weighted versions of these measures outperformed the original ones.

 It was reasoned that if a k-mer is present/absent in a small fraction or most of the sequences, it does not markedly contribute to distinguishing the different sequences. Therefore, weighting the different k-mers according to the frequency of being present/absent in the sequences of interest can increase our understanding of the relationships among the sequences <cit.>. For a k-mer 𝐰, let F_𝐰 be the fraction of the sequences with 𝐰 present. The entropy of the word is defined as

    H_𝐰 = - ( F_𝐰log_2 (F_𝐰)  + (1 - F_𝐰) log_2 (1  - F_𝐰)).

The weighted similarity measure between sequence i and sequence j was defined as

    K_ij = ∑_𝐰∈A^k  H_𝐰 f_i𝐰 f_j𝐰,

 and then normalized using

    K^'_ij = K_ij/√(K_ii K_jj).

Finally, the dissimilarity between the two sequences was defined as d_ij = √(2(1 - K^'_ij)) <cit.>. To speed up computational time as well as to save memory, Murray et al. <cit.> bin the k-mers into different groups so that a group contains multiple k-mers. The authors showed that this weighted version outperformed the traditional d_2 statistic and the Mash program <cit.>.












§ DETERMINATION OF WORD SIZE K

In alignment-free sequence comparison using word counts, an important yet challenging problem is the length of word patterns. Although many studies are available, there are still no definitive answers to the optimal choice of word length. The optimal word length depends on the statistical measures for comparing the sequences, and the background models, the lengths, and the diversity of the sequences to be compared. For example, if the sequences are short, the optimal word length may be short since the sequences do not contain a large number of distinct words. Otherwise, the sequences may rarely share common word patterns.
However, short word patterns do not have high power to discriminate closely related sequences. If the sequences to be compared are highly similar, we expect that the optimal word length should be long as short word patterns will not be able to distinguish them. On the other hand, if the sequences to be compared are diverse, relatively short word patterns  may suffice to distinguish the sequences.

Recently, Bai et al. <cit.> investigated the optimal word length when comparing two Markovian sequences using the χ^2-statistic in <cit.>.
Bai et al. <cit.> framed sequence comparison as a hypothesis testing problem of evaluating if the two sequences come from two different Markov chains and used power under the alternative hypothesis as an optimality criterion.  They showed both theoretically and by simulations  that the optimal word length equals the maximum of the Markov orders of the two sequences plus one. This conclusion also holds for NGS data. Using the estimated Markov orders resulted in minimal loss of power when comparing two sequences. Applications to real sequences to find homologs of the human protein HSLIPAS and the cis-regulatory modules (CRM) in four mouse tissues (forebrain, heart, limb and midbrain) confirmed the theoretical results. Preliminary simulation results showed that this k-mer length may also be optimal for other measures including  CVTree <cit.>, d_2^* and d_2^S <cit.>. However, we could not prove this claim theoretically.


In a series of papers, Kim and colleagues <cit.> investigated the optimal word length when using the Jensen–Shannon (JS) divergence between the word frequency vectors to measure the dissimilarity between two sequences. The lower limit of the word length was suggested as log_L(n), where n is the average length of the sequences to be compared and L is the alphabet size. To obtain an upper bound, they defined cumulative relative entropy (CRE) as follows. Let F_k = (f_𝐰, 𝐰∈A^k) be the frequency vector of all the words of length k and F̂_k  = (f̂_𝐰, 𝐰∈𝒜^k) be the corresponding expected frequency under the k-2-th order Markov chain.  f̂_𝐰 =f̂_𝐰/(n-k+1).  For a given word 𝐰 = w_1 w_2 ⋯ w_k, define ^-𝐰 = w_2 ⋯ w_k,  𝐰^- = w_1 w_2 ⋯ w_k-1, ^-𝐰^- = w_2 w_3 ⋯ w_k-1. Then

    f̂_𝐰 = f_^-𝐰 f_𝐰^-/f_^-𝐰^-.

The CRE function is defined by

    CRE(t) = ∑_k = t^∞ KL(F̂_k, F_k),

where KL is the Kullback-Leibler divergence. The upper bound of the optimal k is the value of t such that CRE(t) is close to zero. In practice, they used the t such that CRE(t) is less than 10% of the maximum CRE.  
For the pairwise comparison among a set of sequences, if the lengths of the sequences to be compared are not highly different,

the above approach will give similar lower and upper bounds for the optimal word length.
The final k-mer length can be chosen within the overlapping ranges of the optimal word length among the sequences.
If the sequences have highly different lengths, the authors suggested to divide the large genomes into blocks of equal length so that the sequences to be compared have similar length. They applied the method to investigate the relationships among the Escherichia coli/Shigella group <cit.>, prokaryotes <cit.>, and dsDNA viruses <cit.>.
Recently, Zhang et al. <cit.> used the approach to investigate the relationship among close to 4,000 viruses with very different lengths. They divided the virus genomes into four groups according to their lengths. First, based on the CRE criterion in <cit.>, the upper limit for the optimal word length was determined to be between 9 to 13. Secondly, the authors defined the average common features (ACF) of a sequence as the average shared feature of this sequence with the other sequences. 












Another commonly used weighting method is to weight each k-mer by the TF-IDF (Frequency-Inverse Document Frequency) <cit.>. TF-IDF was first used for the classification of documents <cit.> and was recently used to study the relationship among molecular sequences. 






§ INTEGRATED SOFTWARE FOR ALIGNMENT-FREE SEQUENCE COMPARISON


As reviewed in the above sections, a large number of alignment-free sequence comparison approaches have been developed and most of the individual studies have accompanying software tools available. To facilitate the use of the different alignment-free methods, a general-purpose alignment-free platform is desirable, which is expected to include the support of both assembled genome sequences and unassembled NGS shotgun reads as input, integration of exhaustive alignment-free sequence comparison measures, and visualization of results.

CAFE <cit.> is a stand-alone alignment-free sequence comparison platform for studying the relationships among genomes and metagenomes through a user-friendly graphical user interface. Overall, CAFE integrates 28 distinct alignment-free measures, including 10 conventional measures based on k-mer counts (e.g., Euclidean, Manhattan, d_2, Jensen-Shannon divergence <cit.>, feature frequency profiles (FFP) <cit.>, Co-phylog <cit.>, etc.), 15 measures based on presence/absence of k-mers (e.g., Jaccard, Hamming, etc.), and 3 measures based on background adjusted k-mer counts (CVTree <cit.>, d_2^* <cit.>, and d_2^S <cit.>). All measures have been evaluated using whole primate and vertebrate genomes, whole microbial genomes, and NGS short reads from mammalian gut metagenomic samples. CAFE significantly speeds up the calculation of the background-adjusted measures such as CVTree, d_2^*, and d_2^S, with reduced memory requirements. Moreover, the resulting pairwise dissimilarities among the sequences form a symmetric distance matrix, which can be directly saved in a standard PHYLIP format (http://evolution.genetics.washington.edu/phylip/credits.html). CAFE also provides four types of built-in downstream visualized analyses, including clustering the sequences into dendrograms using the UPGMA algorithm, heatmap visualization of the matrix, projecting the matrix to a two-dimensional space using principal coordinate analysis (PCoA), and network display. A screenshot of CAFE is shown in Figure  <ref>.


Alternatively, Alfree <cit.> provides a publicly accessible web-based sequence comparison platform for studying the relationships among nucleotide and protein sequences. Alfree integrates 38 popular alignment-free measures, including 25 word-based measures (e.g., Euclidean, Minkowski, FFP, Jaccard, Hamming, etc.), 8 Information-theoretic measures (e.g., Lempel-Ziv complexity <cit.>, normalized compression distance <cit.>, etc.), 3 graph-based measures <cit.>, and 2 hybrid measures (i.e., Kullback–Leibler divergence <cit.> and W-metric <cit.>). The majority of measures have been evaluated using simulated DNA sequences, primate mitochondrial genomes, prokaryotic genomes and proteomes, plant genomes, etc. Moreover, the resulting dissimilarities among the sequences are reported as phylogenetic trees, heat maps, and tables. 








With the advances of efficient and affordable sequencing technologies, the high volumes of sequence data have brought computational challenges even for alignment-free sequence comparison. This concern is alleviated by Mash <cit.> that uses the MinHash dimensionality-reduction technique to reduce large amount of sequences to compressed sketch representations. Generally, Mash estimates the Jaccard distance between pairwise k-mer vectors in terms of compressed sketch representations, with moderate memory and computation overhead.
Similarly, kWIP <cit.> counts k-mers, hashes them into a compressed sketch, and introduces an information-theoretic weighting to elevate the relevant k-mers against irrelevant ones. Finally, it computes the similarity as inner products of weighted frequency vectors, normalized by Shannon entropy.
In addition, Benoit et al. <cit.> developed a program, Simka, for fast calculation of various distance measures between sequences for k-mers up to 30 bps long.

















§ DISCUSSION AND CONCLUSIONS

With the development of NGS technologies, huge amounts of sequencing data can be generated efficiently and economically. Sequence comparison plays crucial roles to analyze the large amount of sequence data and to extract biological knowledge from them. Although alignment-based sequence comparison will continue to dominate molecular sequence analysis, alternative alignment-free sequence comparison has become increasingly important due to its efficiency in analyzing huge amount of sequence data as well as its comparable performance with alignment-based methods. In recent years, there is a surge of interest in using alignment-free sequence comparison approaches for investigating a variety of different problems including the study of evolutionary relationships of whole genome sequences and gene regulatory regions, comparison of metagenomes and metatranscriptomes, binning of contigs, detection of horizontal gene transfer, and virus-host infectious associations based on NGS data. Among the many types of alignment-free sequence comparison approaches, word-count based approaches are most popular due to their easy adaption to NGS data.

Most word-count based alignment-free approaches use the absolute word frequencies for sequence comparison. These approaches have the advantage of being simple, easy to calculate, and using less memory. On the other hand, relative word frequency based alignment-free methods that were originally developed by Karlin's group <cit.> and Hao's group <cit.> and were recently revitalized by our group <cit.>  outperformed absolute word-count based approaches in all the applications we have investigated including the comparison of genomes <cit.>, gene regulatory regions <cit.>,  metagenomes <cit.>, and metatranscriptomics <cit.>. They have also been used to improve the binning of contigs in metagenomes  <cit.> and to predict virus-host interactions <cit.>. By subtracting the expected word counts based on the background MC model from the observed word counts, the words distinguishing the sequences are strengthened while the weights of the irrelevant words are minimized resulting in the excellent performance of the background adjusted methods. However, the calculation of the background adjusted measures such as CVTree, d_2^* and d_2^S adds extra burdens in memory and computational speed. Further improvements to speed up the computation of these measures and to reduce memory are needed.

Although there have been some studies on the optimal choice of word length for some measures such as χ^2-statistic <cit.> and Jensen-Shannon entropy <cit.>, the optimal word length for many other measures is not known. In these studies, the optimal word length was determined by the individual sequences, not by the relationship among the sequences. We expect that for the comparison of highly divergent sequences, short word length should suffice, while for the study of closely related sequences, long word patterns are needed. However, no studies are available on the optimal word length considering the divergency among the sequences. A few recent studies <cit.> used long words of length up to 30 bps and absolute word frequencies to compare genome sequences with excellent results and fast computation speed. It will be interesting to compare the performance of these approaches with the background adjusted measures with relatively short k-mers under realistic assumptions on sequencing errors and NGS data.


With the large number of alignment-free sequence comparison measures available, it is time to establish some benchmark data sets to evaluate the pros and cons of the different measures. Zielezinski et al. <cit.> built a benchmark data set of protein structures and evaluated a variety of different alignment-free sequence comparison measures and the Smith-Waterman algorithm.
Following up from their data set, there is a need for a collection of community-agreed data sets for the comparison of genomes, gene regulation regions, and metagenomes.

In summary, alignment-free sequence comparison methods have shown great promise for NGS data analysis as shown by many applications. They are generally computationally fast and use less memory compared to alignment based methods. Further studies on the choice of length of k-mers, differential weighting of the k-mers, and benchmark data sets are needed to explore the full potential of alignment-free methods.






§ DISCLOSURE STATEMENT

The authors are not aware of any affiliations, memberships, funding, or financial holdings that might be perceived as affecting the objectivity of this review.




§ ACKNOWLEDGMENTS

We thank Drs. Nathan A. Ahlgren, David Chew, Minghua Deng, Jed A. Fuhrman, Bai Jiang, Kai Song, Lin Wan, Michael S. Waterman, Xuegong Zhang, Ms. Weinan Liao and Ms. Kun Wang for collaborations on the investigation of alignment-free sequence comparison and applications.   The preparation of the manuscript was supported by US NSF National Science Foundation (NSF) [DMS-1518001] and National Institutes of Health [R01GM120624]. Dr. Ying Wang was supported by National Natural Science Foundation of
China (61673324, 61503314),  China Scholarship Council (201606315011)
and Natural Science Foundation of Fujian (2016 J01316).














ar-style3













