---
title: "phylotypr"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{phylotypr}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
```
In the microbiome field, 16S rRNA gene sequencing has been a popular method of characterizing the types of bacteria found in a community. The most frequently used tool for classifying these sequences has been the [naive Bayesian classifier](https://doi.org/10.1128/AEM.00062-07) that used to be accessible through the Ribosomal Database Project and is currently found in tools like mothur and QIIME. For the first time, the `{phylotypr}` package provides an R-based implementation of that tool. There are minimal dependencies and the easy to use code executes rapidly.
## Getting started with phylotypr
To get started you'll need to install `{phylotypr}` using either of these commands:
```{r, eval = FALSE}
# devtools::install_github("mothur/phylotypr") #install development version
install.packages("phylotypr") # install stable version from CRAN
```
To get access to the packages function and a demonstration version of the training data, you'll need to run the `library()` function.
```{r setup}
library(phylotypr)
```
## Building the database
Before we can classify a sequence, we need to generate the database. This can be done using the `build_kmer_database()` function and takes several seconds to execute. `{phylotypr}` includes two databases: `rdp_trainset9` and `rdp_trainset9_pds`. If you want others, you are encouraged to install the `{phylotyprrefdata}` package from GitHub. You can also create your own reference set using `create_database()`. For more details, see the "Alternative databases" section at the end of this vignette. We'll use `rdp_trainset9_pds` for this tutorial.
```{r}
db <- build_kmer_database(rdp_trainset9_pds)
```
By default, `build_kmer_database()` will use 8-nt long kmers to calculate its probabilities. This size can be altered by the user if needed, but 8-nt performs very well and shouldn't need to be adjusted.
## Classifying a single unknown sequence
To classify a sequence we need the DNA sequence to be a character object like `unknown` is below.
```{r}
unknown <- "TACGGAGGATGCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGATCGTTAAGTCAGTGGTCAAATTGAGGGGCTCAACCCCTTCCCGCCATTGAAACTGGCGATCTTGAGTGGAAGAGAAGTATGCGGAATGCGTGGTGTAGCGGTGAAATGCATAGATATCACGCAGAACCCCGATTGCGAAGGCAGCATGCCGGCTTCCTACTGACGCTCATGCACGAAAGTGTGGGTAACGAACAGG"
```
The `unknown` object can then be classified using `classify_sequences()` along with the database object from above (i.e., `db`). Because the classification algorithm makes use of a random number generator, I strongly encourage you to set the seed to get consistent results. This needs to be set in `classify_sequences()` and not with `set.seed()`. The default seed is 1.
```{r, echo = TRUE, eval = FALSE}
consensus <- classify_sequences(sequence = unknown, database = db, seed = 123)
```
```{r, echo = FALSE}
consensus <- classify_sequences(sequence = unknown, database = db, seed = 123,
num_processors = 2)
```
The output of `classify_sequences` is a data frame with taxonomic levels from kingdom (or domain) down to the genus and the confidence scores for that sequence in each of those taxonomic levels
```{r}
consensus
```
In this example, you'll notice that the family and genus assignment have confidence scores below 80%. This is the value commonly used in mothur and which was used on the RDP website. To filter the classification to that level of confidence, you can use the `filter_taxonomy()` function:
```{r}
filtered <- filter_taxonomy(consensus)
```
You should see that the "Porphyromonadaceae" and "Barnesiella" names are removed along with its confidence score.
```{r}
filtered
```
To convert that list format to a string, you can use the `print_taxonomy()` function:
```{r}
print_taxonomy(filtered)
```
## Classifying multiple unknown sequences
Most likely you will have multiple sequences that you would like to classify. The `{phylotypr}` package installs with an example FASTA-formatted file, "miseq_sop.fasta". The code below shows you can read in a FASTA-formatted file and classify the sequences. Note that for classifying your own data you would replace `phylotypr_example(...)` with the path to your own data. By default `classify_sequences()` will use as many processors as your computer will allow. We will show how to limit the number of processors later in this vignette:
```{r, echo = TRUE, eval = FALSE}
read_fasta(phylotypr_example("miseq_sop.fasta.gz")) |>
classify_sequences(database = db, seed = 123) |>
print_taxonomy() |>
head()
```
```{r, echo = FALSE}
read_fasta(phylotypr_example("miseq_sop.fasta.gz")) |>
classify_sequences(database = db, seed = 123, num_processors = 2) |>
print_taxonomy() |>
head()
```
Alternatively, you may have your data stored as a `{strollur}` object you can do the following
```{r, echo = TRUE, eval = FALSE}
strollur::new_dataset(dataset_name = "sop") |>
strollur::add(
table = phylotypr::read_fasta(phylotypr_example("miseq_sop.fasta.gz"))
) |>
classify_sequences(db, seed = 123) |>
print_taxonomy() |>
head()
```
```{r, echo = FALSE}
strollur::new_dataset(dataset_name = "sop") |>
strollur::add(
table = phylotypr::read_fasta(phylotypr_example("miseq_sop.fasta.gz"))
) |>
classify_sequences(database = db, seed = 123, num_processors = 2) |>
print_taxonomy() |>
head()
```
## Classifying multiple unknown sequences in parallel
The `{phylotypr}` package allows you to run `classify_sequences()` in parallel. So, if your computer has 8 processors, `classify_sequences()` will split your data into 8 chunks and then classify each chunk separately before pulling it all together. To see how many processors you have available, you can run the following:
```{r}
parallelly::availableCores()
```
By default, `classify_sequences()` will use all of those processors. If you need to limit the number of processors, you can use the `num_processors` argument
```{r}
read_fasta(phylotypr_example("miseq_sop.fasta.gz")) |>
classify_sequences(database = db, seed = 123, num_processors = 2) |>
print_taxonomy() |>
head()
```
## Run-to-run variation
Consider the results of classifying `unknown` three times with three different seeds:
```{r}
classify_sequences(sequence = unknown, database = db, seed = 1) |>
print_taxonomy()
classify_sequences(sequence = unknown, database = db, seed = 2) |>
print_taxonomy()
classify_sequences(sequence = unknown, database = db, seed = 3) |>
print_taxonomy()
```
You'll notice that the classification of the three replicates varies slightly. Again, this is due to the use of the random number generator and getting confidence scores that bounce around 80%, which results in different levels of filtering. To get a more stable classification, you should be setting your random number generator seed with a consistent value of `seed`.
To get a more precise confidence score, you could set `num_bootstraps = 1000` in `classify_sequences()`. Be forewarned that the function will take 10-times longer to run with 10-times the number of bootstraps!
```{r}
classify_sequences(sequence = unknown, database = db, num_bootstraps = 1000,
seed = 123) |>
print_taxonomy()
```
## Alternative databases
The `{phylotypr}` package ships with the RDP's v.9 of their training data. This is relatively small and old (2010) relative to their latest versions. You are encouraged to install newer versions of the RDP, greengenes, and SILVA databases from the `{phylotyprrefdata}` package on GitHub. Note that installing the package will take about 20 minutes to install. If it sits at "moving datasets to lazyload DB" for a long time, this is normal :)
```{r, eval = FALSE}
devtools::install_github("mothur/phylotyprrefdata")
library(phylotyprrefdata)
```
The following will list the references that are available in `{phylotyprrefdata}`:
```{r, eval = FALSE}
data(package = "phylotyprrefdata")
```