Fast Approximation of Frequent $k$-mers and Applications to Metagenomics

TL;DR

SAKEIMA is a sampling-based method that efficiently approximates frequent $k$-mers in large sequencing datasets with rigorous guarantees, enabling faster biological analyses without processing entire datasets.

Contribution

This work introduces SAKEIMA, a novel sampling approach that uses VC dimension analysis to accurately estimate frequent $k$-mers with fewer data processed.

Findings

SAKEIMA achieves significant speed-ups in large dataset analysis.

It provides rigorous approximation guarantees for $k$-mer frequencies.

Accurate $k$-mer based distances are computed using only a fraction of the data.

Abstract

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 SAKEIMA, 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. SAKEIMA 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 SAKEIMA allows to rigorously approximate frequent $k$-mers by processing only a fraction of a dataset and that the frequencies estimated by SAKEIMA lead to accurate estimates of $k$-mer based distances between high-throughput sequencing datasets. Overall, SAKEIMA is an efficient and rigorous tool to estimate $k$-mers abundances providing significant speed-ups in the analysis of large sequencing datasets.

Tables

Table 1. Table 1: Datasets for our experimental evaluation. For each dataset 𝒟 𝒟 \mathcal{D} the table shows: the dataset name and site ( (s) for stool , (t) for tongue dorsum ); the total number t 𝒟 , k subscript 𝑡 𝒟 𝑘 t_{\mathcal{D},k} of k 𝑘 k -mers ( k = 31 𝑘 31 k=31 ) in 𝒟 𝒟 \mathcal{D} ; the number | 𝒟 | 𝒟 |\mathcal{D}| of reads it contains; the maximum read length max n i = max i ⁡ { n i | r i ∈ 𝒟 } subscript subscript 𝑛 𝑖 subscript 𝑖 conditional subscript 𝑛 𝑖 subscript 𝑟 𝑖 𝒟 \max_{n_{i}}=\max_{i}\{n_{i}|r_{i}\in\mathcal{D}\} ; the average read length avg = n i ∑ i = 0 n − 1 n i / n {}_{n_{i}}=\sum_{i=0}^{n-1}n_{i}/n .

dataset	$t_{𝒟,k}$	$\left|𝒟\right|$	${\max }_{n_i}$	avg${}_{n_i}$
SRS024388(s)	$7.92\cdot {10}^9$	$1.20\cdot {10}^8$	102	97.21
SRS011239(s)	$8.13\cdot {10}^9$	$1.24\cdot {10}^8$	102	96.69
SRS024075(s)	$8.82\cdot {10}^9$	$1.38\cdot {10}^8$	96	94.88
SRS075404(t)	$7.75\cdot {10}^9$	$1.22\cdot {10}^8$	102	94.51
SRS062761(t)	$8.26\cdot {10}^9$	$1.18\cdot {10}^8$	101	101.00
SRS043663(t)	$9.15\cdot {10}^9$	$1.31\cdot {10}^8$	101	101.00

Equations

$$FK(\mathcal{D},k,\theta )=\{(A,f_{\mathcal{D}}(A)):f_{\mathcal{D}}(A)\ge \theta \}.$$

$${\hat{f}}_{P_{\ell }}(A)=\frac{1}{m}\sum _{I_{\ell ,i}\in P_{\ell }}\frac{{\hat{\varphi }}_A(I_{\ell ,i})}{\ell }.$$

$$m\ge \frac{2}{(\ell \epsilon {)}^2}\left(\lfloor {\log }_2\min (2\ell ,{\sigma }^k)\rfloor +\ln \left(\frac{1}{\delta }\right)\right).$$

Full text

In this work, we develop, analyze, and test, a sampling-based approach, called SAKEIMA, 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. SAKEIMA 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 SAKEIMA allows to rigorously approximate frequent $k$-mers by processing only a fraction of a dataset and that the frequencies estimated by SAKEIMA lead to accurate estimates of $k$-mer based distances between high-throughput sequencing datasets. Overall, SAKEIMA is an efficient and rigorous tool to estimate $k$-mers abundances providing significant speed-ups in the analysis of large sequencing datasets.

1 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 [21, 31], error correction [9, 24], repeat detection [11], genome comparison [25], digital normalization [3], RNA-seq quantification [20, 33], metagenomic reads classification [30] and binning [7], fast search-by-sequence over large high-throughput sequencing repositories [27]. 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 [13, 15]. 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 [13], DSK [22], KMC 3 [10] and Squeakr-exact [19]. 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 [14], khmer [32], Kmerlight [26] and ntCard [17] 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 [4] 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 [19] relies on a probabilistic data structure to approximate the counts of individual $k$-mers. Turtle [23] 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 [1, 5, 6] to assess the distance between the corresponding datasets. In contrast to presence-based distances [18] (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 (ii) if the abundances of frequent $k$-mers can be computed more efficiently than the counts of all $k$-mers. Recently, preliminary work [8] has shown that, for the cosine distance and $k=12$, the answer to the first question is positive, and in Section 4 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

(ii) 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 $\theta$ 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 $\varepsilon$. Second, we propose a new method, Sampling Algorithm for $K$-mErs approxIMAtion (SAKEIMA), to obtain an approximation to the set of frequent $k$-mers using sampling. SAKEIMA 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 SAKEIMA 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 SAKEIMA 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 SAKEIMA, 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.

2 Preliminaries

Let a dataset $\mathcal{D}$ be a bag of $n$ reads $\mathcal{D}=\{r_{0},\dots,r_{n-1}\}$, where each read $r_{i}$, $0\leq i\leq n-1$, is a string of length $n_{i}$ from an alphabet $\Sigma$ of cardinality $|\Sigma|=\sigma$. For $j\in\{0,\dots,n_{i}-1\}$, let $r_{i}[j]$ be the $j$-th character of $r_{i}$. For a given integer $k\leq\min_{i}\{n_{i}:r_{i}\in\mathcal{D}\}$, we define a $k$-mer $A$ as a string of length $k$ from $\Sigma$, that is $A\in\Sigma^{k}$. We say that a $k$-mer $A$ appears in $r_{i}$ at position $j\in\{0,\dots,n_{i}-k\}$ if $r_{i}[j+h]=A[h],\forall h\in\{0,\dots,k-1\}$. For every $i,0\leq i\leq n-1$, and every $j\in\{0,\dots,n_{i}-k\}$, we define the indicator function $\phi_{r_{i},A}(j)$ that is $1$ if the $k$-mer $A$ appears in $r_{i}$ at position $j$, while $\phi_{r_{i},A}(j)=0$ otherwise. The total number of $k$-mers in $\mathcal{D}$ is $t_{\mathcal{D},k}=\sum_{i=0}^{n-1}(n_{i}-k+1)$. We define the support $o_{\mathcal{D}}(A)$ of a $k$-mer $A$ as the number of distinct positions in $\mathcal{D}$ where $A$ appears: $o_{\mathcal{D}}(A)=\sum_{i=0}^{n-1}\sum_{j=0}^{n_{i}-k}\phi_{r_{i},A}(j)$. We define the frequency $f_{\mathcal{D}}(A)$ of $A$ in $\mathcal{D}$ as the ratio between the number of distinct positions where $A$ appears in $\mathcal{D}$ and the total number of $k$-mers in $\mathcal{D}$: $f_{\mathcal{D}}(A)=o_{\mathcal{D}}(A)/t_{\mathcal{D},k}$.

2.1 Frequent $k$-mers and Approximations

We are interested in obtaining the set $FK(\mathcal{D},k,\theta)$ of frequent $k$-mers in a dataset $\mathcal{D}$ with respect to a minimum frequency threshold $\theta$, defined as follows.

Definition 1

Given a dataset $\mathcal{D}$, an integer $k>0$, and a frequency threshold $\theta\in(0,1]$, the set $FK(\mathcal{D},k,\theta)$ of Frequent $k$-Mers in $\mathcal{D}$ w.r.t. $\theta$ is the collection of all $k$-mers with frequency at least $\theta$ in $\mathcal{D}$ and of their corresponding frequencies in $\mathcal{D}$:

[TABLE]

$FK(\mathcal{D},k,\theta)$ can be computed with a single scan of all the $k$-mers occurrences in $\mathcal{D}$ maintaining the $k$-mers supports in an appropriate data structure; however, when $\mathcal{D}$ is extremely large and $k$ is not small, the exact computation of $FK(\mathcal{D},k,\theta)$ 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 $FK(\mathcal{D},k,\theta)$ may be preferable, provided it ensures rigorous guarantees on its quality. In this work, we focus on the following approximation.

Definition 2

Given a dataset $\mathcal{D}$, an integer $k>0$, a frequency threshold $\theta\in(0,1]$, and a constant $\varepsilon\in(0,\theta)$, an $\varepsilon$-approximation of $FK(\mathcal{D},k,\theta)$ is a collection $C=\{(A,f_{A}):f_{A}\in(0,1]\}$ such that:

• •

for any $(A,f_{\mathcal{D}}(A))\in FK(\mathcal{D},k,\theta)$ there is a pair $(A,f_{A})\in C$;

• •

for any $(A,f_{A})\in C$ it holds that $f_{\mathcal{D}}(A)\geq\theta-\varepsilon$;

• •

for any $(A,f_{A})\in C$ it holds that $|f_{\mathcal{D}}(A)-f_{A}|\leq\varepsilon/2$.

The definition above guarantees that every frequent $k$-mer of $\mathcal{D}$ is in the approximation and that no $k$-mer with frequency $<\theta-\varepsilon$ is in the approximation. The third condition guarantees that the estimated frequency $f_{A}$ of $A$ in the approximation is close (i.e, within $\varepsilon/2$) to the frequency $f_{\mathcal{D}}(A)$ of $A$ in $\mathcal{D}$. It is easy to show that obtaining a $\varepsilon$-approximation of $FK(\mathcal{D},k,\theta)$ with absolute certainty requires to process all $k$-mers in $\mathcal{D}$.

2.2 Simple Sampling-Based Algorithms and Bounds

We aim to provide an approximation to $FK(\mathcal{D},k,\theta)$ with sampling, by processing only randomly selected portions of $\mathcal{D}$. 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_{\mathcal{D},k}=\{(i,j):i\in[0,n-1],j\in[0,n_{i}-k]\}$ (note that $|P_{\mathcal{D},k}|=t_{\mathcal{D},k}$) of all positions where $k$-mers occurs in the dataset $\mathcal{D}$, 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 $\theta\in(0,1]$ one can define the set of frequent $k$-mers (and their frequencies) in the sample $P$ as $FK(P,k,\theta)=\{(A,f_{P}(A)):f_{P}(A)\geq\theta\}$, where $f_{P}(A)$ is the frequency of $k$-mer $A$ in the sample.

Obtaining a $\varepsilon$-approximation from a random sample with absolute certainty is impossible, thus we focus on obtaining a $\varepsilon$-approximation with probability $1-\delta>0$, where $\delta\in(0,1)$ is a confidence parameter, whose value is provided by the user. Intuitively, the set $FK(\mathcal{D},k,\theta)$ 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 $\varepsilon$-approximation. By using Hoeffding’s inequality [16] to bound the deviation of the frequency of a $k$-mer $A$ in the sample from $f_{\mathcal{D}}(A)$ and a union bound on the maximum number $\sigma^{k}$ of $k$-mers, where $\sigma=|\Sigma|$, we have the following result that provides a first such bound, and a corresponding first algorithm to obtain a $\varepsilon$-approximation to $FK(\mathcal{D},k,\theta)$. (Due to space constraints proofs are omitted and will be provided in the full version of this extended abstract.)

Proposition 1

Consider a sample $P$ of size $m$ of $\mathcal{D}$. If $m\geq\frac{2}{\varepsilon^{2}}\left(\ln\left(2\sigma^{k}\right)+\ln\left(\frac{1}{\delta}\right)\right)$ for fixed $\varepsilon\in(0,\theta),\delta\in(0,1)$, then, with probability $\geq 1-\delta$, $FK(P,k,\theta-\varepsilon/2)$ is a $\varepsilon$-approximation of $FK(\mathcal{D},k,\theta)$.

In addition, by using known results in statistical learning theory [29, 16] relating the VC dimension (see Section 3 for its definition) of a family of functions and a novelly derived bound on the family of functions $\{f_{\mathcal{D}}(A)\}$, we obtain the following improved bound and algorithm. (The derivation will be provided in the full version.)

Proposition 2

Let $P$ be a sample of size $m$ of $\mathcal{D}$. For fixed $\varepsilon\in(0,\theta),\delta\in(0,1)$, if $m\geq\frac{2}{\varepsilon^{2}}\left(1+\ln\left(\frac{1}{\delta}\right)\right)$ then $FK(P,k,\theta-\varepsilon/2)$ is an $\varepsilon$-approximation for $FK(\mathcal{D},k,\theta)$ with probability $\geq 1-\delta$.

3 Advanced and Practical Bounds and Algorithms for $k$-mer Approximations

While the bound of Proposition 2 significantly improves the simple bounds of Section 1, since the factor $\ln(2\sigma^{k})$ has been reduced to $1$, it still has an inverse quadratic dependency with respect to the accuracy parameter $\varepsilon$, that is problematic when the quantities to estimate are small. In these cases, one needs a small $\varepsilon$ to produce a meaningful approximation (since $\varepsilon<\theta$), and the inverse quadratic dependence of the sample size from $\varepsilon$ 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 $\sum_{A\in\Sigma^{k}}f_{\mathcal{D}}(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 $\theta$ (and therefore $\varepsilon$) typically needs to be set to a very low value. For example, a typical dataset from the Human Microbiome Project (HMP) has $n\approx 10^{8}$ reads of (average) length $\approx 100$: therefore if we are interested in $k$-mers for $k=31$, by setting $\delta=0.05$ the bound of Section 2.2 gives $\varepsilon\approx 10^{-5}$, that is only $k$-mers with frequency $\geq 10^{-5}$ could be reliably reported by sampling. However, in datasets we considered, no or a very small number ($\leq 30$) of $k$-mers have frequency $\geq 10^{-5}$, therefore according to the result from Section 2.2 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.

3.1 Sampling Bags of Positions and VC dimension Bound.

We propose a method to provide an efficiently computable approximation to $FK(\mathcal{D},k,\theta)$ when the minimum frequency $\theta$ is low, by properly defining samples so that any $k$-mer $A$ will appear in a sample with probability higher than $f_{\mathcal{D}}(A)$, thus lessening the the dependence of the sample size from $1/\varepsilon^{2}$.For this to be achievable, we need to relax the notion of approximation defined in Section 2. In particular, the guarantees, provided by our method, in such relaxed approximation are that all $k$-mers with frequency above $\theta^{\prime}$, with $\theta^{\prime}$ slightly higher than $\theta$, are reported in output, and that no $k$-mer having frequency below $\theta-\varepsilon$ is reported in output. (See Proposition 5 for the definition of $\theta^{\prime}$.) Our experiments show that the fraction of $k$-mers having frequency $\in[\theta,\theta^{\prime})$ 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 $\ell$ positions chosen independently at random from the set $P_{\mathcal{D},k}$ of all positions.

Let $I_{\ell}=\{(i_{1},j_{1}),(i_{2},j_{2}),\dots,(i_{\ell},j_{\ell})\}$ be a bag of $\ell$ positions for $k$-mers in $\mathcal{D}$, chosen uniformly at random from the set $P_{\mathcal{D},k}$. We define the indicator functions $\hat{\phi}_{A}(I_{\ell})$ that, for a given bag $I_{\ell}$ of $\ell$ positions, is equal to $1$ if $k$-mer $A$ appears in at least one of the $\ell$ positions in $I_{\ell}$ and is equal to [math] otherwise. That is $\hat{\phi}_{A}(I_{\ell})=\min\left\{1,\sum_{(i,j)\in I_{\ell}}\phi_{r_{i},A}(j)\right\}.$ We define the $\ell$-positions sample $P_{\ell}$ as a bag of $m$ bags $\{I_{\ell,0},I_{\ell,1},\dots,I_{\ell,m-1}\}$, where each $I_{\ell,j},0\leq j\leq m-1$ is a bag of $\ell$ positions, sampled independently, and

[TABLE]

Intuitively, $\hat{f}_{P_{\ell}}(A)$ is the biased version of the unbiased estimator $f_{P_{\ell}}(A)=\frac{1}{m}\sum_{I_{\ell,i}\in P_{\ell}}\frac{\sum_{(i,j)\in I_{\ell,i}}\phi_{r_{i},A}(j)}{\ell}$ of $f_{\mathcal{D}}(A)$, where the bias arises from considering a value of $1$ every time $\sum_{(i,j)\in I_{\ell,i}}\phi_{r_{i},A}(j)>1$.

In our analysis we use the Vapnik-Chervonenkis (VC) dimension [28, 29], 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\subset X$, the projection of $R_{X}$ on $D$ is defined as $proj_{R_{X}}(D)=\{r\cap D:r\in 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)=\infty$.

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.

Definition 3

Let $Q=(X,R_{X})$ be a range space and let $D$ be a finite subset of $X$. For $\varepsilon\in(0,1]$, a subset $B$ of $D$ is an $\varepsilon$-approximation of $D$ if for all $r\in R_{X}$ we have: $\left|\frac{|D\ \cap\ r|}{|D|}-\frac{|B\ \cap\ r|}{|B|}\right|\leq\varepsilon/2.$

The following result [16] relates $\varepsilon$ and the probability that a random sample of size $m$ is an $\varepsilon$-approximation for a range space of VC dimension at most $v$.

Proposition 3 ([16])

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<\varepsilon$, $\delta<1$, then a random subset $B\subset D$ of cardinality $m$ with $m\geq\frac{4c}{\varepsilon^{2}}\left(v+\ln\left(\frac{1}{\delta}\right)\right)$ is a $\varepsilon$-approximation of $D$ with probability at least $1-\delta$.

The universal constant $c$ has been experimentally estimated to be at most $0.5$ [12].

We now prove an upper bound to the VC dimension $VC(Q)$ of the range space $Q$ associated to the class of functions $\hat{\phi}_{A}$ that grows sub-linearly with respect to $\ell$. To this aim, we first define the range space associated to bags of $\ell$ positions of $k$-mers.

Definition 4

Let $\mathcal{D}$ be a dataset of $n$ reads and let $k$ and $\ell$ be two integers $\geq 1$. We define $Q=(X_{\mathcal{D},k,\ell},R_{\mathcal{D},k,\ell})$ to be the following range space:

• •

$X_{\mathcal{D},k,\ell}$ is the set of all bags of $\ell$ positions of $k$-mers in $\mathcal{D}$, that is the set of all possible subsets, with repetitions, of size $\ell$ from from $P_{\mathcal{D},k}$;

• •

$R_{\mathcal{D},k,\ell}=\{P_{\mathcal{D},\ell}(A)|A\in\Sigma^{k}\}$ is the family of sets of starting positions of $k$-mers, such that for each $k$-mer $A$, the set $P_{\mathcal{D},\ell}(A)$ is the set of all bags of $\ell$ starting positions in $\mathcal{D}$ where $A$ appears at least once.

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

Proposition 4

Let $Q$ the range space from Definition 4. Then: $VC(Q)\leq\lfloor\log_{2}(\ell)\rfloor+1$.

Using the result above, we prove the following.

Proposition 5

Let $\ell\geq 1$ be an integer and $P_{\ell}$ be a bag of $m$ bags of $\ell$ positions of $\mathcal{D}$ with

[TABLE]

Then, with probability at least $1-\delta$:

• •

for any $k$-mer $A\in FK(\mathcal{D},k,\theta)$ such that $f_{\mathcal{D}}(A)\geq\theta^{\prime}=1-(1-\ell\theta)^{1/\ell}$ it holds $\hat{f}_{P_{\ell}}(A)\geq\theta-\varepsilon/2$;

• •

for any $k$-mer $A$ with $\hat{f}_{P_{\ell}}(A)\geq\theta-\varepsilon/2$ it holds $f_{\mathcal{D}}(A)\geq\theta-\varepsilon$;

• •

for any $k$-mer $A\in FK(\mathcal{D},k,\theta)$ it holds $f_{\mathcal{D}}(A)\geq\hat{f}_{P_{\ell}}(A)-\varepsilon/2$;

• •

for any $k$-mer $A$ with $\hat{f}_{P_{\ell}}(A)-\varepsilon/2\geq 0$, it holds $f_{\mathcal{D}}(A)\geq 1-(1-\ell(\hat{f}_{P_{\ell}}(A)-\varepsilon/2))^{1/\ell}$;

• •

for any $k$-mer $A$ with $\ell(\hat{f}_{P_{\ell}}(A)+\varepsilon/2)\leq 1$ it holds $f_{\mathcal{D}}(A)\leq 1-(1-\ell(\hat{f}_{P_{\ell}}(A)+\varepsilon/2))^{1/\ell}$.

Note that from Proposition 5 the set $\{(A,f_{P_{\ell}}(A)):\hat{f}_{P_{\ell}}(A)\geq\theta-\varepsilon/2\}$ is almost a $\varepsilon$-approximation to $FK(\mathcal{D},k,\theta)$: in particular, there may be $k$-mers $A$ for which $\mathbb{E}[\hat{f}_{P_{\ell}}(A)]=(1-\left(1-f_{\mathcal{D}}(A)\right)^{\ell})/\ell<\theta$ while $f_{\mathcal{D}}(A)=\mathbb{E}[f_{P_{\ell}}(A)]\geq\theta$ and such that for the given sample $P_{\ell}$ we have $\hat{f}_{P_{\ell}}(A)\approx\mathbb{E}[\hat{f}_{P_{\ell}}(A)]-\varepsilon/2$. While this can happen, we can limit the probability of this happening by appropriately choosing $\ell$, and still enjoy the reduction in sample size of the order of $\frac{\log_{2}\ell}{\ell^{2}}$ w.r.t. Proposition 2 obtained by considering bags of bags of $\ell$ positions. In particular, this result allows the user to set $\theta$, $\varepsilon$, $\delta$, and $\ell$ to effectively find, with probability at least $1-\delta$, all frequent $k$-mers $A$ for which $f_{\mathcal{D}}(A)\geq\theta^{\prime}$ and do not report any $k$-mer with frequency below $\theta-\varepsilon$, while still being able to report in output almost all $k$-mers with frequency $\in[\theta,\theta^{\prime})$. Our experimental analysis (Section 4) shows that in practice choosing $\ell$ close from below to $1/\theta$ is very effective to obtain such result. Then, the third, fourth, and fifth guarantees from Proposition 5 state that we can use the biased estimates $\hat{f}_{P_{\ell}}(A)$ to derive guaranteed upper and lower bounds to $f_{\mathcal{D}}(A)$ that will be much tighter than the one obtained using the bounds of Section 2.2. We will show how to obtain further improved upper and lower bounds to $f_{\mathcal{D}}(A)$ in Section 3.3. Such lower bounds $\ell b_{A}$ can be used, for example, to prove that the set $\{(A,f_{P_{\ell}}(A)):\ell b_{A}\geq\theta-\varepsilon\}$ enjoys the same last four guarantees from Proposition 5 while the first one holds for a $\theta^{\prime}<1-(1-\ell\theta)^{1/\ell}$; therefore, when false negatives are problematic, the set $\{(A,f_{P_{\ell}}(A)):\ell b_{A}\geq\theta-\varepsilon\}$ can be used to obtain a different approximation of $FK(\mathcal{D},k,\theta)$ with fewer false negatives.

3.2 SAKEIMA: An Efficient Algorithm to Approximate Frequent $k$-mers

We now present our Sampling Algorithm for K-mErs approxIMAtion (SAKEIMA), that builds on Proposition 5 and efficiently samples a bag $P_{\ell}$ of bags of $\ell$-positions from $\mathcal{D}$ to obtain an approximation of the set $FK(\mathcal{D},k,\theta)$ with probability $1-\delta$, where $\delta$ is a parameter provided by the user.

SAKEIMA is described in Algorithm 1. SAKEIMA performs a pass on the stream of $k$-mers appearing in $\mathcal{D}$, and for each position in the stream it samples the number $a$ of times that the position appears in the sample $P_{\ell}$ independently at random from the Poisson distribution $Poisson(\lambda)$ of parameter $\lambda=m\ell/t_{\mathcal{D},k}$. SAKEIMA 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_{\ell}$. (Note that $t_{\mathcal{D},k}$, that can be computed with a very quick linear scan of the dataset, where $n_{i}$ is computed for every $r_{i}\in\mathcal{D}$ without extracting and processing (e.g., inserting or updating information for) $k$-mers; in alternative a lower bound to $t_{\mathcal{D},k}$ can be used, simply resulting in a number of samples higher than needed). Then, such occurrences are partitioned into the $m$ bags $I_{\ell,0},\dots,I_{\ell,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 $\leq 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_{\ell}$ (line 11) divided by the total number of positions in the sample, while the biased estimate $\hat{f}_{A}$ is computed as the number $|\mathcal{P}_{A}|$ of distinct bags of $P_{\ell}$ where $A$ appears at least once divided by the number $m$ of bags (lines 12-13). Then SAKEIMA flags $A$ as frequent if $\hat{f}_{A}\geq\theta-\varepsilon/2$ (line 14) and, in this case, the couple $(A,f_{A})$ is added to the output set $\mathcal{O}$ (line 14), since $f_{A}$ is the best (and unbiased) estimate to $f_{\mathcal{D}}(A)$. Note that bags for different values of $\ell$ (on the same sampled positions) can be obtained by maintaining a table $T_{\ell}$ and a set $\mathcal{P}_{A,\ell}$ for each value $\ell$ of interest.

Note that SAKEIMA does not sample $m$ bags of exactly $\ell$ positions each, since the number of occurrences of each position in $\mathcal{D}$ in the sample $P_{\ell}$ is sampled independently from a Poisson distribution, even if the expected number of total occurrences sampled from the algorithm is $m\ell$. However, the independent Poisson distributions used by SAKEIMA provide an accurate approximation of the random sampling of exactly $m\ell$ positions used in the analysis of Section 3.1. In particular, this holds when one focuses on the events of interests for our approximation of Section 3.1 (e.g., the event “there exists a $k$-mer $A$ such that $|\mathbb{E}[\hat{f}_{P_{\ell}}(A)]-\hat{f}_{P_{\ell}}(A)|>\varepsilon/2$”). In fact, a simple adaptation of a known result (Corollary 5.11 of [16]) on the relation between sampling with replacement and the use of independent Poisson distributions gives the following.

Proposition 6

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 $\mathcal{O}$ features the guarantees of Proposition 5 with probability $\geq 1-\delta^{\prime}$, with $\delta^{\prime}=2\delta$.

3.3 Improved Lower and Upper Bounds to $k$-mers Frequencies

Note that Proposition 5 guarantees that we can obtain upper and lower bounds to $f_{\mathcal{D}}(A)$ for every $A\in FK(\mathcal{D},k,\theta)$ from the sample of bags of $\ell$ positions. These bounds are meaningful only in specific ranges of the frequencies; for example, the lower bound from the third guarantee in Proposition 5 is meaningful when the frequency of $A$ is fairly low, i.e $f_{\mathcal{D}}(A)\approx 1/\ell$, while for very frequent $k$-mers they could be a multiplicative factor $1/\ell$ away from than the correct value. For example, if a $k$-mer is very frequent and appears in all bags of $\ell$ $k$-mers in a sample $S$, its corresponding lower bound is still only $1/\ell-\varepsilon/2$.

However, Proposition 5 can be generalized to obtain tighter upper and lower bounds to the frequency of all $k$-mers. For given $\ell$, $\varepsilon$, and $\delta$, let $m$ as given in Proposition 5. Note that the total number of $k$-mer’s positions in the sample $P_{\ell}$ is $m\ell$. Let $\mathcal{L}$ be a set of integer values $\mathcal{L}=\{\ell_{i}\}$ with $\ell_{i}\in[1,m\ell],\forall i=0,\dots,|\mathcal{L}|-1$. Now, for every $\ell_{i}\in\mathcal{L}$, we can partition the same $m\ell$ $k$-mers that are in $P_{\ell}$ into $m_{i}=m\ell/\ell_{i}$ partitions having size $\ell_{i}$. Let $P_{\ell_{i}}$ be such a random partition of such positions into $m_{i}$ bags of $\ell_{i}$ positions each. Note that each $P_{\ell_{i}}$ is a “valid” sample (i.e., a sample of independent bags of positions, each obtained by uniform sampling with replacement) for Proposition 5, even if the $P_{\ell_{i}}$’s are not independent. From each $P_{\ell_{i}}$, we define a maximum deviation $\varepsilon_{i}$ from Proposition 5 as $\varepsilon_{i}=\frac{1}{\ell_{i}}\sqrt{\frac{2}{m_{i}}\left(\lfloor\log_{2}(\min(2\ell_{i},\sigma^{k}))\rfloor+\ln\left(|\mathcal{L}|/\delta\right)\right)}$. We have the following result.

Proposition 7

With probability at least $1-\delta$, for all $k$-mers $A$ simultaneously and for all the random partitions induced by $\mathcal{L}$ it holds

• •

$f_{\mathcal{D}}(A)\geq\max\{\hat{f}_{P_{\ell_{i}}}(A)-\varepsilon_{i}/2:i=0,\dots,|\mathcal{L}|-1\}$;

• •

$f_{\mathcal{D}}(A)\geq\max\{1-(1-\ell(\hat{f}_{P_{\ell_{i}}}(A)-\varepsilon_{i}/2))^{1/\ell}:i=0,\dots,|\mathcal{L}|-1$* and $\hat{f}_{P_{\ell}}(A)-\varepsilon_{i}/2\geq 0\}$;*

• •

$f_{\mathcal{D}}(A)\leq\min\{1-(1-\ell(\hat{f}_{P_{\ell_{i}}}(A)+\varepsilon_{i}/2))^{1/\ell}:i=0,\dots,|\mathcal{L}|-1$* and $\hat{f}_{P_{\ell}}(A)+\varepsilon_{i}/2\leq 1/\ell\}$.*

In our experiments, we use $\mathcal{L}=\{\ell_{i}\}$ with $\ell_{i}=\ell/2^{i},\forall i\in[0,\lfloor\log_{2}\ell\rfloor-1]$; in this case, note that $P_{\ell_{0}}=P_{\ell}$. 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 $FK(\mathcal{D},k,\theta)$ with different guarantees. For example, by reporting all $k$-mers (and their frequencies) that have an upper bound $\geq\theta$, we have an approximation that guarantees that all $k$-mers $A$ with $f_{\mathcal{D}}(A)\geq\theta$ are in the approximation.

4 Experimental Results

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

4.1 Datasets and Implementation

We considered six datasets from the Human Microbiome Project (HMP)222https://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 stool and the three largest of tongue dorsum (Table 1). 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 SAKEIMA in C++ as a modification of Jellyfish [13] (the version we used is 2.2.10333https://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 SAKEIMA with the exact counting performed by Jellyfish, since they share the underlying characteristics, allowing us to evaluate the impact of SAKEIMA sampling strategy. We did not include the time to compute $t_{\mathcal{D},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 C++ 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 g++ 4.9.4. SAKEIMA can be used in combination with more efficient algorithms and implementations for the computation of these (and other) distances [1], resulting in speed-ups analogous to the ones we present below. For all the experiments of SAKEIMA, given $\theta$ and a dataset $\mathcal{D}$, we fixed the parameters $\delta=0.1$, $\varepsilon=\theta-2/t_{\mathcal{D},k}$, $m=100$, and we fix $\ell$ to the minimum value satisfying the $\varepsilon$-approximation. For all the experiments we have $\ell$ close from below to $1/\theta$. For all the metrics we considered, we report the results for one random run.

4.2 Approximation of the Frequent $k$-mers

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

(i) running time444Every instance of SAKEIMA and Jellyfish was executed with $1$ worker, i.e., sequentially. Note that the Poisson approximation employed by SAKEIMA allows multiple workers to independently process the input $k$-mers, therefore SAKEIMA 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;

(ii) memory requirement.

We also assessed the accuracy of the output of SAKEIMA.

Figure 2 shows the running times and the peak memory as function of $\theta$. Note that for the exact counting algorithm these metrics do not depend on $\theta$, since it always counts all $k$-mers. SAKEIMA is always faster than the exact counting, with a difference that increases when $\theta$ increases and a speed-up around $2$ even for $\theta=2\cdot 10^{-8}$. The memory requirement of SAKEIMA reduces when $\theta$ increases, and for $\theta=2\cdot 10^{-8}$ it is half of the memory required by the exact counting. This is due to SAKEIMA’s sample size being much smaller than the dataset size (Figure 2(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 2(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 SAKEIMA satisfied the guarantees given by Proposition 5 for all runs of our experiments, therefore with probability higher than $1-\delta$. While SAKEIMA may incur in false negatives, its false negative ratio (i.e., the fraction of $k$-mers in $FK(\mathcal{D},k,\theta)$ not reported by SAKEIMA) is always $\leq 3\cdot 10^{-4}$ (Figure 3(a)), even if the sampling technique of Section 3.1 does not provide rigorous guarantees on such quantity. Therefore SAKEIMA is very effective in reporting almost all frequent $k$-mers. As mentioned in Section 3.3, SAKEIMA 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 $\theta-\varepsilon$ and $\theta$. In addition, the estimated frequencies $f_{A}$ reported by SAKEIMA are always close to the true values $f_{\mathcal{D}}(A)$, with a small maximum deviation $|f_{A}-f_{\mathcal{D}}(A)|$ (Figure 3(b)), and an even smaller average deviation (Figure 3(c)). In addition, the upper and lower bounds computed as in Section 3.3 provide small confidence intervals always containing the value $f_{\mathcal{D}}(A)$ (e.g., Figure 3(d) for dataset SRS062761), and could be used to obtain sets of $k$-mers with various guarantees from the sample used by SAKEIMA.

4.3 Application to Metagenomics: Computation of Ecological Distances

We evaluate the use of SAKEIMA to speed up the computation of commonly used $k$-mer based ecological distances [1] 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 $\theta$) instead that the full spectrum of $k$-mers appearing in the data. Therefore, given a pair of datasets $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ and $\theta$, we computed the sets $\mathcal{O}_{1}=FK(\mathcal{D}_{1},k,\theta)$ and $\mathcal{O}_{2}=FK(\mathcal{D}_{2},k,\theta)$ 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(\mathcal{D}_{1},\mathcal{D}_{2},\mathcal{O}_{1},\mathcal{O}_{2})=1-2\frac{\sum_{A\in\mathcal{O}_{1}\cap\mathcal{O}_{2}}\min\{o_{\mathcal{D}_{1}}(A),o_{\mathcal{D}_{2}}(A)\}}{\sum_{A\in\mathcal{O}_{1}}o_{\mathcal{D}_{1}}(A)+\sum_{A\in\mathcal{O}_{2}}o_{\mathcal{D}_{2}}(A)}$.

Note that when $\theta\leq 10^{-10}$ then $FK(\mathcal{D},k,\theta)$ coincides with the set of all $k$-mers, for any of the datasets we tested. The results (Figure 4(a)) show that for $\theta$ up to $5\times 10^{-8}$ the values of the distances are fairly stable and therefore one can use only frequent $k$-mers for such values of $\theta$ to compute the distances, and that for $\theta$ 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 $FK(\mathcal{D},k,\theta)$. We then used SAKEIMA to extract approximations (of $k$-mers and their frequencies) of $FK(\mathcal{D}_{1},k,\theta)$ and $FK(\mathcal{D}_{2},k,\theta)$ and used such approximations to compute the distances among datasets (Figure 4(b)). Strikingly, the distances computed from the output of SAKEIMA are very close to their exact variant (Figure 4(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 $\theta$, the total running time required to compute the approximations of the frequent $k$-mers using SAKEIMA for all datasets in Table 1 and all distances among such datasets using SAKEIMA approximations with the running time required when the exact counting algorithm is used for the same tasks. SAKEIMA reduces the computing time by more than $75\%$ (Figure 4(d)). This result comes from both the efficiency of SAKEIMA 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 SAKEIMA can be used for a very fast comparison of metagenomic datasets while preserving the ability of distinguishing similar datasets from different ones.

5 Conclusion

We presented SAKEIMA, a sampling-based algorithm to approximate frequent $k$-mers and their frequencies with rigorous guarantees on the quality of the approximation. We show that SAKEIMA 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 SAKEIMA allows to compute accurate approximations also for presence-based distances (e.g., the Jaccard distance), even if for such distances other, potentially faster, tools [18] are available. SAKEIMA 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 [1]. While we presented results for $k$-mers from datasets of short reads, SAKEIMA can also be used for the analysis of spaced seeds [2], large datasets of long reads, and whole genome sequences.