BLAST-based validation of metagenomic sequence assignments

Related work

Use of a “first-pass” taxonomic classifier

We selected Kraken (Wood & Salzberg, 2014) (version 1.1) as the “first-pass” taxonomic classification program to be used in this study, primarily because of its widespread use in the bioinformatics community at large (the bioforensics community being no exception). Kraken was run in paired-end mode with default parameters and used standard Kraken databases for bacteria, archaea, viruses, plasmids, and human sequences.

Read simulation

For read simulation we used ART (Huang et al., 2011) (version 2016-06-05). To ensure thorough sampling, all experiments used simulated reads equivalent in total to 10× coverage of the source genome. For 150-bp reads, we used the built-in HS25 quality profile with an insert size of 200 ± 10 bp (mean ± standard deviation). For 250-bp reads, we used a custom Illumina HiSeq quality profile that we generated from recent runs of our HiSeq 2500, with an insert size of 868 ± 408 bp determined from recent library preparations. We supplied this information so that the simulated reads would have characteristics that closely matched what we would expect to obtain from a real HiSeq run in our laboratory, thus ensuring that the simulation results would be maximally useful to us. We recommend that others who emulate our procedures customize the attributes of their simulated reads to correspond to the real data they anticipate analyzing.

Similarity searches

We used BLASTN from the BLAST+ package (Camacho et al., 2009) (version 2.2.25+) together with the NCBI non-redundant nucleotide database (NCBI Resource Coordinators, 2016) (downloaded February 2017) for all classification experiments. Default parameters were used, except when excluding taxa from the reference database, in which case the -negative_gilist option was added. The BLAST computation was distributed over many cluster nodes to complete the analyses in a timely manner.

BLAST result filters

The output from BLAST includes a number of statistics that can potentially be used to filter the results, including alignment length, alignment percent identity, E-value (the number of similar scoring alignments one can “expect” to see by chance in a database of the size being searched), and bit score (a database size-independent measure of alignment quality).

We developed two basic ways of filtering BLAST results. The first we term an “absolute” filter, which simply removes BLAST hits that do not meet a particular criterion. Various possible criteria include minimum alignment length, minimum alignment percent identity, or maximum E-value. Of these three filters, this study only uses the E-value filter (abbreviated E), as E-value is fundamentally a composite of alignment length and alignment similarity. (Our software supports the use of all three filters, however, either individually or in combination.) If the best BLAST hit matches the target taxon after application of the absolute filters, it is then possible to apply a “relative” filter by computing the difference in E-value or bit score between the best hit and the best hit to a non-target taxon (should the latter exist). As very small E-values are typically rounded to zero, our software uses relative bit scores in this context for maximum applicability; we call this quantity the “bit score difference,” abbreviated b. If b is greater than or equal to a threshold determined via read simulation experiments, then we have validated the assignment of the read to the target taxon. Examples of the application of the bit score difference filter are given in Figs. 3 and 4.

Evaluation of classification performance

The two main metrics used in this study to evaluate classification performance are sensitivity and precision. The formula we use for sensitivity is slightly different from the standard one, though, and we also use a non-standard formula for precision in part of the analysis. Thus, here we explain the derivation of our sensitivity and precision formulas in detail.

BLAST result parsing, final taxonomic assignment, and calculation of statistics

BLAST result parsing and final taxonomic assignment of each read was performed with a custom Perl script capable of querying the NCBI taxonomy database (NCBI Resource Coordinators, 2016). If a target taxon is supplied as an argument to the script, assignments to the target taxon lineage that are more specific than the target taxon are simply reassigned to the target taxon. BLAST hits that do not meet the criteria specified by the absolute filters (minimum alignment length, minimum alignment percent identity, or maximum E-value) are removed, as are hits to the “other sequences” clade (NCBI taxon ID 28384), which are presumed to be erroneous. To make the final taxonomic assignment for each read, the lowest common ancestor (LCA) algorithm (Huson et al., 2007) is applied to the remaining hits that have a difference in bit score from the best hit less than a specified amount. If multiple parameter values are supplied for one or more filters, the script parses the BLAST results once for each possible combination of parameter values and writes the results to separate “LCA files,” thus enabling the user to efficiently perform parameter sweeps. The ultimate output from the script is one or more LCA files, each containing the final taxonomic assignment of each read for a particular combination of filter parameter values. Counts of TPs, FNs, IAs, and FPs (from which sensitivity and precision were calculated) were obtained using a separate Perl script that parses the LCA files produced by the BLAST result parser.

Determination of optimal BLAST filter parameter values

When deciding how the absolute and relative BLAST filters should be parameterized, an optimality criterion is needed. In the execution of bioforensic casework, it is important that any assignments made are correct. Thus, we first chose filter parameter values that maximized precision (i.e., minimized incorrect read assignments). In the event that multiple combinations of parameter values yielded exactly the same maximum precision value, we chose from among these the combination that maximized sensitivity (i.e., maximized detection of the target taxon). In the event that multiple combinations of parameter values yielded exactly the same maximum precision and sensitivity values, we reported the strictest combination.

In this study, we present examples aimed at detecting a variety of pathogenic target taxa including B. anthracis, Clostridium botulinum, pathogenic Escherichia coli, and Yersinia pestis. Due to inherent variation in the degree of interrelatedness among genomes from different taxonomic clades, optimal filter parameter values need to be set differently for each target taxon. To determine optimal filter settings, we must know the true origin of our test sequences; thus, we simulate reads from each target taxon genome, BLAST them against nt, and evaluate classification performance under different combinations of filter parameter settings. In each read simulation experiment, 81 different combinations of filter parameter values were tested—i.e., all combinations of maximum E-value (E) = {10⁰, 10⁻¹, 10⁻², 10⁻⁴, 10⁻⁸, 10⁻¹⁶, 10⁻³², 10⁻⁶⁴, 10⁻¹²⁸} and bit score difference (b) = {0, 1, 2, 4, 8, 16, 32, 64, 128}. The parameter optimization workflow is shown in Fig. 5.

Selection of near neighbor and alternate representative genomes

For each target taxon, we used the “Genome neighbor report” feature of the NCBI Genome database (NCBI Resource Coordinators, 2016) to select the most closely related complete genome of a different species or strain, as appropriate, to be used as the “near neighbor.” For species-level target taxa, we used the Genome neighbor report to select the complete genome of the same species that was most distantly related to the original representative genome, which we call the “alternate representative genome” (Table 1).

Clade-level exclusion

In the final read simulation experiment, clade-level exclusion (Brady & Salzberg, 2009) was performed to assess classification performance in the situation where the taxon for which one has sequence data is not represented in the reference database. In these tests, we simulated 250-bp reads from the taxon hypothetically missing from reference database, excluded this taxon from the reference database when performing BLAST searches, and then obtained optimal filter parameter values for classification of the target taxon, which in this case was the taxon immediately above the excluded taxon in taxonomic rank.

Evaluation of a “first-pass” taxonomic classifier

To demonstrate typical use of a first-pass taxonomic classification program, we analyzed all simulated reads from the B. cereus JEM-2 genome (Venkateswaran et al., 2017a, 2017b) with Kraken. The majority of the reads (79%) were assigned to the B. cereus group; of these, only 32% of the reads were assigned more specifically to B. cereus. Worryingly, however, a relatively small number of reads were assigned incorrectly to other B. cereus group species, including B. anthracis, B. cytotoxicus, B. mycoides, B. thuringiensis, and B. weihenstephanensis. Had this benign strain of B. cereus (JEM-2) been the sole representative of the B. cereus group in a metagenomic sample, an analyst using Kraken might have erroneously declared that a variety of B. cereus group species were present in the sample, including pathogenic B. anthracis. As false-positive assignments are relatively commonplace with first-pass classification programs, we were motivated to develop a procedure to validate the assignments of interest made by such classifiers.

Taxon selection

To demonstrate the BLAST-based validation procedure, we selected four target taxa, all of which are biological agents that could conceivably be of interest in a bioforensic setting. The first is B. anthracis, the bacterium that causes anthrax. The second is C. botulinum, a bacterium capable of producing the lethal botulinum neurotoxin. The third is a pathogenic strain of E. coli, E. coli O157:H7 str. Sakai, a bacterium that has been associated with major outbreaks of foodborne illness. The fourth and final target taxon is Y. pestis, the bacterium that causes bubonic plague. Thus, three out of the four target taxa represent particular species to be identified (B. anthracis, C. botulinum, and Y. pestis), whereas one target taxon represents a particular strain to be identified (E. coli O157:H7 str. Sakai). Species-level evaluations were performed using the representative strains indicated in Table 1. In two of the three evaluations, the genome chosen was a reference genome for the species (C. botulinum A str. ATCC 3502 and Y. pestis CO92). For the B. anthracis evaluation, the genome of the Ames Ancestor strain was used to ensure that the pXO plasmids were included, as presence of the pXO plasmids is normally required for B. anthracis to be fully virulent (Okinaka et al., 1999a, 1999b; Pannucci et al., 2002). To evaluate the implications of representative genome choice, an alternate representative genome was selected for each species. Additional information about the target taxa and near neighbors is provided in Table 1.

Simulated read experiments

A total of four simulated read experiments were performed to determine optimal BLAST filter parameter values for the identification of various target taxa. A comparison of sensitivity across experiments on a per-taxon basis is available in Figs. S1–S4, online.

Calculating precision: “near neighbor” vs. “false negatives”

Our results show that it was possible to achieve perfect precision in 3/4 simulated read experiments when using either the near neighbor or the FNs approach (the exception being the clade-level exclusion experiment). In 7/11 cases, the filter parameter values recommended by the two approaches were identical; in the cases where they differed, the FNs approach uniformly recommended more stringent filter parameter values than the near neighbor approach, resulting in reduced sensitivity (Tables 2–4). As mentioned previously, deciding which of the two approaches to calculating precision is superior would require a comparison to a provably optimal approach, which we currently deem computationally intractable. Each heuristic makes assumptions that may not always hold: the near neighbor approach assumes that a single genome that is closely related to the target taxon is sufficient to serve as a proxy for all other non-target taxa in the database, and the FNs approach assumes that the sequencing error process is symmetric. When seeking to avoid an erroneous claim that a particular biological agent is present in a sample, one may wish to use the more conservative set of parameter values recommended by the two approaches.

Practical application of the BLAST-based validation procedure

To demonstrate the practical application of the BLAST-based validation procedure, we downloaded a subset of metagenomic data collected from the New York City subway system (NCBI SRA ID SRR1748708), which the original study indicated might contain some reads from B. anthracis (Afshinnekoo et al., 2015). Indeed, analysis of this data with Kraken, our first-pass classifier, assigned 676 reads to B. anthracis (≈0.04% of reads). However, BLAST-validation of these 676 reads using the most conservative parameters recommended by our study (E = 10⁻⁶⁴ and b = 128; Table 2) resulted in zero reads assigned to B. anthracis. Even after significantly relaxing the minimum required bit score difference (setting b = 8), which was shown in Experiment 1 to significantly increase sensitivity (Table 2), still zero reads were assigned to B. anthracis. Thus, we would conclude that the 676 reads that Kraken assigned to B. anthracis were in fact false-positive assignments, which agrees with other follow-up studies that have been performed on the New York City subway data (Gonzalez et al., 2016).