Release 3 8

DESCRIPTION

@@ -1,5 +1,5 @@
 Package: bumphunter
-Version: 1.11.6
+Version: 1.24.5
 Title: Bump Hunter
 Description: Tools for finding bumps in genomic data
 Authors@R: c(person(c("Rafael", "A."), "Irizarry", role = c("cre","aut"), email
@@ -9,7 +9,7 @@ Authors@R: c(person(c("Rafael", "A."), "Irizarry", role = c("cre","aut"), email
     "ctb"), person(c("Andrew", "E."), "Jaffe", role = "ctb"), person("Harris",
     "Jaffee", role = "ctb"), person("Leonardo", "Collado-Torres", role="ctb"))
 Depends:
-    R (>= 2.10),
+    R (>= 3.5),
     S4Vectors (>= 0.9.25),
     IRanges (>= 2.3.23),
     GenomeInfoDb,
@@ -32,9 +32,11 @@ Imports:
     BiocGenerics,
     utils,
     GenomicFeatures,
-    AnnotationDbi
+    AnnotationDbi,
+    stats
 License: Artistic-2.0
 LazyData: yes
 URL: https://github.com/ririzarr/bumphunter
 biocViews: DNAMethylation, Epigenetics, Infrastructure,
-    MultipleComparison
+    MultipleComparison, ImmunoOncology
+RoxygenNote: 6.1.0

NAMESPACE

@@ -8,6 +8,7 @@ importClassesFrom(methods, matrix)
 import(matrixStats)
 import(BiocGenerics)
 import(S4Vectors)
+importFrom(IRanges,PartitioningByEnd)
 importMethodsFrom(IRanges, as.list, as.matrix,
                   diff, distance, distanceToNearest,
                   end, mean,
@@ -24,6 +25,7 @@ importFrom(utils, data, packageDescription)
 importFrom(parallel, mclapply, stopCluster)
 importFrom(GenomicFeatures, exonsBy, cdsBy, transcriptsBy)
 importFrom(AnnotationDbi, mappedkeys)
+importMethodsFrom(AnnotationDbi, mapIds)
 
 export(bumphunter, bumphunterEngine, clusterMaker, boundedClusterMaker,
        dummyData, getSegments, 

R/bumphunter.R

@@ -98,7 +98,7 @@ bumphunterEngine<-function(mat, design, chr = NULL, pos,
         tmp <- .getEstimate(mat = mat, design = design,
                             coef = coef, full = TRUE)
         rawBeta <- tmp$coef
-        weights <- tmp$sigma
+        weights <- 1/tmp$sigma
         rm(tmp)
     } else {
         rawBeta <- .getEstimate(mat = mat, design = design, coef = coef,

R/matchGenes.R

@@ -13,7 +13,7 @@ annotateNearest <- function(x, subject, annotate=TRUE, ...) {
 	stop("select=\"all\" not supported (yet)")
 
     # find nearests
-    NN <- nearest(x, subject, ...)
+    NN <- smarterNearest(x, subject, ...)
     nearest_exists <- !is.na(NN)
 
     # full return value
@@ -89,23 +89,23 @@ annotateNearest <- function(x, subject, annotate=TRUE, ...) {
 #######################################################
 #######################################################
 
-annotateTranscripts <-function(txdb, annotationPackage=NULL, by=c("tx","gene"),codingOnly=FALSE,verbose=TRUE,requireAnnotation=FALSE){
+annotateTranscripts <-function(txdb, annotationPackage=NULL, by=c("tx","gene"),codingOnly=FALSE,verbose=TRUE,requireAnnotation=FALSE, mappingInfo = NULL, simplifyGeneID = FALSE){
 
     if( class(txdb)!="TxDb") stop("txdb must be of class TxDb")
 
     if(is.null(annotationPackage)){
-        species <- species(txdb)
-        species <- strsplit(species," ")[[1]]
-        species <- paste0(substr(species[1],1,1),
-                          tolower(substr(species[2],1,1)))
-        annotationPackage <-  paste("org",species,"eg.db",sep=".")
+        organism <- organism(txdb)
+        organism <- strsplit(organism," ")[[1]]
+        organism <- paste0(substr(organism[1],1,1),
+                           tolower(substr(organism[2],1,1)))
+        annotationPackage <-  paste("org",organism,"eg.db",sep=".")
         if(verbose)
             message(sprintf("No annotationPackage supplied. Trying %s.",annotationPackage))
     }  
 
     if(!require(annotationPackage,character.only=TRUE)){
         if(requireAnnotation){
-            stop("Can't load ",annotationPackage,".\nMake sure library is installed.\nAnd make sure the species argument follows the convention here http://www.bioconductor.org/packages/release/data/annotation/.\nFor example for human use Hs")}else{
+            stop("Can't load ",annotationPackage,".\nMake sure library is installed.\nAnd make sure the organism argument follows the convention here http://www.bioconductor.org/packages/release/data/annotation/.\nFor example for human use Hs")}else{
                 message("Could not load ",annotationPackage,". Will continue without annotation")
                 annotationPackage <- NULL
             }
@@ -136,7 +136,9 @@ annotateTranscripts <-function(txdb, annotationPackage=NULL, by=c("tx","gene"),c
     strand <- unlist(CharacterList(strand(tt)))
 
     RR <- ranges(tt)
-    entrez <- names(RR)
+    ## Useful for switching Gencode IDs to Ensembl IDs
+    if(simplifyGeneID) names(RR) <- gsub('\\..*', '', names(RR))
+    geneid <- names(RR)
     ## starts / ends (left/right endpoints):
     TSS <- unlist(start(RR))
     TSE <- unlist(end(RR))
@@ -175,24 +177,30 @@ annotateTranscripts <-function(txdb, annotationPackage=NULL, by=c("tx","gene"),c
 
     ##now annotate genes
     if(by=="tx") TT <- elementNROWS(tt) else TT <- rep(1,length(tt))
-    Entrez <- Rle(entrez, TT)
+    Geneid <- Rle(geneid, TT)
 
     Tx <- txx
 
     if(!is.null(annotationPackage)){
         if(verbose) message("Annotating genes.")
+
+        if(!is.null(mappingInfo)) {
+            stopifnot(all(c('column', 'keytype', 'multiVals') %in% names(mappingInfo)))
+            geneid <- mapIds(get(annotationPackage), keys = geneid, column= mappingInfo$column, keytype = mappingInfo$keytype, multiVals = mappingInfo$multiVals)
+        }
+
 
         ##Annotate transcrtipt
         ##cant use select cause map is not 1-1 
-        map <- get(gsub("\\.db", "SYMBOL",annotationPackage))
+        map <- get(gsub("\\.db", "SYMBOL", annotationPackage))
         which <- mappedkeys(map)
         symbols <- sapply(as.list(map[which]), paste, collapse=" ")
-        genes <- symbols[entrez]	
+        genes <- symbols[geneid]	
 
-        map <- get(gsub("\\.db", "REFSEQ",annotationPackage))
+        map <- get(gsub("\\.db", "REFSEQ", annotationPackage))
         which <- mappedkeys(map)
         symbols <- sapply(as.list(map[which]), paste, collapse=" ")
-        refseq <- symbols[entrez]	# a handful of NA
+        refseq <- symbols[geneid]	# a handful of NA
 
         Gene <- Rle(genes, TT)
         Refseq <- Rle(refseq, TT)
@@ -204,13 +212,12 @@ annotateTranscripts <-function(txdb, annotationPackage=NULL, by=c("tx","gene"),c
     transcripts=GRanges(ranges=IRanges(start=TSS, end=TSE),
         seqnames=chr, strand=strand,
         seqinfo=seqinfo,
-        CSS, CSE, Tx, Entrez, Gene, Refseq, Nexons, Exons)
+        CSS, CSE, Tx, Geneid, Gene, Refseq, Nexons, Exons)
     if(codingOnly) transcripts <- transcripts[!is.na(values(transcripts)$CSS),]
     attributes(transcripts)$description <- "annotatedTranscripts"
     return(transcripts)
 }
 
-
 matchGenes <- function(x,subject, type=c("any","fiveprime"),
                        promoterDist=2500,
                        skipExons=FALSE,
@@ -225,9 +232,9 @@ matchGenes <- function(x,subject, type=c("any","fiveprime"),
     type=match.arg(type)
 
     if(type=="fiveprime"){
-        map <- nearest(x,resize(subject,width=1))
+        map <- smarterNearest(x,resize(subject,width=1))
     } else{
-        map <- nearest(x,subject)
+        map <- smarterNearest(x,subject)
     }
 
     ind <- which(!is.na(map))
@@ -252,13 +259,13 @@ matchGenes <- function(x,subject, type=c("any","fiveprime"),
     subdist=rep(0,length(x))
     genenames=rep("",length(x))
     geneannotation=rep("",length(x))
-    entrez=rep("",length(x))
+    geneid=rep("",length(x))
     strands=rep("",length(x))
 
 
     genenames[ind]<-as.character(values(subject[map[ind],])$Gene)
     geneannotation[ind]<-as.character(values(subject[map[ind],])$Refseq)
-    entrez[ind]<-as.character(values(subject[map[ind],])$Entrez)
+    geneid[ind]<-as.character(values(subject[map[ind],])$Geneid)
     nexons[ind] <- values(subject[map[ind],])$Nexons
     strands[ind] <- as.character(strand(subject[map[ind],]))
 
@@ -339,7 +346,7 @@ matchGenes <- function(x,subject, type=c("any","fiveprime"),
         if( type[j]%in%c("overlaps 5'","overlaps 3'","inside") & !skipExons ){
 
             ir=IRanges(start=c(S,E),width=1)
-            map2 = nearest(ir, exons)
+            map2 = smarterNearest(ir, exons)
             dist2<-distance(ir,exons[map2])
             pos <- start(ir) ##start end the same
             dist2 <- dist2*ifelse(
@@ -476,9 +483,25 @@ matchGenes <- function(x,subject, type=c("any","fiveprime"),
                       strand=strands,
                       geneL=geneL,
                       codingL=codingL,
-                      Entrez=entrez,
+                      Geneid=geneid,
                       subjectHits=map)
     return(tmp)
+}                               
+
+poverlapWidth <- function(query, subject)  {
+  overlap_score <- pmin(end(query), end(subject)) - 
+    pmax(start(query), start(subject)) + 1L
+
+  pmax(overlap_score, 0L)
+}
+
+smarterNearest <- function(x, subject)  {
+  hits <- nearest(x, subject, select="all")
+  mcols(hits)$ovwidth <- poverlapWidth(x[queryHits(hits)],subject[subjectHits(hits)])
+
+  partitioning <- PartitioningByEnd(queryHits(hits),NG=queryLength(hits))
+  ovwidth <- relist(mcols(hits)$ovwidth, partitioning)
+  hits <- hits[unlist(ovwidth==max(ovwidth))]
+  selectHits(hits, select="arbitrary")
 }
 
-

README.md

@@ -1,4 +1,4 @@
-This is the developer version of the [bumphunter](http://bioconductor.org/packages/devel/bioc/html/bumphunter.html) R/Bioconductor package.
+This is the developer version of the [bumphunter](http://bioconductor.org/packages/devel/bioc/html/bumphunter.html) R/Bioconductor package. 
 
 [Bioconductor build/check report](http://master.bioconductor.org/checkResults/devel/bioc-LATEST/bumphunter/)
 

man/annotateTranscripts.Rd

@@ -7,7 +7,7 @@ Annotate transcripts
 Annotate transcripts
 }
 \usage{
-annotateTranscripts(txdb, annotationPackage = NULL, by = c("tx","gene"), codingOnly=FALSE, verbose = TRUE, requireAnnotation = FALSE)
+annotateTranscripts(txdb, annotationPackage = NULL, by = c("tx","gene"), codingOnly=FALSE, verbose = TRUE, requireAnnotation = FALSE, mappingInfo = NULL, simplifyGeneID = FALSE)
 }
 \arguments{
   \item{txdb}{
@@ -16,7 +16,7 @@ annotateTranscripts(txdb, annotationPackage = NULL, by = c("tx","gene"), codingO
   \item{annotationPackage}{
     An annotation data package from which to obtain gene/transcript
     annotation. For example \code{org.Hs.eg.db}. If none is provided the
-    function tries to infer it from \code{species(txdb)} and if it can't it
+    function tries to infer it from \code{organism(txdb)} and if it can't it
     proceeds without annotation unless \code{ requireAnnotation = TRUE}.
   }
   \item{by}{
@@ -34,6 +34,17 @@ annotateTranscripts(txdb, annotationPackage = NULL, by = c("tx","gene"), codingO
     logical value. If 'TRUE' function will stop if no annotation package
     is successfully loaded.
   }
+  \item{mappingInfo}{
+    a named list with elements 'column', 'keytype' and 'multiVals'. If specified
+    this information will be used with \link[AnnotationDbi]{mapIds} when mapping
+    the gene ids using \code{annotationPackage}. This is useful when working
+    with a \code{txdb} object from ENSEMBL or GENCODE among other databases.
+  }
+  \item{simplifyGeneID}{
+    logical value. If 'TRUE', gene ids will be shortened to before a dot is
+    present in the id. This is useful for changing GENCODE gene ids to
+    ENSEMBL ids.
+  }
 }
 \details{
   This function prepares a \code{GRanges} for the \code{\link{matchGenes}}
@@ -51,7 +62,8 @@ annotateTranscripts(txdb, annotationPackage = NULL, by = c("tx","gene"), codingO
   number of exons, \code{Exons} is an \code{IRanges} with the exon information.
 }
 \author{
-  Harris Jaffee and Rafael A. Irizarry
+  Harris Jaffee and Rafael A. Irizarry.
+  'mappingInfo' and 'simplifyGeneID' contributed by Leonardo Collado-Torres.
 }
 
 \seealso{
@@ -66,4 +78,4 @@ annotateTranscripts(txdb, annotationPackage = NULL, by = c("tx","gene"), codingO
     genes <- annotateTranscripts(TxDb.Hsapiens.UCSC.hg19.knownGene,annotation="org.Hs.eg.db")
 
 }
-}
+}