Segal’s Law, 16S rRNA gene sequencing, and the perils of foodborne pathogen detection within the American Gut Project

American Gut Project data

American Gut Project data as of August 1, 2016 were analyzed. Given that our original objective was to determine the percentage of samples that contained either Salmonella or Listeria we queried the QIIME output files (i.e.,  ftp://ftp.microbio.me/AmericanGut//ag-July-29-2016/07-taxa/notrim/ag-cleaned_L6.txt) provided by AGP for those taxa. The AGP data contained 10,294 fecal, skin, or oral samples, and yielded 188 samples positive for Listeria and 1,495 positive for Salmonella; 31 samples were positive for both. This yielded 1,652 samples that were carried forward and analyzed here to determine the robustness of these classifications across multiple taxonomy assigners. The sequences from each of the 1,652 were extracted from the fasta file available from AGP (ftp://ftp.microbio.me/AmericanGut//ag-July-29-2016/02-filtered/sequences-notrim.fna), which resulted in a 48,312,131 sequences to be classified. No additional filtering of reads was done.

Taxonomic assigners

For each of the 1,652 samples, we classified the reads associated with them to taxonomy using four different classifiers—RDP (Wang et al., 2007), UCLUST (Edgar, 2010), BLAST (Altschul et al., 1990), and SortMeRNA (Kopylova et al., 2014) (Table 1). SortMeRNA was the method used to produce the AGP results. The program QIIME v1.9.1 (Caporaso et al., 2010) was used to run each of these classifiers. These assigners do not represent an exhaustive list of those available and were chosen as they are among the most frequently used; additional classifiers are likely to only increase the variability that we observe and, therefore, not change our findings. Additionally, we make no statements about the accuracy of the methods as we do not know the true taxonomy of the reads but rather are estimating the variability in classification across the methods, which illustrates the degree to which conclusions about pathogen presence is dependent on the classifier used.

Within each classifier a range of cutoffs were used for the parameter that measures confidence of taxonomic assignment (Table 1). The cutoff settings for these parameters are specific to each assigner and were determined through exploratory analyses that evaluated which cutoff values result in a substantial difference in number of reads classified. An example is UCLUST, which has a min_consensus_fraction parameter, which determines the fraction of hits that must agree in the number of assignments considered (max_accepts). With max_accepts set to 3 (the default), all that matters is whether min_consensus_fraction is more or less than 2/3. Thus, we evaluated two settings (1) min_consensus_fraction = 0.51 (the default, at least 2/3 agreement is required) and (2) min_consensus_fraction = 0.68 (complete agreement is required). The latter threshold is more conservative (higher specificity) but also runs the risk of a higher false-negative rate (lower sensitivity).

Reference databases

Analyses were performed against two different reference databases. The first database was the Greengenes (gg_13_8) database (DeSantis et al., 2006). This database was used to produce the AGP results and we primarily focus on the results using this database. Within the Greengenes database there are 99,322 sequences of which four are classified to Salmonella (one to only genus and three to species enterica) and five sequences classified to Listeria (one to each of the following species monocytogenes, seeligeri, grayi, fleischmannii, weihenstephanensis).

The second database was the human microbiome specific reference database HITdb described in Ritari et al. (2015). HITdb includes many fewer sequences than Greengenes (n = 2,472) of which a single sequence is classified to Salmonella enterica and zero are classified to the genus Listeria. The fewer number of representatives at this taxonomic level for these pathogens, and likely others, may be a desirable property as it could avoid false-positives. However, such a smaller database relative to Greengenes may run the risk of being too restrictive where a read could be reliably classified at one of those taxa given a sufficient algorithm. The taxonomy of within HITdb begins at the phylum level and, therefore, we do not present congruence among assigners at that kingdom level when using HITdb.

Congruence between classifiers

The results from each of the classifiers with the various parameter settings were combined in a pairwise fashion. To evaluate the congruence between the classifiers, we calculated the percent of total reads within each of five different comparison categories (1) match: both classifiers assigned the read to the same taxon, (2) mismatch: the classifiers differed in the taxonomy assigned to a read, (3 & 4) method 1 unclassified and method 2 unclassified: cases in which one classifier did not assign the read to a taxon but the other method did, and (5) noMatch: not classified by either classifier. The small fraction of the reads in the noMatch category (∼0.8% for the default parameter settings using Greengenes) were not analyzed. As these categories suggest, we include reads classified to the species level and reads partially classified (e.g., the lowest taxonomic rank assigned to a read was only phylum or family). We also explored more explicit methods for determining congruence between classifiers (i.e., inter-rater agreement; e.g., Cohen’s κ (Cohen, 1960)). Such statistics attempt to account for the fact that two classifiers could agree by chance, which percent of reads that match (like that used here) does not. However, preliminary analyses showed these statistics to be difficult to interpret as a result of the large number of observations and at times uneven distribution of classifications across a small number of categories (e.g., at the kingdom or phylum level). For example, at higher taxonomic ranks there is clearly better agreement between the classifiers (see Results) but prevalence and bias adjusted Cohen’s κ (PABAK) was substantially less than that observed at the species level (Supplemental Information 1).

As an example of how the incongruence among classifiers impacts the detection of a specific taxon, we focused on the number of reads under each method that were classified to the genera Salmonella and Listeria and quantified differences between classifiers. This was done for only the results produced with Greengenes. The congruence between methods in which reads were classified to Salmonella and Listeria was evaluated using the Jaccard Similarity index (0 meaning two methods classified completely different reads to a given taxon; 1 denotes complete overlap in the reads classified to a given taxon).

Default settings and Greengenes database

The initial investigation into taxonomy assignments provided by the AGP indicated 188 samples positive for Listeria, 1,495 samples positive for Salmonella, and 31 samples positive for both, yielding 1,652 samples. Taxonomic resolution was only to genus so the AGP results do not include whether pathogenic Listeria (e.g., Listeria monocytogenes) was present. To determine the robustness of these classifications we analyzed the sequences from those samples using default parameter settings and the Greengenes database (similar to the AGP analyses) and found little congruence among the four classifiers investigated in the assignment of reads to lower taxonomic ranks. For example, in all but the highest taxonomic rank (kingdom), there was an average difference of approximately 30% in the number of unique taxa detected (Fig. 1). The pattern of difference in number of taxa detected was also not consistent across the ranks where, for instance, the RDP classifier detected more taxa than SortMeRNA at the genus and species level but not at the phylum, class, order or family level (Fig. 1).

Focusing on congruence in the taxonomic identity assigned to reads, we found that congruence among the classifiers decreased with decreasing taxonomic rank (Fig. 2). This suggests that it is more difficult and there is greater error associated with classifying reads to lower taxonomic ranks, which has also been observed in other studies (Mizrahi-Man, Davenport & Gilad, 2013; Wang et al., 2007). This incongruence was primarily driven by a read being classified to a taxonomic group with one method but not the other rather than each method classifying the same read to different taxonomic groups (Fig. 2). Certain methods also seemed to result in greater incongruence. For example, RDP classified many more reads at the genus and species levels than the other methods (Fig. 2).

When we focused on the congruence across the classifiers in assigning reads to Salmonella or Listeria, we found that 447 and 1,521 samples contained reads classified to Listeria or Salmonella, respectively. The number of Listeria ‘positive’ samples reported by AGP was 188 and for Salmonella it was 1,495 (Table 2). However, the different methods varied substantially in the number of reads classified to each genus with BLAST consistently being the most liberal (i.e., classified the most reads) followed by UCLUST; RDP was the most conservative (i.e., did not classify reads) (Table 2). Perhaps of more importance is that there was little congruence between the methods in which reads were classified to each taxonomic group as indicated by the low Jaccard similarity indices between the pairwise comparisons of the methods (Fig. 3). The metric suggests there is little consistency among the methods when classifying Listeria or Salmonella (genus or species) to the reads. For example, only 22% of the reads classified by UCLUST or BLAST to a species within Salmonella were actually classified to the same species—some of those reads may not have been classified at all by one method while others could have been classified to a different species (Fig. 3).

Some of these differences in classifications could likely be reconciled by changing the confidence parameter. For example, with the default settings, read 10317.000001040_2356 was clustered by UCLUST and BLAST to the genus Salmonella, by RDP to the species Salmonella enterica, and by SortMeRNA to the family Enterobacteriaceae. If one reduced the confidence setting for UCLUST, BLAST and SortMeRNA then the methods might converge on the classification of Salmonella enterica. Alternatively, the confidence settings for UCLUST, BLAST, and RDP could be increased so the methods would converge on the same classification at the family level. However, looking at read 10317.000001127_16840, it was classified to the genus Enterococcus within the order Lactobacillales by SortMeRNA and UCLUST, only to the class level (Bacilli) with RDP (even with the low default confidence threshold of 0.5), and Listeria weihenstephanensis within the order Bacillales by BLAST. This read, which is representative of many of the reads classified to Listeria by at least one method, cannot be congruently classified below the rank of class. Although these examples illustrate the ways in which differences in classification may be reconciled, this would be problematic for an entire set of reads (e.g., each read may require its own threshold settings and reads may be ‘under’ classified if one takes a lowest common ancestor approach).

Impact of varying confidence threshold and reference database

The level of congruence among the classifiers varied substantially depending on the parameter settings that each used to determine whether a read was classified to a taxon. However, no parameter combinations produced complete congruence between two classifiers and, in general, the default settings resulted in the highest congruence among the classifiers (fewest number of reads classified by one method but not the other and highest number of reads being classified to the same taxon (i.e., match; Figs. 2 and 4; Supplemental Information 1). An example of classification at the family rank with the RDP classifier and the Greengenes database illustrates this. As the cutoff for the confidence parameter for RDP is increased from the default, a greater number of reads were not classified; thus, fewer reads were classified to the same taxon (Figs. 4A & 4B). This suggests that a higher threshold setting within RDP has the unwanted consequence of producing many false negative assignments (a read is not classified to a taxon when it should be).

The pattern that different threshold parameter cutoffs affected the degree of congruence among the assigners and that no combination of threshold settings resulted in complete congruence was also observed when classifying reads based on the HITdb database (Fig. 4; Supplemental Information 1).