Multiple comparative metagenomics using multiset k-mer counting

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. Figure 2A 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, Lavenier & Chikhi, 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.

Furthermore, to limit disk bandwidth and avoid I/O bottleneck, partitions are compressed. A dictionary-based approach, such as the one provided in zlib (Deutsch & Gailly, 1996), is used. This type of compression is very well suited here since it efficiently packs the high redundancy of sorted k-mers.

Merging Count.

Here, the data partitioning introduced in the previous step is advantageously used to generate abundance vectors. The N files associated to a partition P_i, are taken as input of a merging process. These files contain k-mer counts sorted in lexicographical order. A Merge-Sort algorithm can thus be efficiently applied to directly generate abundance vectors.

In that scheme, P processes can be run independently, resulting in the generation of P abundance vectors in parallel, allowing to compute simultaneously P contributions of the ecological distance. Note that the abundance vectors do not need to be stored. They are only used as input streams for the next step. Figure 2B illustrates the Merging Count phase.

k-mer abundance filter.

Distinct k-mers with very low abundance usually come from sequencing errors. As a matter of fact, a single sequencing error creates up to k erroneous distinct k-mers. Filtering out these k-mers speeds-up the Simka process, as it greatly reduces the overall number of distinct k-mers, but may also impact the content of the distance matrix. This point is evaluated and discussed in the result section.

This filter is activated during the count process. Only k-mers whose abundance is equal to or greater than a given abundance threshold are kept. By default the threshold is set to 2. The k-mers that pass the filter are called “solid k-mers.”

Computing the Bray–Curtis distance.

The Bray–Curtis distance is given by the following equation: (1)${\displaystyle BrayCurti{s}_{Ab}\left(S_i,S_j\right)=1-2\frac{\sum _{w\in S_i\cap S_j}min\left(N_{S_i}\left(w\right),N_{S_j}\left(w\right)\right)}{\sum _{w\in S_i}{N}_{S_i}\left(w\right)+\sum _{w\in S_j}{N}_{S_j}\left(w\right)}}$

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∈SiN_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∩Sjmin(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–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²).

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 & 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, 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 & De Cáceres (2013), can be expressed for two datasets S_i and S_j as: (2)${\displaystyle Distance\left(S_i,S_j\right)=g\left(\sum _{w\in S_i\cup S_j}f\left(N_{S_i}\left(w\right),N_{S_j}\left(w\right),C_{S_i},C_{S_j}\right)\right)}$ 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∈SiN_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. Those distances easily fit in the previous framework as a = ∑_w∈Si∩Sj1_{{NSi(w)NSj(w)>0}}, C_Si = ∑_w∈Si1_{NSi(w)>0} = a + b and similarly C_Sj = a + c. Therefore, a is a crossed term and b and c can be deduced from a and the marginal terms.

In the same vein, Chao et al. (2006) introduced variations of presence-absence distances incorporating abundance information to account for unobserved species. The main idea is to replace “hard” quantities such as a∕(a + b), the fraction of distinct k-mers from S_i shared with S_j, by probabilistic “soft” ones: here the probability U ∈ [0, 1] that a k-mer from S_i is also found in S_j. Similarly, the “hard” fraction a∕(a + c) of distinct k-mers from S_j shared with S_i is replaced by the “soft” probability V that a k-mer from S_j is also found in S_i. U and V play the same role as a, b and c do in qualitative distances and are sufficient to compute the variants named AB-Jaccard, AB-Ochiai and AB-Sorensen. However and unlike the quantities a, b c, which can be observed from the data, U and V are not known in practice and must be estimated from the data. Chao et al. (2006) proposed several estimates for U and V. The most elaborate ones attempt to correct for differences in sampling depths and unobserved species by considering the complete k-mer counts vector of a sample. Those estimates are unfortunately untractable in our case as we stream only a few k-mer counts at a time. Instead we resort to the simplest estimates presented in Chao et al. (2006), which lend themselves well to the additive and distributed nature of Simka: U = Y_SiSj∕C_Si and V = Y_SjSi∕C_Sj where Y_SiSj = ∑_w∈Si∩SjN_Si(w)1_{NSj(w)>0} and C_Si = ∑_w∈SiN_Si(w). Note that Y_SiSj corresponds to crossed terms and is asymmetric, i.e., Y_SiSj ≠ Y_SjSi. Intuitively, U is the fraction of k-mers (not distinct anymore) from S_i also found in S_j and therefore gives more weights to abundant k-mers that its qualitative counterpart a∕(a + b).

Table 1 gives the definitions of the collection of distances computed by Simka while replacing species counts by k-mer counts. These are qualitative, quantitative and abundance-based variants of qualitative ecological distances. The table also provides their expression in terms of C_i, f and g, adopting the notations of Eq. (2).

Finally, note that the additive nature of the computed distances over k-mers is instrumental in achieving a linear time complexity (in W_s, the amount of distinct solid k-mers) and in the parallel nature of the algorithm. The algorithm is therefore not amenable to other, more complex classes of distances that account for ecological similarities between species (Pavoine et al., 2011), or edit distances between k-mers as those complex distances require all versus all k-mer comparisons.

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 (Fig. 3A). 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. Both simple and complex distance computation algorithms have theoretical worst case quadratic time complexity relatively to N (the number of samples). The difference of execution time comes then from the amount of operations required, in practice, to calculate the crossed terms of the distances. For a given abundance vector, the simple distances only need to be updated for each pair (S_i, S_j) such that N_Si > 0 and N_Sj > 0 whereas complex distances need to be updated for each pair such that N_Si > 0 or N_Sj > 0, entailing a lot more update operations. It is noteworthy that among all distances listed in Table 1, all distances are simple, except the Whittaker, Jensen–Shannon and Canberra distances.

When compared to other state of the art tools, namely Commet, Metafast and Mash, we parameterized Simka to compute only the Bray–Curtis distance, since all other tools compute only one such simple distance. The Figs. 3B–3D shows respectively the CPU time, the memory footprint and the disk usage of each tool with respect to an increasing number of samples N. Mash has definitely the best scalability but limitations of its computed distance are shown in the next section. Commet is the only one to show a quadratic time behaviour with N. For N = 40, Simka is 6 times faster than Metafast and 22 times faster than Commet. All tools, except Metafast, have a constant maximal memory footprint with respect to N. For metafast, we could not use its max memory usage option since it often created “out of memory” errors. The disk usage of the four tools increases linearly with N. The linear coefficient is greater for Simka and MetaFast, but it remains reasonable in the case of Simka, as it is close to half of the input data size, which was 11 GB for N = 40.

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 h 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 h (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 h 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. Notably, this does not degrade distance quality, at least for the HMP dataset, as shown in the next section. Additional tests on the impact of k on the performances show that the disk usage increases sub-linearly with k whereas the computation time and the memory usage stay constant (see Fig. S1).

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 Fig. 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.989).

Similarly, clear correlations (r > 0.89) are also observed between the percentage of matched k-mers and the percentage of similar reads as defined by BLAT alignments (Fig. 5). Interestingly, the correlation depends on the k-mer size and the identity threshold used for BLAT: larger k-mer sizes correlate better with higher identity thresholds and vice versa. The highest correlation is 0.987, obtained for Simka with k = 21 compared to BLAT results with 95% identity.

These results demonstrate that we can safely replace read-based metrics by a kmer-based one, and this enables to save huge amounts of time when working on large metagenomics projects. Moreover, the k-mer size parameter seems to be the counterpart of the identity threshold of alignment-based methods if one wants to tune the taxonomic precision level of the comparisons.

Correlation with taxonomic distances on the gut sample.

A traditional way of comparing metagenomics samples rely on so called taxonomic 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.89, see Fig. 6). On this figure, one may also notice that Simka measures are robust with the whole range of distances. On the other hand, Mash distances correlate badly with taxonomic ones (r = 0.51, see Fig. S3 and the comparison protocol in Article S1). This is probably due to the fact that gut samples differ more in terms of relative abundances of microbes than in terms of composition (see next section). As Mash can only output a qualitative distance, it is ill equipped to deal with such a case. Additionally, as shown in Fig. S3, this conclusion stands for the HMP samples from other body sites for which one disposes of high quality taxonomic distances.

Interestingly, these Simka results are robust with the k-mer filtering option and the k-mer size, as long as k is larger than 15 and with an optimal correlation obtained with k = 21 (see Fig. S4). Notably, with very low values of k (k < 15), the correlation drops (r = 0.5 for k = 12). This completes previous results suggesting that the larger the k the better the correlation, that were limited to k values smaller than 11 (Dubinkina et al., 2016).

Visualizing the structure of the HMP samples.

We propose to visualize the structure of the HMP samples and see if Simka is able to reproduce known biological results. To easily visualize those structures, we used the Principal Coordinate Analysis (PCoA) (Borg & Groenen, 2013) to get a 2-D representation of the distance matrix and of the samples: distances in the 2-D plane optimally preserve values of the distance matrix.

Figure 7 shows the PCoA of the quantitative Ochiai distance computed by Simka on the full HMP samples. We can see that the samples are clearly segregated by body sites. This is in line with results from studies of the HMP consortium (Human Microbiome Project Consortium, 2012a; Costello et al., 2009; Koren et al., 2013). Moreover, one may notice that different distances can lead to different distributions of the samples, with some clusters being more or less discriminated (see Fig. S5). This confirms the fact that it is important to conduct analyses using several distances as suggested in (Koren et al., 2013; Legendre & De Cáceres, 2013) as different distances may capture different features of the samples.

We conduct the same experiment on the 138 gut samples from the HMP project. Arumugam et al. (2011) showed that the gut samples are organized in three groups, known as enterotypes, and characterized by the abundance of a few genera: Bacteroides, Prevotella and genera from the Ruminococcaceae family. The original enterotypes were built from Jensen–Shannon distances on taxonomic profiles. The Fig. 8 shows the PCoA of the Jensen–Shannon distances obtained with Simka. Mapping the relative abundance of those genera in each sample, as provided by the HMP consortium, on the 2-D representation reveals a clear gradient in the PCoA space. Simka distances therefore recover biological features it had no direct access to: here, the fact that gut samples are structured along gradients of Bacteroides, Prevotella and Ruminococcaceae. The fact that Simka is able to capture such subtle signal raises hope of drawing new interesting biological insights from the data, in particular for those metagenomics project lacking good references (soil, seawater for instance).