r - 通过基因ID获取SNP列表的最佳方法？

我有一个很长的基因数据框以及各种形式的id（例如OMIM，Ensembl，Genatlas）。我想获取与每个基因相关的所有SNP的列表。（这是this question的反向。）

到目前为止，我发现的最佳解决方案是使用biomaRt package (bioconductor)。有一个我需要进行here的查找类型的示例。适合我的目的，这是我的代码：

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
)

输出的数据帧如下所示：

  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行）

该代码可以工作，但是运行时间非常长。对于上述情况，大约需要45秒。我以为这可能与等位基因频率有关，服务器可能会即时计算这些等位基因频率。但是，仅查找SNP rs id的最低要求大约需要25秒。我有数千个基因，因此这将花费一整天的时间（假设没有超时或其他错误）。这是不对的。我的互联网连接速度并不慢（20-30兆位）。

我尝试在每个查询中查找更多基因。这确实起到了帮助作用。一次查找10个基因的速度大约是查找单个基因的速度的10倍。

获得与基因ID载体相关的SNP载体的最佳方法是什么？

如果我可以下载两个表，一个表包含基因和它们的位置，一个表包含SNP和它们的位置，那么我可以使用dplyr（或者也许是data.table）轻松解决此问题。我还没有找到这样的表。

最佳答案

由于您使用的是R，因此这里有一个使用包rentrez的想法。它利用NCBI的Entrez数据库系统，尤其是eutils函数elink。您将不得不为此编写一些代码，并可能需要调整参数，但这可能是一个不错的开始。

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'

我建议您手动仔细检查SNP的数量是否与您对感兴趣基因的期望相符-您可能需要进一步细化并按转录本进行限制等。

对于多个基因ID：

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'

结果用by_id=TRUE按基因分组

关于r - 通过基因ID获取SNP列表的最佳方法？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/41490657/