dataesr/scientific-paragraphs-categorization huggingface dataset:
Actually, k-mers are a natural unit for comparing communities: (1) sufficiently long k-mers are usually specific of a genome  (Fofanov et al. , 2004) , (2) k-mer frequency is linearly related to genome's abundance  (Wu and Ye, 2011) , (3) k-mer aggregates organisms with very similar k-mer composition (e. g. related strains from the same bacterial species) without need for a classification of those organisms  (Teeling et al. , 2004). Dubinkina et al. (2016)  conducted an extensive comparison between k-mer-based distances and taxonomic ones (ie. based on taxonomic assignation against a reference database) for several large scale metagenomic projects. They demonstrate that k-mer-based distances are well correlated to taxonomic ones, and are therefore accurate enough to recover known biological structure, but also to uncover previously unknown biological features that were missed by reference-based approaches due to incompleteness of reference databases.
In this paper, we present Simka. Simka compares N metagenomic datasets based on their k-mers counts. It computes a large collection of distances classically used in ecology to compare communities. Computation is performed by replacing species counts by k-mer counts, for a large range of kmer sizes, including large ones (up to 30). Simka is, to our knowledge, the first method able to rapidly compute a full range of distances enabling the comparison of any number of datasets. This is performed by processing data on-the-fly (i. e. without storage of large temporary results). With the exception of Mash that is, thanks to sub-sampling, approximately two to five time faster, Simka outperforms state-of-the-art read comparison methods in terms of computational needs. For instance, Simka ran on 690 samples from the Human Microbiome Project (HMP) (Human Microbiome Project Consortium, 2012a) (totalling 32 billion reads) in less than 10 hours and using no more than 70 GB RAM.
The contributions of this manuscript are three-fold. First we propose a new method for efficiently counting k-mers from a large number of metagenomic samples. The usefulness of such counting is not limited to comparative metagenomics and may have applications in many other fields. Second, we show how to derive a large number of ecological distances from k-mer counts. And third, we show on real datasets that k-mer-based distances are highly correlated to taxonomic distances: they therefore capture the same underlying structure and lead to the same conclusions.
Overview Given N metagenomic datasets, denoted as S 1 , S 2 , S i ,. S N , the objective is to provide a N × N distance matrix D where D i,j represents an ecological distance between datasets S i and S j. Such possible distances are listed in Table  1. The computation of the distance matrix can be theoretically decomposed into two distinct steps: 1. k-mer count. Each dataset is represented as a set of discriminant features, in our case, k-mer counts. More precisely, a k-mer count matrix KC of size W × N is computed. W is the number of distinct k-mer among all the datasets. KC i,j represents the number of times a k-mer i is present in the dataset S j. 2. distance computation. Based on the k-mer count information, the distance matrix D is computed. Actually, many ecological distances (cf Table  1 ) can be derived from matrix KC when replacing species counts by k-mer counts.
Actually, Simka does not require to have the full KC matrix to start the distance computation. But for sake of simplicity, we will first consider this matrix to be available. The k-mer count step splits all the reads of the datasets into k-mers and performs a global count. This can be done by counting individually k-mers in each dataset, then merging the overall k-mer counts. The output is the matrix KC (of size W × N ). Efficient algorithms, such as KMC2  (Deorowicz et al. , 2015) , have recently been developed to count all the occurrences of distinct k-mers in a read dataset, allowing the computation to be executed in a reasonable amount of time and memory even on very large datasets. However, the main drawback of this approach is the huge main memory space it requires which is computed as follow: M em KC = W s * (8 + 4N ) bytes, with W s the number of distinct k-mers, N the number of samples, and 8 and 4 the number of bytes required to store respectively 31-mers and a k-mer count.
For example, experiments on the HMP (Human Microbiome Project Consortium, 2012a) datasets (690 datasets containing on average 45 millions of reads each) would require a storage space of 260T B for the matrix KC. However, a careful look at the definition of ecological distances (Table  1 ) shows that, up to some final transformation, they are all additive over the k-mers. Independent contributions to the distance can thus be computed in parallel from disjoint sets of k-mers and aggregated later on to construct the final distance matrix. Furthermore, each independent contribution can itself be constructed in an iterative way by receiving lines of the KC matrix, called abundance vectors, one at a time. The abundance vector of a specific k-mer simply consists of its N counts in the N datasets.
To sum up, instead of computing the complete k-mer count matrix KC, the alternative computation scheme we propose is to generate successive abundance vectors from which independent contributions to the distances can be iteratively updated in parallel. The great advantage is that the huge k-mer count matrix KC does not need to be stored anymore. However, this approach requires a new strategy to generate abundance vectors. We propose and describe below a new efficient multiset k-mer counting algorithm (called MKC) that can be highly parallelized on large computing resources infrastructures. As illustrated Fig. 1 , Simka uses abundance vectors generated by MKC for computing ecological distances. The first step takes as input N datasets and generates multiple streams of abundance vector from disjoint sets of k-mers.
Multiset k-mer Counting Starting from N datasets of reads, the aim is to generate abundance vectors that will feed the ecological distance computation step. This task is divided into two phases: 1. Sorting Count, 2. Merging Count. Sorting Count Each k-mer of a dataset is extracted and its canonical representation is stored (the canonical representation of a k-mer is the smallest lexicographic value between the k-mer and its reverse complement). Canonical k-mers are then sorted in lexicographical order. Distinct k-mers can thus be identified and their number of occurrences computed. As the number of distinct k-mers is generally huge, the sorting step is divided into two sub-tasks and proceeds as follows: the k-mers are first separated into P partitions, each stored on disk. After this preliminary task, each partition is sorted and counted independently, and stored again on disk.
Conceptually, at the end of the sorting count process, we dispose of N × P sorted partitions. As the same distribution function is applied to all datasets, a partition P i contains a specific subset of k-mers common to all datasets. Fig. 2 -A illustrates the Sorting Count phase. The Sorting Count phase has a high parallelism potential. A first parallelism level is given by the independent counts of each dataset. N processes can thus be run in parallel, each one dealing with a specific dataset. A second level is given by the fine grained parallelism implemented in software such as DSK  (Rizk et al. , 2013)  or KMC2  (Deorowicz et al. , 2015)  that intensively exploit today multicore processor capabilities. Thus, the overall Sorting Count process is especially suited for grid infrastructures made of hundred of nodes, and where each node implements 8 or 16-core systems.
Ecological distance computation Simka computes a collection of distances for all pairs of datasets. As detailed in the previous section, abundance vectors are used as input data. For the sake of simplicity, we first explain the computations of the Bray-Curtis distance. All other distances, presented later on, can be computed in the same way, with only small adaptations. Computing the Bray-Curtis distance The Bray-Curtis distance is given by the following equation: BrayCurtis Ab (S i , S j ) = 1 -2 w∈Si∩Sj min(N Si (w), N Sj (w)) w∈Si N Si (w) + w∈Sj N Sj (w) (1) where w is a k-mer and N Si (w) is the abundance of w in the dataset S i. We consider here that w ∈ S i ∩ S j if N Si (w) > 0 and N Sj (w) > 0. The equation involves marginal (or dataset specific) terms (i. e. w∈Si N Si (w) is the total amount of k-mers in dataset S i ) acting as normalizing constants and crossed terms that capture the (dis)similarity between datasets (i.
e. w∈Si∩Sj min(N Si (w), N Sj (w)) is the total amount of k-mers in the intersection of the datasets S i and S j ). Marginal and crossed terms are then combined to compute the final distance. Algorithm 1 shows that it is straightforward to compute the distance matrix between N datasets from the abundance vectors. Inputs of this algorithm are provided by the Multiple k-mer Counting algorithm (MKC). These are the P streams of abundance vectors and the marginal terms of the distance, i. e. the number of k-mers in each dataset, determined during the first step of the MKC which counts the k-mers. A matrix, denoted M ∩ , of dimension N × N is initialized (step 1) to record the final value of the crossed terms of each pair of datasets. P independent processes are run (step 2) to compute P partial crossed term matrices, denoted M ∩part (step 3), in parallel.
Each process iterates over its abundance vector stream (step 4). For each abundance vector, we loop over each possible pair of datasets (steps 5-6). The matrix M ∩part is updated (step 8) if the k-mer is shared, meaning that it has positive abundance in both datasets S i and S j (step 7). Since a distance matrix is symmetric with null diagonal, we limit the computation to the upper triangular part of the matrix M ∩part. The current abundance vector is then released. Each process writes its matrix M ∩part on the disk when its stream is done (step 9). When all streams are done, the algorithm reads each written M ∩part and accumulates it to M ∩ (step 10-11). The last loop (steps 13 to 16) computes the Bray-Curtis distance for each pair of datasets and fills the distance matrix reported by Simka.
The amount of abundance vectors streamed by the MKC is equal to W s , which is also the total amount of distinct solid k-mers in the N datasets. This algorithm has thus a time complexity of O(W s × N 2 ). Other ecological distances The distance introduced in Eq. 1 is a single example of ecological distance. There exists numerous other ecological distances that can be broadly classified into two categories (see  Legendre and De Cáceres (2013)  for a finer classification): distances based on presence-absence data (hereafter called qualitative) and distances based on proper abundance data (hereafter called quantitative). Qualitative distances are more sensitive to factors that affect presence-absence of organisms (such as pH, salinity, depth, Algorithm 1: Compute the Bray-Curtis distance (equation 1) between N datasets Input: -V s : vector of size P representing the abundance vector streams -V ∪ : vector of size N containing the number of k-mers in each dataset Output: a distance matrix Dist 1 M ∩ ← empty square matrix of size N // number of k-mers in each dataset intersection 2 In parallel: foreach abundance vector stream S in V s do 3 M ∩part ← empty squared matrix of size N // part of M ∩ 4 foreach abundance vector v in S do 5 for i ← 0 to N -1 do 6 for j ← i + 1 to N -1 do 7 if v[i] > 0 and v[j] > 0 then 8 M ∩part [i, j] ← M ∩part [i, j] + min(v[i], v[j]) 9 Write M ∩part to disk 10 foreach each written matrix M ∩part do 11 M ∩ ← M ∩ + M ∩part 12 Dist ← empty squared matrix of size N // final distance matrix 13 for i ← 0 to N -1 do 14 for j ← i + 1 to N -1 do 15 Dist[i, j] = 1 -2 * M ∩ [i, j] / (V ∪ [i] + V ∪ [j]) 16 Dist[j, i] = 1 -2 * M ∩ [i, j] / (V ∪ [i] + V ∪ [j]) 17 return Dist humidity, absence of light, etc) and therefore useful to study bioregions.
Quantitative distances focus on factors that affect relative changes (seasonal changes, nutrient availability, concentration of oxygen, depth, diet, disease, etc) and are therefore useful to monitor communities over time or along an environmental gradient. Note that some factors, such as pH, are likely to affect both presence-absence (for large changes in pH) and relative abundances (for small changes in pH). Algorithmically, most ecological distances, including most of those mentioned in  Legendre and De Cáceres (2013) , can be expressed for two datasets S i and S j as: Distance(S i , S j ) = g   w∈Si∪Sj f N Si (w), N Sj (w), C Si , C Sj   (2) where g and f are simple functions, and C Si is a marginal (i. e. dataset-specific) term of dataset S i , usually of size 1 (i. e. a scalar). In most distances, C Si is simply the total number of k-mers in S i.
By contrast, the value of f corresponds to crossed terms and requires knowledge of both N Si (w) and N Sj (w) (and potentially C Si and C Sj as well). For instance, for the abundance-based Bray-Curtis distance of Eq. 1, we have C Si = w∈Si N Si (w), g(x) = 1 -2x and f (x, y, X, Y ) = min(x, y)/(X + Y ). Those distances can be computed in a single pass over the data using a slightly modified variant of Algorithm 1. The marginal terms C Si are computed during the first step of the MKC which counts the k-mers of each dataset. The crossed terms involving f are computed and summed in steps 7-8 (but exact instructions depend on the nature of f ). Finally, the actual distances are computed in steps 15-16 and depend on both f and g. Qualitative distances form a special case of ecological distances: they can all be expressed in terms of quantities a, b and c where a is the number of distinct k-mers shared between datasets S i and S j , b is the number of distinct k-mers specific to dataset S i and c is the number of distinct k-mers specific to dataset S j.
Implementation Simka is based on the GATB library  (Drezen et al. , 2014) , a C++ library optimized to handle very large sets of k-mers. It includes a powerful implementation of the sorting count algorithm with the latest improvements in terms of k-mer counting introduced by  Deorowicz et al. (2015). Simka is usable on standard computers and has also been entirely parallelized for grid or cloud systems. It automatically splits the process into jobs according to the available number of nodes and cores. These jobs are sent to the job scheduling system, while the overall synchronization is performed at the data level. Simka is an open source software, distributed under GNU affero GPL License, available for download at https://gatb. inria. fr/software/simka/.
Results First, Simka performances are evaluated in terms of computation time, memory footprint and disk usage and compared to those of other state of the art methods. Then, the Simka distances are evaluated with respect to de novo and reference-based distances and with respect to known biological results. We conduct our numerical experiments on data from the Human Microbiome Project (HMP) (Human Microbiome Project Consortium, 2012a) which is currently one of the largest publicly available metagenomic datasets: 690 samples gathered from different human body sites (http://www. hmpdacc. org/HMASM/). The whole dataset contains 2*16 billions of Illumina paired reads distributed non uniformly across the 690 samples. One advantage of this dataset is that it has been extensively studied, in particular the microbial communities are relatively well represented in reference databases (Human Microbiome Project Consortium, 2012a,b) (see http://hmpdacc.
Performance Evaluation Performances on small datasets The scalability of Simka was first evaluated on small subsets of the HMP project, where the number of compared samples varied from 2 to 40. When computing a simple distance, such as Bray-Curtis for instance, Simka running time shows a linear behavior with the number of compared samples (Figure  3-A ). As expected, counting the kmers for each sample (MKC-count) consumes most of the time. This task has a theoretical time complexity linear with the number of kmers, and thus the number of samples, and this explains the observed linear behavior of the overall program. In fact, most steps of Simka, namely MKC-count, MKC-merge and simple distance computation, show a linear behavior between running time and the number of compared samples. The only exception is the computation of complex distances, where the time devoted to this task increases quadratically.
In summary, Simka and Mash seems to be the only tools able to deal with very large metagenomics datasets, such as the full HMP project. Performances on the full HMP samples Remarkably, on the full dataset of the HMP project (690 samples), the overall computation time of Simka is about 14 hours with very low memory requirements (see Table  2 ). By comparison, Metafast ran out of memory (it also ran out of memory while considering only a sub-sample composed of the 138 HMP gut samples) and Commet took several days to compute one 1-vs-all distance matrix and therefore would require years of computation to achieve the N × N distance matrix. Conversely, Mash ran in less than 5 hours (255 min) and is faster than Simka. This was expected since Mash outputs an approximation of a simple qualitative distance, based on a sub-sample of 10,000 k-mers.
By comparison, Simka computes numerous distances, including quantitative ones, over 15 billion distinct k-mers (see Table  2 ). Note that Simka is also designed for coarse-grain parallelism, and such computation took less than 10 hours on a 200-CPU platform. These results were obtained with default parameters, namely filtering out k-mers seen only once. On this dataset, this filter removes only 5 % of the data: solid k-mers (k-mers seen at least twice) account for 95% of all base pairs of the whole dataset (see Table  2 ). But interestingly, when speaking in terms of distinct k-mers, solid distinct k-mers represent less than half of all distinct k-mers before merging across all samples and only 15% after merging. Consequently, Simka performances are greatly improved, both in terms of computation time and disk usage when considering only solid k-mers.
Evaluation of the distances We evaluate the quality of the distances computed by Simka answering two questions. First, are they similar to distances between read sets computed using other approaches? Second, do they recover the known biological structure of HMP samples? For the first evaluation, two types of other approaches are considered, either de novo ones (similar to Simka but based on read comparisons), or taxonomic distances, e. g. approaches based on a reference database. Correlation with read-based approaches In this section, we focus on comparing Simka k-mer-based distance to two read-based approaches: Commet  (Maillet et al. , 2014)  and an alignment-based method using BLAT  (Kent, 2002). Both these read-based approaches define and use a read similarity notion. They derive the percentage of reads from one sample similar to at least one read from the other sample as a quantitative similarity measure between samples.
Commet considers that two reads are similar if they share at least t non-overlapping k-mers (here t = 2, k = 33). For BLAT alignments, similarity was defined based on several identity thresholds: two reads were considered similar if their alignment spanned at least 70 nucleotides and had a percentage of identity higher than 92%, 95% or 98%. For ease of comparison, Simka distance was transformed to a similarity measure, such as the percentage of shared kmers (see Article S1 for details of transformation). Looking at the correlation with Commet is interesting because this tool uses a heuristic based on shared k-mers but its final distance is expressed in terms of read counts. As shown in Figure  4 , on a dataset of 50 samples from the HMP project, Simka and Commet similarity measures are extremely well correlated (Spearman correlation coefficient r = 0.
875 distances that are based on sequence assignation to taxons by mapping to reference databases. To compare Simka to such traditional reference-based method, we used the HMP gut samples, which is a well studied dataset comprising 138 samples. The HMP consortium provides a quantitative taxonomic profile for each sample on its website. These profiles were obtained by mapping the reads on a reference genome catalog at 80% of identity. From these profiles, we computed the Bray-Curtis distance, latter used as a reference. The complete protocol to obtain taxonomic distances is given in Article S1. Only Mash and Simka results have been considered for this experiment. As previously mentioned, Commet and MetaFast could not scale this dataset. Simka k-mer-based distance appears very well correlated to the traditional taxonomic distance (r = 0.
Discussions In this article, we introduced Simka, a new method for computing a collection of ecological distances, based on k-mer composition, between many large metagenomic datasets. This was made possible thanks to the Multiple k-mer Count algorithm (MKC), a new strategy that counts k-mers with state-of-the-art time, memory and disk performances. The novelty of this strategy is that it counts simultaneously k-mers from any number of datasets, and that it represents results as a stream of data, providing counts in each dataset, k-mer per k-mer. The distance computation has a time complexity in O(W × N 2 ), with W is the number of considered distinct k-mers and N is the number of input samples. N is usually limited to a few dozens or hundreds and can not be reduced. However, W may range in the hundreds of billions.
73%) PC2 (12. 37%) provides large speed improvement without affecting the results, at least on the tests performed on the HMP datasets. However, the HMP dataset is not representative of all metagenomics projects and, in some cases, this filter may not be desired. For instance, in the case of samples with low coverage or when performing qualitative studies where the rare species have more impact. As a matter of fact, it is notable that Simka is able to scale large datasets even with the solid filter disabled as shown in the performance section. Interestingly, when applied on a low coverage dataset, namely the Global Ocean Sampling  (Yooseph et al. , 2007) , Simka was able to capture the essential underlying biological structure with or without the k-mer solid filter (see Fig. S6 ). However, an important incoming challenge is to precisely measure the impact of applied thresholds together with the choice of k, depending on the input dataset features such as community complexity and sequencing effort.
On the other hand, genome size heterogeneity and k-mer composition variation along a microbe genome will bias the k-mer-based distance upwards. However, species composition based approaches are not feasible for large read sets from complex ecosystems (soil, seawater) due to the lack of good references and/or mapping scaling limitations. Moreover, our proposal has the advantage of being a de novo approach, unbiased by reference banks inconsistency and incompleteness. Finally, numerical experiments on the HMP datasets show that k-mer based and taxonomic distances are well correlated (r > 0. 8 for k ≥ 21) and consequently that Simka recovers the same biological structure as taxonomic studies do. There is nevertheless room for improving Simka distances. For instance, recently,  Břinda et al. (2015)  showed that spaced seeds can improve the k-mer-based metagenomic classification obtained with the popular tool Kraken  (Wood and Salzberg, 2014).
Figure1. Simka strategy. The first step takes as input N datasets and generates multiple streams of abundance vector from disjoint sets of k-mers. The abundance vector of a k-mer consists of its N counts in the N datasets. These abundance vectors are taken as input by the second step to iteratively update independent contributions to the ecological distance in parallel. Once an abundance vector has been processed, there is no need to keep it on record. The final step aggregates each contribution and computes the final distance matrix.