Merge pull request #119 from charles-plessy/annot_from_TxDB

annotateCTSS method for TxDB objects

Closes #30, #53

DESCRIPTION

@@ -18,6 +18,7 @@ Imports:
     formula.tools,
     GenomeInfoDb,
     GenomicAlignments,
+    GenomicFeatures,
     GenomicRanges (>= 1.37.16),
     ggplot2 (>= 2.2.0),
     gtools,

NAMESPACE

@@ -105,6 +105,10 @@ importFrom(GenomeInfoDb,seqnames)
 importFrom(GenomeInfoDb,sortSeqlevels)
 importFrom(GenomicAlignments,qwidth)
 importFrom(GenomicAlignments,readGAlignments)
+importFrom(GenomicFeatures,exons)
+importFrom(GenomicFeatures,genes)
+importFrom(GenomicFeatures,promoters)
+importFrom(GenomicFeatures,transcripts)
 importFrom(GenomicRanges,GPos)
 importFrom(GenomicRanges,GRanges)
 importFrom(GenomicRanges,GRangesList)

R/Annotations.R

@@ -440,6 +440,8 @@ msScope_annotation <- function(libs) {
 #' annotateCTSS(exampleCAGEexp, exampleZv9_annot)
 #' colData(exampleCAGEexp)
 #' 
+#' @importFrom GenomicFeatures genes
+#' 
 #' @export
 
 setGeneric("annotateCTSS", function(object, ranges, upstream=500, downstream=500)
@@ -459,6 +461,22 @@ setMethod("annotateCTSS", c("CAGEexp", "GRanges"), function (object, ranges, ups
   object
 })
 
+#' @rdname annotateCTSS
+
+setMethod("annotateCTSS", c("CAGEexp", "TxDb"), function (object, ranges){
+  g <- genes(ranges)
+  g$gene_name <- g$gene_id
+  annotateCTSS(object, g)
+  CTSScoordinatesGR(object)$genes      <- ranges2genes(CTSScoordinatesGR(object), g)
+  CTSScoordinatesGR(object)$annotation <- txdb2annot(CTSScoordinatesGR(object), ranges)
+
+  annot <- sapply( CTSStagCountDF(object)
+                   , function(X) tapply(X, CTSScoordinatesGR(object)$annotation, sum))
+  colData(object)[levels(CTSScoordinatesGR(object)$annotation)] <- DataFrame(t(annot))
+  validObject(object)
+  object
+})
+
 #' `annotateTagClusters` annotates the _tag clusters_ of a _CAGEr_ object.
 #'   
 #' @return `annotateTagClusters` returns the input object with the same
@@ -619,6 +637,33 @@ ranges2annot <- function(ranges, annot, upstream=500, downstream=500) {
   Rle(annot)
 }
 
+#' @name txdb2annot
+#' 
+#' @rdname ranges2annot
+#' 
+#' @importFrom GenomicFeatures promoters exons transcripts
+
+txdb2annot <- function(ranges, annot) {
+  findOverlapsBool <- function(A, B) {
+    overlap <- findOverlaps(A, B)
+    overlap <- as(overlap, "List")
+    any(overlap)
+  }
+
+  classes <- c("promoter", "exon", "intron", "unknown")
+  p <- findOverlapsBool(ranges, trim(suppressWarnings(promoters(annot, 500, 500))))
+  e <- findOverlapsBool(ranges, exons(annot))
+  t <- findOverlapsBool(ranges, transcripts(annot))
+  annot <- sapply( 1:length(ranges), function(i) {
+    if      (p[i]) {classes[1]}
+    else if (e[i]) {classes[2]}
+    else if (t[i]) {classes[3]}
+    else           {classes[4]}
+  })
+
+  annot <- factor(annot, levels = classes)
+  Rle(annot)
+}
 
 #' ranges2genes
 #'