Best way to get list of SNPs by gene id?

Question

I have a long data frame of genes and various forms of ids for them (e.g. OMIM, Ensembl, Genatlas). I want to get the list of all SNPs that are associated with each gene. (This is the reverse of this question.)

So far, the best solution I have found is using the biomaRt package (bioconductor). There is an example of the kind of lookup I need to do here. Fitted for my purposes, here is my code:

library(biomaRt)

#load the human variation data
variation = useEnsembl(biomart="snp", dataset="hsapiens_snp")

#look up a single gene and get SNP data
getBM(attributes = c(
  "ensembl_gene_stable_id",
  'refsnp_id',
  'chr_name',
  'chrom_start',
  'chrom_end',
  'minor_allele',
  'minor_allele_freq'),
  filters = 'ensembl_gene',
  values ="ENSG00000166813",
  mart = variation
)

This outputs a data frame that begins like this:

  ensembl_gene_stable_id  refsnp_id chr_name chrom_start chrom_end minor_allele minor_allele_freq
1        ENSG00000166813  rs8179065       15    89652777  89652777            T          0.242412
2        ENSG00000166813  rs8179066       15    89652736  89652736            C          0.139776
3        ENSG00000166813 rs12899599       15    89629243  89629243            A          0.121006
4        ENSG00000166813 rs12899845       15    89621954  89621954            C          0.421126
5        ENSG00000166813 rs12900185       15    89631884  89631884            A          0.449681
6        ENSG00000166813 rs12900805       15    89631593  89631593            T          0.439297

(4612 rows)

The code works, but the running time is extremely long. For the above, it takes about 45 seconds. I thought maybe this was related to the allele frequencies, which the server perhaps calculated on the fly. But looking up the bare minimum of only the SNPs rs ids takes something like 25 seconds. I have a few thousand genes, so this would take an entire day (assuming no timeouts or other errors). This can't be right. My internet connection is not slow (20-30 mbit).

I tried looking up more genes per query. This did dot help. Looking up 10 genes at once is roughly 10 times as slow as looking up a single gene.

What is the best way to get a vector of SNPs that associated with a vector of gene ids?

If I could just download two tables, one with genes and their positions and one with SNPs and their positions, then I could easily solve this problem using dplyr (or maybe data.table). I haven't been able to find such tables.

Answer 1

Since you're using R, here's an idea that uses the package rentrez. It utilizes NCBI's Entrez database system and in particular the eutils function, elink. You'll have to write some code around this and probably tweak parameters, but could be a good start.

library(rentrez)
# for converting gene name -> gene id
gene_search <- entrez_search(db="gene", term="(PTEN[Gene Name]) AND Homo sapiens[Organism]", retmax=1)
geneId <- gene_search$ids 
# elink function
snp_links <- entrez_link(dbfrom='gene', id=geneId, db='snp')

# access results with $links
length(snp_links$links$gene_snp)
5779

head(snp_links$links$gene_snp)
'864622690' '864622594' '864622518' '864622451' '864622387' '864622341' 

I suggest you manually double-check that the number of SNPs is about what you'd expect for your genes of interest -- you may need to drill down further and limit by transcript, etc...

For multiple gene ids:

multi_snp_links <- entrez_link(dbfrom='gene', id=c("5728", "374654"), db='snp', by_id=TRUE)
lapply(multi_snp_links, function(x) head(x$links$gene_snp))
1.    '864622690' '864622594' '864622518' '864622451' '864622387' '864622341' 
2.    '797045093' '797044466' '797044465' '797044464' '797044463' '797016353' 

The results are grouped by gene with by_id=TRUE