annotateRanges with TxDb object

Author: ischeller

DESCRIPTION

@@ -60,6 +60,8 @@ Suggests:
     rmarkdown,
     testthat,
     covr,
+    TxDb.Hsapiens.UCSC.hg19.knownGene,
+    org.Hs.eg.db,
 LinkingTo:
     Rcpp,
     RcppArmadillo

NAMESPACE

@@ -20,6 +20,7 @@ export(N)
 export(aberrant)
 export(addCountsToFraseRDataSet)
 export(annotateRanges)
+export(annotateRangesWithTxDb)
 export(bamFile)
 export(bestQ)
 export(calculatePSIValues)
@@ -134,6 +135,7 @@ importFrom(GenomeInfoDb,seqnames)
 importFrom(GenomicAlignments,junctions)
 importFrom(GenomicAlignments,readGAlignments)
 importFrom(GenomicAlignments,summarizeJunctions)
+importFrom(GenomicFeatures,genes)
 importFrom(GenomicFeatures,intronsByTranscript)
 importFrom(GenomicFeatures,makeTxDbFromGFF)
 importFrom(GenomicRanges,GRanges)

R/FRASER-package.R

@@ -18,7 +18,7 @@
 #'
 ### GRange/Experiment/bamFile packages
 #' @importFrom BiocGenerics updateObject counts counts<- strand strand<-
-#' @importFrom GenomicFeatures makeTxDbFromGFF intronsByTranscript
+#' @importFrom GenomicFeatures makeTxDbFromGFF intronsByTranscript genes
 #' @importFrom GenomicAlignments junctions readGAlignments summarizeJunctions
 #' @importFrom SummarizedExperiment assay assay<- assays assays<- assayNames
 #'          colData colData<- rowData rowRanges rowRanges<- SummarizedExperiment

R/annotationOfRanges.R

@@ -14,22 +14,32 @@
 #' @param featureName Name of the feature in the FraseRDataSet mcols.
 #' @param biotype The biotype.
 #' @param ensembl The ensembl that should be used. If NULL, the default one is 
-#' used (hsapiens_gene_ensembl, GRCh38).
-#' @param GRCh GRCh version to connect to if not the current GRCh38, 
-#' currently this can only be 37
+#' used (hsapiens_gene_ensembl, GRCh37).
+#' @param GRCh GRCh version to connect to. If this is NULL, then the current 
+#' GRCh38 is used. Otherwise, this can only be 37 (default) at the moment 
+#' (see \code{\link{useEnsebml}}).
+#' @param txdb A TxDb object (default: TxDb.Hsapiens.UCSC.hg19.knownGene).
+#' @param orgDb An orgDb object (default: org.Hs.eg.db).
 #' 
 #' @return FraseRDataSet
 #' 
 #' @examples
 #'
 #' fds <- makeExampleFraseRDataSet()
-#' fds <- annotateRanges(fds)
+#' fds <- annotateRanges(fds, GRCh=NULL)
+#' 
+#' require(TxDb.Hsapiens.UCSC.hg19.knownGene)
+#' txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
+#' require(org.Hs.eg.db)
+#' orgDb <- org.Hs.eg.db
+#' fds <- annotateRangesWithTxDb(fds, txdb=txdb, orgDb=orgDb)
 #'
 #' rowRanges(fds, type="psi5")[,"hgnc_symbol"]
 #'
+#' @rdname annotateRanges
 #' @export
 annotateRanges <- function(fds, feature="hgnc_symbol", featureName=feature,
-            biotype=list("protein_coding"), ensembl=NULL, GRCh=NULL){
+            biotype=list("protein_coding"), ensembl=NULL, GRCh=37){
 
     # check input
     stopifnot(is(fds, "FraseRDataSet"))
@@ -73,6 +83,39 @@ annotateRanges <- function(fds, feature="hgnc_symbol", featureName=feature,
     return(fds)
 }
 
+#' @rdname annotateRanges
+#' @export
+annotateRangesWithTxDb <- function(fds, feature="SYMBOL", 
+                                    featureName="hgnc_symbol",
+                                    txdb=TxDb.Hsapiens.UCSC.hg19.knownGene, 
+                                    orgDb=org.Hs.eg.db){
+ 
+    # check input
+    stopifnot(is(fds, "FraseRDataSet"))
+    if(length(fds) == 0) return(fds)
+    
+    for(i in c("psi3", "psiSite")){
+        # get GRanges object with the split reads which should be annotated
+        gr <- rowRanges(fds, type=i)
+        
+        # get the annotation to compare to
+        anno <- genes(txdb)
+        mcols(anno)[[featureName]] <- 
+            select(orgDb, keys=mcols(anno)[,"gene_id"], columns=feature, 
+                   keytype="ENTREZID")[,feature]
+        anno <- anno[!is.na(mcols(anno)[,featureName]),]
+        anno <- anno[mcols(anno)[,featureName] != "",]
+        if(any(strand(gr) == "*")){
+            strand(anno) <- "*"
+        }
+        
+        # retrieve the feature of interest for the split reads
+        mcols(fds, type=i)[[featureName]] <- 
+            getAnnotationFeature(gr, featureName, anno)
+    }
+    
+    return(fds)   
+}
 
 #'
 #' use biomart to extract the current feature annotation

R/methods-makeExampleObj.R

@@ -2,7 +2,7 @@
 #' Create example data sets for FraseR
 #'
 #' Creates an example data set from files with example counts. 
-#' \Rfunction{makeFittedExampleFraseRDataSet} additionally runs the FRASER 
+#' \code{makeFittedExampleFraseRDataSet} additionally runs the FRASER 
 #' pipeline.
 #' 
 #' @return FraseRDataSet

man/annotateRanges.Rd

@@ -204,7 +204,12 @@ fds <- calculatePSIValues(fds)
 # use PCA to speed up the tutorial
 fds <- FraseR(fds, q=2, correction="PCA")
 
-fds <- annotateRanges(fds)
+# fds <- annotateRanges(fds) # uses biomart
+require(TxDb.Hsapiens.UCSC.hg19.knownGene)
+txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
+require(org.Hs.eg.db)
+orgDb <- org.Hs.eg.db
+fds <- annotateRangesWithTxDb(fds, txdb=txdb, orgDb=orgDb)
 
 # get results: usually, using p-value cutoff is recommended, but on the small 
 # example dataset, this would give no results, so we use z-scores here instead
@@ -477,8 +482,13 @@ cutoff or a cutoff on the $\Delta\Psi$ values between the observed and expected
 $\Psi$ values or both.
 
 <<result table, echo=TRUE>>=
-# annotate
-fds <- annotateRanges(fds)
+# annotate introns with the hgnc symbols of the corresponding gene
+require(TxDb.Hsapiens.UCSC.hg19.knownGene)
+txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
+require(org.Hs.eg.db)
+orgDb <- org.Hs.eg.db
+fds <- annotateRangesWithTxDb(fds, txdb=txdb, orgDb=orgDb)
+# fds <- annotateRanges(fds) # alternative way using biomart
 
 # retrive results
 res <- results(fds, padjCutoff=0.05, zScoreCutoff=NA, deltaPsiCutoff=0.3)