I'm attempting to use bumphunter::annotateNearest() with Canis familiaris 3.1 (see https://github.com/ririzarr/bumphunter/issues/1). annotateNearest() appears to accept an object with a set of transcripts and annotate against that, even if the transcripts are from a non-human genome, but there isn't specific documentation about how to do this. My attempt is below along with the issue I've bumped up against.

I built my transcripts object like so:

library(TxDb.Cfamiliaris.UCSC.canFam3.xenoRefGene)
library(RMySQL)
library(bumphunter)
tx <- transcripts(TxDb.Cfamiliaris.UCSC.canFam3.xenoRefGene)
names(tx) <- unlist(mcols(tx)[,2])

con <- dbConnect("MySQL", host="genome-mysql.cse.ucsc.edu", user="genome", dbname="canFam3")
refGene <- dbGetQuery(con, "select * from xenoRefGene;")
refGene <- refGene[match(names(tx), refGene$name),]
Nexons = sapply(strsplit(refGene$exonStarts, ","), length)
Exons <- IRangesList(start = lapply(strsplit(refGene$exonStarts, ","), function(x) x[x != "" & !is.na(x)]), end = lapply(strsplit(refGene$exonEnds, ","), function(x) x[x != "" & !is.na(x)]))
mcols(tx) <- DataFrame(CSS = refGene$cdsStart, CSE = refGene$cdsEnd, Tx = refGene$name, Gene = refGene$name2, Nexons = Nexons, Exons = Exons)

I annotated like so:

annotation <- annotateNearest(regions$regions, tx)

(where regions$regions is derfinder output -- I am happy to provide more info there if anyone is curious).

The results are somewhat annotated:

> head(annotation) dist matchIndex type amountOverlap insideDist size1 size2 1 0 15454 inside NA 59 125 5819 2 0 10018 inside NA -1811 141 15559 3 0 10017 inside NA -34060 118 68867 4 0 13154 inside NA -22676 93 169910 5 0 876 inside NA -30307 59 157801 6 0 13167 inside NA 2095 114 12769

I am surprised that the annotation does not include the gene name or RefSeq ID, as the docs suggest it will. My suspicion is that I did not hand enough information along in my tx object, but the object does contain seqnames, and those weren't passed back in the annotation object.

> tx GRanges object with 254017 ranges and 6 metadata columns: seqnames ranges strand | CSS CSE <Rle> <IRanges> <Rle> | <numeric> <numeric> NM_011767 chr1 [363561, 368894] + | 75038343 75122798 NM_001011177 chr1 [363601, 367239] + | 363606 367239 NM_001090507 chr1 [363601, 367239] + | 363606 367239 NM_199558 chr1 [363604, 367325] + | 363671 367325 NM_001044283 chr1 [440183, 441303] + | 65252835 65294349 ... ... ... ... ... ... ... NM_001090227 chrX [123325049, 123390403] - | 123325048 123390403 NM_001102830 chrX [123325049, 123390403] - | 123325048 123390403 NM_001135068 chrX [123325049, 123390403] - | 123325048 123390403 NM_001012575 chrX [123325049, 123408056] - | 123325048 123408056 NM_001184797 chrX [123326032, 123438981] - | 123326060 123408176 Tx Gene Nexons <character> <character> <integer> NM_011767 NM_011767 Zfr 30 NM_001011177 NM_001011177 zfr 16 NM_001090507 NM_001090507 zfr 17 NM_199558 NM_199558 zfr 15 NM_001044283 NM_001044283 Snx3 10 ... ... ... ... NM_001090227 NM_001090227 MGC68497 6 NM_001102830 NM_001102830 tmlhe 7 NM_001135068 NM_001135068 tmlhe 7 NM_001012575 NM_001012575 TMLHE 7 NM_001184797 NM_001184797 TMLHE 7 Exons <IRangesList> NM_011767 [75038283, 75038380] [75038682, 75038775] [75056588, 75056878] ... NM_001011177 [363600, 363660] [363682, 363706] [363729, 364308] ... NM_001090507 [363600, 363660] [363682, 363706] [363729, 364308] ... NM_199558 [363603, 363703] [363723, 364293] [364310, 364412] ... NM_001044283 [65252221, 65252345] [65252377, 65252423] [65252439, 65252532] ... ... ... ------- seqinfo: 3268 sequences (1 circular) from canFam3 genome

I would love advice about how to proceed -- is bumphunter behaving as expected (and it's my problem) or is this a bug?

session_info to follow.

Thanks,
Jessica

Jessica, I've made some major changes to bumphunter to accommodate this. If you install the latest devel version 1.7.3 (still being tested) the matchGenes and annotateNearest will no longer assume hg19. Here is an example:

library(bumphunter)
islands=makeGRangesFromDataFrame(read.delim("http://rafalab.jhsph.edu/CGI/model-based-cpg-islands-hg19.txt")[1:100,])

library("TxDb.Hsapiens.UCSC.hg19.knownGene")

genes <- annotateTranscripts(TxDb.Hsapiens.UCSC.hg19.knownGene, annotationPackage="org.Hs.eg.db")
tab<- matchGenes(islands,genes)

Note that if you don't specify the annotationPackage argument the function will try to figure out what library to load using the species method. If it can't find an annotation package it will continue and return an object with no gene annotation but all the other components. You should be able to substitute the human TxDb for TxDb.Cfamiliaris.UCSC.canFam3.xenoRefGene

Please try it out and report back if you run into trouble.

Note we are not making this back compatible. Apologies if this breaks code, but the previous way had to many problems to keep around. 

ps - You might need to install from gitbub like this devtools::install_github("ririzarr/bumphunter")

Hi Jessica,

You can now use matchGenes() (similar to what annotateTranscripts used to do) with non-human genomes as Rafa already mentioned. Basically, use:

library('GenomicFeatures')
library('bumphunter')
can_txdb <- makeTranscriptDbFromUCSC("canFam3", "refGene")
can_genes <- annotateTranscripts(can_txdb, 'org.Cf.eg.db')
annotation <- matchGenes(x, can_genes)

where 'x' is a data.frame or GRanges. For example:

> matchGenes(head(can_genes), can_genes)
          name                                          annotation
1        KRT12                           NM_001082421 NP_001075890
2     CEACAM23                           NM_001097552 NP_001091021
3     CEACAM23                           NM_001097552 NP_001091021
4     CEACAM28 NM_001097546 NP_001091015 XM_005615767 XP_005615824
5     CEACAM24                           NM_001097554 NP_001091023
6 LOC100049001                           NM_001097555 NP_001091024
     description region distance      subregion insideDistance
1         covers covers        0 covers exon(s)             NA
2         covers covers        0 covers exon(s)             NA
3         covers covers        0 covers exon(s)             NA
4 covers exon(s) inside        5 covers exon(s)              0
5         covers covers        0 covers exon(s)             NA
6         covers covers        0 covers exon(s)             NA
  exonnumber nexons            UTR strand geneL codingL    Entrez
1         NA      8           <NA>      +  5581    5190 100034663
2         NA      9           <NA>      + 16264   16245 100048998
3         NA      9           <NA>      - 32617   16245 100048998
4          1      9 overlaps 3'UTR      + 22076   22029    484473
5         NA      6           <NA>      + 10598   10434 100049000
6         NA      5           <NA>      +  1424    1266 100049001
  subjectHits
1           1
2           2
3           3
4        1740
5           5
6           6
> 

Note that I proposed some fixes for some minor bugs at https://github.com/ririzarr/bumphunter/pull/3

Cheers,

Leo

> devtools::session_info()
Session info------------------------------------------------------------
 setting  value                                             
 version  R Under development (unstable) (2014-11-01 r66923)
 system   x86_64, darwin10.8.0                              
 ui       AQUA                                              
 language (EN)                                              
 collate  en_US.UTF-8                                       
 tz       America/New_York                                  

Packages----------------------------------------------------------------
 package           * version  date       source        
 AnnotationDbi     * 1.29.13  2015-01-14 Bioconductor  
 base64enc           0.1.2    2014-06-26 CRAN (R 3.2.0)
 BatchJobs           1.4      2014-09-24 CRAN (R 3.2.0)
 BBmisc              1.7      2014-06-21 CRAN (R 3.2.0)
 Biobase           * 2.27.1   2014-12-20 Bioconductor  
 BiocGenerics      * 0.13.4   2014-12-31 Bioconductor  
 BiocParallel        1.1.11   2015-01-10 Bioconductor  
 biomaRt             2.23.5   2014-11-22 Bioconductor  
 Biostrings          2.35.7   2014-12-13 Bioconductor  
 bitops              1.0.6    2013-08-17 CRAN (R 3.2.0)
 brew                1.0.6    2011-04-13 CRAN (R 3.2.0)
 bumphunter        * 1.7.4    2015-01-15 Bioconductor  
 checkmate           1.5.0    2014-10-19 CRAN (R 3.2.0)
 codetools           0.2.9    2014-08-21 CRAN (R 3.2.0)
 DBI               * 0.3.1    2014-09-24 CRAN (R 3.2.0)
 devtools            1.6.1    2014-10-07 CRAN (R 3.2.0)
 digest              0.6.4    2013-12-03 CRAN (R 3.2.0)
 doRNG               1.6      2014-03-07 CRAN (R 3.2.0)
 fail                1.2      2013-09-19 CRAN (R 3.2.0)
 foreach           * 1.4.2    2014-04-11 CRAN (R 3.2.0)
 GenomeInfoDb      * 1.3.12   2014-12-22 Bioconductor  
 GenomicAlignments   1.3.24   2015-01-14 Bioconductor  
 GenomicFeatures   * 1.19.10  2015-01-09 Bioconductor  
 GenomicRanges     * 1.19.33  2015-01-13 Bioconductor  
 IRanges           * 2.1.35   2015-01-07 Bioconductor  
 iterators         * 1.0.7    2014-04-11 CRAN (R 3.2.0)
 lattice             0.20.29  2014-04-04 CRAN (R 3.2.0)
 locfit            * 1.5.9.1  2013-04-20 CRAN (R 3.2.0)
 matrixStats         0.10.3   2014-10-15 CRAN (R 3.2.0)
 org.Cf.eg.db      * 3.0.0    2014-09-26 Bioconductor  
 pkgmaker            0.22     2014-05-14 CRAN (R 3.2.0)
 R.methodsS3         1.6.1    2014-01-05 CRAN (R 3.2.0)
 RCurl               1.95.4.3 2014-07-29 CRAN (R 3.2.0)
 registry            0.2      2012-01-24 CRAN (R 3.2.0)
 rngtools            1.2.4    2014-03-06 CRAN (R 3.2.0)
 Rsamtools           1.19.26  2015-01-13 Bioconductor  
 RSQLite           * 1.0.0    2014-10-25 CRAN (R 3.2.0)
 rstudioapi          0.1      2014-03-27 CRAN (R 3.2.0)
 rtracklayer         1.27.7   2015-01-14 Bioconductor  
 S4Vectors         * 0.5.16   2015-01-07 Bioconductor  
 sendmailR           1.2.1    2014-09-21 CRAN (R 3.2.0)
 stringr             0.6.2    2012-12-06 CRAN (R 3.2.0)
 XML                 3.98.1.1 2013-06-20 CRAN (R 3.2.0)
 xtable              1.7.4    2014-09-12 CRAN (R 3.2.0)
 XVector           * 0.7.3    2014-11-24 Bioconductor  
 zlibbioc            1.13.0   2014-10-14 Bioconductor 

Thank you so much! I am in the middle of final exams at the moment but I will check it out in the next few days and report back.

I really appreciate the work.

Jessica