Alignment-free phylogeny of whole genomes using underlying subwords

Background With the progress of modern sequencing technologies a large number of complete genomes are now available. Traditionally the comparison of two related genomes is carried out by sequence alignment. There are cases where these techniques cannot be applied, for example if two genomes do not share the same set of genes, or if they are not alignable to each other due to low sequence similarity, rearrangements and inversions, or more specifically to their lengths when the organisms belong to different species. For these cases the comparison of complete genomes can be carried out only with ad hoc methods that are usually called alignment-free methods. Methods In this paper we propose a distance function based on subword compositions called Underlying Approach (UA). We prove that the matching statistics, a popular concept in the field of string algorithms able to capture the statistics of common words between two sequences, can be derived from a small set of “independent” subwords, namely the irredundant common subwords. We define a distance-like measure based on these subwords, such that each region of genomes contributes only once, thus avoiding to count shared subwords a multiple number of times. In a nutshell, this filter discards subwords occurring in regions covered by other more significant subwords. Results The Underlying Approach (UA) builds a scoring function based on this set of patterns, called underlying. We prove that this set is by construction linear in the size of input, without overlaps, and can be efficiently constructed. Results show the validity of our method in the reconstruction of phylogenetic trees, where the Underlying Approach outperforms the current state of the art methods. Moreover, we show that the accuracy of UA is achieved with a very small number of subwords, which in some cases carry meaningful biological information. Availability http://www.dei.unipd.it/∼ciompin/main/underlying.html


Background
The global spread of low-cost high-throughput sequencing technologies has made a large number of complete genomes publicly available, and this number is still growing rapidly. In contrast, only few computational methods can really handle as input entire chromosomes, or entire genomes.
Traditionally the comparison of related genomes is carried out by sequence alignment. Popular methods extract gene-specific sequences from all species under examination and build a multiple sequence alignment for each gene [1]. Then all multiple sequence alignments are merged to form the final phylogeny. Other methods http://www.almob.org/content/7/1/34 prohibitive task even for supercomputers, hence simply infeasible. To overcome these obstacles, in the last ten years a variety of alignment-free methods have been proposed. In principle they are all based on counting procedures that characterize a sequence based on its constituents, e.g., k-mers [3,4].
An important aspect in phylogeny reconstruction is the fact that gene-based methods strictly focus on comparing the coding regions of genomes, which can account for as little as 1% of the genomic sequence in humans [5]. Whereas it is known that the use of whole genomes may provide more robust information when comparing different organisms [6]. Also most alignment-free methods in the literature use only a portion of complete genomes [7]. For instance, there are approaches that use only genic regions [3] or mitochondria; other methods filter out regions that are highly repetitive or with low complexity [4]. Recently, it has been shown that the evolutionary information is also carried by non-genic regions [8]. For certain viruses, we are not even able to estimate a complete phylogeny by analyzing their genes, since these organisms may share a very limited genetic material [7].
Sims et al. recently applied the Feature Frequency Profiles method (FFP) presented in [4] to compute a wholegenome phylogeny of mammals [8]-i.e., large eukaryotic genomes including the human genome -and of bacteria. This method needs to estimate the parameter k in order to compute a feature vector for each sequence, where the vector represents the frequency of each possible k-mer. Each feature vector is then normalized by the total number of k-mers found (i.e., by the sequence length minus k-1), obtaining a probability distribution vector, or feature frequency profile, for each genome. FFP finally computes the distance matrix between all pairs of genomes by applying the Jensen-Shannon [9] divergence to their frequency profiles.
This general characterization of strings based on their subsequence composition closely resembles some of the information theory problems, and is tightly related with the compression of strings. In fact, compositional methods can be viewed as the reinterpretation of data compression methods, well known in the literature [10,11]. For a comprehensive survey on the importance and implications of data compression methods in computational biology, we refer the reader to [12].
When comparing entire genomes we want to avoid that large non-coding regions, which by nature tend to be highly repetitive, may contribute to our scoring function a multiple number of times, thus misleading the final similarity score. In fact, while analyzing massive genomes, the number of repeated patterns is very high, particularly in the non-genic regions. Furthermore if we allow mismatches the number of patterns can grows exponentially [13][14][15]. In this paper we will address this problem by controlling the distribution of subwords over the sequences under consideration, so that their contribution will not be overcounted.
Moreover, when comparing genomes it is well known that different evolutionary mechanisms can take place. In this framework, two closely related species are expected to share larger portions of DNA than two distant ones, whereby also other large complements and reversecomplements, or inversions, may occur [16]. In this work we will take into account all these symmetries, in order to define a measure of similarity between whole genomes.

Matching statistics
Among the many distance measures proposed in the literature, which in most cases are dealing with k-mers, an effective and particularly elegant method is the Average Common Subword approach (ACS), introduced by Ulitsky et al. [7]. They use a popular concept in the field of string algorithms, known as matching statistics [17]. In short, given two sequences s 1 and s 2 , where s 1 is the reference sequence, the matching statistic is a vector l such that l[i] is the length of the longest subword starting at position i of s 1 that is also a subword of s 2 , for every possible position i of s 1 (see Table 1).
A popular measure of similarity between strings is the average of this vector. In fact the general form of ACS is: We can notice the similarity with the cross entropy of two probability distributions P and Q: where p(x) log q(x) measures the number of bits needed to code an event x from P if a different coding scheme based on Q is used, averaged over all possible events x.
From the theoretical prospective it can be shown [7] that the ACS approach mimics the cross entropy estimated between two large sequences generated by a finite-state Markov process. In practice, this is closely related to the Kullback-Leibler information divergence, and represents the minimum number of bits needed to describe one string, given the other: D KL (P Q) = H(P, Q) − H(P).  This is perhaps the most frequently used informationtheoretic similarity measure. The advantage of using the matching statistics is that it is not based on fixed-length subwords, but it can capture also variable length matches, in contrast to most methods that are based on fixed sets of k-mers. In fact, with the latter the choice of the parameter k is critical, and every method needs to estimate k from the data under examination, typically using empirical measurements [4].
For this reason ACS proved to be useful for reconstructing whole-genome phylogenies of viruses, bacteria, and eukaryotes, outperforming in most cases the state-of-theart methods [7].

Methods
In this section we propose a distance measure between entire genomes based on the notion of underlying subwords. In order to build a sound similarity measure between genomes, we need first to study the properties of the matching statistics. Our first contribution is the characterization of the subwords that are needed to compute the matching statistics. A second contribution is the selection of these subwords so that the resulting similarity measure does not contain overcounts. Our main idea is to avoid overlaps between selected subwords, more precisely by discarding common subwords occurring in regions covered by other more significant subwords.

Irredundant common subwords
In the literature, the values l[i] used by the ACS approach are called the matching statistics, as described in detail by Gusfield [17]. Our first contribution is to characterize the matching statistics in order to identify which subwords are essentials.
It is well known that the total number of distinct subwords of any length found in a sequence of length n can be at most (n 2 ). Remarkably a notable family of fewer than 2n subwords exist that is maximal in the host sequence, in the sense that it is impossible to extend a word in this class by appending one or more characters to it without losing some of its occurrences [18]. It has been shown that the matching statistics can be derived from this set of maximal subwords [19]. Here we will further tighten this bound by showing that to compute the matching statistics it is enough to consider a subset of the maximal subwords, called irredundant common subwords.
The notion of irredundancy was introduced in [20] and later modified for the problem of protein comparison [21,22]. It proved useful in different contexts from data compression [23] to the identification of transcription factors [24]. In this paper we introduce the concept of irredundant common subwords (i.e., without mismatches/wildcards). This ensures that there exists a close correspondence between the irredundant common subwords and the matching statistics.

Definition 1. (Irredundant/Redundant common subword) A common subword w is irredundant if and only
if at least an occurrence of w in s 1 or s 2 is not covered by other common subwords. A common subword that does not satisfy this condition is called a redundant common subword.
We observe that the number of irredundant common subwords I s 1 ,s 2 is bounded by m + n, where |s 1 | = n and |s 2 | = m, since it is a subset of the set of maximal common subwords (see [19,25] for a more complete treatment of this topic). Proof. To show that the vector l s 1 (i) can be derived from the irredundant common subwords, we define a new vector of scores l w for each subword w, where l w [ j] = |w| − j + 1 represents the length of each suffix j of w, with j = 1, . . . , |w|. Then, for each subword w in I s 1 ,s 2 we superimpose the vector l w on all the occurrences of w in s 1 . For each position i, in s 1 , l s 1 (i) is the maximum value of the scores max w (l w [ j] ), such that k + j = i and k is an occurrence of w.
To complete the proof we have to show that every occurrence of a common subword of s 1 and s 2 is covered by some occurrence of a subword in I s 1 ,s 2 . By definition of irredundant common subword, any occurrence of a subword corresponds to an irredundant common subwords or is covered by some subword in I s 1 ,s 2 . Moreover every irredundant common subword w has at least an occurrence i that is not covered by other subwords. Thus, l s 1 (i) corresponds exactly to |w| and the subword w is necessary to compute the matching statistics. In conclusion, by using the method described above for l s 1 (i), we can compute for each position the length of the maximum common subword starting in that location, which corresponds to the matching statistics.
In summary, the notion of irredundant common subwords is useful to decompose the information provided by the matching statistics into several patterns. Unfortunately these subwords can still overlap in some position. This might lead to an overcount in the matching statistics, in which the same region of the string contributes more than once. Our aim is to remove the possibility of overcount by filtering the most representative common subwords for each region of the sequences s 1 and s 2 , which will also remove all overlaps. http://www.almob.org/content/7/1/34

Underlying subwords
When comparing entire genomes we want to avoid that large non-coding regions, which by nature tend to be highly repetitive, may contribute to the similarity score a multiple number of times, thus misleading the final score. In fact, while analyzing massive genomes, the number of repeated patterns is very high, particularly in the nongenic regions. Therefore we need to filter out part of this information, and select the "best" common subword, by some measure, for each region of the sequences.
In this regard, we must recall the definition of pattern priority and of underlying pattern, adapted from [26] to the case of pairwise sequence comparison. We will take as input the irredundant common subwords and the underlying quorum u = 2, i.e. they must appears at least twice. Let now w and w be two distinct subwords. We say that w has priority over w , or w → w , if and only if either |w| ≥ |w |, or |w| = |w | and the first occurrence of w appears before the first occurrence of w . In this case, every subword can be defined just by its length and one of its starting positions in the sequences, meaning that any set of subwords is totally ordered with respect to the priority rule. We say that an occurrence l of w is tied to an occurrence l of a subword w , if the two occurrences over- Otherwise, we say that l is untied from l .

s 2 is said to be underlying if and only if:
(i) every subword w in U s 1 ,s 2 , called an underlying subword, has at least two occurrences, one in s 1 and the other in s 2 , that are untied from all the untied occurrences of other subwords in U s 1 ,s 2 \ w, and (ii) there does not exist a subword w ∈ I s 1 ,s 2 \ U s 1 ,s 2 such that w has at least two untied occurrences, one per sequence, from all the untied occurrences of subwords in U s 1 ,s 2 .
This subset of I s 1 ,s 2 is composed only by those subwords that rank higher with our priority rule with respect to s 1 . In fact, if there are overlaps between subwords that are in I s 1 ,s 2 , we will select only the subwords with the highest priority. Similarly to the score ACS(s 1 , s 2 ), the set U s 1 ,s 2 is asymmetric and depends on the order of the two sequences; we will address this issue in Section "A distance-like measure based on underlying subwords". As for the underlying patterns [26], one can show that the set of underlying subwords exists, and is unique. As a corollary we know that the untied occurrences of the underlying subwords can be mapped into the sequences s 1 and s 2 without overlaps. Moreover, by definition, the total length of the untied occurrences cannot exceed the length of the sequences. These two properties are crucial when building a similarity measure, because any similarity that is based on these subwords will count the contribution of a region of the sequence only once.

Efficient computation of underlying subwords
To summarize we select the irredundant common subwords that best fit each region of s 1 and s 2 , employing a technique that we call Underlying Approach or, in short, UA. This technique is based on a simple pipeline. We first select the irredundant common subwords and subsequently filter out the subwords that are not underlying. From a different perspective, we start from the smallest set of subwords that captures the matching statistics and remove the overlaps by applying our priority rule. In the following we show how to compute the irredundant common subwords and the matching statistics, and then we present an approach for the selection of the underlying subwords among these subwords. The general structure of the Underlying Approach (UA) is the following: • 1) Compute the set of the irredundant common subwords I s 1 ,s 2 • 2) Rank all subwords in I s 1 ,s 2 according to the priority and initialize U to an empty set. • 3) Iteratively select a subwords p from I s 1 ,s 2 following the ranking. • 4a) If p has at least two untied occurrences: add p to U and update the corresponding regions of (see next) in which p occurs; • 4b) otherwise, discard p and return to (3).

Discovery of the irredundant common subwords
In step (1) we construct the generalized suffix tree T s 1 ,s 2 of s 1 and s 2 . We recall that an occurrence of a subword is (left)right-maximal if it cannot be covered from the (left)right by some other common subword. The first step consists in making a depth-first traversal of all nodes of T s 1 ,s 2 , and coloring each internal node with the colors of its leaves (each color corresponds to an input sequence). In this traversal, for each leaf i of T s 1 ,s 2 , we capture the lowest ancestor of i having both the colors c 1 and c 2 , say the node w. Then, w is a common subword, and i is one of its right-maximal occurrences (in s 1 or in s 2 ); we select all subwords having at least one right-maximal occurrence. The resulting set will be linear in the size of the sequences, that is O(m + n). This is only a superset of the irredundant common subwords, since the occurrences of these subwords could be not left-maximal.
In a second phase, we map the length of each rightmaximal occurrence i into l s 1 (i), and, using Proposition 1, we check which occurrences i have length greater than or equal to the length stored in the location i − 1 (for locations i ≥ 2). These occurrences are also left-maximal, http://www.almob.org/content/7/1/34 since they cannot be covered by a subword appearing at position i − 1. Finally we can retain all subwords that have at least an occurrence that is both right-and left-maximal, i.e, the set of irredundant common subwords I s 1 ,s 2 . Note that, by employing the above technique, we are able to directly discover the irredundant common subwords and the matching statistics l s 1 (i).
The construction of the generalized suffix tree T s 1 ,s 2 and the subsequent extraction of the irredundant common subwords I s 1 ,s 2 can be completed in time and space linear in the size of sequences.

Selection of the underlying subwords
In this section we describe, given the set of the irredundant common subwords I s 1 ,s 2 , how to filter out the subwords that are not underlying, obtaining the set of underlying subwords U s 1 ,s 2 .
The extraction of underlying subwords takes as input the set I s 1 ,s 2 and the tree T s 1 ,s 2 from the previous section. First we need to sort all subwords in I s 1 ,s 2 according to the priority rule (step 2). Then, starting from the top subword, we analyze iteratively all subwords by checking their untied occurrences (step 3). If the subword passes a validity test we select it as underlying (step 4a), otherwise we move on with the next subword (step 4b).
The two key steps of this algorithm are: sorting the subwords (step 2) and checking for their untied occurrences (step 4a).
Step 2 is implemented as follows. For all subwords we retrieve their lengths and first occurrences in s 1 from the tree T s 1 ,s 2 . Then each subword is characterized by its length and the first occurrence. Since these are integers in the range [ 0, n] we can apply radix sort [27], first by length and then by occurrence. This step can be done in linear time.
In order to implement step 4a we need to define the vector of n booleans, representing the locations of s 1 . If [ i] is TRUE, then the location i is covered by some untied occurrence. We also preprocess the input tree and add a link for all nodes v to the closest irredundant ancestor, say prec(v). This can be done by traversing the tree in preorder. During the visit of a the node v if it is not irredundant we transmit to the children prec(v) otherwise if v is irredundant we transmit v. This preprocess can be implemented is linear time and space.
For each subword w in I s 1 ,s 2 we consider the list L w of occurrences to be checked. All L w are initialized in the following way. Every leaf v, that represent a position i, send its value i to the location list of the closest irredundant ancestor using the link prec(v). Again this preprocess takes linear time and space since all positions appear in exactly one location list. We will updated these lists L w only with the occurrences to be checked, i.e. that are not covered by some underlying subword already discovered.
We start analyzing the top subword w and for this case L w is composed by all the occurrences of w.
For Now, to select w , the longest prefix of w with |w | ≤ d, we employ an algorithm proposed by Kopelowitz and Lewenstein [28] for solving the weighted ancestor problem, where weights correspond to the length of words spelled in the path from the root to each node, in case of a suffix tree. In the weighted ancestor problem one preprocesses a weighted tree to support fast predecessor queries on the path from a query node to the root. That is, with a linear preprocessing on a tree of height n, using the above algorithm it is possible to locate any ancestor node w that has a weight less than d in time O(log log n). In our case, the maximum length for an irredundant subword is min{m, n}, thus we can find a suitable ancestor w of w in time O(log log min{m, n}), with O(m+n) preprocessing of the tree T s 1 ,s 2 .
At the end of the process, if the subword w has at least one untied occurrence per sequence, then we mark w as http://www.almob.org/content/7/1/34 underlying subword. Otherwise, all the occurrences of w that are not covered are sent to its ancestors, using the previous procedure.
To analyze the overall complexity we need to compute how many times the same location i is evaluated. Suppose, for example, that i belongs to L w of the subword w. The location i is evaluated again for somew, and inserted into the list Lw, only if [ i] is FALSE and [ i + |w| − 1] is TRUE. Note that the locations not already covered are in the range [ i, i+|w|−d−1], with d > 0. Then, the subword w is the longest prefix of w that is an irredundant common subword and that lives completely in the locations [ i, i + |w| − d − 1]; howeverw may not cover the entire interval. Now, the occurrence i will be evaluated again only if there exists another subword w that overlaps withw, and that has a higher priority with respect tow. The worst case is when w ends exactly at position i+|w|−d−1 and overlaps withw by only one location. Since w must be evaluated beforew, then |w | ≥ |w|. Thus the worst case is when the two subwords have about the same length. In this settings the length of the subwordw can be at most (|w|−d)/2. We can iterate this argument at most O(log |w|) times for the same position i. Therefore any location can be evaluated at most O(log min{m, n}) times. In conclusion, our approach requires O((m + n) log min{m, n} log log min{m, n}) time and O(m + n) space to discover the set of all underlying subwords U s 1 ,s 2 .

Extension to inversions and complements
In this section we discuss the extension of the algorithmic structure discussed above to accommodate also inversion and complement matches.
A simple idea is to concatenate each sequence with its inverse and its complement, while keeping separate the occurrences coming from direct matches, inversions, and complements. In brief, we first definex as the concatenation of a string x with its inverse, followed by its complement, in this exact order. Then, we compute the irredundant common subwords, I s 1 ,ŝ 2 , on the sequences s 1 andŝ 2 . We subsequently select the underlying subwords by ranking all the irredundant common subwords in the set I s 1 ,ŝ 2 . Using the same algorithm described above we compute the set U s 1 ,ŝ 2 , and then we map each subword occurrence to the reference sequences s 1 . This will include also inversions and complements of s 2 that are shared by s 1 . In this way, we can store all the untied occurrences and consider all possible matches for each region of s 1 .
In this framework, we choose to take into account all these symmetries, and thus the experiments presented will use this extended approach. We will also measure the contribution of inversions and complements to our similarity measure.

A distance-like measure based on underlying subwords
In the following we report the basic steps of our distancelike measure. Let us assume that we have computed U s 1 ,s 2 , and the other specular set U s 2 ,s 1 . For every subword w ∈ U s 1 ,s 2 we sum up the score h s 1 s 2 ), where h s 1 w is the number of its untied occurrences in s 1 , similarly to ACS [7]. Then, we average UA(s 1 , s 2 ) over the length of the first sequence, s 1 , yielding This is a similarity score that is large when two sequences are similar, therefore we take its inverse. Moreover, for a fixed sequence s 1 this score can also grow with the length of s 2 , since the probability of having a match in s 1 increases with the length of s 2 . For this reason, we consider the measure log 4 (|s 2 |)/UA(s 1 , s 2 ); we use a base-4 logarithm since DNA sequences have four bases. Another issue with the above formula is the fact that it is not equal to zero for s 1 = s 2 ; thus we subtract the correction term log 4 (|s 1 |)/UA(s 1 , s 1 ), which ensures that this condition is always satisfied. Since U s 1 ,s 1 contains only one subword, the sequence s 1 itself, which trivially has only one untied occurrence in s 1 , this yields to UA(s 1 , s 1 ) = |s 1 |(|s 1 |+1)/(2|s 1 |) = (|s 1 |+1)/2. The following formulas accommodate all of these observations in a symmetrical distance-like measure d UA (s 1 , s 2 ) between the sequences s 1 and s 2 : We can easily see that the correction term rapidly converges to zero as |s 1 | increases. Moreover, we notice that d UA (s 1 , s 2 ) grows as the two sequences s 1 and s 2 diverge. From now we will simply refer to the measure d UA (s 1 , s 2 ) as the Underlying Approach measure, or UA.

Genome datasets and reference taxonomies
We assess the effectiveness of the Underlying Approach on the estimation of whole-genome phylogenies of different organisms. We tested our distance function on three types of datasets: viruses, prokaryotes, and unicellular eukaryotes.
In the first dataset we selected 54 virus isolates of the 2009 human pandemic Influenza A -subtype H1N1, also called the "Swine Flu. " The Influenza A virion has eight segments of viral RNA with different functions. These http://www.almob.org/content/7/1/34 RNAs are directly extracted from infected host cells, and synthesized into complementary DNA by reverse transcription reaction, where a specific gene amplification is performed for each segment [29]. We concatenate these segments according to their conventional order given by the literature [30]; this step, in general, does not affect the final phylogeny computed by our algorithm, and is used to sort subwords by location. The resulting sequences are very similar to each other, and have lengths in the order of 13,200 nucleotides each, accounting for a total of 714,402 b. To compute a reference taxonomic tree, we perform multiple sequence alignment using the ClustalW2 [31] tool a as suggested by many scientific articles on the 2009 Swine Flu [29,30]. Then, we compute the tree using the Dnaml tool from the PHYLIP [32] software package, b which implements the maximum likelihood method for aligned DNA sequences. In Dnaml we used the parameters suggested in [29,30], which consider empirical base frequencies, constant rate variation among sites (with no weights), a transition ratio of 2.0, and best tree search based on proper searching heuristics.
In the second dataset we selected 18 prokaryotic organisms among the species used in [7] for a DNA phylogenomic inference. We chose the species whose phylogenomic tree can be inferred by well-established methods in the literature (see Table 2). The organisms come from both the major prokaryotic domains: Archaea, 8 organisms in total, and Bacteria, 10 organisms in total. The sequences in question have lengths ranging from 0.6 Mbp to 8 Mbp, accounting for a total 48 Mbp. We compute their tree-of-life by using genes that code for the 16S RNA, the main RNA molecule inside the small ribosomal subunit characterizing prokaryotes and widely used to reconstruct their phylogeny; the considered sequences are called 16S rDNA. We can extract a multiple alignment of 16S rDNA sequences of the selected organisms from the Ribosomal Database Project [33]; c our experiments are based on the release 8.1. Next, we perform a maximum likelihood estimation on the aligned set of sequences, employing Dnaml from PHYLIP with standard parameters, in order to compute a reference tree based on the resulting estimation.
In the third dataset we selected five eukaryotic organisms of the protozoan genus Plasmodium whose genomes have been completely sequenced (Table 3). Plasmodium are unicellular eukaryotic parasites best known as the etiological agents of malaria infectious disease. The sequences have lengths ranging from 18 Mbp to 24 Mbp, accounting for a total 106 Mbp. We used as reference tree the taxonomy computed by Martinsen et al. [34], as suggested by the Tree of Life Project.

Whole-genome phylogeny reconstruction
We exploited the above datasets to compare our method, the Underlying Approach (UA), with other efficient stateof-the-art approaches in the whole-genome phylogeny reconstruction challenge: ACS [7], FFP [4]  We reconstruct the phylogenomic trees from the distance matrices using the Neighbor-joining method as implemented in the PHYLIP package. We compare the resulting topologies with the respective reference trees using the symmetric difference of Robinson and Foulds (R-F) and the triplet distance. For two unrooted binary trees with n ≥ 3 leaves, the R-F score is in the range [0, 2n − 6]. A score equal to 0 means that the two trees are isomorphic, while 2n − 6 means that all non-trivial bipartitions are different. The R-F difference between two or more trees can be computed using the TreeDist tool from the PHYLIP package. We ran FFP and FFP RY for different values of k (the fixed subword length) as suggested by [4], retaining the best results in agreement with the reference trees. Table 4 compares our method with the other state-of-the-art approaches, by showing the R-F difference with respect to the reference taxonomy tree.
Our method, UA, achieves good performance in every test considering the R-F difference with the reference taxonomy tree, and very good performance if we further analyze the resulting phylogenies, as in Figures 1, 2, and 3. For every dataset the best results are shown in bold. We can observe that UA is constantly the best performing method, and that this advantage becomes more evident for large dataset, where sequences share large parts, such as the Influenza A (H1N1) viruses.
The Robinson and Foulds distance is a standard method to evaluate topological discordance between trees. However when dealing with large trees it is known that small variations can generate very large R-F scores (typically, already for n ∼ 10). For this reason we conducted a second series of experiments using the triplet distance [35]. The triplet distance is a more refined measure that does not suffer this problem. Moreover, to better compare all taxonomies, we report the triplet distance between all trees. Tables 5, 6 and 7 show the triplet distance between all trees for all datasets. This more refined measure confirms the applicability of UA with respect to FFP and ACS.
In more detail, Figure 1 shows that our approach can distinguish the two main clades of the 2009 Influenza A-H1N1 (in green and red), which have been outlined in [30]. The origin of the flu could reside in the Mexican isolate (Mexico/4108, in green), from which all other green isolates may have ensued. Two sub-clades for the U.S. states of California and Texas are highlighted within the red clade, most probably corresponding to the first major evolutions of the viral disease.
Similar results are obtained for the second dataset, as shown in Figure 2. UA can easily distinguish the Archaea domain, in red, from the Bacteria domain, in green, and also other sub-clades with respect to the reference tree (these sub-clades are highlighted in the figure with different colors). The organisms in black do not form a clade with other organisms in the reference tree. For the third dataset (Figure 3), the whole-genome phylogeny of the genus Plasmodium generated by UA corresponds exactly to the taxonomy found in the literature. The accuracy results are promising, but we believe that of equal interest are the patterns used for the classification. Our approach, by construction, uses only a very small number of patterns. For this reason we report in Table 8 some statistics for the underlying subwords selected, averaged over all experiments. We can notice that the number of irredundant patterns is in general smaller than the length of the genomes, and this is a first   form of information filtering. Moreover we can observe that only a few underlying subwords are selected on average among the irredundant common subwords. This number is always very small when compared with all possible irredundant subwords, and much smaller than the length of the sequences. Similar considerations can be drawn for the underlying subwords length. On average they can be very long, especially with respect to FFP that uses only k-mers with k in the range [ 5,10]. Furthermore, each underlying subword occurs only a few times per sequence, and in general about one occurrence per sequence. Removing the highfrequency subwords, we can notice that the underlying subwords typically have length ≥ log 4 min{m, n}, and in the case of viruses they can be very large, capturing more information than FFP. The longest underlying subwords appear in the virus dataset, and they are on the order of a thousand bases. We checked if these subwords may have some biological meaning and we found that in some cases they correspond to whole viral segments that are shared between two genomes. This confirms that, in some cases, the underlying subwords used for classification can capture some biological insight.
Another interesting aspect is the contribution of inversions and complements in our similarity measure, with respect to the classical notion of match. We compute the average number of occurrences used in our scoring function that is due to inversions and complements. The contribution of inversions and complements is about 28-33% and 19-20%, respectively. This fact may be due to the nature of the sequences considered, but we believe that this topic deserves more attention.

Conclusion
In conclusion, we have shown that the underlying subwords can be used for the reconstruction of phylogenetic trees. Preliminary experiments have shown very good performance in the identification of major clusters for viruses, prokaryotes, and unicellular eukaryotes. An important observation that distinguishes our methods from the others is that only a small number of