The insect R package

Introduction

Welcome to the insect R package, a tool for assigning taxon IDs to amplicon libraries using informatic sequence classification trees. The insect learning algorithm takes a set of reference sequences (obtained from GenBank and/or other sources) to build a classification tree, which is then used to assign taxonomic IDs to a set of query sequences (e.g. those generated from an NGS platform such as Illumina MiSeq). The learning and classification functions are best suited to computing environments with multiple processors and access to large amounts of memory; however, most modest-sized datasets can be processed on standard personal computers if time is available. While not a prerequisite, insect is designed to be used in conjunction with the ape package (Paradis et al., 2004; Paradis, 2012), which features memory-efficient binary formats for DNA and amino acid sequences (“DNAbin” and “AAbin” object types), and the dada2 package (Callahan et al., 2016) which contains essential functions for de-noising high-throughput sequencing data and other important pre-processing steps.

The insect package can be used to analyze environmental DNA (eDNA) meta-barcode libraries as well as single-source NGS/Sanger amplicon sequences.

The most time-consuming and memory-intensive stage of the insect work-flow generally involves training the classifier. For example, the COI classifier used in this tutorial, which was built from the MIDORI UNIQUE reference trainingset consisting of nearly a million COI barcode sequences, took around three days to run on a 24 x multithread. The insect classification trees are amplicon specific, so a unique tree is generally required for each primer set. However, trees are already available for some of the more commonly used barcoding primers here:

The insect package also includes functions for downloading, trimming and filtering reference datasets (including a “virtual PCR” tool and an annotation quality filter), building a hierarchical taxonomy database, and training the classifier; however, these methods are beyond the scope of this introductory tutorial. New classification trees and updates are frequently added to the collection, so please feel free to suggest a barcoding primer set with which to train a classifier and we will endeavor to add it to the list.

The INSECT learning and classification algorithms

To learn a classification tree, a reference sequence dataset is first obtained from GenBank and/or other sources from which barcode sequences with accurate taxon IDs are available. These sequences are filtered to remove any with obvious taxonomic labeling issues, and trimmed to retain only the region of interest using the virtualPCR function (sequences that do not span the entire amplicon region are removed). the learn function then splits the training sequences into two subsets, maximizing the dissimilarity between the two groups. A profile hidden Markov model is then derived for each group (see Durbin et al. (1998) for a detailed description of these models). The partitioning and model training procedure then continues recursively, splitting the reference sequences into smaller and smaller subsets while adding new nodes and models to the classification tree.

Once the classifier has been trained, query sequences obtained from the specified primer set can be assigned taxonomic IDs along with probabilistic confidence values. The classification algorithm works as follows: starting from the root node of the classification tree, the likelihood of the query sequence (the full probability of the sequence given a particular model) is computed for each of the models at the two immediate child nodes using the forward algorithm (see Durbin et al. (1998)). The competing likelihood values are then compared by computing their Akaike weights (see Johnson and Omland, 2004). If one model is overwhelmingly more likely to have produced the sequence than the other, that child node is selected and the classification is updated to reflect the lowest common taxonomic rank of the training sequences belonging to the node.

This procedure is repeated recursively, descending down the tree until either an inconclusive result is returned from a model comparison test (i.e. the Akaike weight is lower than a pre-defined threshold, e.g. 0.9), or a terminal leaf node is reached, at which point a species-level ID is generally returned. The classify function outputs the taxon name, rank and ID number (i.e. NCBI taxon ID, WoRMS aphia ID, or other identifier depending on the taxonomy database used in the training step), along with the final Akaike weight value, which can be interpreted as a confidence score (between 0 and 1, with values close to 1 indicating high confidence). Note that the default behavior is for the Akaike weight to ‘decay’ as it moves down the tree, by computing the cumulative product of all preceding Akaike weight values. This is perhaps an overly conservative approach, but it minimizes the chance of mis-classifying or over-classifying the query sequences.

A worked example

This tutorial demonstrates the insect work-flow using an example dataset of COI sequences derived from autonomous reef monitoring structures (ARMS) in Ofu, American Samoa, amplified using the metazoan COI barcoding primers mlCOIintF and jgHCO2198 (GGWACWGGWTGAACWGTWTAYCCYCC and TAIACYTCIGGRTGICCRAARAAYCA, respectively; Leray et al. (2013)).

The dataset was first de-noised and tabulated using the DADA2 pipeline, to produce a table of chimera-free amplicon sequence variants (ASVs). The most abundant 16 ASVs are included in the insect package as an example dataset (see below).

If using other tools for de-noising, trimming, merging, etc, the data can be read in using either the readFASTA or readFASTQ functions to produce a “DNAbin” object compatible with the insect classifier.

devtools::install_github("shaunpwilkinson/insect")
library(insect)

download.file("https://www.dropbox.com/s/dvnrhnfmo727774/classifier.rds?dl=1", 
              destfile = "classifier.rds", mode = "wb")

## read in the example seqtab.nochim ASV table
data(samoa)
## get sequences from table column names
x <- char2dna(colnames(samoa))
## name the sequences sequentially
names(x) <- paste0("ASV", seq_along(x))
## optionally remove column names that can flood the console when printed
colnames(samoa) <- NULL 

classifier <- readRDS("classifier.rds")
classifier
names(attributes(classifier))

#> 'dendrogram' with 2 branches and 113833 members total, at height 70
#>  [1] "k"           "height"      "midpoint"    "members"     "class"      
#>  [6] "taxonomy"    "clade"       "frame"       "remainder"   "sequences"  
#> [11] "minscore"    "seqlengths"  "pointers"    "key"         "kmers"      
#> [16] "minlength"   "maxlength"   "model"       "trainingset" "alternative"
#> [21] "nunique"     "ntotal"      "taxID"       "numcode"

longDF <- classify(x, classifier, threshold = 0.8)

longDF <- cbind(longDF, t(samoa))

taxa <- aggregate(longDF[3:12], longDF["taxID"], head, 1)
counts <- aggregate(longDF[13:ncol(longDF)], longDF["taxID"], sum)
shortDF <- merge(taxa, counts, by = "taxID")

Acknowledgements

This software was developed with funding from a Rutherford Foundation Postdoctoral Research Fellowship from the Royal Society of New Zealand. Unpublished COI data care of Molly Timmers (NOAA).

References

Callahan,B.J. et al. (2016) DADA2: High-resolution sample inference from illumina amplicon data. Nature Methods, 13, 581–583.

Durbin,R. et al. (1998) Biological Sequence Analysis: Probabilistic Models of Proteins and Nucleic Acids. Cambridge University Press, Cambridge.

Johnson,J.B. and Omland,K.S. (2004) Model selection in ecology and evolution. Trends in Ecology and Evolution, 19, 101–108.

Leray,M. et al. (2013) A new versatile primer set targeting a short fragment of the mitochondrial COI region for metabarcoding metazoan diversity: application for characterizing coral reef fish gut contents. Frontiers in Zoology, 10, 34.

Paradis,E. (2012) Analysis of Phylogenetics and Evolution with R. Second Edition. Springer, New York.

Paradis,E. et al. (2004) APE: analyses of phylogenetics and evolution in R language. Bioinformatics, 20, 289–290.