This package uses the Naive Bayes classifier (from e1071) to assign probability values to putative polyadenylation sites (pA sites) based on training data from zebrafish. This will allow the user to separate true, biologically relevant pA sites from false, oligodT primed pA sites.
cleanUpdTSeq 1.16.0
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ï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.
A paper has been submitted to Bioinformatics and currently under revision [3].
Here is a step-by-step guide on using cleanUpdTSeq to classify a list of putative polyadenylation sites
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, cbind, colMeans, colSums, colnames, 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## 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: e1071testFile <- 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 ACAAAAACTACATCAAATTTGAGGACAAAAThe 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.
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 AAAGAGGAAAACAGCATTGTTCCTCCTGGATMeyer, D., et al., e1071: Misc Functions of the Department of Statistics (e1071), TU Wien. 2012.
Pages, H., BSgenome: Infrastructure for Biostrings-based genome data packages.
Sarah Sheppard, Nathan D. Lawson, and Lihua Julie Zhu. 2013. Accurate identification of polyadenylation sites from 3’ end deep sequencing using a naïve Bayes classifier. Bioinformatics. Under revision
sessionInfo()R version 3.4.2 (2017-09-28) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Ubuntu 16.04.3 LTS
Matrix products: default BLAS: /home/biocbuild/bbs-3.6-bioc/R/lib/libRblas.so LAPACK: /home/biocbuild/bbs-3.6-bioc/R/lib/libRlapack.so
locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages: [1] stats4 parallel stats graphics grDevices utils datasets [8] methods base
other attached packages: [1] cleanUpdTSeq_1.16.0 e1071_1.6-8
[3] seqinr_3.4-5 BSgenome.Drerio.UCSC.danRer7_1.4.0 [5] BSgenome_1.46.0 rtracklayer_1.38.0
[7] Biostrings_2.46.0 XVector_0.18.0
[9] GenomicRanges_1.30.0 GenomeInfoDb_1.14.0
[11] IRanges_2.12.0 S4Vectors_0.16.0
[13] BiocGenerics_0.24.0 BiocStyle_2.6.0
loaded via a namespace (and not attached): [1] Rcpp_0.12.13 compiler_3.4.2
[3] class_7.3-14 bitops_1.0-6
[5] tools_3.4.2 zlibbioc_1.24.0
[7] digest_0.6.12 evaluate_0.10.1
[9] lattice_0.20-35 Matrix_1.2-11
[11] DelayedArray_0.4.0 yaml_2.1.14
[13] GenomeInfoDbData_0.99.1 stringr_1.2.0
[15] knitr_1.17 ade4_1.7-8
[17] rprojroot_1.2 grid_3.4.2
[19] Biobase_2.38.0 XML_3.98-1.9
[21] BiocParallel_1.12.0 rmarkdown_1.6
[23] bookdown_0.5 magrittr_1.5
[25] backports_1.1.1 Rsamtools_1.30.0
[27] htmltools_0.3.6 matrixStats_0.52.2
[29] GenomicAlignments_1.14.0 SummarizedExperiment_1.8.0 [31] stringi_1.1.5 RCurl_1.95-4.8