Bioconductor系列之biomaRt

Posted on by

biomaRt是一个超级网络资源库，里面的信息非常之多，就是网页版的biomaRt的R语言接口。谷歌搜索 the biomart user’s guide filetype:pdf 这个关键词，就看到关于这个包的详细介绍以及例子，我这里简单总结一下它的用法。

这个包一直在更新，函数总是变化，而且需要联网才能使用，基本上等于废物一个，希望大家不要使用它~

包的安装

Bioconductor系列包的安装方法都一样

source("http://bioconductor.org/biocLite.R")
biocLite(“biomaRt”)

选择数据库

然后我们就可以加载这个包来看看它的一些性质啦，大家也可以根据这个包的说明书里面的代码一行行代码敲进去看看效果，这样学习最快也最容易记住。

library("biomaRt")
listMarts() #看看有多少数据库资源
ensembl=useMart("ensembl")
listDatasets(ensembl) #看看选择的数据库里面有多少数据表，这个跟物种相关
ensembl = useDataset("hsapiens_gene_ensembl",mart=ensembl)
ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl")
# 这是一步法选择人类的ensembl数据库代码。

三个主要函数getBM,getSequence,getLDS

我们选择好了数据库就要开始干活啦，这个数据库的检索主要是三个函数getBM,getSequence,getLDS, 其中getBM这个函数可以部分用select语句替代。

首先我们讲讲getBM函数，它就四个参数。

  • filter来控制根据什么东西来过滤，可是不同数据库的ID，也可以是染色体定位系统坐标
  • Attributes来控制我们想获得什么，一般是不同数据库的ID
  • Values是我们用来检索的关键词向量
  • Mart是我们前面选择好的数据库，一般都是ensembl = useMart(“ensembl”,dataset=”hsapiens_gene_ensembl”)

getBM函数唯一的用处，就是做各种ID转换。

首先我们可以查看filter和attribute有哪些东西。

filters = listFilters(ensembl)
filters[1:5,]
attributes = listAttributes(ensembl)
attributes[1:5,]

然后我们简单讲几个例子咯：
Ps：这些都是在线注释，所以都是要网络的，网速慢的会非常坑

几个实用的例子

一．对几个芯片探针的ID号，注释它所捕获的基因的entrezID

ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl")
affyids=c("202763_at","209310_s_at","207500_at")
getBM(attributes=c('affy_hg_u133_plus_2', 'entrezgene'), filters = 'affy_hg_u133_plus_2', values = affyids, mart = ensembl)

结果如下：

affy_hg_u133_plus_2 entrezgene
1 202763_at 836
2 207500_at 838
3 209310_s_at 837

成功的把affymetrix的芯片探针ID，转为了对应的基因的entrez ID

二．对刚才的那三个探针ID号进行多个内容注释，每个探针都对应着基因名已经染色体及起始终止坐标。

affyids=c("202763_at","209310_s_at","207500_at")
getBM(attributes=c('affy_hg_u133_plus_2', 'hgnc_symbol', 
'chromosome_name','start_position','end_position', 'band'),
filters = 'affy_hg_u133_plus_2', values = affyids, mart = ensembl
)

affy_hg_u133_plus_2 hgnc_symbol chromosome_name start_position end_position band
1 202763_at CASP3 4 184627696 184649509 q35.1
2 207500_at CASP5 11 104994235 105023168 q22.3
3 209310_s_at CASP4 11 104942866 104969436 q22.3

三．对给定的基因ID号进行GO注释

library("biomaRt")
ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl")
entrez=c("673","837")
goids = getBM(attributes=c('entrezgene','go_id'), filters='entrezgene', values=entrez, mart=ensembl)
head(goids)

entrezgene go_id
1 673 GO:0004674
2 673 GO:0005524
3 673 GO:0006468
4 673 GO:0006464
5 673 GO:0008150
6 673 GO:0003674

四．通过染色体及起始终止坐标来挑选基因

getBM(c('affy_hg_u133_plus_2','ensembl_gene_id'), filters = c('chromosome_name','start','end'),
 values=list(16,1100000,1250000), mart=ensembl)

affy_hg_u133_plus_2 ensembl_gene_id
1 ENSG00000260702
2 215502_at ENSG00000260532
3 ENSG00000273551
4 205845_at ENSG00000196557
5 ENSG00000196557
6 ENSG00000260403
7 ENSG00000259910
8 ENSG00000261294
9 220339_s_at ENSG00000116176
10 ENSG00000277010
11 217023_x_at ENSG00000197253
12 210084_x_at ENSG00000197253
13 215382_x_at ENSG00000197253
14 216474_x_at ENSG00000197253
15 207134_x_at ENSG00000197253
16 205683_x_at ENSG00000197253
17 217023_x_at ENSG00000172236
18 210084_x_at ENSG00000172236
19 215382_x_at ENSG00000172236
20 207741_x_at ENSG00000172236
21 216474_x_at ENSG00000172236
22 207134_x_at ENSG00000172236
23 205683_x_at ENSG00000172236

五．对特定的GO ID号来查询该go通路上面的基因是哪些。

getBM(c('entrezgene','hgnc_symbol'), filters='go_id', values='GO:0004707', mart=ensembl)

entrezgene hgnc_symbol
1 5596 MAPK4
2 225689 MAPK15
3 5601 MAPK9
4 6300 MAPK12
5 5594 MAPK1
6 51701 NLK
7 5597 MAPK6
8 5599 MAPK8
9 1432 MAPK14
10 5603 MAPK13
11 5595 MAPK3
12 5600 MAPK11
13 5602 MAPK10
14 5598 MAPK7

六．根据refseq数据库的NM系列ID号来获取信息

refseqids = c("NM_005359","NM_000546")
ipro = getBM(attributes=c("refseq_dna","interpro","interpro_description"), filters="refseq_dna",values=refseqids, mart=ensembl)

这个例子的代码有错误，因为refseq的信息没有refseq_dna
Error in getBM(attributes = c(“refseq_dna”, “interpro”, “interpro_description”), :
Invalid attribute(s): refseq_dna
Please use the function ‘listAttributes’ to get valid attribute names
我简单检查了一下，发现需要更正为refseq_mrna

tmp=listAttributes(ensembl)
grep("refseq",tmp[,1])
# [1] 81 82 83 84 85 86
tmp[grep("refseq",tmp[,1]),]

name description
81 refseq_mrna RefSeq mRNA [e.g. NM_001195597]
82 refseq_mrna_predicted RefSeq mRNA predicted [e.g. XM_001125684]
83 refseq_ncrna RefSeq ncRNA [e.g. NR_002834]
84 refseq_ncrna_predicted RefSeq ncRNA predicted [e.g. XR_108264]
85 refseq_peptide RefSeq Protein ID [e.g. NP_001005353]
86 refseq_peptide_predicted RefSeq Predicted Protein ID [e.g. XP_001720922]

然后我用了新的代码

refseqids = c("NM_005359","NM_000546")
ipro = getBM(attributes=c("refseq_mrna","interpro","interpro_description"), filters="refseq_mrna",values=refseqids, mart=ensembl)
ipro

refseq_mrna interpro interpro_description
1 NM_000546 IPR011615 p53, DNA-binding domain
2 NM_000546 IPR010991 p53, tetramerisation domain
3 NM_000546 IPR013872 p53 transactivation domain
4 NM_000546 IPR002117 p53 tumour suppressor family
5 NM_000546 IPR008967 p53-like transcription factor, DNA-binding
6 NM_005359 IPR001132 SMAD domain, Dwarfin-type
7 NM_005359 IPR003619 MAD homology 1, Dwarfin-type
8 NM_005359 IPR013019 MAD homology, MH1
9 NM_005359 IPR008984 SMAD/FHA domain

七．根据基因的entrez ID号来挑选该基因的指定上下游区域信息或者蛋白序列

entrez=c("673","7157","837")
getSequence(id = entrez, type="entrezgene",seqType="coding_gene_flank",upstream=100, mart=ensembl)

这样就得到了三条序列，是给定基因的上游100bp序列 coding_gene_flank entrezgene

1 CCTCCGCCTCCGCCTCCGCCTCCGCCTCCCCCAGCTCTCCGCCTCCCTTCCCCCTCCCCGCCCGACAGCGGCCGCTCGGGCCCCGGCTCTCGGTTATAAG 673
2 TCCTTCTCTGCAGGCCCAGGTGACCCAGGGTTGGAAGTGTCTCATGCTGGATCCCCACTTTTCCTCTTGCAGCAGCCAGACTGCCTTCCGGGTCACTGCC 7157
3 CACGTTTCCGCCCTTTGCAATAAGGAAATACATAGTTTACTTTCATTTTTGACTCTGAGGCTCTTTCCAACGCTGTAAAAAAGGACAGAGGCTGTTCCCT 837

其实这个getSequence函数还有非常多的用法，当然主要的变化在其readme上面可以看到，主要是seqType可以有多个选择。

utr5 = getSequence(chromosome=3, start=185514033, end=185535839,
 type="entrezgene",seqType="5utr", mart=ensembl)
utr5
#这样就查询到了 chromosome=3, start=185514033, end=185535839,这个定位里面的所有utr5信息。
protein = getSequence(id=c(100, 5728),type="entrezgene",
 seqType="peptide", mart=ensembl)
protein
#这样就查到了这两个基因转录本的所有蛋白序列，100这个基因有5条序列，5728这个基因有四条序列。

八，选择其它数据库来进行查询，比如snp数据库

当然还有一些数据库的小技巧，第一个是参数 archive = TRUE，设置只用能获取的数据库
ensembl = useMart(“ensembl_mart_46”, dataset=”hsapiens_gene_ensembl”, archive = TRUE)
然后是设置特定选取hg19对应的信息。

listMarts(host='may2009.archive.ensembl.org')
ensembl54=useMart(host='may2009.archive.ensembl.org', biomart='ENSEMBL_MART_ENSEMBL')
ensembl54=useMart(host='may2009.archive.ensembl.org', biomart='ENSEMBL_MART_ENSEMBL', dataset='hsapiens_gene_ensembl')

还可以选取其它物种，秀丽隐杆线虫

wormbase=useMart("WS220",dataset="wormbase_gene")
listFilters(wormbase)
listAttributes(wormbase)
getBM(attributes = c("public_name","rnai","rnai_phenotype_phenotype_label"),
filters="gene_name", values=c("unc-26","his-33"),
mart=wormbase)

写在最后

最后我简单提一下select函数是如何部分替代getBM函数的，因为biomaRt是在线数据库，本来只能用它自己的getBM系列函数，但是为了对接其它bioconductor系列包，也可以用select函数来操作这个在线数据库。

mart<-useMart(dataset="hsapiens_gene_ensembl",biomart='ensembl') head(keytypes(mart), n=3) head(columns(mart), n=3) k = keys(mart, keytype="chromosome_name") head(k, n=3) k = keys(mart, keytype="chromosome_name", pattern="LRG") head(k, n=3) affy=c("202763_at","209310_s_at","207500_at") select(mart, keys=affy, columns=c('affy_hg_u133_plus_2','entrezgene'), keytype='affy_hg_u133_plus_2')
​

 

Comments are closed.