Deconvolute individual genomes from metagenome sequences through short read clustering

Clustering strategies

An overview of the two clustering strategies is shown in Fig. 1. The local clustering strategy in the original SpaRC has been described in Shi et al. (2018). In brief, during the local clustering step, we cluster reads by their overlap. The read clusters are further clustered into bigger clusters by the global clustering strategy, which we will describe in detail below.

Local clustering improvement with minimizers

In SpaRC, the number of shared k-mers is used to estimate similarity between reads (Shi et al., 2018). As it takes 100–200 times more space after reads are transformed into k-mers and edges, SpaRC is neither space nor time-efficient. To improve computing efficiency, we implemented a new function to use minimizers (Roberts et al., 2004) instead of k-mers to estimate similarity between reads. As many adjacent shared k-mers can be represented by a single minimizer without losing sensitivity, in theory, the minimizer-based method should greatly reduce the memory requirement in SpaRC (as fewer k-mers and edges will be produced). In practice we did observe a 3.2-fold memory usage reduction, and 3.3-fold speed-up (Fig. S1). It is worth noting that minimizers may not be applicable to uncorrected long reads from PacBio and Nanopore sequencing technologies due to their high error rate.

Global clustering

Reads from a genome can form many read clusters after the local clustering step, leading to low clustering completeness. The ultimate goal of global clustering is to predict all the read clusters originated from the same genome. It does this based on the assumption that the sequencing coverage of each region of a genome, defined by the read clusters, closely resembles the sequencing coverage of the same genome across different metagenomic samples. In other words, if two clusters, c1 and c2, belong to the same genome g. After c1 and c2 are assembled into contigs C1 and C2, the coverage of C1 and C2 in sample S, in theory, should be equal to the coverage of g.

Datasets

Computing environments

Read clustering experiments were performed on Amazon Web Service (AWS)’s Elastic MapReduce (EMR, emr-5.17.0). Depending on the size of the dataset, a number of r4.2×large instances were used to form a cluster, and the configuration details are shown in Table 2.

Global clustering greatly improves short read clustering performance

In order to test whether multiple samples derived from the same microbial community could be leveraged to improve short read clustering performance, we designed a control dataset by taking 10% of the reads from 50 samples from the CAMI2 synthetic metagenome dataset (Materials and Methods). The smaller dataset allowed us to reduce computation cost and obtain results faster. We clustered the reads using two clustering methods: in the “local clustering” method, we combined all the reads from the 50 samples for clustering and selected clusters with 50 reads or larger; in the “global clustering” method, we further applied the global clustering module to these clusters to form big clusters (Material and Methods). This labeled synthetic dataset enabled us to systematically compare the clustering performance between these two methods for cluster size, purity, completeness. The results are shown in Fig. 2. We used the same purity and completeness metrics as in Shi et al. (2018). Briefly, the purity of a cluster is defined as the percentage of reads from the predominant genome within the cluster, while the completeness of a cluster is defined as the percentage of all the reads from the predominant genome that are captured by the cluster. Because almost identical strains from the same species were engineered in the dataset, both of the two metrics, especially purity, likely underestimate species-level clustering performance. For example, if a species has two closely related strains with equal number of reads, clusters derived from this species will have a purity of 50%. We therefore also measured cluster purity at the species level. In this experiment, the parameters were k = 41, m = 22, min_shared_kmers=2, max_degree=25, representative k-mer count = 100, and cosine threshold=0.925.

The local clustering step resulted in 78.29% of the reads clustered into many small clusters (n = 378,829), with the largest cluster having only 10,151 reads. The majority of the clustered reads 99.08% are distributed in clusters with 1,000 reads or less (Fig. 2A). In contrast, after the global clustering step the number of clusters is significantly reduced (n = 10,083), 78.49% of reads are in clusters with 1,000,000 or more reads, and the largest cluster contains 6,836,687 reads. Consequently, the median completeness of the clusters from the global clustering is also 19.58 times larger than that from the local clustering (19.97% vs 1.02%, Fig. 2B). The increase in completeness came at the expense of some purity loss (median purity from 97.95% to 93.25% at the genome level%), but most of the clusters are still pure , especially at the species level (median purity decreased slightly, from 100% to 99.78%) (Fig. 2B).

As the success of genome clustering should heavily depend on its sequencing coverage, we next explored the relationship of cluster completeness as a function of genome sequencing coverage (Fig. 2C). For local clustering, higher sequencing coverage seems to have little effect on cluster completeness. In contrast, higher sequencing coverage does translate into higher completeness, suggesting global clustering effectively leverages multiple sample statistics for read clustering. After the sequencing coverage reaches a sufficient threshold (100×, Fig. 2C), the completeness of most genomes exceeds 80%.

The performance of short read clustering can be improved by increasing the number of samples

We next explored the relationship between the number of samples in a dataset and the global clustering performance. Intuitively, more samples should enable a more robust estimation of the similarity between clusters and lead to better clustering performance. By limiting the total size of the datasets to 25 Gb to reduce the computation cost, we made several datasets with varying number of samples (5, 10, 20, 50) randomly selected from the CAMI2 synthetic metagenome dataset (Materials and Methods). We obtained clusters from these datasets by running SpaRC with the same parameters (k = 41, m = 22, min_shared_kmers=2, min_read_per_cluster=50) for local clustering, and representative_kmer_count=100 for global clustering.

As in the previous section, we evaluated the purity and completeness of the resulting clusters. As shown in Fig. 3, the median purity of the clusters at the genome level slightly dropped as the number of samples increases, from the highest 96.96% (at n = 5) to the lowest 88.89% (at n = 20, Fig. 3A). This is likely caused by the fact that more samples may contain more strain variation, and currently SpaRC can not distinguish very similar strains. Consistent with this notion, the purity measured at the species level remains largely unchanged (Table S1). In contrast, the completeness continuously increases with increasing number of samples (median completeness rises from 7.69% at n = 5 to 19.97% at n = 50, Fig. 3B). This result supports the hypothesis that more samples enhance global clustering performance, likely due to better estimation of cluster similarity.

The ultimate goal of read clustering is to recover the complete set of any genome without any contamination from other genomes. In order to measure how many genomes can be recovered by read clustering, here we define “a clustered genome” as a read cluster that simultaneously satisfies two criteria: purity >95% and completeness >80%. It is worth noting that these criteria are very strict. As strain-level heterogeneity can greatly reduce purity, and one small region larger than a read with no sequencing coverage, either due to statistical sampling or systematic sequencing biases, will greatly decrease completeness.

The clustered genomes from 5, 10, 20 and 50 samples are 0, 1, 5 and 7, respectively (Fig. 3C). Since the ability to obtain a clustered genome depends on the sequencing coverage, and there are 20 genomes with at least 100× coverage, this translates to a recovery rate of 35% with 50 samples. Other genomes with lower sequencing coverage also benefit from more samples included in clustering, as Fig. 3C shows the extent of recovery for all genomes with a sequence coverage >10×. The details of these 7 recovered genome can be found in Table S1.

Parameters that may impact clustering performance

There are 10 parameters in the local clustering algorithm (Shi et al., 2018). In the global clustering step two more parameters are added, where the number of representative k-mers used to estimate cluster abundance (rp) and cosine similarity threshold (cs) to control graph complexity. While some of these parameters only affect computing efficiency, there are four parameters that in theory may affect clustering accuracy: k-mer(k)/minimizer(m) length, min_shared_kmers among reads (minsk), rp and cs. In theory higher k, minsk, and cs all lead to smaller clusters with lower completeness but higher purity, and vice versa. Higher rp should make the estimation of cluster abundance more accurate with a small cost in computing efficiency. To explore the effect of these parameters on clustering accuracy, we ran SpaRC with different sets of parameters on the 25 Gb dataset with 50 samples. In each of the parameter sets, only one parameter varies while the other three were held constant. The rest of the parameters were used as their default values. We used several metrics including number of reads clustered (bigger is better), number of clusters formed (bigger is worse), median cluster purity, median cluster completeness and number of clustered genomes to measure clustering performance. The results are shown in Table 3.

Overall, there is no single best parameter that can maximize all metrics. In general, the most abundant genomes (number of clustered genomes) are less likely to be affected by these parameters except a few extreme cases (where k-mers are set too small or cluster similarity thresholds are too low). In those extreme cases over-clustering happened, as only a small number of clusters formed with very low median purity. These parameters, except the number of representative k-mers, can greatly affect the median purity and completeness. Longer k-mers and requiring more shared k-mers among reads increases median purity and completeness. These improvements are likely driven by better clusters from the genomes with medium to high sequencing coverage. If these parameters became very large, the number of reads clustered dramatically decreases (under-clustering), but these genomes do not seem to be affected. For example, the number of clustered reads drops from 60,953,377 at minsk=1 to only 26,643,712 at minsk=5, but the median completeness and purity reaches their peak. The un-clustered reads at high minsk presumably are derived from a lot of genomes with low sequencing coverage, as bigger k-mer size decreases the correct k-mers present in the dataset and makes their counts more noisy (Chikhi & Medvedev, 2013).

The number of representative k-mers used for cluster abundance estimation seems to have a very minor effect on clustering performance. If the number is too low, then the performance is slightly lower.

As expected, very low cluster similarity thresholds cause over-clustering, and very high ones lead to under-clustering. Same as the other parameters, the best parameter choice should be determined by the underlying scientific requirements to balance sensitivity and specificity.

Assembly quality comparison between SpaRC-based and the classic approach

Compared to the current metagenome assembly strategy without clustering, that is, assemble the entire dataset followed by binning (hereafter referred as “the classical approach”), clustering the reads into individual clusters by SpaRC followed by assembly and binning (hereafter referred as “SpaRC-based approach”) may produce better results. To test this hypothesis, we carried out both approaches on the above CAMI2 testing dataset (Methods) and compared their assembly results.

As shown above, the effectiveness of clustering depends on sequencing coverage, we therefore evaluated the assembly performance for genomes with coverage of 100× or above, 50–100× and 30–50×, respectively. We observed a comparable number of mis-assemblies from both approaches (Table S2), so we focused on the following four metrics: genome coverage (percent of reference covered by the assembled contigs), contamination (percent of contigs not belonging to the reference at the species level), N50 and L50 (measuring contiguity of the assembly). Specially, we only reported genome bins with at least 95% purity, or less than 5% contigs from other species. We evaluated their assembly accuracy in terms of near-complete genomes (95% genome coverage and above) and fairly complete genomes (80–95% genome coverage). The results are shown in Fig. 4.

At high sequencing coverage (100× or above), the SpaRC-based clustering approach was able to recover more near-complete genomes than the classical approach, 16 for SpaRC vs 11 for Classic (Fig. 4A). Among them, the SpaRC approach assembled 6 genomes that were not completely assembled by the classical approach, while missed only one assembled by the classical approach (Fig. 4B). The number of fairly complete genomes assembled by the two approaches are comparable at this sequencing coverage. Besides, 11 genomes assembled by the classical approach have smaller L50s and larger N50s while 6 genomes assembled by SpaRC approach do (Table S2). As sequencing coverage is getting lower, the classical approach has an advantage over the SpaRC-based approach, especially when coverage drops to below 50× (Fig. 4A). These genomes tend to spread across many small pure clusters. This result is consistent with the above clustering performance analyses, suggesting there is still an “under-clustering” problem for global clustering.

We added a binning step using MetaBAT after the assembly of each of the clusters because some of the large ones are mixtures of a few genomes. Clusters with more than 100,000 reads contain 3.125 genomes on average, suggesting the species complexity of these clusters are much more reduced comparing to the original metagenome. The metrics presented in Fig. 4 were calculated after the binning step.

To test whether these conclusions can be generalized to real-world datasets, we applied the two approaches to a human microbiome dataset (MetaHIT, “Methods”). As there is no ground truth for this dataset, we mapped the metagenomic short reads to NCBI non-redundant reference genome database, and used the 68 genomes with 30× or more sequencing coverage as a reference set.

Similar to the previous experiment on CAMI2 dataset, the SpaRC-based approach can cluster genomes with 100× sequencing coverage from MetaHIT dataset, and this ability decreases dramatically as coverage decreases. The SpaRC-based clustering approach was able to recover 9 near complete genomes with 100× or more coverage (total 16), while the classic method recovered 10 (Fig. 4C). Among them, SpaRC approach recovered 3 genomes that were not completely recovered by the classical approach, while missed 4 genomes (Fig. 4D). For a complete list of the recovered genomes from these two approaches on the two datasets, please refer to Table S3.

The above results suggest that adding a clustering step can complement the common metagenome assembly strategy, at least for the species with high sequencing coverage.

Execution time of SpaRC on the MetaHIT dataset

Assembling large metagenome datasets using MEGAHIT requires a large amount of RAM on the server. It took 24.63 h on a single node with 64 CPU cores and 488 GB RAM (r4.16×large, a Memory-Optimized instance of AWS Elastic Compute Cloud) to assemble the MetaHIT dataset. As for SpaRC-based approach, we were able to distribute the read clustering step on an AWS Elastic MapReduce (EMR) cluster with 350 nodes, each with 8 CPU cores and 61 GB of RAM (Table 1). The clustering step took 8.9 h to complete. These results suggest that SpaRC is not as good as the classic approach in terms of costs and computational efficiency, instead it offers an advantage to scale up to bigger datasets (over assemblers only run on single nodes), and overall shorter computational execution time.