Low-density locality-sensitive hashing boosts metagenomic binning

Institute for Interdisciplinary Information Sciences, Tsinghua University, Beijing, China
Computer Science and Artificial Intelligence Laboratory, MIT, Cambridge, MA, USA
Department of Mathematics, MIT, Cambridge, MA, USA
Department of Computer Science, University of Illinois at Urbana-Champaign, Urbana, IL, USA
Corresponding authors: bab@mit.edu and jianpeng@illinois.edu

1 Introduction

Metagenomics techniques enable researchers to analyze the functional and genetic composition of microbial communities from environmental samples. Amplicon-based sequencing methods, which focus on the diversity of given marker genes, e.g. the 16S rRNA gene, provide efficient phylogenetic and functional diversity surveys of microbial communities. Due to the cost effectiveness of 16S rRNA sequencing, marker gene based analysis has frequently been used for studies involving large sample sets. With recent advances in next-generation sequencing (NGS) technologies, whole genome- or fragment-based metagenomics provides much richer information on and broader functional characterisics of microbial communities in the samples. For instance, novel hypotheses of microbial functions and potential enzymes have been identified through such metagenomic analysis [1].

During the past several years, high-throughput metagenomic sequencing has been extensively applied; however, the inherent complexity of metegenomic sequencing data poses a number of computational and statistical challenges for data analysis. Normally, DNA fragments, such as sequence reads or contigs, need to be first assigned to their organisms of origin (also called “binning”), since genes are typically sequenced from multiple, diverse organisms. After sequence fragments are assigned to taxonomic origins, downstream data analysis can be applied to elucidate the structure of microbial populations and assign functional annotations [1]. Note that in this work, we focus on the whole-genome metagenomic DNA sequencing, instead of the marker-gene based or gene-centric methods that only analyze the protein-coding regions for which other protein search algorithms have been proposed [2, 3].

Arguably the most popular metagenomic binning approaches are alignment-based methods. A sequence fragment is searched against a reference database with full genomes of organisms, and the highest scoring organism is assigned as the taxonomic origin. Although efficient sequence alignment algorithms, including BWA-MEM [4], Bowtie2 [5] and (mega)BLAST [6], can be readily used for this purpose, the computational cost of alignment-based methods becomes prohibitive as the size of the sequence dataset dramatically grows, which is often the case in recent studies.

Another completely different binning approach is based on genomic sequence composition. Codon usage, oligonucleotide frequencies and GC content often are distinct in different genomes. Computational classification methods have exploited such differences to identify sequences with similar compositional features. Typically, a supervised classifier, such as a naive Bayesian classifier, a neural network or a support vector machine (SVM), is trained on a set of reference genome sequences to classify the origins of metagenomic fragments [7, 8, 9, 10]. Since the lengths of metagenomic fragments can vary from 200 to 10,000 basepairs, sequence compositional features are often designed to be within a fixed dimensionality. Short -mers, contiguous nucleotide fragments with basepairs, have been shown to be both efficient and effective for metagenomic binning. For example, PhyloPythia [11] uses an ensemble of SVM models trained on contiguous 6-mers and demonstrates good performance on large datasets. Its successor, PhyloPythiaS [8], further improves the binning accuracy by tweaking the SVM model and simultaneously including -mers of multiple sizes () as compositional features. Since compositional methods need to compute only the -mer profiles for query sequence fragments, these methods are significantly faster than alignment-based methods on large datasets, although without providing alignment resolution and often suffering a moderate loss of the robustness. While longer -mers, which capture compositional dependency within larger contexts, could potentially lead to higher binning accuracy, they are more prone to noise and errors if used in the supervised setting. Moreover, incorporating long -mers as features increases computational cost exponentially and requires significantly larger training datasets. Note that there are existing methods use mid-size -mers (e.g. ) but they are mainly used for fast indexing and nearest (or exact) search [12, 13, 14, 15] but not in the supervised manner.

Here we overcome these bottlenecks of handling long -mers for classification, enabling fast, accurate and robust metagenomic binning. We introduce a novel compositional metagenomic binning algorithm, Opal, which efficiently encodes long -mers using low-dimensional profiles. To use long -mers as features in an SVM, we would need up to dimensions, which becomes practically infeasible if . Inspired by the low-density parity-check codes (also known as Gallager codes) from coding theory [16, 17], we propose the use of a set of low-density locality-sensitive hash (LSH) functions [18] to represent long -mers or sequence fragments. We have two major conceptual advances in this work. First, although LSH has been previously used for fast sequence alignment and assembly [19, 20], to the best of our knowledge, it is the first time that the idea of LSH has been proposed for compositional-based metagenomic binning. Second, we have developed LSH functions based on the Gallager design for very long -mers (e.g. ), which makes LSH practically applicable for this problem. Gallager codes were initially designed for error correction but we use it to design highly efficient LSH functions for fast binning of metagenomic fragments. We have also observed error tolerance for binning, which is partially due to the design of Gallager codes. Methodologically, starting from a Gallager design matrix with row weight , we construct hash functions to encode high-order sequence compositions within a -mer. In contrast to the complexity it would take to represent contiguous -mers, our proposed Gallager LSH adaptation requires only time. For very long -mers, we can construct the Gallager LSH functions in a hierarchical way to further capture compositional dependencies from both local and global contexts. To evaluate the performance of the Gallager LSH method, we trained an SVM model with features generated by the Gallager LSH method. When tested on a large dataset with 50 microbial species, Opal achieved better binning accuracy than the traditional method that uses contiguous -mer profiles as features [9, 8, 7]. Moreover, our method is more robust to mutations and sequencing errors, compared to the method with the contiguous -mer representation. We also compared Opal with BWA-MEM [4], the state-of-the-art alignment-based method. Remarkably, we achieved up to two orders of magnitude improvement in binning speed on large datasets; our method is also substantially more accurate than BWA-MEM when the rate of sequencing errors is high (e.g., 10-15%). It is remarkable to show that a compositional binning approach can be as robust as or even more robust than alignment-based approaches, in the presence of high sequencing errors or mutations in metagenomic sequence data. Finally, we demonstrate that it is possible to combine both compositional and alignment-based methods, by applying the compositional SVM with the Gallager LSH coding as a “coarse-search” procedure to reduce the taxonomic space for a subsequent alignment-based BWA-MEM “fine search”, to enable both efficient and accurate metagenomic binning, with improved binning accuracy, metagenomic alignment and near 20 times speedup. Previously, a similar “coarse search” approach has been proposed to speed up the metagenomic mapping of protein-coding sequences and speed up Diamond, an earlier state-of-the-art method, by 10 times [3, 2]. Note that with Opal, “coarse search” is performed by a supervised method instead of a unsupervised clustering approach, thus potentially better encoding the dependency within the data and leading to a larger speedup.

2 Metagenomic binning revisited

Metagenomic sequencing techniques produce a large data sets of DNA fragments (e.g. reads or contigs) from environmental samples. To understand the microbial communities and functional structures within the samples, we need to first assign or bin these sequence fragments with the taxonomic origins from which they were derived to facilitate downstream analyses. A straightforward approach for metagenomic binning is through sequence assembly. Since the DNA fragments are sampled from chromosomes of some unknown species, we should be able to identify the original species if we can reconstruct the chromosomes from the sequence fragments. However, it is often not feasible to generate accurate assemblies from metagenomic sequence fragments, due to potential undersampled organisms, ambiguity among closely related species and the limited capability and complexity of existing assembly algorithms. To deal with large datasets, efficient and accurate metagenomic binning algorithms are thus a pressing need.

2.1 Alignment-based methods

Possibly the most widely used binning methods are based on sequence alignment. Metagenomic fragments are binned according to their sequence similarity to a reference database consisting of genomes with taxonomical annotations. Binning tools, including MEGAN [21], incorporating sequence alignment programs such as (mega)BLAST and BWA-MEM and assigning taxonomic groups or organisms, have been successfully applied in many studies. Although alignment-based methods can provide high accuracy and high resolution, the demanding requirement of computational cost makes them prohibitive for large metagenomic sequence datasets, as one must align each fragment to every genome in the reference database.

2.2 Compositional-based methods

Instead of the time-consuming sequence alignment, sequence compositional-based binning methods exploit the sequence characteristics of metagenomic fragments and apply machine learning classification algorithms to assign putative taxnomic origins for all fragments. Since classifiers, such as support vector machines, are trained on the whole reference genome sequences beforehand, compositional methods normally are substantially faster than alignment-based methods on large datasets. The rationale of compositional-based binning methods is based on the fact that different genomes have different conserved sequence composition patterns, such as GC content, codon usage or a particular abundance distribution of consecutive nucleotide -mers. To design a good compositional-based algorithm, we need to extract informative and discriminative features from the reference genomes. Most existing methods, including PhyloPythia(S) [8, 11], use the -mer frequency to represent sequence fragments. In the rest of this section, we will give a brief review of compositional-based methods. Here we also want make clear that there are existing methods which utilize mid-size -mers (e.g. ) for fast indexing and nearest (or exact) search but not in the supervised manner [12, 13, 14, 15]. This work is focused on comparisons of fragment feature representation for supervised binning, so we will leave comparisons to these methods in the future work.

-mer profile. We assume that a sequence fragment , where , contains nucleotides. A -mer, with , is a short word of contiguous nucleotides. We define the -mer profile of in a vector representation . If we index each -mer as a binary string with length , then we have a one-to-one mapping between any -mer and an integer from to . In the rest of the paper, we will not distinguish the -mer string with its integer presentation for notational simplicity. Each coordinate in the -mer profile stores the frequency of -mer in the sequence fragment . For instance, for a fragment , its -mer profile has non-zero entries: , , and . In this way, instead of representing a -nucleotide fragment in , we can use -mer profile to represent it in . Many previous studies have shown that a small , e.g. , works reasonably well in practice, although longer can improve the binning accuracy but model training becomes a serious issue because of the high dimensionality which grows exponentially in . A recent study [9] has found that even with a highly tuned indexing technique, we cannot easily handle -mers with in the RAM.

Classification. After the -mer profile has been constructed, we can use supervised machine learning classification algorithms, such as logistic regression, naive Bayes classifier and support vector machines, to train a binning model. The training data can be generated by sampling -nucleotide fragments from the reference genomes with taxonomic annotations. Since metagenomic fragment can have different lengths depending on the applied sequencing technologies, it is possible to construct a number of binning models, each corresponding to a particular fragment length. Because the binning classifier often only involves vector multiplication, the speed of compositional-based binning algorithms is much faster than that of alignment-based methods, thus more suitable for large datasets. On the other hand, due to the fact that -mer profile can only capture the local patterns within a fragment, existing compositional binning algorithms usually have lower binning accuracy than the alignment-based methods which compare fragments and references in a global way. In addition, compositional-based classification methods are generally more sensitive to mutations or sequencing errors, partially due to the way -mer profile is constructed.

3 Opal: Gallager locality-sensitive hashing for fragment binning

In this work, we introduce Opal, a novel compositional-based metagenomic binning algorithm, that robustly represents long -mers (e.g. ) in a compact way to better capture the long-range compositional dependency in a fragment. The key idea of our algorithm is built on locality-sensitive hashing, a dimensionality reduction technique that hashes input high-dimensional data into low-dimensional buckets, with the goal to maximize the probability of collisions for similar input data. LSH has been widely used in bioinformatics for fast indexing for sequence alignment and assembly [19, 20]. To the best of our knowledge, it is the first time that LSH functions have been applied for compositional-based metagenomic binning. We propose to use them for represent metagenomic fragments compactly and subsequently for machine learning classification algorithms to train metagenomic binning models. Since metagenomic fragments can be very long, sometimes from hundreds of bp to tens of thousands of bp, we hope to construct compositional profiles to encode long-range dependency in long -mers. To handle large , we develop string LSH functions to compactly encode the global dependency with -mers in a low-dimensional feature vector, as oppose to directly using a -length -mer profile vector. Although LSH functions are usually constructed in a uniformly random way, we propose a new and efficient design of LSH functions based on the idea of the low-density parity-check (LDPC) code invented by Robert G Gallager for noisy message transmission [16, 17]. A key observation is that Gallager’s LDPC design not only leads to a family of LSH functions but also makes them efficient such that even a small number of random LSH functions can well encode the long fragments. Different from uniformly random LSH functions, the Gallager LSH functions are constructed structurally and hierarchically to ensure the compactness of the feature representation and the robustness when sequencing noise appears in the data.

3.1 Locality-sensitive hashing

LSH is a family of hash functions that have the property that two similar objects are mapped to the same hash value [18]. For the metagenomic binning problem, we are only interested in strings of length . Then a family of LSH functions can be defined as functions which map -mers into a -dimensional Euclidean space. Assume that we consider Hamming distances between -mers, if we choose randomly and for two -mers and with at most different positions, holds with probability at least . For two -mers and with more than different positions, holds with probability at least . With the construction of a LSH family, we can amplify or by sampling multiple hash functions from the family. Compared with the straightforward -mer indexing representation, the LSH scheme can be more compact and more robust. For example, we can construct LSH functions such that . Moreover, when a small number of sequencing errors or mutations appear in the -mer, LSH can still map the noisy -mer into a feature representation that is very similar to original -mer. This observation is highly significant since mutations or sequencing errors are generally inevitable in the data, and we hope to develop compositional-based methods less sensitive to such noises.

One way to construct LSH functions on strings under Hamming distance is to construct index functions by uniformly sampling a subset of positions from the -mer. Specifically, given a string of length over , we choose indices uniformly at random from without replacement. Then, the spaced -mer can be generated according to and these indices. More formally, we can define a random hash function to generate a spaced -mer explicitly:


The hash value can also be seen as a dimensional binary vector with only the string ’s corresponding coordinate set to and otherwise . It is not hard to see that such LSH function has the property that it maps two similar -mers to the same hash value with high probability. For example, consider two similar -mers and that differ by at most nucleotides, then the probability that they are mapped to the same value is given by


For two -mers and that differ at least nucleotides, the probability that they are mapped to different value is given by


With the family of LSH functions, we randomly sample a set of LSH functions and concatenate them together as the feature vector for a long -mer. Note that the complexity of the LSH-based feature vector is only , much smaller compared to that is the complexity of the complete -mer profile. More importantly, the LSH-based feature vector is not sensitive to errors or mutations in the -mer if and are well chosen, but for the traditional -mer profile, even one nucleotide change can change the feature vector completely. To compute the feature vector for a metagenomic fragment with length , we first extract all -mers by sliding a window of length over the sequence, and then apply on each -mer to generate LSH-based feature vectors and then normalize the sum of the feature vectors by . In this way, one can easily show that similar fragments can also be mapped to similar LSH-based feature vectors. After the feature vectors are generated for fragments with taxonomic annotations, we train a linear classifier for metagenomic binning. It is also fairly straightforward to show that similar fragments have similar classification responses if the coefficients of the linear classification function are bounded. One may expect that the complexity of linear classification with -mer profiles would be lower since there are at most different -mers in a fragment and can be computed easily using sparse vector multiplications, but we find that the LSH-based feature vector is also sparse in practice and the indexing overhead is much smaller when constructing the feature vectors, since the LSH-based method can have much smaller dimensionality. In practice, the LSH-based methods can sometimes be even faster if and are not too large.

3.2 Gallager low-density locality-sensitive hashing

Despite that the random LSH function family described above has a lot of nice theoretical properties, uniformly sampled LSH functions are usually not optimal in practice. Theoretical properties of LSH functions hold probabilistically, which means that we need to sample a large number of random LSH functions to make sure the bounds are tight. However, practically, we simply cannot use a very large number of random LSH functions to build feature vectors for metagenomic fragments, given the limited computational resources. Thus it would be ideal if we could construct a small number of random LSH functions that are sufficiently discriminative and informative to represent long -mers. Here we take inspiration from the Gallager code or low-density parity-check code that has been widely used for noisy communication. The idea behind the Gallager code is similar to our LSH family but with a different purpose, namely error correction. The goal of the LDPC code is to generate a small number of extra bits when transmitting a binary string via a noisy channel [16, 17]. These extra bits are constructed to capture the long-range dependency in the binary string before the transmission. After the message string and these extra bits have been received, a decoder can perform error correction by performing probabilistic inference to compare the differences between the message string and these code bits to infer the correct message string. In the same spirit, we here adopt the idea behind the design of the LDPC code to construct a compact set of LSH functions for metagenomic binning.

To construct compact LSH functions, we hope to not waste coding capacity on any particular position in the -mer. While, under expectation, uniformly sampled spaced -mers on average cover each position equally, with a small number of random LSH functions, it is likely that we will see imbalanced coverage among positions since the probability of a position being chosen is binomially distributed. The Gallager’s design of LDPC, on the other hand, generates a subset of positions not uniformly random but make sure to equally cover each position [16]. So we can use the Gallager’s design to generate spaced -mers. The Gallager’s LDPC matrix is a binary matrix with dimension , and has exactly 1’s in each rows and 1’s in each column. The matrix can be divided into blocks with rows in each block. We first define the first block of rows as an matrix :

where each row of matrix has exactly consecutive 1’s from left to right across the columns. Every other block of rows is a random column permutation of the first set, and the LDPC matrix is given by:

where is a uniform random permutation matrix for . An example with is shown in Figure 1. An equivalent bipartite graph with the Gallager design matrix as the adjacency matrix also is shown. The algorithm for constructing the LDPC design matrix is shown in Algorithm S1 in Supplementary Information.

Figure 1: An illustration of Gallager LSH method. Left: an example of Gallager LDPC matrix . Right: The bipartite graph corresponding to . Each cycle node corresponds to a position in a -mer, and each square node corresponds to a row in , which generates a spaced -mer.

We use each row of to extract a spaced -mer to construct an LSH function. Note that the first set of gives contiguous -mers. With Gallager LSH functions, we can see that each position in a -mer is equally covered times, while the same uniformly sampled LSH function is very likely to have very imbalanced coverage times for different positions because of the high variance (). To further improve the efficiency, we construct random LSH functions with minimal overlap using a modified Gallager design algorithm. The idea is to avoid the “4-cycles” in the bipartite graph representation, as we hope not to encode two positions together in two “redundant” LSH functions [17]. An algorithm which finds “4-cycles” and removes them is shown in Algorithm S2 in Supplementary Information. For very long -mers, we can use a hierarchical approach to generate low-dimensional LSH functions for very long-range compositional dependency in -mers. We first generate a number of intermediate spaced -mers using the Gallager’s design matrix. Then from these -mers, we again apply the Gallager’s design to generate -mers to construct the hierarchical LSH functions. Moreover, it is not hard to see that the hash functions generated from the Gallager design are also a family of LSH functions. Finally, in this work, we use a simple linear SVM to train a classification model with the Gallager LSH features. We will also test other sophisticated classifiers, such as the structured classifier that considers the taxonomic hierarchy during training, in the future. Evaluating on a small dataset with species, we found that LSH functions generated by hierarchical Gallager design is more robust and more accurate than the uniformly sampled LSH functions, indicating its efficiency for long -mers (see Supplementary Figure S2).

3.3 Compositional-based binning as “coarse search” in compressive metagenomics

After we train the Gallager LSH-based binning classifier, we can use it as a “coarse search” procedure in the compressive genomics manner to reduce the search space of alignment-based methods [3]. For example, if we want to perform binning against reference organisms, instead of comparing a fragment to all reference genomes, we first apply the compositional-based binning classifier to identify a very small subset or group of putative taxonomic origins that are ranked very highly by the classifier. Then we perform sequence alignment between the fragment and the reference genomes of the top-ranked organisms. This natural combination of compositional-based and alignment-based methods provides metagenomic binning with high scalability, high accuracy and high-resolution alignments.

4 Results

4.1 Experimental setting

We downloaded 50 complete bacterial genomes from NCBI database as suggested by [9] (see Supplementary Table S1), and simulated metagenomic samples by generating fragments from these reference genome sequences. We set the coverage , and generate fragments of length bp and bp from the reference genome sequences. The number of fragments sampled from a genome sequence of a species, , is determined by a coverage number , which is defined by , where is then length of the whole genome sequence of the species. This coverage value is chosen since we found that larger will not further improve the performance of the classifier but it may vary if we include more species in the dataset. In addition, with a very large , it would be difficult to design fair experiments, because there would be a lot substantially overlapped fragments in training and test sets. For , we sampled 71,259 training fragments and for , we sampled 35,631 training fragments. To assess the robustness of our method, we also simulated mutation and sequencing error in a genome sequence with rates in {5%, 10%, 15%}. In each mutated position, the original nucleotide is substituted by one of the other three type of nucleotides with equal probability. We have also experimented on noisy data with 1% indels and observed similar results. In this work, we use this setting to show the effectiveness of the Gallager LSH method as a proof-of-concept. We plan to further investigate other values of parameters in the future, for example, fragments with or nucleotides, other noise models and a much larger set of reference genomes.

To evaluate the performance of our method, we simulated a test set of fragments that are not used and with less than 50% overlapped positions with any training fragment. The ratio between the size of training data and that of test data is 5:1. For each of the -mer profile method and different Opal settings, we trained a multiclass SVM with a inner-loop cross-validation on training data only for hyperparameter selection. Here we selected SVM because it works better than several other classification methods, including naive Bayes classifier and logistic regression in our local test. Then we measured the binning performance by first computing the portion of misclassified fragments within each species, and taking the mean error rate across all species. This evaluation indicator is less biased when there are over-represented species or species with genomes of extremely imbalanced sizes in the data.

We compared the classification errors of -mer profile method, alignment-based method, uniformly random LSH-based method and our Opal method based on the Gallager LSH. For the -mer profile method, we constructed the 12-mer profile, which is an optimal -mer profile that can be loaded into memory, as shown in a previous work [9] and also in our in-house experiment (see Supplementary Figure S1). For the alignment-based method, we chose BWA-MEM with default settings as suggested in [9]. In a in-house experiment, we also find BWA-MEM outperforms megaBLAST in terms of both speed and accuracy on our dataset. For the uniformly random LSH-based method, we sampled a set of spaced -mers to construct LSH functions, denoted as LSH(64,8). For Opal, we randomly generated 32 Gallager LSH functions, denoted as Opal(64,8). We also constructed hierarchical Gallager LSH functions with the first layer generated by the Gallager design of and the second layer with , denoted as Opal(64,32,8). For all LSH-based methods, we randomly drawn hash functions to construct the compositional feature vectors. Note that the dimensionality of LSH-based feature vector under this setting is , much smaller than , the dimensionality of -mer profile. It is expected that we see better performance if more hashing functions are sampled. For Opal, we used the hierarchical (64,32,8) as the default setting.

4.2 Comparison to previous compositional-based methods

Figure 2: Performance comparison of metagenomic binning methods.

Opal and LSH-based methods outperforms traditional -mer profile. Due to the limit of computational resources, the -mer profile method fails to handle very large , while the LSH-based method can capture long range information in a very long -mer. We compared the performance between -mer profile and LSH-based methods. The results are shown in Figure 2. Even with a smaller dimensionality, LSH-based methods, including LSH(64,8), Opal(64,8), Opal(64,32,8), achieve almost the same performance at mutation/sequencing error rate 0. As the mutation/sequencing error rate increases, the performance of the -mer profile method drops dramatically. The -mer profile method suffers severely from the mutations/sequencing errors and its binning error increases to 21.14% and 9.77% with mutation/sequencing error rate 15%, for and , respectively. The uniformly random LSH-based method, however, shows the robustness to the mutations/sequencing errors. Its misclassification rate increases slightly, with only around 5% for and 3% for , while the misclassification rate of the -mer profile method increases by 20% and 7%, respectively.

The Gallager LSH is more efficient than the uniformly random LSH. Next, we showed the comparison between the Gallager low-density LSH and the uniformly random LSH. The Gallager LSH has an advantage over random LSH in that each position in the contiguous -mer gets equal coverage when the spaced -mers are generated, hence more efficient when the number of hash functions sampled is practically manageable (see Supplementary Figure S2). We demonstrated this advantage by comparing the classification error of spaced -mer generated by LDPC and random LSH, respectively. We observed that the classification error of Opal(64,8) is consistently lower than that of LSH(64,8). For example, when and mutation/sequecing error is 15%, spaced -mer generated by random LSH has a classification error of 8.31%, while Opal method benefits from the equal coverage and gives an classification error of 7.81%. In addition, the hierarchical -LDPC further reduces the classification error to 7.07%, due a better structured manner of generating LSHs for very long -mers. All these comparisons are statistically significant by paired t-test.

These results demonstrate that Opal, based on the Gallager LSH method, is able to capture very long-range compositional dependency for long metagenomic fragments and shows substantial improvements over traditional -mer profile methods. With the efficient design of the Gallager LDPC code which was designed for error correction in noisy communication, we observe a similar effect in metagenomic binning, that is high robustness to the errors in the data. These observations indicate the practical applicability of Opal to large-scale and noisy metagenomic sequencing datasets from environmental samples.

4.3 Comparison to the alignment-based method

We also compared Opal with BWA-MEM, a fully optimized read mapper. For , we observed that Opal works remarkably better than BWA-MEM, especially when the mutation/sequencing error rate is high. For , the comparison is similar but the gap between Opal and BWA-MEM is smaller because the alignment-based method would benefit from the longer-range dependency. We have also compared BWA-MEM and Opal on smaller datasets with fewer species. Their performance are comparable as there is less ambiguity for BWA-MEM.

Figure 3: Comparison of classification speed. These figures show the time (in microseconds) required to classify each fragment using BWA-MEM, Opal, and Opal+BWA-MEM.

For scalability, we compared the classification speed of our method and BWA-MEM which is highly optimized for mapping. The experiments are conducted on a machine equipped with an Intel Xeon E5-2650 CPU with 8 cores running at 2.00GHz and 32G of RAM. The classification speed of the Opal shows no variation across different mutation/sequecing error rates. The mean classification time per fragment of Opal is 13.1 and 33.8 microseconds for fragment length and , respectively. However, the classification speed of BWA-MEM increases with the mutation rates. For example, BWA needs 697 microseconds per fragment at mutation/sequecing error rate 0% and 1211 microseconds per fragment at mutation/sequecing error rate 15%, for fragment length . We observed that Opal performs up to near 100 times faster than BWA-MEM. In addition, since Opal only stores the weights for the classifier, its memory usage is only about of BWA-MEM’s usage. We believe that these scalability improvements would be more significant on datasets with a more complete set of reference species and longer fragments.

4.4 Compositional-based binning as “coarse search”

While our compositional-based method has numerous advantages over alignment-based methods, it is not an aligner so that it cannot provide the high-resolution mapping results for users, which could be critical for certain downstream applications, such as comparative analysis. Given that the binning accuracy of Opal is outstanding, we can combine it with the alignment-based method in the compressive genomics manner by first applying it as a “coarse search” step to identify a very small group of putative taxonomic origins that are ranked very top by Opal and then performing the alignment-based method as a “fine search” to find the mapping locations and hopefully further improve the ranking predicted by Opal. In the experiment, we only picked the top 2 predicted species by Opal and see whether we can improve by using BWA-MEM subsequently. Remarkably, we found that BWA-MEM is able to further improve Opal by performing sequence alignment against the top 2 ranked reference genomes. The accuracy is very nearly optimal given the top 2 accuracy by Opal, suggesting that BWA-MEM picks the best possible taxonomic origin from the two putative species. Moreover, since we greatly reduced the search space for BWA-MEM, this integrated approach is almost 20 times faster than original BWA-MEM and also with substantially improved binning performance on noisy data. Moreover it also provides high-resolution mapping results that the compositional-based methods cannot generate. These results indicate that this natural combination of compositional-based and alignment-based methods provides metagenomic binning with high scalability and high accuracy along with high-resolution alignments.

mutation 0% 5% 10% 15%
Opal top 1 0.0399 0.0497 0.0542 0.0707
Opal top 2 0.0305 0.0417 0.0445 0.0621
Opal+BWA-MEM 0.0325 0.0426 0.0466 0.0645
Table 1: Compositional-based binning as “coarse search” for alignment-based method further improves binning accuracy.

5 Discussion

We have presented Opal, a novel compositional-based method for metagenomic binning. By drawing ideas from the Gallager LDPC code from coding theory, we designed a family of efficient and discriminative LSH functions to construct compositional features that capture the long-range dependencies within metagenomic fragments. Our method can also be seen as a dimensionality reduction approach for genomic sequence data, which extends the previously-used -mer profile.

By comparing the Gallager LSH method with the traditional -mer based binning method, we have demonstrated substantial improvement on large metagenomic datasets with high sequencing errors or mutations presented, which is mainly because of the theoretical properties of LSH and the efficient design of Gallager code. Clearly, when only contains and and the LSH functions maps the binary strings to vectors of , our Gallager LSH method degenerates to the original LDPC code. Compared to BWA-MEM, a state-of-the-art alignment-based method, Opal achieves comparable performance when the fragments have a small number of sequencing errors or mutations and performs much better and more robust in the presence of high sequencing errors or mutations. Moreover, our method is up to two orders of magnitude faster than the alignment-based method, indicating its practical advantage on very large metagenomic datasets. Finally, we also have demonstrated that by using Opal as a “coarse search” step to identify a small candidate set of taxonomic origins for a subsequent alignment-based method, we are able to provide metagenomic binning with high scalability and high accuracy along with high-resolution alignments. Overall, Opal enables us to perform accurate metagenomic analysis for very large metagenomic studies with greatly reduced computational cost.

In the future, we plan to further explore the improvements in metagenomic binning with Opal. For example, here we only used the simplest linear multiclass SVM, which is agnostic to the structure of taxonomy and can only provide predictions on the species level. We believe that with a structured SVM or other hierarchical classification algorithms, we will be able to perform binning on different taxonomic levels phylogenetically and even provide insights for new species or clades [14, 15, 10, 8]. In addition, we will compare our method to the recent indexing or nearest search-based methods, such as [14, 15, 10] on larger real-world datasets and expect to see further application of the Gallager LSH method for fast similarity search-based binning. Finally, we hope to find a better way to integrate the compositional-based and alignment-based binning methods. For example, in a “compressive genomics” manner, we can further devise the compositional-based method to handle compressed metagenomic fragments with high sequence similarity, as the LSH functions are theoretically capable for efficient coding for similar sequences. Also we hope to investigate principled guidances on how to use Opal as the “coarse search” to better suit the subsequent alignment-based method.


  • [1] Victor Kunin, Alex Copeland, Alla Lapidus, Konstantinos Mavromatis, and Philip Hugenholtz. A bioinformatician’s guide to metagenomics. Microbiology and molecular biology reviews, 2008.
  • [2] Benjamin Buchfink, Chao Xie, and Daniel H Huson. Fast and sensitive protein alignment using diamond. Nature methods, 12:59–60, 2015.
  • [3] Y. William Yu, Noah M. Daniels, David Christian Danko, and Bonnie Berger. Entropy-scaling search of massive biological data. Cell Systems, 2:130–140, 2015.
  • [4] Heng Li. Aligning sequence reads, clone sequences and assembly contigs with bwa-mem. arXiv preprints, arXiv:1303.3997, 2013.
  • [5] Ben Langmead and Steven Salzberg. Fast gapped-read alignment with bowtie 2. Nature methods, 9:357–359, 2012.
  • [6] Stephen F Altschul, Warren Gish, Webb Miller, Eugene W Myers, and David J Lipman. Basic local alignment search tool. Journal of molecular biology, 215:403–410, 1990.
  • [7] Qiong Wang, George M Garrity, James M Tiedje, and James R Cole. Naive bayesian classifier for rapid assignment of rrna sequences into the new bacterial taxonomy. Applied and environmental microbiology, 73(16):5261–5267, 2007.
  • [8] Kaustubh R Patil, Peter Haider, Phillip B Pope, Peter J Turnbaugh, Mark Morrison, Tobias Scheffer, and Alice C McHardy. Taxonomic metagenome sequence assignment with structured output models. Nature methods, 8(3):191–192, 2011.
  • [9] Kévin Vervier, Pierre Mahé, Maud Tournoud, Jean-Baptiste Veyrieras, and Jean-Philippe Vert. Large-scale machine learning for metagenomics sequence classification. arXiv preprint arXiv:1505.06915, 2015.
  • [10] Arthur Brady and Steven L Salzberg. Phymm and phymmbl: metagenomic phylogenetic classification with interpolated markov models. Nature methods, 6:673–676, 2009.
  • [11] Alice Carolyn McHardy, Héctor García Martín, Aristotelis Tsirigos, Philip Hugenholtz, and Isidore Rigoutsos. Accurate phylogenetic classification of variable-length dna fragments. Nature methods, 4(1):63–72, 2007.
  • [12] Rachid Ounit, Steve Wanamaker, Timothy J Close, and Stefano Lonardi. CLARK: fast and accurate classification of metagenomic and genomic sequences using discriminative k-mers. BMC Genomics, 16:236, 2015.
  • [13] Karel Břinda, Maciej Sykulski, and Gregory Kucherov. Spaced seeds improve k-mer-based metagenomic classification. Bioinformatics, 2015.
  • [14] Sasha K. Ames, David A. Hysom, Shea N. Gardner, G. Scott Lloyd, Maya B. Gokhale, and Jonathan E. Allen. Scalable metagenomic taxonomy classification using a reference genome database. Bioinformatics, 2014.
  • [15] Derrick E Wood and Steven L Salzberg. Kraken: ultrafast metagenomic sequence classification using exact alignments. Genome Biology, 2014.
  • [16] Robert Gallager. Low-density parity-check codes. IEEE Transactions on Information Theory, 8(1):21–28, 1962.
  • [17] David MacKay and Radford Neal. Near shannon limit performance of low density parity check codes. Electronics Letters, 32:1645–1646, 1996.
  • [18] Alexandr Andoni and Piotr Indyk. Near-optimal hashing algorithms for approximate nearest neighbor in high dimension. Foundations of Computer Science, 2006.
  • [19] Jeremy Buhler. Efficient large-scale sequence comparison by locality-sensitive hashing. Bioinformatics, 17(5):419–429, 2001.
  • [20] Konstantin Berlin, Sergey Koren, Chen-Shan Chin, James P Drake, Jane M Landolin, and Adam M Phillippy. Assembling large genomes with single-molecule sequencing and locality-sensitive hashing. Nature Biotechnology, 33:623–630, 2015.
  • [21] Daniel H. Huson, Suparna Mitra, Hans-Joachim Ruscheweyh, Nico Weber, and Stephan C. Schuster. Integrative analysis of environmental sequences using megan4. Genome research, 21(9):1552–1560, 2011.

Supplementary Information

Acetobacter pasteurianus Listeria monocytogenes
Acinetobacter baumannii Methylobacterium extorquens
Bacillus amyloliquefaciens Mycobacterium tuberculosis
Bacillus anthracis Mycoplasma fermentans
Bacillus subtilis Mycoplasma genitalium
Bacillus thuringiensis Mycoplasma mycoides
Bifidobacterium bifidum Mycoplasma pneumoniae
Bifidobacterium longum Neisseria gonorrhoeae
Borrelia burgdorferi Propionibacterium acnes
Brucella abortus Pseudomonas aeruginosa
Brucella melitensis Pseudomonas stutzeri
Buchnera aphidicola Ralstonia solanacearum
Burkholderia mallei Rickettsia rickettsii
Burkholderia pseudomallei Shigella flexneri
Campylobacter jejuni Staphylococcus aureus
Corynebacterium pseudotuberculosis Streptococcus agalactiae
Corynebacterium ulcerans Streptococcus equi
Coxiella burnetii Streptococcus mutans
Desulfovibrio vulgaris Streptococcus pneumoniae
Enterobacter cloacae Streptococcus thermophilus
Escherichia coli Thermus thermophilus
Francisella tularensis Treponema pallidum
Helicobacter pylori Yersinia enterocolitica
Legionella pneumophila Yersinia pestis
Leptospira interrogans Yersinia pseudotuberculosis
Table S1: List of names of the 50 microbial species used in this work.
1:Input: , ,
2: all zero matrix
3:for  to  do
4:     for  to  do
6:     end for
7:end for
8:choose uniform random permutation matrix , for .
10:Output: Gallager’s LDPC Matrix
Algorithm S1 Gallager’s LDPC Matrix
1:Input: Gallager’s LDPC Matrix
3:     for  to  do
4:         for  to  do
5:              if  then check if 4-cycle exists
6:                   row index of the first same element in and .
8:                  swap the elements of and that belong to the -th block.
9:              end if
10:         end for
11:     end for
12:until no 4-cycle
13:Output: 4-Cycle-free Gallager’s LDPC Matrix
Algorithm S2 Removing 4-Cycles
Figure S1: Performance of -mer profiles with different values on the 50-species dataset with .
Figure S2: Comparison of position coverage between LSH and Gallager LSH. The left figure shows how many times a position being covered when spaced -mers are generated by uniformly sampled LSH and Gallager LSH functions, respectively. The right figure compares the relationship between the number of hash functions and the classification performance. The experiment is done on a 20-species dataset.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description