THAPBI PICT—a fast, cautious, and accurate metabarcoding analysis pipeline

Introduction

Metabarcoding of DNA is a sensitive and powerful method to detect, identify, and potentially quantify the diversity of biological taxa present in any given environmental sample. It is based on PCR amplification of a “barcode” region diagnostic for the groups of organisms of interest followed by high-throughput sequencing of the amplimers, and is often applied to environmental DNA (eDNA) samples (Deiner et al., 2017). This method is revolutionising areas of research including wildlife conservation, ecological processes and microbiology, by highly-sensitive detection of biodiversity across many taxa simultaneously (Arulandhu et al., 2017). Metabarcoding enables early detection of invasive threats to plant and human health in support of biosecurity (Batovska et al., 2021; Green et al., 2021), and is applicable to many complex and intractable systems, such as soil (Ahmed et al., 2019), in which standard methods of microbial isolation and characterisation are impractical or costly.

Our motivating use case is metabarcoding in which multiple environmental samples are multiplexed for high-throughput sequencing on the Illumina platform using paired-end reads, and for which the expected PCR amplification product is short enough to be fully covered by the overlapping paired reads. Each sample is expected to yield taxon-specific marker sequences that can be matched to a high-quality database of marker sequences with known taxonomic identity, to give a taxonomic breakdown reflecting the community composition. One of our goals was to minimise false positive identification of the presence of any taxon on the basis of small or disputable quantities of physical evidence. Metabarcoding is prone to generation of artefactual sequence variation and sufficiently highly sensitive to register low-abundance sample reads at the same level as such sequences, and sequences originating from cross-sample contamination and “splashover” in even a careful laboratory. We therefore chose to prioritise accurate reporting of taxonomic assignment for high abundance sequences over sensitive detection of low-abundance marker sequences.

This manuscript was initially released as a preprint (Cock et al., 2023). We describe THAPBI PICT v1.0.0, a metabarcoding tool developed as part of the UKRI-funded Tree Health and Plant Biosecurity Initiative (THAPBI) Phyto-Threats project, which focused on identifying Phytophthora species in commercial forestry and horticultural plant nurseries (Green et al., 2021). Phytophthora (from Greek meaning plant-destroyer) is an economically important genus of oomycete plant pathogens that causes severe losses and damage to plants in agricultural, forest and natural ecosystems. The Phyto-Threats project’s metabarcoding protocol used nested PCR primers designed to target the Internal Transcribed Spacer 1 marker sequence (ITS1; a genomic region located between 18S and 5.8S rRNA genes in eukaryotes) of Phytophthora and related plant pathogenic oomycetes (Scibetta et al., 2012). This approach is the current de facto standard within the oomycete community (Robideau et al., 2011), and these primers have been used in conjunction with THAPBI PICT in recent Phytophthora surveys (Vélez et al., 2020; La Spada et al., 2022). PICT was short for Phytophthora ITS1 Classification Tool.

We describe the implementation, operation, performance and output of THAPBI PICT using datasets from the Phyto-Threats project, and public metabarcoding datasets. Although originally designed as a Phytophthora ITS1 Classification Tool (PICT), we show that with appropriate primer settings and a custom database of genus/species distinguishing markers, THAPBI PICT is an effective tool for analysis of short read amplicon sequencing data with barcode marker sequences from other organisms.

Workflow overview

The THAPBI PICT core workflow comprises (i) data reduction to unique marker sequences, often called amplicon sequence variants (ASVs) (ii) discard of low abundance sequences to remove noise and artifacts (rather than attempting error correction by default), and (iii) classification using a curated reference database. This approach contrasts with commonly-used operational taxonomic unit (OTU) clustering approaches (as implemented, for example, in QIIME (Caporaso et al., 2010), UPARSE (Edgar, 2013), and MOTHUR (Schloss et al., 2009)), which can be sensitive to changes in the input data resulting in unpredictable clustering behaviour (Callahan, McMurdie & Holmes, 2017) and overestimate population diversity (Nearing et al., 2018).

THAPBI PICT’s approach of reducing amplicons to ASVs is similar to that of DADA2 (Callahan et al., 2016) but, by contrast, THAPBI PICT does not by default attempt to correct sequencing errors with a denoising model. Our approach is instead to discard low-abundance sequences because we consider that they are likely not to represent meaningful biological information in the sequenced sample. We observe, using synthetic control sequences, that the abundance of such controls accidentally transferred between samples tends to exceed by no small margin the abundance of amplicons whose sequence variation might constitute “noise” in the amplicon sequence data. We consider the observed abundance of (e.g., synthetic) control sequences, which could not have been present in the biological sample, to be a lower bound for the abundance of reads we can confidently claim derive from that sample. Consequently, ASVs with much lower total abundance cannot confidently be determined to derive from the analysed sample, and so are discarded. In general, we consider that proper use of negative and synthetic controls, to account for alternative sources of experimental error, such as accidental transfer or “splashover” from one well to another, should be considered best practice in metabarcoding.

Figure 1 gives an overview of the workflow. Paired raw Illumina FASTQ files for each sample are merged by overlap, trimmed to remove primers, and reduced to a list of observed unique marker sequences (labelled by MD5 checksum) with abundance counts. Discarding low abundance sequences further reduces the data volume; unique reads alone may represent half the data (and 90% of the ASVs), but may not derive from the sequenced sample. The remaining higher abundance sequences are then classified by matching them to a curated database. By default a species-level assignment is made when a database entry is identical or different by at most one base pair (1 bp; algorithm onebp) to the query. The matching algorithm can be chosen to adjust sensitivity for taxonomic classification (Table 1).

Following read preparation and ASV classification, the pipeline generates two tables describing (i) taxon presence/absence for each sample, and (ii) ASV presence/absence for each sample (Figs. 2A and 2B respectively), in both plain text and Excel format. If the user provides suitably formatted sample metadata, cross-referenced by the filename stem, this can be incorporated into the report to make for easier interpretation. Additionally, an edit-graph showing the distances between the ASVs recorded in the sample can be exported (e.g., Fig. 3).

Classifiers and databases

Classification assessment

In assessing classification performance, it is the combination of both classification method (algorithm) and marker database which matters. Settings like the abundance threshold are also important, and the tool default settings partly reflect one of the original project goals being to avoid false positives.

To objectively assess a metabarcoding classifier we require sequenced samples of known composition, which generally means single isolates (where a single marker sequence is typically expected), or mock communities of known species (the bulk of our examples). Carefully controlled environmental samples may also be used, such as Muri et al. (2020) in our worked examples. Here a lake was drained to collect and identify all the individual fish, but this is problematic as the lake was large enough that DNA from each fish could not be expected at all the sampling points, giving an inflated false negative count.

Our tool includes a presence/absence based assessment framework based on supplying expected species lists for control samples, from which the standard true positive (TP), false positive (FP), true negative (TN), and false negative (FN) counts can be computed for each species. These are the basis of standard metrics like sensitivity (recall), specificity, precision, F-score (F-measure, or F1), and Hamming Loss. It is simple but not overly helpful to apply metrics like this to each species, rather the overall performance is more informative.

However, some scores like the Hamming Loss are fragile with regards to the TN count when comparing databases. The Hamming Loss is given by the total number of mis-predicted class entries divided by the number of class-level predictions, thus $(FP+FN)/(TP+FP+FN+TN)$. Consider a mock community of 10 species, where the classifier made 11 predictions which break down as $9$ TP and $2$ FP, meaning $10-9=1$ FN. Suppose the database had a hundred species (including all ten in the mock community), that leaves $100-9-1-2=88$ TN, and a Hamming Loss of $3/100=0.03$. Now suppose the database was extended with additional references not present in this mock community, perhaps expanding from European Phytophthora species to include distinct entries for tropical species, or a sister group like Peronospora. The denominator would increase, reducing the Hamming Loss, but intuitively the classifier performance on this mock community has not changed. To address this, the classifier assessment also includes a modified ad hoc loss metric calculated as the total number of mis-predicted class entries divided by the number of class-level predictions ignoring TN, or $(FP+FN)/(TP+FP+FN)$ which in this example would give $3/12=0.25$ regardless of the number of species in the database. This is an intuitive measure weighting FP and FN equally (smaller is better, zero is perfect), a potential complement to the F-score.

Note that the assessment framework only considers species level predictions, ignoring genus only predictions and unknowns, and thus will not distinguish between the default onebp classifier and variants like 1s3g (see Table 1).

As a benchmark of the default classifier and Phytophthora focused database, we used the 10 and 15 species mixes in Riddell et al. (2019), see Table 2. This was originally analysed with the metapy pipeline with a high stringency classifier using bowtie to find perfect alignments, and a more relaxed classifier using swarm for clustering. In both samples and both classifiers, Phytophthora boehmeriae was not found, and this was attributed to uncompetitive amplification in a mixed DNA sample due to poor PCR primer binding. That being so, the best classifier results would be either 14 TP and 9 TP respectively, with 0 FP if the markers were unique, and 1 FN. However, note that not all the markers are unique, both mixes contain species known to share their ITS1 marker with other species, giving unavoidable technical FP, also discussed in Riddell et al. (2019).

Using F1 score or our ad hoc loss ranking, THAPBI PICT is clearly performing best on the 10 species mix (and better than metapy did). However, with default settings it does worse on the 15 species mix, due the high FN count where the default ASV abundance threshold of 100 reads is excluding expected species. In this MiSeq dataset the Illumina Nano Kit was used giving lower yields, making the default 100 read threshold overly harsh. Optimising on maximising the F1 score and minimising ad-hoc-loss, and weighting the two communities equally, suggests running THAPBI PICT with an ASV read abundance threshold of around 50 reads performs best overall here. This is a fundamental problem however, low abundance community members can be indistinguishable from background noise/contamination, meaning without controls the best threshold is arbitrary.

Reporting

The pipeline produces two tabular reports (which can also be requested directly with the summary command), output as both tab-separated plain text, and Excel format with colouring and conditional formatting (Fig. 2). These include information on read counts from the preparation stage (as used in Fig. 4), information on the abundance thresholds, and foremost the species classification from the chosen method. The user may provide a table of metadata cross referenced by the sample FASTQ filename stem, which will be used for sorting the samples and if possible colouring inferred sample groupings (e.g., sample source, or replicates) to ease interpretation. This allows quick visual comparison of replicates as adjacent rows/columns.

The read report by default sorts the ASVs by their taxonomic classification, and then by abundance. This makes it easy to identify the most common unknowns or genus-only predictions for manual review (using the ASV sequence). This sorting also means that when the thresholds are low enough to let through noise, the grey halo effect shown in the edit graph (see Fig. 3) is also visually distinct as highly abundant rows followed by less abundance variants. This read report can also be exported in BIOM format.

For many of the worked examples the sample metadata on the NCBI Short Read Archive (SRA) or European Nucleotide Archive (ENA) had to be supplemented by information in the associated publication. Providing such metadata to the archives using an approved ontology based checklist is non-trivial, but adds greatly to the reuse potential (Tedersoo et al., 2015). We provide an ena-submit command which facilitates using the interactive ENA upload step for matching FASTQ filenames to previously entered sample information.

The tool’s repository includes a number of helper scripts, including a pooling script written for the Phyto-Threats project for preparing plant nursery specific summary reports. This simplifies the sample report by combining replicate samples into a single row, and can either use the read count sum, or just “Y” (present) or “N” (absent).

The other noteworthy report from the tool is an edit graph, invoked via the edit-graph command, as shown in Fig. 3. By default this outputs the edit graph in XGMML format which can then be visualised in a tool like Cytoscape (Shannon et al., 2003), with a choice of node layouts and representations (e.g., customising node size by sample count, or colour by genus). The graph can help guide the choice of minimum abundance threshold (as discussed above), and the choice of classifier. In the example shown with a 3 bp maximum edit-distance shown, the cliques formed are for the most part clearly distinct species, with a single central node. With the default ITS1 marker used for Phytophthora we find greater sequence variation and therefore more diverse non-simple clusters for species like Phytophthora nicotianae and P. gonapodyides, but most species show a single central ITS1 sequence.

Development practices

THAPBI PICT is released as open source software under the MIT licence. It is written in Python, a free open source language available on all major operating systems. Version control using git hosted publicly on GitHub at https://github.com/peterjc/thapbi-pict/ is used for the source code, documentation, and database builds including tracking the hand curated reference set of Phytophthora related ITS1 sequences. Continuous integration of the test suite is currently run on both CircleCI (Linux) and AppVeyor (Windows). Software releases are to the Python Packaging Index (PyPI) as standard for the Python ecosystem, and additionally packaged for Conda via the BioConda channel (Grüning et al., 2018). This offers simple installation of the tool itself and all the command line dependencies on Linux or macOS. Installation on Windows requires manual installation of some dependencies. The documentation is currently hosted on Read The Docs, updated automatically from the GitHub repository.

Worked examples

In this section we briefly discuss the application of THAPBI PICT to public data sets from several published articles, covering a range of organisms and markers. The selection has prioritised examples including mock communities and negative controls, and have been included in the tool documentation as worked examples. These worked examples generally are highly concordant with the published analyses, with differences largely down to the exact choice of thresholds.

The example scripts first-run times range from a few minutes with under 1 GB of raw FASTQ data (Bakker, 2018; Riddell et al., 2019; Walker et al., 2019; Muri et al., 2020), to a few hours with the larger datasets like Ahmed et al. (2019) with 12 GB of input. These times are dominated by the merging the paired reads during read preparation stage, and as discussed earlier, running the read preparation stage in parallel across a cluster is advised on larger projects.

The first worked example is a simple one using the provided Phytophthora ITS1 database we have generated for this work to reexamine Riddell et al. (2019). This example does not include the synthetic controls introduced later, but does have blanks as negative controls and simple mock communities as DNA mixes (discussed above for classifier assessment). The second example uses Redekar, Eberhart & Parke (2019) but focuses on how to build a database, including how species names can optionally be validated against the NCBI taxonomy.

The example based on Muri et al. (2020) is a single 12S marker for fish, with a custom database including numerous off-target matches like humans and sheep. In this case the lake contents were determined by draining the lake and collecting the fish, but this did not determine which of the sampling sites any given fish might have visited, complicating interpretation compared to an artificial mock community. Another single marker example based on Walker et al. (2019) uses COI in simple mock communities of bats, and shows the importance of the database content with the default classifier. The most interesting single marker example is based on Palmer et al. (2018), fungal ITS2 markers with mock biological fungal communities and synthetic control sequences. This has been discussed above in the context of setting abundance thresholds.

There are examples with multiple markers which were sequenced separately in Klymus, Marshall & Stepien (2017), two different 16S mitochondrial markers with mock communities, and Ahmed et al. (2019), four different markers in mock nematode communities. The example in Batovska et al. (2021) uses three markers together, while Arulandhu et al. (2017) sequences over a dozen markers together. Here the primer sequences themselves are non-overlapping and so serve to separate out the amplicons for each sample, allowing them to be matched to the relevant reference set. Note currently a primer cocktail as used for the COI example in this data set is not supported. This article is also noteworthy as an inter-laboratory replication study of metabarcoding.

Datasets from some potentially useful publications could not be used directly, generally for technical reasons. Many used custom multiplexing (Elbrecht & Leese, 2015; Port et al., 2016; Elbrecht et al., 2016, 2019; Elbrecht, Peinert & Leese, 2017; Braukmann et al., 2019), and thus would require separate de-multiplexing before use. Some like Braukmann et al. (2019) and Duke & Burton (2020), use an amplicon too long to span with overlapping Illumina MiSeq paired reads. Sometimes, however, articles did not provide the raw FASTQ files. For instance, Blanckenhorn et al. (2016) did not share the raw FASTQ files at all, while Hänfling et al. (2016) and Zaiko et al. (2022) provided primer trimmed FASTQ files. Some older articles (also) used the Roche 454 or Ion Torrent platforms, which would require re-engineering mainly around the different error profile, which is potentially unsuited for our default strict classifier.

Discussion

Here we present a novel and flexible pipeline for the objective analysis of metabarcode data, not just of single markers but also pooled markers where the amplicons can be separated via their primer sequences. Some of the design choices and default settings reflect the initial use case being Phytophthora ITS1 sequence markers in a context where specificity was favoured over sensitivity. In general, appropriate abundance thresholds and classifier algorithm will be experiment and/or marker specific, with the quality of the reference database a key factor. All amplicon barcoding experiments should be designed with suitable controls to assess the limits of quantification vs presence/absence (Lamb et al., 2019), including the effects of the PCR (Thielecke et al., 2017) and contamination (Thalinger et al., 2021).

By design, the tool currently reports lists of genus/species for each ASV, without attempting anything like a most recent common ancestor analysis. This limitation can be a handicap with some use-cases where the markers may not readily resolve at species level, and/or an ASV is often shared between genera. See for example, the Brassicaceae discussed in Arulandhu et al. (2017), and fish examples in Muri et al. (2020). Moreover, it makes the tool unsuited to markers like regions of the bacterial 16S rRNA gene which are typically used at phylum level with environmental datasets (Straub et al., 2020). Rather it is appropriate for comprehensive analyses of better defined taxonomic markers such as the plant pathogenic oomycete ITS1 marker used primarily for Phytophthora, where it is proving valuable for the ongoing characterisation of a comprehensive set of several thousand samples from plant nurseries in the Phyto-Threats project (Green et al., 2021) and in natural ecosystems (Riddell et al., 2019).

Our pipeline supports using negative or synthetic spike-in controls to set an abundance threshold on groups of samples (such as each sequencing run). Rather than ASV-specific thresholds as in Palmer et al. (2018), THAPBI PICT takes the more cautious approach of interpreting the worst case as indicative of cross contamination rates for that sample batch in general.

The pipeline does not currently explicitly attempt to find and remove chimera sequences beyond the use of abundance filters. As discussed in Edgar (2016), chimeras which are also only 1 bp away from a reference sequence cannot be distinguished from a point error, and would be matched to that reference by all but our strictest identity classifier. Apart from this special case, any high abundance chimera would likely appear in our reports as an unknown, and would most likely be only in a single sample. Regular manual inspection of the high abundance unknown reads appearing in multiple samples was part of the ongoing quality assurance during the Phyto-Threats project, primarily to identify any gaps in the database coverage. The only clear chimeras identified were from our synthetic controls, where part of our non-biological sequence was fused to some unexpected sequence. Potentially more complex mock communities of synthetic sequences could be used to generate a gold standard for identifying chimeras which might serve as a benchmark dataset for testing chimera algorithms.

Another important difference from other ASV based tools like DADA2 (Callahan et al., 2016), obitools (Boyer et al., 2016) and UNOISE2 (Edgar, 2016), is THAPBI PICT does not by default attempt read correction. From the halo pattern of PCR induced variants seen from synthetic inputs as viewed on an edit-graph, like Fig. 4B, there is usually a natural central node to which a variant can be attributed. However, the situation is not always clear cut, with some species like Phytophthora gonapodyides showing a range of known ITS1 sequences. Rather our approach is to exclude most PCR noise through the abundance filters, and allow a modest amount of variation when matching the higher abundance sequences to the reference database. As an option, however, the pipeline can apply our re-implementation of the original UNOISE2 Levenstein distance based read-correction described in Edgar (2016), or invoke the UNOISE3 algorithm implemented in the USEARCH or VSEARCH tools. Read-correction seems most appropriate where the reference sequences are well separated, unlike our default Phytophthora ITS1 amplicon where a single base pair can distinguish known species, and thus read correction can mask lower abundance species under their more abundant neighbours.

Examination of mock community samples of our synthetic spike-in sequences showed ASV abundance to be at best semi-quantitative, as found in other work (Palmer et al., 2018; Lamb et al., 2019). For the Phyto-Threats project reports sent to plant nursery owners, we therefore only use species presence/absence (above or below the abundance threshold, and pooled replicates). However, the raw abundances are in the main tool reports, and can be used for plots or a quantitative interpretation where appropriate. The nested primer protocol with two rounds of PCR may be a factor in undermining quantitative interpretation, and increasing the risk of cross-sample or other sample contamination.

Conclusion

Here we present a novel and flexible pipeline for the objective analysis of metabarcode data, with user friendly reports including ASV read counts enabling custom graphs, as well as summary species lists per sample. While initially designed for Phytophthora ITS1 sequence markers, the THAPBI PICT tool can be applied more generally, including to samples containing multiple marker regions. It is best suited to markers where ASV are at least genus specific. Care should be taken picking appropriate abundance thresholds, which can be set using negative and/or synthetic controls, and in applying read-correction for de-noising. While high-throughput amplicon sequencing does give read counts per species (or per ASV), we and others caution against treating these as quantitative (Palmer et al., 2018; Lamb et al., 2019). The most suitable classifier algorithm will be marker specific, with the quality and coverage of the reference database a key factor. Including mock communities in your experiment allows the performance of classifier and database to be evaluated objectively.