按基因 id 获取 SNP 列表的最佳方法？答案

【问题标题】：Best way to get list of SNPs by gene id?按基因 id 获取 SNP 列表的最佳方法？
【发布时间】：2017-05-20 07:52:39
【问题描述】：

我有一个很长的基因数据框和各种形式的 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 mbit）。

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

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

如果我可以只下载两张表，一张包含基因及其位置，另一张包含 SNP 及其位置，那么我可以使用 dplyr（或者可能是 data.table ）。我一直找不到这样的表。

【问题讨论】：

由于您的查询非常简单，长时间运行是 Biomart 数据库的限制，不太可能有更快的方法。也就是说，这非常慢，通常应该运行得更快。在我的电脑上，它运行不到 5 秒。但 Biomart 服务器以反复无常着称。
理想情况下，我只需在某处下载两个 CSV 格式的表格并进行自己的过滤。但我找不到任何这样的表。根据this 5 year old question，显然曾经有一些。
您仍然可以下载 dbSNP 批量数据，例如：ftp.ncbi.nih.gov/snp/organisms/human_9606_b146_GRCh38p2
是的，但是格式很糟糕。我花了几个小时试图打开这些文件，但没有运气。为什么他们不能以 csv 等普通格式提供数据？
@Deleet 这是行业标准格式。它有大量的解析器，包括用于 R/Bioconductor 的解析器，例如‹VariantAnnotation›。事实上，我相信 Bioconductor 可能也有预制的 dbSNP 数据库，尽管您必须自己进行调查。

标签： r bioinformatics bioconductor biomart

【解决方案1】：

由于您使用的是 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

【讨论】：