[BioC] Gene-snp matching

James W. MacDonald jmacdon at med.umich.edu
Mon Sep 14 15:42:38 CEST 2009


Hi Mohamed,

First, you don't need to ask the same question twice with different 
subject lines; the list traffic is typically low over the weekend, so 
you could have waited until today for a response. In addition, you will 
either get a response or not - asking again probably reduces your 
chances of getting an answer.

Mohamed Lajnef wrote:
> Hi  All,
> 
> I have a list of genes and I want the snps matching of these genes (the 
> other way too e.g: Snps ----> genes)
> 
> Is there any R package  to do this ?

Several. You can do this with either an org.Xx.eg.db package in concert 
with a SNPloc package, using biomaRt, or RMySQL with a direct query of 
the UCSC genome browser public MySQL server.

1.)

library(org.Hs.eg.db)
library("SNPlocs.Hsapiens.dbSNP.20080617")
gns <- c("VEGFA","TP53", "BRCA1")
egs <- unlist(mget(gns, org.Hs.egALIAS2EG))
beg <- mget(egs, org.Hs.egCHRLOC)
finish <- mget(egs, org.Hs.egCHRLOCEND)
chr6 <- getSNPlocs("chr6")
chr6[chr6$loc > beg[[1]] & chr6$loc < finish[[1]],]

        RefSNP_id alleles_as_ambig      loc
215667   2010963                S 43846328
215668   6928818                R 43846352
215669     25648                Y 43846955
215670  58159269                Y 43847424
215671  61085192                R 43847452
<snip>

2.)

mart <- useMart("snp","hsapiens_snp")
embl <- unlist(mget(egs, org.Hs.egENSEMBL))
snps <- getBM(c("ensembl_gene_stable_id", "refsnp_id"),"ensembl_gene", 
embl , mart)
  head(snps)
   ensembl_gene_stable_id  refsnp_id
1        ENSG00000012048  rs2355954
2        ENSG00000012048 rs35578914
3        ENSG00000012048  rs3074670
4        ENSG00000012048 rs34130132
5        ENSG00000012048   rs762533
6        ENSG00000012048   rs762532
<snip>

3.)

library(RMySQL)
con <- dbConnect("MySQL", user = "genome", host = 
"genome-mysql.cse.ucsc.edu", dbname = "hg18")
sql <- paste("select name from snp128 where chrom='chr6' and chromStart 
between", beg[[1]], "and", finish[[1]], ";")
dbGetQuery(con, sql)
  dbGetQuery(con, sql)
           name
1    rs2010963
2    rs6928818
3      rs25648
4   rs35440538
5     rs833067
6    rs1885657
7     rs943070
8    rs1413711
<snip>

I leave it as an exercise to figure out how to go the other direction.

Best,

Jim




> 
> Regards
> M
> 

-- 
James W. MacDonald, M.S.
Biostatistician
Douglas Lab
University of Michigan
Department of Human Genetics
5912 Buhl
1241 E. Catherine St.
Ann Arbor MI 48109-5618
734-615-7826



More information about the Bioconductor mailing list