Fast Approximation of Frequent $k$-mers and Applications to Metagenomics
Leonardo Pellegrina Cinzia Pizzi Fabio Vandin















L. Pellegrina et al.



Department of Information Engineering, University of Padova,
Padova (Italy)
{pellegri,cinzia.pizzi}@dei.unipd.it;

fabio.vandin@unipd.it


Fast Approximation of Frequent k-mers
 and Applications to MetagenomicsThis work is supported, in part, by the University of Padova grants SID2017 and STARS: Algorithms for Inferential Data Mining.

    Leonardo Pellegrina
Cinzia PizziFabio Vandin
    December 30, 2023
======================================================================================================================================================================================================




Estimating the abundances of all k-mers in a set of biological sequences is a fundamental and challenging problem with many applications in biological analysis. While several methods have been designed for the exact or approximate solution of this problem, they all require to process the entire dataset, that can be extremely expensive for high-throughput sequencing datasets. While in some applications it is crucial to estimate all k-mers and their abundances, in other situations reporting only frequent k-mers, that appear with relatively high frequency in a dataset, may suffice. This is the case, for example, in the computation of k-mers' abundance-based distances among datasets of reads, commonly used in metagenomic analyses.

In this work, we develop, analyze, and test, a sampling-based approach, called , to approximate the frequent k-mers and their frequencies in a high-throughput sequencing dataset while providing rigorous guarantees on the quality of the approximation.  employs an advanced sampling scheme and we show how the characterization of the VC dimension, a core concept from statistical learning theory, of a properly defined set of functions leads to practical bounds on the sample size required for a rigorous approximation. Our experimental evaluation shows that  allows to rigorously approximate frequent k-mers by processing only a fraction of a dataset and that the frequencies estimated by  lead to accurate estimates of k-mer based distances between high-throughput sequencing datasets. Overall,  is an efficient and rigorous tool to estimate k-mers abundances providing significant speed-ups in the analysis of large sequencing datasets.








 

§ INTRODUCTION

 

The analysis of substrings of length k, called k-mers, is ubiquitous in biological sequence analysis and is among the first steps of processing pipelines for a wide spectrum of applications, including: de novo assembly <cit.>, error correction <cit.>, repeat detection <cit.>, genome comparison <cit.>, digital normalization <cit.>, RNA-seq quantification <cit.>, metagenomic reads classification <cit.> and binning <cit.>, fast search-by-sequence over large high-throughput sequencing repositories <cit.>. A fundamental task in k-mer analysis is to compute the frequency of all k-mers, with the goal to distinguish frequent k-mers from infrequent k-mers <cit.>. For example, this task is relevant in the analysis of high-throughput sequencing data, since infrequent k-mers are often assumed to result from sequencing errors. For several applications, the computation of k-mers frequencies is among the most computationally demanding steps of the analysis.
 
Many algorithms have been proposed for computing the exact frequency of all k-mers, such as Jellyfish <cit.>, DSK <cit.>, KMC 3 <cit.> and Squeakr-exact <cit.>. These methods typically perform a linear scan of the sequence to analyze, and use a combination of parallelism and efficient data structures (such as Bloom filters and Hash tables) to maintain membership and counting information associated to all k-mers.

Since the computation of exact k-mer frequencies is computationally demanding, in particular for large sequence analysis or for high-throughput sequence datasets, recent methods have focused on providing approximate solution to the problem,
improving the time and memory requirements. KmerStream <cit.>, khmer <cit.>, Kmerlight <cit.> and ntCard <cit.> proposed streaming approaches for the approximation of the k-mer frequencies histogram. Of these, only Kmerlight and ntCard provide analytical bounds on their accuracy guarantee. KmerGenie <cit.> performs a linear scan of the input to compute the frequencies of a (random) subset of the k-mers that appear in the input, and uses these frequencies to approximate the abundance histogram. The recently proposed Squeakr <cit.> relies on a probabilistic data structure to approximate the counts of individual k-mers. Turtle <cit.> focuses on finding k-mers that appear at least twice in the dataset, but still processes all the k-mer occurrences in the input dataset, as all the other aforementioned methods do.

All the methods cited above try to estimate the frequency of all k-mers or of all k-mers that appear at least few times (e.g., twice) in the dataset. While this is crucial in some applications (e.g., in genome assembly k-mers that occur exactly once often represents sequencing errors and it is therefore important to estimate the count of all observed k-mers), in other applications this is less justified. For example, in the comparison of high-throughput sequencing metagenomic datasets, abundance-based distances or dissimilarities (e.g., the Bray-Curtis dissimilarity) between k-mer counts of two datasets are often used <cit.> to assess the distance between the corresponding datasets.  In contrast to presence-based distances <cit.> (e.g., Jaccard distance), abundance-based distances take into account the frequency of each k-mer, with frequent k-mers contributing more to the distance than k-mers that appear with low frequency, but still more than a handful of times, in the dataset. 
Thus, two natural questions are [(i)] 
  * whether the results obtained considering all k-mers can be estimated by considering the abundances of frequent k-mers  only, and 
  * if the abundances of frequent k-mers can be computed more efficiently than the counts of all k-mers. Recently, preliminary work <cit.> has shown that, for the cosine distance and k=12, the answer to the first question is positive, and in Section <ref> we show that this indeed the case for larger values of k and other abundance-based distances as well as presence-based distances (e.g., the Jaccard distance). To the best of our knowledge, the second question is hitherto unexplored.  In addition, considering only frequent k-mers allows to focus on the most reliable information in a metagenomic dataset, since a high stochastic variability in low frequency k-mers is to be expected due to the sampling process inherent in sequencing.

A natural approach to reduce time and memory requirements for frequency estimation problems is to process only a portion of the data, for example by sampling some parts of a dataset. Sampling approaches are appealing because infrequent k-mers naturally tend to appear with lower probability in a sample, allowing to directly focus on frequent k-mers in subsequent steps. However, major challenges in sampling approaches are [(i)]

  * to  provide rigorous guarantees relating the results obtained by processing the sample and the results that would be obtained from the whole dataset, and

  * to provide effective bounds on the size of the sample required to achieve such guarantees.

The application of sampling to k-mers is even more challenging than in other scenarios since, for values of k in the typical range of interest to applications (e.g., 20-60), even the most frequent k-mers have relatively low frequency in the data.
To the best of our knowledge, no approach based on sampling a portion of the input dataset has been proposed to approximate frequent k-mers and their frequencies while providing rigorous guarantees.




Our Contribution. We study the problem of approximating frequent k-mers, i.e., k-mers that appear with frequency above a user-defined threshold θ in a high-throughput sequencing dataset. In these regards, our contributions are fourfold. First, we define a rigorous definition of approximation, governed by an accuracy parameter ε. Second, we propose a new method, Sampling Algorithm for K-mErs approxIMAtion (), to obtain an approximation to the set of frequent k-mers using sampling.  is based on a sampling scheme that goes beyond naïve sampling of k-mers and allows to estimate low frequency k-mers considering only a fraction of all k-mers occurrences in the dataset. Third, we provide analytical bounds to the sample size needed to obtain rigorous guarantees on the accuracy of the estimated k-mer frequencies, with respect to the ones measured on the entire dataset. Our bounds are based on the notion of VC dimension, a fundamental concept from statistical learning theory. To our knowledge, ours is the first method that applies concepts from statistical learning to provide a rigorous approximation of the k-mers frequencies. Fourth, we use  to extract frequent k-mers from metagenomic datasets from the Human Microbiome Project (HMP) and to approximate abundance-based and presence-based distances among such datasets, showing that  allows to accurately estimate such distances by analyzing only a fraction of the entire dataset, resulting in a significant speed-up.


Our approach is orthogonal to previous work: any exact or approximate algorithm can be applied to the sample extracted by , that can therefore be used before applying previously proposed methods, thus reducing their computational requirements while providing rigorous guarantees on the results w.r.t. to the entire dataset. 
While we present our methodology in the case of finding frequent k-mers from a set of sequences representing a high-throughput sequencing dataset of short reads, our results can be applied to datasets of long reads and to whole-genome sequences as well.



§ PRELIMINARIES



Let a dataset  be a bag of n reads ={r_0 , … , r_n-1}, where
each read r_i, 0 ≤ i ≤ n-1,  is a string of length n_i from an alphabet Σ of cardinality |Σ|=σ. For j ∈{0,…,n_i-1}, let r_i[j] be the j-th character of r_i. For a given integer k ≤min_i{n_i: r_i ∈}, we define a k-mer A as a string of length k from Σ, that is A ∈Σ^k. We say that a k-mer A appears in r_i at position j ∈{0,…, n_i-k} if r_i[j+h]=A[h], ∀ h ∈{0,…,k-1}.  For every i, 0≤ i ≤ n-1, and every j ∈{0,… , n_i - k}, we define the indicator function ϕ_r_i,A(j) that is 1 if the k-mer A appears in r_i at position j, while  ϕ_r_i,A(j)=0 otherwise. The total number of k-mers in  is t_,k = ∑_i=0^n-1 (n_i - k + 1). We define the support  o_(A) of a k-mer A as the number of distinct positions in  where A appears: o_(A) = ∑_i=0^n-1∑_j=0^n_i-kϕ_r_i,A(j). We define the frequency f_(A) of A in  as the ratio between the number of distinct positions where A appears in  and the total number of k-mers in : f_(A) = o_(A) / t_,k.



 §.§ Frequent k-mers and Approximations


We are interested in obtaining the set (,k,θ) of frequent k-mers in a dataset  with respect to a minimum frequency threshold θ, defined as follows.


Given a dataset , an integer k>0, and a frequency threshold θ∈ (0,1], the set (,k,θ) of Frequent k-Mers in  w.r.t. θ is the collection of all k-mers with frequency at least θ in  and of their corresponding frequencies in :

    (,k,θ) = { (A,f_(A)) : f_(A) ≥θ}.



(,k,θ) can be computed with a single scan of all the k-mers occurrences in  maintaining the k-mers supports in an appropriate data structure; however, when  is extremely large and k is not small, the exact computation of (,k,θ) is extremely demanding in terms of time and memory, since the number of k-mers grows exponentially with k. In this case, a fast to compute approximation of the set (,k,θ) may be preferable, provided it ensures rigorous guarantees on its quality. In this work, we focus on the following approximation. 


Given a dataset , an integer k>0, a frequency threshold θ∈ (0,1], and a constant ε∈ (0,θ), an ε-approximation of (,k,θ) is a collection C = { (A,f_A) : f_A ∈ (0, 1]  } such that:


  * for any  (A,f_(A)) ∈(,k,θ)  there is a pair (A,f_A) ∈ C;

  * for any (A,f_A) ∈ C it holds that f_(A) ≥θ - ε;

  * for any (A,f_A) ∈ C it holds that |f_(A) - f_A| ≤ε / 2.



The definition  above guarantees that every frequent k-mer of  is in the approximation and that no k-mer with frequency <θ - ε is in the approximation. The third condition guarantees that the estimated frequency f_A of A in the approximation is close (i.e, within ε/2) to the frequency f_(A) of A in . It is easy to show that obtaining a ε-approximation of (,k,θ) with absolute certainty requires to process all k-mers in .



 §.§ Simple Sampling-Based Algorithms and Bounds



We aim to provide an approximation to (,k,θ) with sampling, by processing only randomly selected portions of . The simplest sampling scheme is the one in which a random sample is a bag P of m positions taken uniformly at random, with replacement, from the set P_,k = {(i,j):  i ∈ [0,n-1], j ∈ [0,n_i-k] }  (note that |P_,k| = t_,k) of all positions where k-mers occurs in the dataset , corresponding to m 
occurrences of k-mers (with repetitions) taken uniformly at random. Given such sample P, an integer k>0, and a minimum frequency threshold θ∈ (0,1] one can define the set of frequent k-mers (and their frequencies) in the sample P as
(P,k,θ) = { (A,f_P(A)): f_P(A) ≥θ}, where f_P(A) is the frequency of k-mer A in the sample.

Obtaining a ε-approximation from a random sample with absolute certainty is impossible, thus we focus on obtaining a ε-approximation with probability 1-δ>0, where δ∈ (0,1) is a confidence parameter, whose value is provided by the user. 
Intuitively, the set (,k,θ) of frequent k-mers  is well approximated by the set of frequent k-mers in a random sample P when P is sufficiently large. One natural question regards how many samples are needed to obtain the desired ε-approximation. 
By using Hoeffding's inequality <cit.> to bound the deviation of the frequency of a k-mer A in the sample from f_(A) and a union bound on the maximum number σ^k of k-mers, where σ = |Σ|, we have the following result that provides a first such bound, and a corresponding first algorithm to obtain a ε-approximation to (,k,θ).  (Due to space constraints proofs are omitted and will be provided in the full version of this extended abstract.)



Consider a sample P of size m of . If
m ≥2/ε^2( ln( 2 σ^k)+ln(1/δ) ) for fixed ε∈ (0,θ), δ∈ (0,1),
then, with probability ≥ 1- δ, (P,k,θ-ε/2) is a ε-approximation of (,k,θ). 


In addition, by using known results in statistical learning theory <cit.> relating the VC dimension (see Section <ref> for its definition) of a family of functions and a novelly derived bound on the family of functions {f_(A)}, we obtain  the following improved bound and algorithm. (The derivation will be provided in the full version.)




Let P be a sample of size m of . For fixed ε∈(0,θ), δ∈ (0,1), if 
m ≥2/ε^2(1 + ln(1/δ))
then (P,k,θ - ε / 2) is an ε-approximation for (,k,θ) with probability ≥ 1 - δ.





§ ADVANCED AND PRACTICAL BOUNDS AND ALGORITHMS FOR K-MER APPROXIMATIONS



While the bound of Proposition <ref> significantly improves the simple bounds of Section <ref>, since the factor ln(2 σ^k) has been reduced to 1, it still has an inverse quadratic dependency with respect to the accuracy parameter ε, that is problematic when the quantities to estimate are small. In these cases, one needs a small ε to produce a meaningful approximation (since ε < θ), and the inverse quadratic dependence of the sample size from ε often results in a sample size larger than the entire input, defeating the purpose of sampling. The case of k-mers is particularly challenging, since the sum ∑_A ∈Σ^k f_(A) of all k-mers frequencies is exactly 1. Therefore the higher the number of distinct k-mers appearing in the input, the lower their frequencies will be, with the consequence that θ (and therefore ε) typically needs to be set to a very low value. For example, a typical dataset from the Human Microbiome Project (HMP) has n ≈ 10^8 reads of (average) length ≈ 100: therefore if we are interested in k-mers for k=31, by setting δ=0.05  the bound of Section <ref> gives ε≈ 10^-5, that is only k-mers with frequency ≥ 10^-5 could be reliably reported by sampling. However, in datasets we considered, no or a very small number (≤ 30) of k-mers have frequency ≥ 10^-5,
therefore according to the result from Section <ref> we cannot obtain a meaningful approximation of k-mers and their frequencies.
In the remaining of this section we develop more refined sampling schemes and estimation techniques leading to a practical sampling-based algorithm.



 §.§ Sampling Bags of Positions and VC dimension Bound.



We propose a method to provide an efficiently computable approximation to (,k,θ) when the minimum frequency θ is low, 
by properly defining samples so that any k-mer A will appear in a sample with probability higher than f_(A), thus lessening the the dependence of the sample size from 1/ε^2.
For this to be achievable, we need to relax the notion of approximation defined in Section <ref>. In particular, the guarantees, provided by our method, in such relaxed approximation are that all k-mers with frequency above θ', with θ' slightly higher than θ, are reported in output, and that no k-mer having frequency below θ - ε is reported in output.  (See Proposition <ref> for the definition of θ'.)
 Our experiments show that the fraction of k-mers having frequency ∈ [θ , θ') which are non reported is very small. Our method works by sampling bags of positions instead than single positions. In particular, an element of the sample is now a set of ℓ positions chosen independently at random from the set P_,k of all positions.
 
Let I_ℓ = {(i_1,j_1) , (i_2,j_2) , … , (i_ℓ,j_ℓ) } be a bag of ℓ positions for k-mers in , chosen uniformly at random from the set P_,k.
We define the indicator functions ϕ̂_A(I_ℓ) that, for a given bag I_ℓ of ℓ positions, is equal to 1 if k-mer A appears in at least one of the ℓ positions in I_ℓ and is equal to 0 otherwise. That is
ϕ̂_A(I_ℓ)  = min{ 1 , ∑_(i,j) ∈ I_ℓϕ_r_i,A(j) }.
We define the ℓ-positions sample P_ℓ as a bag of m bags {I_ℓ,0,I_ℓ,1,…,I_ℓ,m-1}, where each I_ℓ,j, 0 ≤ j≤ m-1 is a bag of ℓ positions, sampled independently, and


    f̂_P_ℓ(A) = 1/m∑_I_ℓ,i∈ P_ℓϕ̂_A(I_ℓ,i) /ℓ.


Intuitively, f̂_P_ℓ(A) is the biased version of the unbiased estimator 
f_P_ℓ(A) = 1/m∑_I_ℓ,i∈ P_ℓ∑_(i,j) ∈ I_ℓ,iϕ_r_i,A(j) /ℓ
of f_(A), where the bias arises from considering a value of 1 every time .

In our analysis we use the Vapnik-Chervonenkis (VC) dimension <cit.>, a statistical learning concept that measures the expressivity of a family of binary functions. We define a range space Q as a pair Q = (X,R_X) where X is a finite or infinite set and R_X is a finite or infinite family of subsets of X. The members of R_X are called ranges. Given D ⊂ X, the projection of R_X on D is defined as proj_R_X(D) = { r ∩ D : r ∈ R_X }. We say that D is shattered by R_X if proj_R_X(D) = 2^|D|. The VC dimension of Q, denoted as VC(Q), is the maximum cardinality of a subset of X shattered by R_X. If there are arbitrary large shattered subsets of X shattered by R_X, then VC(Q)=∞.

A finite bound on the VC dimension of a range space Q implies a bound on the number of random samples required to obtain a good approximation of its ranges, defined as follows.



Let Q=(X,R_X) be a range space and let D be a finite subset of X. For ε∈ (0,1], a subset B of D is an ε-approximation of D if for all r ∈ R_X we have:

||D ∩ r|/|D| - |B ∩ r|/|B|| ≤ε/2.



The following result <cit.> relates ε and the probability that a random sample of size m is an ε-approximation for a range space of VC dimension at most v.



There is an absolute positive constant c such that if (X,R_X) is a range-space of VC dimension at most v, D is a finite subset of X, and 0 < ε, δ < 1, then a random subset B ⊂ D of cardinality m with
m ≥4c/ε^2(v + ln(1/δ))
is a ε-approximation of D with probability at least 1 - δ.


The universal constant c has been experimentally estimated to be at most 0.5 <cit.>. 


We now prove an upper bound to the VC dimension VC(Q) of the range space Q associated to the class of functions ϕ̂_A that grows sub-linearly with respect to ℓ. To this aim, we first define the range space associated to bags of ℓ positions of k-mers.


 

Let  be a dataset of n reads and let k and ℓ be two integers ≥ 1. We define Q = (X_,k,ℓ,R_,k,ℓ) to be the following range space:


  * X_,k,ℓ is the set of all bags of ℓ positions of k-mers in , that is the set of all possible subsets, with repetitions, of size ℓ from  from P_,k;

  * R_,k,ℓ = {P_, ℓ(A) | A ∈Σ^k } is the family of sets of starting positions of k-mers, such that for each k-mer A, the set P_, ℓ(A) is the set of all bags of ℓ starting positions in  where A appears at least once.


 
We prove the following results on the VC dimension of the above range space.
 


Let Q the range space from Definition <ref>. Then: VC(Q) ≤⌊log_2(ℓ) ⌋ + 1.


Using the result above, we prove the following.



Let ℓ≥ 1 be an integer and P_ℓ be a bag of m bags of ℓ positions of  with

    m ≥2/(ℓε)^2(⌊log_2min( 2 ℓ , σ^k ) ⌋ + ln(1/δ)).

Then, with probability at least 1 - δ:


  * for any k-mer A ∈(,k,θ) such that  f_(A) ≥θ'= 1 - (1-ℓθ)^1/ℓ it holds f̂_P_ℓ(A)  ≥θ - ε / 2;

  * for any k-mer A with f̂_P_ℓ(A) ≥θ - ε / 2 it holds f_(A) ≥θ - ε;

  * for any k-mer A ∈(,k,θ) it holds  f_(A) ≥f̂_P_ℓ(A) - ε / 2;

  * for any k-mer A with f̂_P_ℓ(A) - ε/2 ≥ 0, it holds f_(A) ≥ 1- ( 1 - ℓ (f̂_P_ℓ(A) - ε/2))^1 / ℓ;

  * for any k-mer A with ℓ (f̂_P_ℓ(A) + ε/2) ≤ 1 it holds f_(A) ≤ 1- ( 1 - ℓ (f̂_P_ℓ(A) + ε/2))^1 / ℓ.



Note that from Proposition <ref> the set {(A,f_P_ℓ(A)): f̂_P_ℓ(A)≥θ - ε / 2 } is almost a ε-approximation to (,k,θ): in particular, there may be k-mers A for which [f̂_P_ℓ(A)]=  (1 - ( 1 - f_(A) )^ℓ) / ℓ < θ while f_(A) = [f_P_ℓ(A)] ≥θ and such that for the given sample P_ℓ we have f̂_P_ℓ(A) ≈[f̂_P_ℓ(A)] - ε/2. 
While this can happen, we can limit the probability of this happening by appropriately choosing ℓ, and still enjoy the reduction in sample size of the order of log_2 ℓ/ℓ^2 w.r.t. Proposition <ref> obtained by considering bags of bags of ℓ positions.
In particular, this result allows the user to set θ, ε, δ, and ℓ to effectively find, with probability at least 1 - δ, all frequent k-mers A for which f_(A) ≥θ' and do not report any k-mer with frequency below θ - ε, while still being able to report in output almost all k-mers with frequency ∈ [θ , θ'). Our experimental analysis (Section <ref>) shows that in practice choosing ℓ close from below to 1/θ is very effective to obtain such result.
Then, the third, fourth, and fifth guarantees from Proposition <ref> state that we can use the biased estimates f̂_P_ℓ(A) to derive guaranteed upper and lower bounds to f_(A) that will be much tighter than the one obtained using the bounds of Section <ref>. We will show how to obtain further improved upper and lower bounds to f_(A) in Section <ref>. Such lower bounds ℓ b_A can be used, for example, to prove that the set {(A,f_P_ℓ(A)): ℓ b_A ≥θ - ε} enjoys the same last four guarantees 
from Proposition <ref> while the first one holds for a θ' < 1 - (1-ℓθ)^1/ℓ; therefore, when false negatives are problematic, the set {(A,f_P_ℓ(A)): ℓ b_A ≥θ - ε} can be used to obtain a different approximation of (,k,θ) with fewer false negatives.




 §.§ : An Efficient Algorithm to Approximate Frequent k-mers


We now present our Sampling Algorithm for K-mErs approxIMAtion (), that builds on Proposition <ref> and efficiently samples  a bag P_ℓ of bags of ℓ-positions from  to obtain an approximation of the set (,k,θ) with probability 1 - δ, where δ is a parameter provided by the user.




 is described in Algorithm <ref>.
 performs a pass on the stream of k-mers appearing in , and for each position in the stream it samples the number a of times that the position appears in the sample P_ℓ independently at random from the Poisson distribution Poisson(λ) of parameter λ = m ℓ/t_,k.  stores such values in a counting structure T (lines 3-7) that keeps, for each k-mer A, the total number of occurrences of A in the sample P_ℓ. 
(Note that t_,k, that can be computed with a very quick linear scan of the dataset, where n_i is computed for every r_i ∈ without extracting and processing (e.g., inserting or updating information for) k-mers; in alternative a lower bound to t_,k can be used, simply resulting in a number of samples higher than needed). Then, such occurrences are partitioned into the m bags I_ℓ,0,…,I_ℓ,m-1  (line 9); this can be efficiently implemented by assigning each occurrence to a random bag while keeping the difference between the final size of the bags ≤ 1. For each k-mer A appearing at least once in the sample (line 10), the unbiased estimate f_A is computed as the number T[A] of occurrences of A in the sample P_ℓ (line 11) divided by the total number of positions in the sample, while the  biased estimate f̂_A is computed
as the number |𝒫_A| of distinct bags  of P_ℓ where A appears at least once divided by the number m of bags (lines 12-13).
Then  flags A as frequent if f̂_A ≥θ - ε/2 (line 14) and, in this case, the couple (A,f_A) is added to the output set 𝒪 (line 14), since f_A is the best (and unbiased) estimate to f_(A). Note that bags for different values of ℓ (on the same sampled positions) can be obtained by maintaining a table T_ℓ and a set 𝒫_A,ℓ for each value ℓ  of interest.


Note that  does not sample m bags of exactly ℓ positions each, since the number of occurrences of each position in  in the sample P_ℓ is sampled independently from a Poisson distribution, even if the expected number of total occurrences sampled from the algorithm is m ℓ. However, the independent Poisson distributions used by  provide an accurate approximation of the random sampling of exactly  m ℓ positions used in the analysis of Section <ref>.  In particular, this holds when one focuses on the events of interests for our approximation of Section <ref> (e.g., the event “there exists a k-mer A such that |[f̂_P_ℓ(A)] - f̂_P_ℓ(A)| > ε / 2"). In fact, a simple adaptation of a known result (Corollary 5.11 of <cit.>) on the relation between sampling with replacement and the use of independent Poisson distributions gives the following.

Let E be an event whose probability is either monotonically increasing or monotonically decreasing in the number of sampled positions. If E has probability p when the independent Poisson distributions are used, then E has probability at most 2p when the sampling with replacement is used.

As a simple corollary, the output 𝒪 features the guarantees of Proposition <ref> with probability ≥ 1-  δ', with δ'=2 δ.



 §.§ Improved Lower and Upper Bounds to k-mers Frequencies


Note that Proposition <ref> guarantees that we can obtain upper and lower bounds to f_(A) for every A ∈(,k,θ) from the sample of bags of ℓ positions. These bounds are meaningful only in specific ranges of the frequencies; for example, the lower bound  from the third guarantee in Proposition <ref> is meaningful when the frequency of A is fairly low, i.e f_(A) ≈ 1 / ℓ, while for very frequent k-mers they could be a multiplicative factor 1/ℓ away from than the correct value. For example, if a k-mer is very frequent and appears in all bags of ℓ k-mers in a sample S, its corresponding lower bound is still only 1 / ℓ - ε/2.

However, Proposition <ref> can be  generalized to obtain tighter upper and lower bounds to the frequency of all k-mers. For given ℓ, ε, and δ, let m as given in Proposition <ref>. Note that the total number of k-mer's positions in the sample P_ℓ is mℓ.
Let ℒ be a set of integer values ℒ = {ℓ_i } with ℓ_i ∈ [1,mℓ], ∀ i=0,…,|ℒ|-1. Now, for every ℓ_i ∈ℒ, we can partition the same mℓ k-mers that are in P_ℓ into m_i=mℓ/ℓ_i partitions having size ℓ_i. Let P_ℓ_i be such a random partition of such positions into m_i bags of ℓ_i positions each. Note that each P_ℓ_i is a “valid” sample (i.e., a sample of independent bags of positions, each obtained by uniform sampling with replacement) for Proposition <ref>, even if the P_ℓ_i's are not independent.  From each P_ℓ_i, we define a maximum deviation ε_i from Proposition <ref> as ε_i  = 1/ℓ_i√(2/m_i( ⌊log_2 (min (2ℓ_i , σ^k)) ⌋ + ln(|ℒ| / δ) )).
We have the following result.



With probability at least 1-δ, for all k-mers A simultaneously and for all the random partitions induced by ℒ it holds 


  * f_(A) ≥max{f̂_P_ℓ_i(A) - ε_i/2: i=0,…,|ℒ|-1 };

  * f_(A) ≥max{1- ( 1 - ℓ (f̂_P_ℓ_i(A) - ε_i/2))^1 / ℓ: i=0,…,|ℒ|-1 and f̂_P_ℓ(A) - ε_i/2 ≥ 0 };

  * f_(A) ≤min{1- ( 1 - ℓ (f̂_P_ℓ_i(A) + ε_i/2))^1 / ℓ: i=0,…,|ℒ|-1 and f̂_P_ℓ(A) + ε_i/2 ≤ 1/ℓ}.





In our experiments, we use ℒ = {ℓ_i } with ℓ_i = ℓ / 2^i, ∀ i ∈ [0 , ⌊log_2ℓ⌋ - 1]; in this case, note that P_ℓ_0=P_ℓ. Using this scheme, we can compute upper and lower bounds for k-mers having frequencies of many different orders of magnitude, but any (application dependent) distribution can be specified by the user. These upper and lower bounds can be used to obtain different approximations of (,k,θ) with different guarantees. For example, by reporting all k-mers (and their frequencies) that have an upper bound ≥θ, we have an approximation that guarantees that all k-mers A with f_(A) ≥θ are in the approximation.




§ EXPERIMENTAL RESULTS


In this section we present the results of our experimental evaluation for . 
Section <ref> describes the datasets, our implementation for [Available at https://github.com/VandinLab/SAKEIMA], and the baseline for comparisons.
In Section <ref>, we report the results for computing the approximation of the frequent k-mers using .
Section <ref> reports the results of using our approximation to compute abundance-based and presence-based distances between metagenomic datasets. 



 §.§ Datasets and Implementation



We considered six datasets from the Human Microbiome Project (HMP)[https://hmpdacc.org/HMASM], one of the largest publicly available collection of metagenomic datasets from high-throughput sequencing. In particular, we selected the three largest datasets of  and the three largest of  (Table <ref>). 
These datasets constitute the most challenging instances, due to their size, and provide a test case with different degrees of similarities among datasets.




We implemented  in  as a modification of Jellyfish <cit.> (the version we used is 2.2.10[https://github.com/gmarcais/Jellyfish]), a very popular and efficient algorithm for exact k-mer counting. Doing so, our algorithm enjoys the succinct counting data structure provided by Jellyfish publicly available implementation. We remark that our sampling-based approach can be used in combination with any other highly tuned method available for exact, approximate, and parallel k-mer counting. For this reason, we only compare  with the exact counting performed by Jellyfish, since they share the underlying characteristics, allowing us to evaluate the impact of  sampling strategy.


We did not include the time to compute t_,k in our experiments since it was always negligible (i.e., less than 2 minutes) w.r.t. the time for counting k-mers.


For the computation of the abundance-based distances
from the k-mer counts of two dataset, we implemented in  a simple algorithm  that loads the counts of one dataset in main memory and then performs one pass on the counts of the other dataset, producing the distances in output.
We executed all our experiments on the same machine with 512 GB of RAM and a 2.30 GHz Intel Xeon CPU, compiling both implementations with .  can be used in combination with more efficient algorithms and implementations for the computation of these (and other) distances <cit.>, resulting in speed-ups analogous to the ones we present below.
For all the experiments of , given θ and a dataset , we fixed the parameters δ = 0.1, ε = θ - 2/t_,k, m = 100, and we fix ℓ to the minimum value satisfying the ε-approximation. For all the experiments we have ℓ close from below to 1/θ. For all the metrics we considered, we report the results for one random run.



 §.§  Approximation of the Frequent k-mers



We fixed k=31, and we compared  with the exact counting of all k-mers (from Jellyfish) in terms of:
[(i)]

  * running time[Every instance of  and Jellyfish was executed with 1 worker, i.e., sequentially. 
Note that the Poisson approximation employed by  allows multiple workers to independently process the input k-mers, therefore  can be used in a parallel scenario. We will investigate the impact of parallelism in the extended version of this work.
], including, for both algorithms, the time required to write the output on disk;

  * memory requirement.

We also assessed the accuracy of the output of .



Figure <ref> shows the running times and the peak memory as function of θ. Note that for the exact counting algorithm these metrics do not depend on θ, since it always counts all k-mers.  is always faster than the exact counting, with a difference that increases when θ increases and a speed-up around 2 even for θ = 2 · 10^-8. The memory requirement of  reduces when θ increases, and for θ = 2 · 10^-8 it is half of the memory required by the exact counting. This is due to 's sample size
being much smaller than the dataset size (Figure <ref>(d)), therefore a large portion of extremely low frequency k-mers are naturally left out from the random sample and do not need to be accounted for in the counting data structure, as confirmed by counting the number of distinct k-mers that are inserted in the counting data structure by the two algorithms (Figure <ref>(c)). (The difference between the memory requirement and the number of distinct k-mers is given by Jellyfish's strategy to doubles the size of the counting data structure when it is full.)







In terms of quality of the approximation, the output of  satisfied the guarantees given by Proposition <ref> for all runs of our experiments, therefore with probability higher than 1 - δ. While  may incur in false negatives, its false negative ratio (i.e., the fraction of k-mers in (,k,θ) not reported by ) is always ≤ 3 · 10^-4 (Figure <ref>(a)),
 even if the sampling technique of Section <ref> does not provide rigorous guarantees on such
quantity. Therefore  is very effective in reporting almost all frequent k-mers. As mentioned in Section <ref>,  can be easily modified so to report all frequent k-mers in output, even if at the cost of reporting also more k-mers with frequency between θ-ε and θ. In addition, the estimated frequencies f_A reported by  are always close to the true values f_(A), with a small maximum deviation |f_A - f_(A)| (Figure <ref>(b)), and an even smaller average deviation (Figure <ref>(c)). In addition, the upper and lower bounds computed as in Section <ref> provide small confidence intervals always containing the value f_(A) (e.g., Figure <ref>(d) for dataset ), and could be used to obtain sets of k-mers with various guarantees from the sample used by .



 §.§ Application to Metagenomics: Computation of Ecological Distances


We evaluate the use of  to speed up the computation of commonly used k-mer based ecological distances <cit.> between datasets of Next-Generation Sequencing (NGS) reads. We present results for the Bray-Curtis distance; analogous results hold for other distances and will be presented in the full version of this extended abstract.


We first investigated how the distances change when those are computed considering only the frequent k-mers (w.r.t. a frequency threshold θ) instead that the full spectrum of k-mers appearing in the data. Therefore, given a pair of datasets _1 and _2 and θ, we computed the sets 𝒪_1 = (_1,k,θ) and 𝒪_2 = (_2,k,θ) using Jellyfish and then computed a generalized version of the distances for all pairs of datasets we used for our experiments. For the Bray-Curtis distance, this generalization is defined as: BC(_1 , _2 , 𝒪_1 , 𝒪_2) = 1 - 2 ∑_A ∈𝒪_1 ∩𝒪_2min{ o__1(A) , o__2(A) }/∑_A ∈𝒪_1  o__1(A) + ∑_A ∈𝒪_2 o__2(A).

Note that when θ≤ 10^-10 then (,k,θ) coincides with the set of all k-mers, for any of the datasets we tested. The results (Figure <ref>(a)) show that for θ up to 5×10^-8 the values of the distances are fairly stable and therefore one can use only frequent k-mers for such values of θ to compute the distances, and that for θ up to 10^-7 the relation between distances of different pairs of datasets are almost always conserved. We underline that the exact counting approach needs to count all the k-mers and only afterwards can filter the infrequent ones before writing them to disk to compute (,k,θ).
We then used  to extract approximations (of k-mers and their frequencies) of (_1,k,θ) and (_2,k,θ) and used such approximations to compute the distances among datasets (Figure <ref>(b)). Strikingly, the distances computed from the output of  are very close to their exact variant (Figure <ref>(c)). Interestingly this holds also for the Jaccard distance,
a presence-based distance that does not depend neither on k-mer abundances nor on k-mer ranking by frequencies (detailed results will be provided in the full version of this extended abstract).





We then compared, for different values of θ, the total running time required to compute the approximations of the frequent k-mers using  for all datasets in Table <ref> and all distances among such datasets using  approximations with the running time required when the exact counting algorithm is used for the same tasks.  reduces the computing time by more than 75% (Figure <ref>(d)). This result comes from both the efficiency of  and from the fact that by focusing on the the most frequent k-mers we greatly reduce the number of distinct k-mers that need to be processed for computing the distances. Therefore  can be used for a very fast comparison of metagenomic datasets while preserving the ability of distinguishing similar datasets from different ones.



§ CONCLUSION


We presented , a sampling-based algorithm to approximate frequent k-mers and their frequencies with rigorous guarantees on the quality of the approximation. We show that  can be used to speed up the analysis of large high-throughput sequencing metagenomic datasets, in particular to compute abundance-based distances among such datasets. Interestingly  allows to compute accurate approximations also for presence-based distances (e.g., the Jaccard distance), even if for such distances other, potentially faster, tools <cit.> are available.  can be combined with any highly optimized method that counts all k-mers in a set of strings, including recent parallel methods designed for comparative metagenomics <cit.>. While we presented results for k-mers from datasets of short reads,  can also be used for the analysis of spaced seeds <cit.>, large datasets of long reads, and whole genome sequences.


plainurl


10

benoit2016multiple
G. Benoit, P. Peterlongo, et al.
Multiple comparative metagenomics using multiset k-mer counting.
PeerJ Computer Science, 2:e94, 2016.

bvrinda2015spaced
K. Břinda, M. Sykulski, and G. Kucherov.
Spaced seeds improve k-mer-based metagenomic classification.
Bioinformatics, 31(22):3584–3592, 2015.

brown2012reference
C.T. Brown, A. Howe, et al.
A reference-free algorithm for computational normalization of shotgun
  sequencing data.
arXiv preprint arXiv:1203.4802, 2012.

chikhi2013informed
R. Chikhi and P. Medvedev.
Informed and automated k-mer size selection for genome assembly.
Bioinformatics, 30(1):31–37, 2013.

danovaro2017submarine
R. Danovaro, M. Canals, et al.
A submarine volcanic eruption leads to a novel microbial habitat.
Nature ecology & evolution, 1(6):0144, 2017.

dickson2017carryover
L. B. Dickson, D. Jiolle, et al.
Carryover effects of larval exposure to different environmental
  bacteria drive adult trait variation in a mosquito vector.
Science advances, 3(8):e1700585, 2017.

girotto2016metaprob
S. Girotto, C. Pizzi, and M. Comin.
Metaprob: accurate metagenomic reads binning based on probabilistic
  sequence signatures.
Bioinformatics, 32(17):i567–i575, 2016.

hrytsenko2018efficient
Y. Hrytsenko, N. M. Daniels, and R. S. Schwartz.
Efficient distance calculations between genomes using mathematical
  approximation.
In Proc. of ACM-BCB, pages
  546–546, 2018.

kelley2010quake
D. R. Kelley, M. C. Schatz, and S. L. Salzberg.
Quake: quality-aware detection and correction of sequencing errors.
Genome biology, 11(11):R116, 2010.

kokot2017kmc
M. Kokot, M. Długosz, and S. Deorowicz.
Kmc 3: counting and manipulating k-mer statistics.
Bioinformatics, 33(17):2759–2761, 2017.

li2003estimating
X. Li and M. S. Waterman.
Estimating the repeat structure and length of dna sequences using
  ℓ-tuples.
Genome research, 13(8):1916–1922, 2003.

loffler2009shape
M. Löffler and J. M. Phillips.
Shape fitting on point sets with probability distributions.
In Proc. of ESA, pages 313–324,
  2009.

marccais2011fast
G. Marçais and C. Kingsford.
A fast, lock-free approach for efficient parallel counting of
  occurrences of k-mers.
Bioinformatics, 27(6):764–770, 2011.

melsted2014kmerstream
P. Melsted and B. V. Halldórsson.
Kmerstream: streaming algorithms for k-mer abundance estimation.
Bioinformatics, 30(24):3541–3547, 2014.

melsted2011efficient
P. Melsted and J. K. Pritchard.
Efficient counting of k-mers in DNA sequences using a Bloom filter.
BMC bioinformatics, 12(1):333, 2011.

mitzenmacher2017probability
M. Mitzenmacher and E. Upfal.
Probability and computing: Randomization and probabilistic
  techniques in algorithms and data analysis.
Cambridge University Press, 2017.

mohamadi2017ntcard
H. Mohamadi, H. Khan, and I. Birol.
ntCard: a streaming algorithm for cardinality estimation in genomics
  data.
Bioinformatics, 33(9):1324–1330, 2017.

ondov2016mash
B. D. Ondov, T. J. Treangen, et al.
Mash: fast genome and metagenome distance estimation using minhash.
Genome biology, 17(1):132, 2016.

pandey2017squeakr
P. Pandey, M. A. Bender, R. Johnson, and R. Patro.
Squeakr: an exact and approximate k-mer counting system.
Bioinformatics, 2017.

patro2014sailfish
R. Patro, S. M. Mount, and C. Kingsford.
Sailfish enables alignment-free isoform quantification from rna-seq
  reads using lightweight algorithms.
Nature biotechnology, 32(5):462, 2014.

pevzner2001eulerian
P. A. Pevzner, H. Tang, and M. S. Waterman.
An eulerian path approach to dna fragment assembly.
Proceedings of the National Academy of Sciences,
  98(17):9748–9753, 2001.

rizk2013dsk
G. Rizk, D. Lavenier, and R. Chikhi.
DSK: k-mer counting with very low memory usage.
Bioinformatics, 29(5):652–653, 2013.

roy2014turtle
R. S. Roy, D. Bhattacharya, and A. Schliep.
Turtle: Identifying frequent k-mers with cache-efficient algorithms.
Bioinformatics, 30(14):1950–1957, 2014.

salmela2016accurate
L. Salmela, R. Walve, E. Rivals, and E. Ukkonen.
Accurate self-correction of errors in long reads using de bruijn
  graphs.
Bioinformatics, 33(6):799–806, 2016.

sims2009alignment
G. E. Sims, S.-R. Jun, G. A. Wu, and S.-H. Kim.
Alignment-free genome comparison with feature frequency profiles
  (ffp) and optimal resolutions.
Proceedings of the National Academy of Sciences,
  106(8):2677–2682, 2009.

sivadasan2016kmerlight
N. Sivadasan, R. Srinivasan, and K. Goyal.
Kmerlight: fast and accurate k-mer abundance estimation.
arXiv preprint arXiv:1609.05626, 2016.

solomon2016fast
B. Solomon and C. Kingsford.
Fast search of thousands of short-read sequencing experiments.
Nature biotechnology, 34(3):300, 2016.

vapnik1998statistical
V. Vapnik.
Statistical learning theory. 1998.
Wiley, New York, 1998.

vapnik1971uniform
V. Vapnik and A. Chervonenkis.
On the Uniform Convergence of Relative Frequencies of Events to Their Probabilities.
Theory of Probability & Its Applications, 16 (2): 264, 1971.

wood2014kraken
D. E. Wood and S. L. Salzberg.
Kraken: ultrafast metagenomic sequence classification using exact
  alignments.
Genome biology, 15(3):R46, 2014.

zerbino2008velvet
D. R. Zerbino and E. Birney.
Velvet: algorithms for de novo short read assembly using de bruijn
  graphs.
Genome research, 18(5):821–829, 2008.

zhang2014these
Q. Zhang, J. Pell, et al.
These are not the k-mers you are looking for: efficient online k-mer
  counting using a probabilistic data structure.
PloS one, 9(7):e101271, 2014.

zhang2014rna
Z. Zhang and W. Wang.
RNA-skim: a rapid method for rna-seq quantification at transcript
  level.
Bioinformatics, 30(12):i283–i292, 2014.



