[BioC] biomaRt - batch query for chromosome location to gene identifier?

Vincent Carey stvjc at channing.harvard.edu
Fri Nov 12 13:18:15 CET 2010


Here is a solution "entirely in R" provided you have acquired the
hg18 transcript database using GenomicFeatures makeTranscriptDbFromUCSC("hg18")
and run saveFeatures on the result to yield file "hg18.txdb.sqlite":

I put your data in space-delimited text, then transformed to GRanges preserving
chromosome and strand, then find overlaps of your regions with transcripts
for hg18.  The names of these transcripts are given in the UCSC 'knownGene'
vocabulary.  After the session info, we use R 2.10 org.Hs.eg.db to resolve these
to (I think) hg18 entrez ids (we only have hg19 mappings UCSCKG<->Entrez for
current R, to the best of my knowledge, and I believe that even these
are deprecated).

> kemdat = read.table("kemdat.txt", sep=" ", h=TRUE,stringsAsFactors=FALSE)
> library(GenomicRanges)
> kem = with(kemdat, GRanges(ranges=IRanges(Cluster_Begin, Cluster_End), strand=
+  Rle(factor(Strand)), seqnames=Rle(factor(Chromosome))))
> kem
GRanges with 5 ranges and 0 elementMetadata values
    seqnames                 ranges strand |
       <Rle>              <IRanges>  <Rle> |
[1]     chr2 [ 74295624,  74295644]      + |
[2]     chr1 [203949843, 203949866]      - |
[3]     chr8 [103464182, 103464207]      - |
[4]    chr17 [ 53272097,  53272117]      - |
[5]    chr12 [ 11150173,  11150194]      - |

seqlengths
  chr1 chr12 chr17  chr2  chr8
    NA    NA    NA    NA    NA
> library(GenomicFeatures)
> hg18.txdb = loadFeatures("hg18.txdb.sqlite")
> tx18 = transcripts(hg18.txdb)
> kg = values(tx18[ findOverlaps(kem,tx18)@matchMatrix[,2] ])$tx_name
> dput(kg)
c("uc002skj.1", "uc002skk.1", "uc010ffb.1", "uc001hdb.2", "uc003ykr.1",
"uc003yks.1", "uc002ivc.1", "uc009zhp.1", "uc001qzb.2", "uc001qzc.2",
"uc001qze.2", "uc001qzf.1", "uc001qzj.2")
> sessionInfo()
R version 2.13.0 Under development (unstable) (2010-10-29 r53474)
Platform: x86_64-apple-darwin10.4.0/x86_64 (64-bit)

locale:
[1] C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] GenomicFeatures_1.3.6 GenomicRanges_1.3.2   IRanges_1.9.6

loaded via a namespace (and not attached):
[1] BSgenome_1.17.7    Biobase_2.10.0     Biostrings_2.17.47 DBI_0.2-5
[5] RCurl_1.4-3        RSQLite_0.9-2      XML_3.1-1          biomaRt_2.5.1
[9] rtracklayer_1.11.3
>

Now use R 2.10.1 and org.Hs.eg.db 2.3.6 (current is 2.4.5 for devel,
and the UCSCKG
mapping is deprecated)

> kgs = c("uc002skj.1", "uc002skk.1", "uc010ffb.1", "uc001hdb.2", "uc003ykr.1",
+ "uc003yks.1", "uc002ivc.1", "uc009zhp.1", "uc001qzb.2", "uc001qzc.2",
+ "uc001qze.2", "uc001qzf.1", "uc001qzj.2")
> unlist(mget(kgs, revmap(org.Hs.egUCSCKG))
+ )
uc002skj.1 uc002skk.1 uc010ffb.1 uc001hdb.2 uc003ykr.1 uc003yks.1 uc002ivc.1
   "10797"    "10797"    "10797"    "64710"    "51366"    "51366"    "51649"
uc009zhp.1 uc001qzb.2 uc001qzc.2 uc001qze.2 uc001qzf.1 uc001qzj.2
   "11272"     "5554"     "5554"     "5554"     "5545"     "5554"
> sessionInfo()
R version 2.10.1 RC (2009-12-10 r50697)
i386-apple-darwin9.8.0

locale:
[1] C

attached base packages:
[1] stats     graphics  grDevices datasets  tools     utils     methods
[8] base

other attached packages:
[1] org.Hs.eg.db_2.3.6  RSQLite_0.7-3       DBI_0.2-5
[4] AnnotationDbi_1.8.2 Biobase_2.6.1       weaver_1.11.1
[7] codetools_0.2-2     digest_0.4.1



On Thu, Nov 11, 2010 at 10:01 PM, Kemal Akat <kakat at mail.rockefeller.edu> wrote:
> Hi all,
>
> I have a list of mapped sequence reads to hg18 for that I have the exact chromosomal location on NCBI build 36.
>
> Cluster ID              Strand  Chromosome      Cluster_Begin   Cluster_End
> slc754_chr2             +               chr2            74295624                74295644
> slc4695_chr1            -               chr1                    203949843               203949866
> slc2213_chr8            -               chr8            103464182               103464207
> slc1866_chr17   -               chr17           53272097                53272117
> slc1642_chr12   -               chr12           11150173                        11150194
> ...
>
> For the downstream analysis I would like to assign each location an identifier (entrez gene id, ensembl gene id and so forth), and the question is simply if I can use the biomaRt package for this at all?
>
> It is easy for a single entry:
>
>> geneid <-getBM(attributes="entrezgene", filters=c("chromosome_name","start","end", "strand"), values = list(2,74295624, 74295644,1), mart=ensembl54) # ensembl54 is using the archived build 54 = NCBI 36
>> geneid
>  entrezgene
> 1      10797
>>
>
> However, so far I have failed to make a batch query out of it.
>
> I imported/created the following 1 column data frame with the localization formatted as necessary
>
>> tdp
>                                 chromosomal_cocation
> slc754_chr2              2,74295624,74295644,1
> slc4695_chr1    1,203949843,203949866,-1
> slc2213_chr8    8,103464182,103464207,-1
> slc1866_chr17           17,53272097,53272117,-1
> slc1642_chr12           12,11150173,11150194,-1
> ...
>
> I have two points where I failed:
>
> 1) I have not found a single filter that replaces the multiple filters above. When I use "chromosomal_region" as single filter and run:
>
>> geneid <-getBM(attributes="entrezgene", filters="chromosomal_region", values = list(tdp$chromosomal_location, mart=ensembl54)
>
> I get 19733 gene ids; my dataset actually has only 161 locations.
>
> 2) If I use multiple filters like I did above in the first example, "values" has to be a vector and the expression "values = list(tdp$chromosomal_location, mart=ensembl54)" yields a "subscript out of bounds" error.
> I tried splitting the localization infos into separate vectors, i.e. chromosome <- c(2,1,8,17,12,...), start <- c(74295624,....), end <- c(...), strand <- c(...) and modified my query:
>
>> geneid <-getBM(attributes="entrezgene", filters=c("chromosome_name","start","end", "strand"), values = list(chromosome, start, end, strand, mart=ensembl54)
>
> But this seems to combine the information in the different vectors as the result is over 20.000 entries.
>
> Finally, I was thinking of a loop to complete the task, but this has been discouraged by another post in the mailing archive!?
>
> Any help/idea appreciated!
>
> Thank you,
> Kemal
>
> Dr. med. Kemal Akat
> Postdoctoral Fellow
> Laboratory of RNA Molecular Biology
> The Rockefeller University
> 1230 York Avenue, Box #186
> New York, NY 10065
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
>



More information about the Bioconductor mailing list