MetaBAT, an efficient tool for accurately reconstructing single genomes from complex microbial communities

An overview of MetaBAT software and its probabilistic models

As a pre-requisite for binning, the user must create BAM files by aligning the reads of each sample separately to the assembled metagenome (Fig. 1 steps from 1 to 3). MetaBAT takes an assembly file (fasta format, required) and sorted bam files (one per sample, optional) as inputs. For each pair of contigs in a metagenome assembly, MetaBAT calculates their probabilistic distances based on tetranucleotide frequency (TNF) and abundance (i.e., mean base coverage), then the two distances are integrated into one composite distance. All the pairwise distances form a matrix, which then is supplied to a modified k-medoid clustering algorithm to bin contigs iteratively and exhaustively into genome bins (Fig. 1).

We use tetranucleotide frequency as sequence composition signatures as it has been previously shown that different microbial genomes have distinct TNF biases (Mrazek, 2009; Pride et al., 2003; Saeed, Tang & Halgamuge, 2012; Teeling et al., 2004a). To empirically derive a distance to discriminate TNFs of different genomes, we calculated the likelihood of inter- and intra-species Euclidean distance (Deza, 2012) by using 1,414 unique, complete genome references from NCBI (Fig. 2A). This empirically derived distance is termed Tetranucleotide frequency Distance Probability (TDP).

To evaluate the effect of contig sizes on inter-species distance, we obtained posterior probability distributions of inter-species distance with several fixed sizes and observed better inter-species separation as contig size increases (Fig. 2B). As contigs in real metagenome assemblies have various sizes, we then modeled TDP between contigs of different sizes by fitting a logistic function to reflect the dynamic nature of the non-linear relationship between Euclidean TNF distance and TDP across different contig sizes. The results (Figs. 2C and 2D) suggest that the values of two parameters of the model, b and c are unstable if the size of either contig is very small (<2 kb) and one should be cautious to allow smaller contigs to be binned.

Although contigs originating from the same genome are expected to have similar sequence coverage, i.e., genome abundance, the coverage of contigs can vary significantly within a library due to biases originated from the current sequencing technology (Benjamini & Speed, 2012; Harismendy et al., 2009; Nakamura et al., 2011; Ross et al., 2013). As illustrated in Fig. 2E, the observed coverage variance derived from data consisting of isolate genome sequencing projects (total 99 from IMG Database (Markowitz et al., 2012), henceforth referred as the IMG dataset) significantly deviate from the theoretical Poisson distribution, consistent with the notion that both the variance and the mean should be modelled (Clark et al., 2013). For computational convenience, we chose the normal distribution as an approximation since it fits the observation much better (Fig. 2E). To compute the abundance distance of two contigs in one sample, we use the area not shared by their inferred normal distributions with given coverage mean and variance (Fig. 2F). A geometric mean of the distances for all samples is used for the final abundance distance probability (ADP) of two contigs. In addition, we applied a progressive weighting mechanism to adjust the relative strength of the information from abundance distance, meaning that we put more weight on abundance distance when it was calculated from many samples (see below).

We then integrate TDP and ADP of each contig pair as the following: ${\displaystyle P\left({\mu }_1,{\sigma }_1^2,{\mu }_2,{\sigma }_2^2\right)=\left\{\begin{array}{ll}max\left(\mathrm{TDP},\mathrm{ADP}\right),\hfill & \text{if }\mathrm{TDP}>0.05\hfill \\ \mathrm{ADP}\cdot w+\mathrm{TDP}\cdot \left(1-w\right),\hfill & \text{otherwise}\hfill \end{array}\right.}$ where w = min[log(n + 1)/log(m + 1), α] ⋅ n, m, and α represent the number of samples, a large number (100 as the default), and the maximum weight of ADP (0.9 as the default), respectively. For instance, in the default setting, the weight would be about 0.5 when the number of samples is 10 and TDP is less than 0.05. The resulting distance matrix is used for binning (see below).

Tetranucleotide frequency probability distance (TDP)

To establish empirical probabilities of intra- and inter-species for tetranucleotide frequency distance, we downloaded 1,414 unique, completed bacterial genomes from the NCBI database and shredded them into fragments ranging from 2.5 kb to 500 kb. Next, we obtained 1 billion random contig pairs from within or between genomes. The empirical posterior probability that two contigs are from different genomes is given as the following: ${\displaystyle P\left(T|D\right)=\frac{P\left(T\right)P\left(D|T\right)}{P\left(T\right)P\left(D|T\right)+P\left(R\right)P\left(D|R\right)}}$ where T or R represent cases where two contigs are from different (inter) or the same (intra) species, respectively. D is the Euclidean TNF distance between two contigs. The same uninformative priors of T and R were chosen. In reality, P(T) is expected to be much bigger than P(R), thus we set P(T) = 10∗P(R) as the default implementation to adjust for the possible under-sampling issue in inter species distance.

The TDP of contig pairs with different sizes is approximated using logistic regression: ${\displaystyle P\left(D_{ij};b_{ij},c_{ij}\right)=\frac{1}{1+e^{-\left(b_{ij}+c_{ij}\ast D_{ij}\right)}}}$ where D_ij represents a Euclidean TNF distance between contig i and j ⋅ b and c, the two parameters for the logistic regression, are estimated from the empirical data.

Abundance distance probability (ADP)

The probabilistic abundance distance was calculated as follows: Suppose two contigs have the mean coverage of μ₁ and μ₂ and the variances of ${\sigma }_1^2$ and ${\sigma }_2^2$, then we defined the abundance distance as the non-shared area of two normal distributions of N(μ₁, ${\sigma }_1^2\left.\right)$ and N(μ₂, ${\sigma }_2^2\left.\right)$: ${\displaystyle P\left({\mu }_1,{\sigma }_1^2,{\mu }_2,{\sigma }_2^2\right)=\frac{1}{2}\int |{\varphi }_{{\mu }_1,{\sigma }_1^2}-{\varphi }_{{\mu }_2,{\sigma }_2^2}|}$ where ϕ represents a normal distribution having two parameters μ and σ². Numerically this can be simplified using cumulative distribution functions as follows assuming ${\sigma }_2^2$ is greater than or equal to ${\sigma }_1^2$: ${\displaystyle P\left({\mu }_1,{\sigma }_1^2,{\mu }_2,{\sigma }_2^2\right)=\left\{\begin{array}{ll}{\Phi }_{{\mu }_1,{\sigma }_1^2}\left(k_0\right)-{\Phi }_{{\mu }_2,{\sigma }_2^2}\left(k_0\right),\hfill & \text{if }{\sigma }_1^2={\sigma }_2^2\hfill \\ {\Phi }_{{\mu }_1,{\sigma }_1^2}\left(k_2\right)-{\Phi }_{{\mu }_1,{\sigma }_1^2}\left(k_1\right)+{\Phi }_{{\mu }_2,{\sigma }_2^2}\left(k_1\right)-{\Phi }_{{\mu }_2,{\sigma }_2^2}\left(k_2\right),\hfill & \text{otherwise}\hfill \end{array}\right.}$ where Φ represents a cumulative normal distribution, and ${\displaystyle k_0=\frac{{\mu }_1+{\mu }_2}{2}}$ ${\displaystyle k_1^{\ast }=\frac{\sqrt{{\sigma }_1^2\cdot {\sigma }_2^2\cdot \left({\left({\mu }_1-{\mu }_2\right)}^2-2\cdot \left({\sigma }_1^2-{\sigma }_2^2\right)\cdot log\left({\sigma }_2/{\sigma }_1\right)\right)}-{\mu }_1\cdot {\sigma }_2^2+{\mu }_2\cdot {\sigma }_1^2}{{\sigma }_1^2-{\sigma }_2^2}}$ ${\displaystyle k_2^{\ast }=\frac{\sqrt{{\sigma }_1^2\cdot {\sigma }_2^2\cdot \left({\left({\mu }_1-{\mu }_2\right)}^2-2\cdot \left({\sigma }_1^2-{\sigma }_2^2\right)\cdot log\left({\sigma }_2/{\sigma }_1\right)\right)}+{\mu }_1\cdot {\sigma }_2^2-{\mu }_2\cdot {\sigma }_1^2}{{\sigma }_1^2-{\sigma }_2^2}}$ ${\displaystyle k_1=min\left(k_1^{\ast },k_2^{\ast }\right)\text{and}{k}_2=max\left(k_1^{\ast },k_2^{\ast }\right).}$

To combine multiple abundance probabilities across different samples, we calculated the geometric mean of probabilities: ${\displaystyle P_{ij}=\sum _n\mathbb{1}\left\{{\mu }_{in}>c\hspace{0.17em}\text{OR}\hspace{0.17em}{\mu }_{jn}>c\right\}\sqrt{\prod _n{P}_{ijn}\left({\mu }_{in},{\sigma }_{in}^2,{\mu }_{jn},{\sigma }_{jn}^2\right)\ast \mathbb{1}\left\{{\mu }_{in}>c\hspace{0.17em}\text{OR}\hspace{0.17em}{\mu }_{jn}>c\right\}}}$ where P_ijn represents the probability calculated from two abundances μ_in and μ_jn, and c represents a cut-off for reasonable minimum abundance for a contig.

As metagenome assemblies contain many small contigs, whether or not to include them is a dilemma; including small contigs will likely improve the genome completeness, at a cost of genome quality because their larger abundance variations make it harder to bin them correctly. We tried to empirically determine a reasonable contig size cut-off by plotting the ratio of mean and variance from the IMG single genome dataset. Although most genome variances are much larger than their means, their ratio becomes stabilized after contig size increases to 2.5 kb (Fig. S1). Therefore, we used 2.5 kb or larger contigs for the initial binning. Smaller contigs can be recruited after the binning, based on their correlation to the bins (Imelfort et al., 2014).

Iterative binning

We modified the k-medoid clustering algorithm (Kaufman & Rousseeuw, 1987) to eliminate the need to input a value for k and to reduce search space for efficient binning. Specifically, the binning algorithm works as follows:

1. Find a seed contig (e.g., having the greatest coverage), and set it as the initial medoid.

2. Recruit all other contigs within a cutoff distance (i.e., parameters p1 and p2) to the seed.

3. Find a new medoid out of all member contigs.

4. Repeat 2–3 until there are no further updates to the medoid. These contigs form a bin.

5. For the rest of the contigs, repeat 1–4 to form more bins until no contigs are left.

6. Keep large bins (e.g., >200 kb), and dissolve all other bins into free contigs.

7. (optional) For dataset with at least 10 samples, recruit additional free contigs to each bin based on their abundance correlation.

Binning performance on “error-free” metagenomic assemblies

A metagenomic dataset (Accession #: ERP000108) from the MetaHIT consortium (henceforth referred to as the MetaHIT dataset) (Qin et al., 2010) was chosen to benchmark MetaBAT because it contains a large number of samples and the community contains many species with reference genomes. To derive a reference genome set, we selected 290 known genomes from NCBI that are present in MetaHIT at >5× mean coverage (Table S1). These reference genomes were then shredded into contigs of random sizes (>2.5 kb) following an exponential distribution modeled to mimic real metagenome assemblies. The abundance of each contig in every sample was also obtained using real data. These “error-free” metagenome contigs, their abundance information, along with their parental reference genomes (as “true answers”), were used in the following analysis to benchmark binning performance. For a full description of the experiment, refer to the MetaBAT wiki page: https://bitbucket.org/berkeleylab/metabat/wiki/Home.

For comparison, we ran several alternative binning tools on the same dataset described above. These software include Canopy (Nielsen et al., 2014), CONCOCT v.0.4.0 (Alneberg et al., 2014), GroopM v.0.3.0 (Imelfort et al., 2014), and MaxBin v.1.4.1 (Wu et al., 2014). Among them, CONCOCT, GroopM, and MaxBin are also fully automated binning tools. An optional manual step in GroopM for improving the quality of bins was excluded. Since MaxBin does not consider multiple samples, we combined multiple samples into one.

We used >90% precision (lack of contamination) and >30% recall (completeness) as the minimum criteria for a bin to be considered “good” which basically means the bin should be composed of one or more strains of a single species (for results of other thresholds, refer to Figs. S3 and S4). Formulas for this calculation are described in the Supplemental Information 1. Among all binning tools, MetaBAT binned the greatest number of genomes at almost every recall threshold (Fig. 3A). CONCOCT is the only tool that produces more genome bins with over 80% or 90% completeness than MetaBAT, but MetaBAT produces many more bins at 70% completeness threshold. Interestingly, we found these tools complement each other in forming genome bins (Fig. 3B, GroopM was omitted since it detects only a few genomes). Among the unique 133 genome bins collectively formed by all tools, MetaBAT binned the most number of genomes (111, 83.5%), with 23 bins (17.2%) that were not found by any other tool.

MetaBAT runs very efficiently in computation; the entire binning process only took 14 min and 4 GB of RAM (Table 1). Multiple simulations produced almost identical performance results and thus were not shown.

Binning performance on real metagenomic assemblies

We next tested the performance of MetaBAT on real metagenomic assemblies. Using the same raw sequence data from the above MetaHIT dataset, we pooled the sequences from all samples and then assembled them using Ray Meta assembler (Boisvert et al., 2012). Because real metagenomic assemblies often contain many small contigs, we lowered the minimum contig size requirement from 2.5 kb to 1.5 kb to include more contigs into our binning experiment. As a result, 118,025 contigs were used for binning. We then ran the above 5 binning tools with their default settings on this dataset. In contrast to the previous simulation experiment, in this experiment we do not have a reference genome for every genome bin; we instead used CheckM (Parks et al., 2015) to calculate the approximate recall (percent of expected single-copy-genes that are binned) and precision (the absence of genes from different genomes) rates. A full description of this experiment is available of MetaBAT wiki page: https://bitbucket.org/berkeleylab/metabat/wiki/Benchmark_MetaHIT.

Similar to the previous “error-free” experiment, MetaBAT again identified the greatest number of unique genome bins having >90% precision (Fig. 4A). In this experiment with real metagenomic contigs, the superior completeness we saw in CONCOCT during the “error-free” experiment was lost. Moreover, the number of genome bins formed by MetaBAT was consistently greater than the others at every completeness threshold. Similarly, different tools produced complementary binning results as before (Fig. 4B). MetaBAT’s contribution appears to be more pronounced this time. It missed 17 bins formed by all other tools combined, but recovered 31 bins that no other tools produced. MetaBAT alone recovered 90.2% (133/144) of genome bins from all tools. These results suggest MetaBAT is very robust when run against a real metagenome assembly. Consistent with the simulation experiment, MetaBAT is computationally very efficient and requires only 4 min to complete this experiment (Table 2).

To test the performance of MetaBAT on large-scale metagenomic data sets, we used a dataset containing 1,704 (with replicates) human gut microbiome samples (Accession #: ERP002061 (Nielsen et al., 2014). The entire dataset was assembled using Ray Meta assembler (Boisvert et al., 2012) and Megahit (Li et al., 2015) resulting in 1,058,952 contigs (>1 kb) that were then used for binning. MetaBAT took less than 2 h to generate 1,634 genome bins (>200 kb) using a single node with 16 CPU cores (32 hyper-threads), and the peak memory consumption was at 17G. In comparison, Canopy took 18 h using 36G memory in the same setting. The other binning tools—CONCOCT, GroopM, and MaxBin—all failed to generate any genome bins for this data set likely due to their inability to scale. For accuracy evaluation, we used CheckM (Parks et al., 2015) and identified 610 high quality bins (>90% precision and >50% completeness) among the bins predicted by MetaBAT, which is 35% more than the published CAG bins (Nielsen et al., 2014) and 11% more than Canopy bins using our assembly. For details on the use of MetaBAT on this large dataset please refer to: https://bitbucket.org/berkeleylab/metabat/wiki/Example_Large_Data.

Post-binning processing to further improve quality

Assembly from pooled samples in the above experiment raises the possibility that similar genomes (e.g., different strains) present in different samples are assembled into chimeric contigs. This level of contamination may not be tolerated in some applications. Based on the assumption that a single sample will contain fewer strains of the same species than all pooled samples, we implemented an optional post-binning process to reduce the strain-level contamination. Briefly, we first selected a single sample with the most reads mapped to a specific bin, and then assembled these reads into a new set of contigs. If the new contigs produces better CheckM results, we subsequently replaced the old contigs in this bin with the new ones.

This post-processing step significantly improved both completeness and precision (lack of contamination) for 61% (992 out of 1,634) of the genome bins. Overall, there were 571 bins with >95% precision and >50% completeness, compared with 375 without post-processing. This improvement was more obvious when we increased the precision threshold to >99%, as the number high quality bins increases from 46 to 186.

By incorporating additional sequencing data and other post-binning optimizations, Nielsen et al. (2014) generated 373 high quality draft genomes (“MGS genomes”). We therefore used these MGS draft genomes as reference for additional quality assessment of the MetaBAT genome bins after post-processing. As shown in Fig. 5A, 31 MGS draft genomes were not well represented by MetaBAT bins, but MetaBAT recovered 55 additional genome bins not reported by MGS draft genomes. For those overlapping bins, most MetaBAT bins closely approximate the MGS draft genomes in accuracy—94% precision and 82% completeness (Fig. 5B).