利用biomaRt包下载HGMD公开版的突变位点

前两面文章介绍了Ensembl的biomart,相信你对biomart应该有所了解了,在此再介绍一种方法,即通过R语言包biomaRt下载HGMD的数据。

HGMD的最新数据是需要购买授权才行,公开版信息不仅滞后,而且不能下载,不能得到基因组位置,在biostart上看到有人说Ensembl整合了HGMD的公开版,心想能获得公开版的数据也不错,于是采用biomaRt包下载。

各位不要高兴,最终的结果是,只得到了所有突变的基因组位置,未能下到具体的突变类型,以及与表型的关系。不过能下载基因组的位置,也算不错,结合对这些位置的注释,能获取不少信息。如果您对这些位置的利用有更多或者更好的想法,欢迎与我讨论。

1,安装biomaRt包

source("http://bioconductor.org/biocLite.R")
chooseCRANmirror()
chooseBioCmirror()
biocLite("biomaRt")

2,显示ensembl的biomart

library(biomaRt)
listEnsembl()
#如果要显示特定版本的,添加version参数
#listEnsembl(version=78)


3,选择snp mart

#HGMD的数据,不用想,一定是属于SNP,所以选择snp
snp = useEnsembl(biomart="snp")
#snp = useEnsembl(biomart="snp",version=78)

4,显示snp mart的数据集

下述命令会显示snp mart的所有数据集,很明显,我们需要的是hsapiens_snp的数据集,相当于选择哪个物种的数据。description是对数据集的描述,version为版本。

listDatasets(snp) 
snp = useEnsembl(biomart="snp",dataset="hsapiens_snp")

5,显示过滤参数

有时候,我们并不是将整个数据集下载下来,有可能只关注一个基因或者一个RSID,就可以根据显示的过滤参数加上我们的条件,这样得到的数据只是我们关心的。因为我们要获得HGMD的数据,后续要选择variation_source。

listFilters(snp)

6,显示属性

如果你看过前两篇文章,应该能明白属性其实对应的是我们想得到的关于某一对象的属性,一个属性对应一列,比如SNP的属性就可以包括chr,pos等。后续我选择了’chr_name’, ‘chrom_start’,’chrom_end’,’refsnp_id’,’associated_gene’,’allele’,’phenotype_description’。我相信,你也关心这些信息,不过最终结果可能没有那么美好,请往下看

listAttributes(snp)

7,通过getBM得到数据

attributes对应你想要获得的属性
filters对应你要过滤的项目,values对应filter的过滤参数
如果你有多个filter,values的参数对应的应该是个list

result <- getBM(attributes=c('chr_name', 'chrom_start','chrom_end','refsnp_id','associated_gene','allele','phenotype_description'), filters = 'variation_source', values ="HGMD-PUBLIC", mart = snp)
write.table(result,sep="t",file="./hgmd.txt",quote=FALSE,row.names=F)
head(result)

在这里只得到了突变的位置和基因,其他的未公开,想想也对,要是都能下载到了,HGMD的生意也会少了很多。不过HGMD官网没有提供基因组位置,这里能获得,再结合相关基因,这些数据还是很有用的。

Ensembl官网提供的另外两个例子,也可以参考下

Example query: Fetch all the Ensembl gene, transcript IDs, HGNC symbols and chromosome locations located on the human chromosome 1

> library(biomaRt)
> ensembl = useEnsembl(biomart="ensembl", dataset="hsapiens_gene_ensembl")
> chr1_genes  head(chr1_gene)


Example query: Fetch Ensembl Gene, Transcript IDs, HGNC IDs and symbols and Uniprot Swissprot accessions mapped to the human Ensembl Gene ID "ENSG00000139618"

> library(biomaRt)
> ensembl = useEnsembl(biomart="ensembl", dataset="hsapiens_gene_ensembl")
> hgnc_swissprot  hgnc_swissprot  

参考:
http://asia.ensembl.org/info/data/biomart/biomart_r_package.html
https://bioconductor.org/packages/release/bioc/html/biomaRt.html

#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
#####################################################################

利用biomaRt包下载HGMD公开版的突变位点》上有2条评论

  1. 小周

    > rm(list=ls())
    > library(“biomaRt”)
    > ensembl=useMart(“ensembl”)
    > listMarts()
    biomart version
    1 ENSEMBL_MART_ENSEMBL Ensembl Genes 89
    2 ENSEMBL_MART_MOUSE Mouse strains 89
    3 ENSEMBL_MART_SNP Ensembl Variation 89
    4 ENSEMBL_MART_FUNCGEN Ensembl Regulation 89
    > ensembl = useMart(“ensembl”,version = 78)
    Error in useMart(“ensembl”, version = 78) :
    Incorrect BioMart name, use the listMarts function to see which BioMart databases are available
    #listMarts()得到的列表其实useMart也没用,反而是直接用的ensembl,而不是函数useMart里面自己说的
    useMart(biomart, dataset, host=”www.ensembl.org”,
    path=”/biomart/martservice”, port=80, archive=FALSE, ssl.verifypeer =
    TRUE, version, verbose = FALSE)
    Arguments

    biomart
    BioMart database name you want to connect to. Possible database names can be retrieved with the functio listMarts
    不知道到底怎么才能正确使用version参数,因为想用同一个版本的

    回复
  2. 王浩

    你好,我想请教一下,我想下载HGMD的突变具体信息,可以您说的方法下载突变位置以后再用其他方法注释突变类型和其他详细信息么?谢谢您!

    回复

发表评论

电子邮件地址不会被公开。 必填项已用*标注

This site uses Akismet to reduce spam. Learn how your comment data is processed.