---
title: "The cleanUpdTSeq user's guide"
author: "Sarah Sheppard, Haibo Liu, Jianhong Ou, Nathan Lawson, Lihua J. Zhu"
date: "`r doc_date()`"
package: "`r pkg_ver('cleanUpdTSeq')`"
abstract: >
   This package implements a published Naive Bayes classifier for accurate
   identification of polyadenylation (pA) sites from oligo(dT)-based 3' end
   sequencing such as PAS-Seq, PolyA-Seq and RNA-Seq. The classifer is highly
   accurate and outperforms other heuristic methods. Even though the classfier
   has been trained using Zebrafish data, it works well with other species as
   demonstrated in several organisms including the human, mouse, rat, dog and
   monkey.
vignette: >
  %\VignetteIndexEntry{cleanUpdTSeq Vignette}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
output:
  BiocStyle::html_document
---

# Introduction
3' ends of transcripts have generally been poorly annotated. With the advent of deep sequencing, many methods have been developed to identify 3' ends. The majority of these methods use an oligo-dT primer, which can bind to internal adenine-rich sequences, and lead to artifactual identification of polyadenylation sites. Heuristic filtering methods rely on a certain number of adenines in the genomic sequence downstream of a putative polyadenylation site to remove internal priming events. We introduce a package to provide a robust method to classify putative polyadenylation sites. cleanUpdTSeq uses a na&iuml;ve Bayes classifier, implemented through the **e1071** [1], and sequence features surrounding the putative polyadenylation sites for classification.

The package includes a training dataset constructed from 6 different Zebrafish sequencing dataset, and functions for fetching surrounding sequences using BSgenome [2], building feature vectors and classifying whether the putative polyadenylations site is a true polyadenylation site or a mis-primed false site.

# Citation
If  you use cleanUpdTSeq, please cite:

> Sheppard, S., Lawson, N.D. and Zhu, L.J., 2013. Accurate identification of polyadenylation sites from 3' end deep sequencing using a naive Bayes classifier. Bioinformatics, 29(20), pp.2564-2571 and 2014, 30(4), pp.596, https://doi.org/10.1093/bioinformatics/btt714

# step-by-step guide

Here is a step-by-step guide on how to use the cleanUpdTSeq package to classify a list of putative polyadenylation sites

## Step 1. Load the  cleanUpdTSeq package, and then use the function BED6WithSeq2GRangesSeq to convert the test dataset in a extended BED6 file to a GRanges object with or without sequence information.

```{r 1}
suppressPackageStartupMessages(library(cleanUpdTSeq))
testFile <- system.file("extdata", "test.bed", 
                        package = "cleanUpdTSeq")
peaks <- BED6WithSeq2GRangesSeq(file = testFile, 
                               skip = 1L, withSeq = FALSE)
```

If test dataset contains sequence information already, then use the following command instead.

```{r 2}
peaks <- BED6WithSeq2GRangesSeq(file = testFile, 
                                skip = 1L, withSeq = TRUE)
```

To work with your own test dataset, please set testFile to the file path that contains the putative pA site information.

Here is how the test dataset look like.

```{r 3}
head(read.delim(testFile, header = TRUE, skip = 0))
```
## Step2. Build feature vectors for the classifier using the function buildFeatureVector.
The zebrafish genome from BSgenome is used in this example for obtaining flanking sequences. For a list of other genomes available through BSgenome, please refer to the BSgenome package documentation [2].

```{r 4}
suppressPackageStartupMessages(library(BSgenome.Drerio.UCSC.danRer7))
testSet.NaiveBayes <- buildFeatureVector(peaks, 
                                         genome = Drerio,
                                         upstream = 40L, 
                                         downstream = 30L, 
                                         wordSize = 6L,
                                         alphabet = c("ACGT"),
                                         sampleType = "unknown", 
                                         replaceNAdistance = 30, 
                                         method = "NaiveBayes",
                                         fetchSeq = TRUE,
                                         return_sequences = TRUE)
```

If sequences are present in the test dataset already, then set fetchSeq = FALSE.

## Step 3. Load the training dataset and classify putative polyadenylation sites.

The output file is a tab-delimited file containing the name of the putative
polyadenylation sites, the probability that the putative polyadenylation site is false/oligodT internally primed, the probability the putative polyadenylation site is true, the predicted class based on the assignment cutoff and the sequence surrounding the putative polyadenylation site.


```{r 5}
data(data.NaiveBayes)
if(interactive()){
   out <- predictTestSet(data.NaiveBayes$Negative,
                         data.NaiveBayes$Positive, 
                         testSet.NaiveBayes = testSet.NaiveBayes, 
                         outputFile = file.path(tempdir(), 
                                          "test-predNaiveBayes.tsv"), 
                        assignmentCutoff = 0.5)
}
```
Alternatively, instead of passing in a positive and a negative training dataset, set the parameter classifier to a pre-built **PASclassifier** to speed up the process. To built a **PASclassifier** using the training dataset, please use function **buildClassifier**. A **PASclassifier** named as **classifier** is included in the package which is generated using the included training dataset with upstream = 40, downstream = 30, and wordSize = 6. Please note that in order to use this pre-built classier, you need to build feature vector using buildFeatureVector from your test dataset with the same setting, i.e., upstream = 40, downstream = 30, and wordSize = 6. 

```{r 6}
data(classifier)
testResults <- predictTestSet(testSet.NaiveBayes = testSet.NaiveBayes,
                              classifier = classifier,
                              outputFile = NULL, 
                              assignmentCutoff = 0.5,
                              return_sequences = TRUE)
head(testResults)
```

# References
1. Meyer, D., et al., e1071: Misc Functions of the Department of Statistics (e1071), TU Wien. 2012.

2. Pages, H., BSgenome: Infrastructure for Biostrings-based genome data packages.

3. Sheppard, S., Lawson, N.D. and Zhu, L.J., 2013. Accurate identification of polyadenylation sites from 3' end deep sequencing using a naive Bayes classifier. Bioinformatics, 29(20), pp.2564-2571.

# Session Info
```{r sessionInfo, results = 'asis'}
sessionInfo()
```

