HAYSTAC: A Bayesian framework for robust and rapid species identification in high-throughput sequencing data

Abstract

Identification of specific species in metagenomic samples is critical for several key applications, yet many tools available require large computational power and are often prone to false positive identifications. Here we describe High-AccuracY and Scalable Taxonomic Assignment of MetagenomiC data (HAYSTAC), which can estimate the probability that a specific taxon is present in a metagenome. HAYSTAC provides a user-friendly tool to construct databases, based on publicly available genomes, that are used for competitive read mapping. It then uses a novel Bayesian framework to infer the abundance and statistical support for each species identification and provide per-read species classification. Unlike other methods, HAYSTAC is specifically designed to efficiently handle both ancient and modern DNA data, as well as incomplete reference databases, making it possible to run highly accurate hypothesis-driven analyses (i.e., assessing the presence of a specific species) on variably sized reference databases while dramatically improving processing speeds. We tested the performance and accuracy of HAYSTAC using simulated Illumina libraries, both with and without ancient DNA damage, and compared the results to other currently available methods (i.e., Kraken2/Bracken, KrakenUniq, MALT/HOPS, and Sigma). HAYSTAC identified fewer false positives than both Kraken2/Bracken, KrakenUniq and MALT in all simulations, and fewer than Sigma in simulations of ancient data. It uses less memory than Kraken2/Bracken, KrakenUniq as well as MALT both during database construction and sample analysis. Lastly, we used HAYSTAC to search for specific pathogens in two published ancient metagenomic datasets, demonstrating how it can be applied to empirical datasets. HAYSTAC is available from https://github.com/antonisdim/HAYSTAC.

Fig 1. HAYSTAC’s workflow.

HAYSTAC consists of three main modules: (i) DATABASE, which builds a database of reference genomes from various user input sources; (ii) SAMPLE, which handles downloading sequencing files from the SRA and pre-processing of samples prior to analysis; and (iii) ANALYSE, which performs an analysis of a sample against a database by applying the mathematical model (see methods) for taxonomic abundance estimation.

Fig 2. Computational performance of HAYSTAC and other methods.

Linear regression of the elapsed time (wall clock) and peak memory usage (maximum resident set size) for a sample of size 1 million reads and reference databases containing either 10, 100 or 500 genomes, each with 5 replicates. When constructing the database, HAYSTAC uses substantially less memory and runs faster than either Kraken2/Bracken, KrakenUniq or MALT for restricted database sizes. When performing analyses, HAYSTAC uses less memory than Kraken2/Bracken, KrakenUniq or MALT, while its runtime was only marginally slower.

Fig 3. Accuracy of HAYSTAC and other methods for a simple simulation.

Bar plot showing the mean count of false positives (red), false negatives (orange), and true detected species (blue) for two versions of the simple simulations dataset, each with two replicates: (A) without human DNA contamination (n = 2); and (B) with human DNA contamination (n = 2). The dotted line shows the number of simulated species in each set of samples (i.e., the maximum true positive; n = 10), and numbers above the error bars indicate the mean species count in each category. For the simulation without human contamination, HAYSTAC outperforms Kraken2/Bracken, KrakenUniq and MALT, and performs equally with Sigma (i.e., no false positives or false negatives). For the simulation with human contamination, HAYSTAC outperforms all four other methods. All of the 5 false positive species identified by HAYSTAC are known to contain human sequences in their reference genomes [37], confounding any analyses which do not explicitly filter human contamination.

Fig 4. Accuracy of HAYSTAC and other methods for an oral microbiome simulation.

Bar plot showing the mean count of false positives (red), false negatives (orange), and true detected species (blue) for two versions of the oral microbiome dataset, each with six replicates: (A) modern simulation, with fixed read lengths (n = 6); and (B) ancient simulation, with variable read lengths and post-mortem damage (n = 6). The dotted line shows the average number of simulated species in each set of samples (i.e., the maximum true positive; n = 178), and numbers above the error bars indicate the mean species count in each category. For the modern simulation, HAYSTAC substantially outperforms Kraken2/Bracken, KrakenUniq and MALT with respect to false positives, and performs equivalently with Sigma. For the ancient simulation, HAYSTAC outperforms all four other methods with respect to false positives. The overall high rates of false negative identifications are due to the absence of many simulated species from the reference database for all four methods. HAYSTAC also outperforms all the other four methods in both the modern and ancient Oral Microbiome datasets by identifying the highest number of true positive species.

Fig 5. Accuracy of HAYSTAC and other methods using a reference database restricted to a single genus.

Bar plot, with a pseudo-log10 transformed y-axis, showing the mean count of false positives (red), false negatives (orange), and true detected species (blue) for nine different genera (Bacteroides, Burkholderia, Campylobacter, Clostridium, Corynebacterium, Desulfitobacterium, Mycobacterium, Solimonas and Streptococcus), each with 20 samples from the general and oral microbiome datasets. The dotted line shows the average number of simulated species in each set of samples (i.e., the maximum true positive; n = 2.0), and numbers above the error bars indicate the mean species count in each category. For the genus specific analysis, HAYSTAC substantially outperforms both Kraken2/Bracken and MALT with respect to false positives and performs better than Sigma.

Fig 6. HAYSTAC inferred posterior abundance levels.

Scatter plot showing the mean posterior abundances across all taxa (n = 362) and samples (n = 20) for either a genus specific database or the entire RefSeq representative database of prokaryotic species. Using a genus specific database has a small positive bias in mean posterior abundance for taxa within that genus (paired t-test p-value < 2.2·10⁻¹⁶, mean of the differences = 3.9·10⁻⁶), nevertheless the overall abundance levels are highly correlated (R² = 0.999). Computational runtime for the genus specific analyses are faster and use less memory, making genus specific analyses suitable for rapid initial screening (e.g. 1 million reads against a Corynebacterium specific database runs approximately 6.15 faster than against a database containing 500 species and uses approximately 3.9 times less memory).

Fig 7. Histogram of the number of assigned pathobiont simulated reads.

Histogram showing the read count frequency of true positive (blue) and false positive (red) pathobiont reads as identified by HAYSTAC and HOPS, after screening 200 spiked iterations of an ancient Oral Microbiome dataset sample (anc200e2repgn). HAYSTAC identifies robustly more pathobiont reads than HOPS, while producing less false positive identifications.

Fig 8. Posterior abundances of Yersinia species in Case Study 1.

Heatmap showing the mean posterior abundances for the seven RISE samples, based on a genus specific analysis of 18 Yersinia species. Yersinia pestis is the species with the highest posterior abundance, followed by Y. pseudotuberculosis, in agreement with the results of (Rasmussen et al., 2015).

Fig 9. Posterior abundances of oral microbiome species in Case Study 2.

Heatmap showing the mean posterior abundances for the 44 dental calculus samples, based on a custom database that combined the prokaryotic representative RefSeq and pathobionts with complete genomes from the following 5 genera: Corynebacterium, Haemophilus, Klebsiella, Streptococcus, and Bordetella. Species from these five genera of interest that naturally colonise the oral cavity can be found in more samples and at higher abundance (e.g. C. matruchotii) compared to pathobionts of the upper respiratory system (e.g. S. pneumoniae).