Xander: employing a novel method for efficient gene-targeted metagenomic assembly

Metagenomics faces scalability challenges stemming from the amount of raw sequencing data necessary to describe complex microbial communities, now often termed the microbiome [1, 2]. Metagenomic assembly has been an area of growing interest in the past decade, with early datasets assembled using single-genome assembly methods that had difficulty with metagenomic samples [3, 4]. The tendency for assemblers to only assemble dominant organisms while producing only fragmented partial assemblies for less dominant organisms, with limited recovery for any individual genes, has been an impetus to develop better metagenomic-specific assembly methods [5, 6].

We propose a gene-targeted assembly approach called Xander for assembling metagenomic datasets. Xander is a de Bruijn graph [7] assembler [8] that uses external information to perform a guided, instead of exhaustive, traversal of the assembly graph. Xander uses profile hidden Markov models (HMMs) [9] to guide graph traversal (HMM-guided assembly). An HMM can be considered as a directed probabilistic graph. Built from the alignment of a set of homologous sequences, an HMM quantifies position-specific conservations of the underlying sequence family. Using an HMM, the paths most likely to code for the target gene can be extended first thus limiting the portion of the assembly graph that must be explored. In addition to limiting the graph traversal, the HMM provides a measure of how likely the resulting assembled contig comes from the supplied model and, by inference, whether the contig actually codes for the gene of interest. Xander enables researchers to target specific genes involved in biologically interesting pathways without using amplicon-based sequencing approaches that are prone to primer limitations and PCR artifacts and that produce sequences too short to more fully characterize the genes.

Gene-targeted assembly is less resource intensive and faster than traditional whole-genome metagenomic assembly. In addition to the de Bruijn graph, only small paths relative to the graph’s size must be kept in memory. In our initial implementation, further reductions in memory usage are achieved by using a probabilistic data structure for holding the graph in memory, a Bloom filter [10–12]. By targeting relatively small segments of the assembly graph using an HMM to guide assembly, the amount of the graph that must be explored during assembly is constrained, providing a speedup over whole-genome approaches.

The targeted assembly of metagenomic datasets has drawn research interest, generating several approaches including EMIRGE, Mira, and SAT-Assembler [13–15]. These targeted methods first identify reads likely to be from targeted genes before assembly, as does the reference-assisted assembler IDBA-Hybrid [16]. Instead, Xander begins with a de Bruijn graph representation based on all reads, as most current bulk assemblers do. A de Bruijn graph approach has also been used for a BLAST-like search algorithm, BlastGraph [17]. BlastGraph uses a graph representation of a reference database that can then be queried with unknowns, whereas Xander queries the graph representation of metagenomic reads using HMMs. Other work used an HMM to search sequences in an indexed tree structure [18], but only short HMMs up to length 12 were tested.

We validated our method using defined community metagenomic data and then applied the tool to 21 rhizosphere soil samples totaling 800 giga base pairs (Gbp) from three different biofuel crops. We targeted three genes: one conserved taxonomic gene and two genes from different parts of the nitrogen cycle that are both critical to plant productivity.

Xander uses a novel data structure combining de Bruijn graphs and HMMs to target assembly of important genes from metagenomic data. Using this data structure allows us to apply powerful graph search techniques to assemble individual genes. The HMM-guided assembly reduces the search space and provides gene annotation. Additional optimizations including caching of best paths from previous searches and heuristic search pruning provide several orders of magnitude of speedup over a naïve approach. Xander requires reference sequences from the targeted gene families for HMM training and start finding and is not meant for discovering completely novel gene families. Although we applied this method to soil rhizosphere metagenomic data, this method should also be applicable to de novo assembly of transcriptome and metatranscriptome data.

We targeted three gene families: rplB, a single-copy core phylogenetic marker gene; nifH, a marker for nitrogen fixation; and nirK, a marker for denitrification. There is considerable interest in using biofuel crops as one way to mitigate global warming. One factor in improving lower cost biofuel crop production on non-food-producing lands is to improve plant-available nitrogen. The Kellogg Biological Station has replicate field plots of corn, a traditionally farmed annual biofuel crop, and two perennial biofuel crops, an exotic, Miscanthus, and a native, switchgrass. We used Xander to assemble the three genes from the pooled rhizosphere samples from each crop. Roughly 17,000 to 20,000 rplB genes assembled for each sample (Table 4). We found about 10 % as many nirK contig clusters and even fewer nifH contig clusters. A previous study at other Midwest sites found higher gene richness in nitrous oxide reductase (nosZ), another marker gene for denitrification, in corn soil communities than in switchgrass and Miscanthus soil communities [36], as did we for nirK. They also found lower abundance of nifH compared to nosZ in all soils, matching our findings. When compared to the available bulk metagenomic assembly using the diginorm/partitioning assembly pipeline used in several recent published works [12, 27], only 125 to 146 rplB gene clusters were found in the available bulk metagenomic assembly (Table 3). In general, Xander assembled 10- to 100-fold more contig clusters for these three genes. The contigs, assembled by Xander, were longer in length and shared higher percent aa identity with known reference sequences.

In ordination analysis, these samples were grouped by crop, with corn being more distinct from the two perennial crops. This is consistent with 16S rRNA amplicon data and whole metagenomic analysis using the same samples that showed corn communities differed from those of perennial crops (Jiarong Guo, personal communication).

Using the closest matching rplB reference sequences, five phyla were present at 5 % or more of the total: Proteobacteria, Actinobacteria, Firmicutes, Acidobacteria, and Bacteroidetes. The distribution is similar to 16S rRNA amplicon sequence data [37] with the exception of the Acidobacteria, a typical major soil phylum. The reference set available for rplB is much smaller than those available for rRNA genes, and for some groups of common soil organisms there are few available reference sequences. This is especially true for the Acidobacteria and Verrucomicrobia where there are few cultured representatives. For the Acidobacteria, there are at least 26 subgroups roughly equivalent to classes [38], yet there are only 8 genome sequences, hence rplB sequences, available and from only 3 of these class-level subgroups. This affects the outcome of Xander at several steps: our HMM model may be biased against sequences with no relatives in the training set; if available references do not share any peptide of, in this case, 15 aa with an unknown, we will fail to identify an appropriate starting kmer; and assembled contigs may be misclassified if there are no related references available. Targeted genome sequencing efforts and single-cell sequencing have been filling in some gaps [39, 40], so this may become less of an issue.

About 10 % of bacteria were estimated to carry the nirK gene for the 21 samples (Table 4), similar to the qPCR estimates of 6 % nirK genes per rRNA gene (a multi-copy gene) in agricultural soil [41]. The nirK gene contigs, as expected, were dominated by matches in the Proteobacteria. However, when extrapolating from nearest gene match to taxonomic affiliation, it is important to keep in mind that some horizontal transfer of nirK has been detected [42]. One reference, B. japonicum USDA 110 (BAC52354), was the most abundant match in all three crop systems. These contigs averaged 91 % aa identity with a maximum of 98.4 % aa identity to the reference. Intriguingly, in bulk soil from Illinois sites growing the same three cropping systems, this same strain was found as having the best match to PCR amplicons targeting nosZ, encoding the terminal step in the denitrification pathway, although these matches were not as abundant [43]. This strain was widely used in soybean inoculants sold throughout the U.S. Midwest.

The HMP-defined community metagenomic data proved to be a good test set for validating Xander. While the community contains members from diverse taxonomic groups, several members are closely related and shared kmers of lengths 30 to 60 in their rplB genes, creating complexity in the graph structure. The overlap between the rplB gene sequences in the community combined with sequencing errors lead to small differences in the assembled sequences. The dataset was also challenging because of the low sequencing depth and regions of zero coverage of many members. Unfortunately, none of the members have nirK genes. Two members have nifH genes, but only one had all kmers of length 30 in the data.

We consider accuracy more important than longer and more contigs when assessing the assembler quality. For the defined community, we can measure the accuracy using the sequence identity of the recovered sequences to the known member sequences. The kmer length is an important factor in determining accuracy. A kmer length of 30 resulted in a high percentage of chimeric contig clusters. A kmer length of 45 assembled contig clusters with fewer mismatches and no chimeras but ended up with shorter contigs due to more regions with zero coverage. None of the rplB sequences from the defined community member are covered by kmers of length 60. For this dataset, length 45 appears to be a good choice to balance the trade-off between length and accuracy. We compared the performance of Xander to SAT-Assembler, another recently published gene-targeted assembly tool. Xander produced more and longer contigs with higher identity for more defined community members using a small dataset (Table 2). In experiments mixing the HMP community data with the corn 1 sample, SAT-Assembler did not return any non-chimeric contigs greater than 450 bp and with greater than 90 % identity to community member rplB genes. In contrast, Xander returned contigs corresponding to three community members.

Chimeras appeared to be a much smaller problem for contigs assembled from the rhizosphere soil data. Bulk assemblers attempt to avoid chimeras during assembly, often by terminating contigs where paths fork. In contrast, Xander relies on appropriate kmer size to reduce chimeras and post-processing to discard chimeric assemblies. On the HMP community data, UCHIME was very accurate and detected all chimeras found. We confirmed these results from UCHIME by manually validating all the contigs. It seems likely that because of the higher diversity, and potentially competitive exclusion, fewer sequences in the rhizosphere data share kmers of length 45.

Because Xander scores paths at the aa level, sequencing errors that translate into less likely amino acids for a given position will not be incorporated into a best path if an error-free read covering the same position exists; however, sequencing errors resulting in synonymous codons or resulting in more likely amino acids may be incorporated. Requiring a kmer count of 2 removed most sequencing errors and had the additional advantage of reducing the number of vertices searched by Xander. The latter is at least partially due to a reduction of starting kmers found in the graph (Table 1). For low-coverage datasets such as these, a minimum count of 2 decreased the length and number of contigs returned. For example, in sample C1, all five high-quality contigs (full-length and sharing 99 to 100 % aa identity with references) have multiple kmers of abundance 1 and so cannot be assembled with count 2. More than 35 % of the kmers in the Xander-assembled nirK and rplB contigs occurred only once in the reads (Fig. 3), while more than 40 % of nirK and rplB contigs identified in bulk assembly had median read coverage of 1 (data not shown). Xander worked well without discarding the reads with low abundance or partitioning the dataset. Such techniques are commonly employed by bulk assemblers but usually cause information loss [4].

When the search reaches a zero coverage region, meaning the gene is not complete, Xander will try to explore all available paths. It will eventually find the least bad path even though it is not related to the target gene. We used a heuristic algorithm to prune non-productive paths that reduces search space by at least tenfold with essentially identical results. This is a very simple heuristic that prunes paths that have not improved the overall path score within the user-specified number of steps. This works well for these target genes, but a more sophisticated heuristic might reduce the search space even further.

To save memory, Xander uses a Bloom filter for the in-memory representation of the de Bruijn graph [12] but is not dependent on this particular de Bruijn representation. In particular, a succinct de Bruijn graph representation has been used in some new assemblers such as MEGAHIT [44] to reduce space. Such a structure could be adopted in Xander to achieve additional space savings. Since only specific genes are targeted, pre-screening of reads could be used before the de Bruijn graph assembly to reduce space requirements even further. For example, SALT can filter out most non-target reads with minimal loss of target-gene reads by using read overlap information determined by scanning the reads through a set of relaxed HMM models [45]. Most bulk assemblers implement heuristics [46] to reduce the search space and remove errors and these could potentially be incorporated into Xander.

Project name: Xander

Project home page: https://github.com/rdpstaff/Xander_assembler

Operating system(s): Platform independent

Programming language: Java [47]

Other requirements: Java 1.6 or above

License: GNU GPLv3

Any restrictions to use by non-academics: None