README.Rmd

---
title: ""
output: 
  md_document
---

```{r, echo = FALSE, message = FALSE, warning = FALSE}
#knitr::opts_chunk$set(out.width='750px', dpi=200)
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```




--------------------------------------------------------------------------------

### Informatic sequence classification trees

`insect` is an R package for taxonomic identification of amplicon 
sequence variants generated by DNA meta-barcoding analysis. 
The learning and classification algorithms implemented in the 
package are based on full probabilistic models (profile hidden Markov models) 
and offer highly accurate taxon IDs, albeit at a relatively high computational cost.

The package also contains functions for searching and downloading reference 
sequences and taxonomic information from NCBI, 
a "virtual PCR" tool for sequence trimming, 
a function for purging erroneously labeled reference sequences, 
and several other tools.  

`insect` is designed to be used in conjunction with the 
[dada2](https://benjjneb.github.io/dada2/index.html) pipeline or other
de-noising tools that produce a list of amplicon sequence variants (ASVs). 
While unfiltered sequences can also be processed with high accuracy, 
the **insect** classification algorithm is relatively slow, 
since it uses a computationally intensive dynamic
programming algorithm to find the likelihood values
of each sequence given the models at each node of the classification tree. 
Hence filtered input datasets are generally be 
much faster to process.


### Installation

To download **insect** from CRAN and load the package, run

```{r, eval = FALSE}
install.packages("insect")
library(insect)
```

To download the latest development version from GitHub, run:

```{r, eval = FALSE}
devtools::install_github("shaunpwilkinson/insect", build_vignettes = TRUE) 
library(insect)
```

```{r, echo = FALSE}
library(insect)
```



### Classifying sequences

Classifiers for some of the more commonly used metabarcoding primer
sets are available here:

<!-- note newlines needed between html tags and code chunk -->

<div class="verysmall">

```{r, echo = FALSE, results='asis'}
tmp <- tempfile(fileext = ".csv")
u <- "https://docs.google.com/spreadsheets/d/1pmTlZBnWIZzxZzS8zm943uGL56qpDPBmoOBht9OzVvQ/export?gid=0&format=csv"
download.file(u, destfile = tmp)
mytab <- read.table(text = readLines(tmp, warn = FALSE), header = TRUE,  sep = ",", stringsAsFactors = FALSE)
mytab <- mytab[order(mytab$Marker, mytab$Target),]
rownames(mytab) <- NULL
knitr::kable(mytab)
```

</div>

To classify a sequence or set of sequences, first read them into R as a "DNAbin" 
list object. FASTA files can be parsed as follows:

```{r, eval = FALSE}
x <- readFASTA("<path-to-file>.fasta")
```

Alternatively users may wish to assign taxon IDs to the output from the 
[DADA2](https://www.nature.com/articles/nmeth.3869)
pipeline, in which case the column names of the ouput
table can be parsed as in the following example:

```{r}
data("samoa") 
x <- char2dna(colnames(samoa))
## name the sequences sequentially
names(x) <- paste0("ASV", seq_along(x))
```

The next step is to download and read in the classifier.
It is important to ensure that the classifier was trained using the
same primer set as that used to generate the query data.
In this example the data were generated from 
autonomous reef monitoring structures in 
American Samoa (ARMS) using the COI metabarcoding primers mlCOIintF and jgHCO2198 
([Leray et al 2013](https://frontiersinzoology.biomedcentral.com/articles/10.1186/1742-9994-10-34)),
and de-noised, filtered and merged following the 
[DADA2 tutorial](https://benjjneb.github.io/dada2/tutorial.html).

The COI classifier was created using the 
MIDORI UNIQUE 20180221 (https://reference-midori.info/download.php)
trainingset, supplemented with around 14,000 non-metazoan 
COI sequences downloaded from GenBank. 

The 140 MB classifier can be downloaded 
and read into R as follows:

```{r}
tmpf <- tempfile()
download.file("https://www.dropbox.com/s/dvnrhnfmo727774/classifier.rds?dl=1", 
              destfile = tmpf, mode = "wb")
classifier <- readRDS(tmpf)
```



There is an option to perform a nearest-neighbor search prior to the 
computationally-expensive recursive model test procedure, which can save 
time and improve resolution ('recall') at lower taxonomic ranks.
Note that this can be a double-edged sword; if multiple species share
an identical or near-identical sequence, and the true taxon of the query sequence
is missing from the trainingset, the algorithm may over-classify the sequence
and return a congeneric taxon.
To perform a nearest-neighbor search with a similarity threshold of 0.99
(meaning any sequence in the trainingset with a similarity greater than or 
equal to 99% is considered a match), set `ping = 0.99`.
To stay on the safe side, we will set `ping = 1` 
(i.e. only sequences with 100% identity are considered matches).


```{r}
out <- classify(x, classifier, threshold = 0.8)
```


<!-- note newlines needed between html tags and code chunk -->

<div class="verysmall">

```{r, echo = FALSE}
knitr::kable(out)
```

</div>



### Further reading

A more detailed overview of the package and its functions can be found
[here](https://rpubs.com/shaunpwilkinson/insect) or by running

```{r, eval = FALSE}
vignette("insect-vignette")
```

### Issues

If you experience a problem using this software please feel free to
raise it as an issue on [GitHub](https://github.com/shaunpwilkinson/insect/issues).


### Acknowledgements

This software was developed at 
[Victoria University of Wellington](https://www.wgtn.ac.nz/) 
with funding from a Rutherford Foundation Postdoctoral Research Fellowship 
award from the Royal Society of New Zealand.