2.1 Step 1. Load the package cleanUpdTSeq, read in the test dataset and then use the function BED2GRangesSeq to convert it to GRanges.

library(cleanUpdTSeq)

## Loading required package: BiocGenerics

## Loading required package: parallel

## 
## Attaching package: 'BiocGenerics'

## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB

## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs

## The following objects are masked from 'package:base':
## 
##     Filter, Find, Map, Position, Reduce, anyDuplicated, append,
##     as.data.frame, basename, cbind, colMeans, colSums, colnames,
##     dirname, do.call, duplicated, eval, evalq, get, grep, grepl,
##     intersect, is.unsorted, lapply, lengths, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, rank, rbind,
##     rowMeans, rowSums, rownames, sapply, setdiff, sort, table,
##     tapply, union, unique, unsplit, which, which.max, which.min

## Loading required package: BSgenome

## Loading required package: S4Vectors

## Loading required package: stats4

## 
## Attaching package: 'S4Vectors'

## The following object is masked from 'package:base':
## 
##     expand.grid

## Loading required package: IRanges

## 
## Attaching package: 'IRanges'

## The following object is masked from 'package:grDevices':
## 
##     windows

## Loading required package: GenomeInfoDb

## Loading required package: GenomicRanges

## Loading required package: Biostrings

## Loading required package: XVector

## 
## Attaching package: 'Biostrings'

## The following object is masked from 'package:base':
## 
##     strsplit

## Loading required package: rtracklayer

## Loading required package: BSgenome.Drerio.UCSC.danRer7

## Loading required package: seqinr

## 
## Attaching package: 'seqinr'

## The following object is masked from 'package:Biostrings':
## 
##     translate

## Loading required package: e1071

testFile <- system.file("extdata", "test.bed", package="cleanUpdTSeq")
testSet <- read.table(testFile, sep="\t", header=TRUE)
peaks <- BED2GRangesSeq(testSet, withSeq=FALSE)

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

peaks <- BED2GRangesSeq(testSet, upstream.seq.ind=7, 
                          downstream.seq.ind=8, withSeq=TRUE)

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

Here is how the test dataset look like.

head(testSet)

##     chr    start     stop        name score strand
## 1 chr10  2965327  2965327 6hpas-22249     1      -
## 2 chr10  2966558  2966558 6hpas-22250     1      -
## 3 chr10  2974251  2974251 6hpas-22251     2      -
## 4 chr10  2978441  2978441 6hpas-22252     1      -
## 5 chr11 16772291 16772291 6hpas-33204     1      -
## 6 chr11 16777848 16777848 6hpas-33205     1      -
##                                   upstream                     downstream
## 1 TCTTCATCATGGTCATCTCGCACCAGAGAGTGTGCCAGGG CAGGAAGTTTTACCTGTCTGTCATTATCGT
## 2 ACCCTGGTGAGGGTATAGAGCTGGTCCAGTGTGCCACGGC AAAGAGGAAAACAGCATTGTTCCTCCTGGA
## 3 TGATTTGTTTGTAACTGATTTTATCTTTTAATAAAAAAGA AAAAAGAAAGTCAAGCCAAGAGGCAAATAC
## 4 GGAGCGCGACCGCATCAACAAAATCTTGCAGGATTATCAG AAGAAAAAGATGGTGAGTTATTATCATTCA
## 5 AGGGAAATAAATACAAAAGAATAAAAATATGATTCATTGT AAGAAAAACACTTTAGCTACAAAAGTCCTT
## 6 ATTTAGTTGGGTATTATTTCAAATAAAGAGAGAGAGAGAC ACAAAAACTACATCAAATTTGAGGACAAAA

2.2 Step2. Build feature vectors for the classifier using the function buildFeatureVector.

The zebrafish genome from BSgenome is used in this example for obtaining surrounding sequences. For a list of other genomes available through BSgenome, please refer to the BSgenome package documentation [2].

library(BSgenome.Drerio.UCSC.danRer7)
testSet.NaiveBayes <- buildFeatureVector(peaks, BSgenomeName=Drerio,
                                         upstream=40, downstream=30, 
                                         wordSize=6, alphabet=c("ACGT"),
                                         sampleType="unknown", 
                                         replaceNAdistance=30, 
                                         method="NaiveBayes",
                                         ZeroBasedIndex=1, fetchSeq=TRUE)

## Not all the sequence names are in the genome.
## Try to convert.
## Not all the sequence names are in the genome.
## Try to convert.

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

2.3 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.

data(data.NaiveBayes)
if(interactive()){
    predictTestSet(data.NaiveBayes$Negative, data.NaiveBayes$Positive, 
                   testSet.NaiveBayes=testSet.NaiveBayes, 
                   outputFile="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.

data(classifier)
testResults <- predictTestSet(testSet.NaiveBayes=testSet.NaiveBayes,
                              classifier=classifier,
                              outputFile=NULL, 
                              assignmentCutoff=0.5)
head(testResults)

##      PeakName prob False/oligodT internally primed    prob True pred.class
## 1 6hpas-78439                         9.999997e-01 3.137403e-07          0
## 2 6hpas-78440                         1.000000e+00 7.431591e-14          0
## 3 6hpas-78441                         4.474675e-06 9.999955e-01          1
## 4 6hpas-78442                         3.085156e-07 9.999997e-01          1
## 5 6hpas-22249                         9.899582e-01 1.004181e-02          0
## 6 6hpas-22250                         1.000000e+00 1.636415e-09          0
##                                UpstreamSeq                   DownstreamSeq
## 1 GGTCATTGTCCTGCAAAATGGACTACTTAACCGAACTGGA GAAGTATAAGAAGTAAGTACATTAAAGCTAC
## 2 TGGATTTAAATAACAAACAAGTTAAATAAAACGATTTGTA AAAAAATAAAACAACTGAAGAAGAAAATGAA
## 3 ATCTGCTTCAAAATGGATGCTCTGTTGAATCCTGAGCTCA GGTAATCTTTCAAGTGCTGCTATTGAGCCAA
## 4 AAATGCTTGCACATAATAAATGTAGGCTTAAAAGATTTCA AAACGTTTGTGAGAGACGGATTTTACTTTGC
## 5 TCTTCATCATGGTCATCTCGCACCAGAGAGTGTGCCAGGG CAGGAAGTTTTACCTGTCTGTCATTATCGTC
## 6 ACCCTGGTGAGGGTATAGAGCTGGTCCAGTGTGCCACGGC AAAGAGGAAAACAGCATTGTTCCTCCTGGAT