可能是个生物信息学数据超市吧

biomaRt这个包很久以前我就给它写过教程(点击阅读),但是排版不好,可读性很差,所以我用R Markdown重新来一个。 当然了,它本身有官方的英文版教程(点击阅读),我在翻译的基础上面,加入了自己的理解, 下面是正文:

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

包的安装

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

source("http://bioconductor.org/biocLite.R")
biocLite(“biomaRt”)
install.packages('DT')

选择数据库

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

library("biomaRt")
library(DT)
listMarts() #看看有多少数据库资源
##                biomart               version
## 1 ENSEMBL_MART_ENSEMBL      Ensembl Genes 91
## 2   ENSEMBL_MART_MOUSE      Mouse strains 91
## 3     ENSEMBL_MART_SNP  Ensembl Variation 91
## 4 ENSEMBL_MART_FUNCGEN Ensembl Regulation 91
ensembl=useMart("ENSEMBL_MART_ENSEMBL")
all_datasets <- listDatasets(ensembl) #看看选择的数据库里面有多少数据表,这个跟物种相关
library(DT)
datatable(all_datasets, options = list(
  searching = FALSE,
  pageLength = 5,
  lengthMenu = c(5, 10, 15, 20)
))
ensembl = useDataset("hsapiens_gene_ensembl",mart=ensembl)
# 这是一步法选择人类的ensembl数据库代码.
ensembl = useMart("ENSEMBL_MART_ENSEMBL",dataset="hsapiens_gene_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,]
library(DT)
datatable(filters, options = list(
  searching = FALSE,
  pageLength = 5,
  lengthMenu = c(5, 10, 15, 20)
))
attributes = listAttributes(ensembl)
#attributes[1:5,]
datatable(attributes, options = list(
  searching = FALSE,
  pageLength = 5,
  lengthMenu = c(5, 10, 15, 20)
))

getSequence函数大同小异,建议大家活用help函数去查看每个函数的帮助文档。

简单讲几个例子咯: 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
)

可以看到结果里面已经成功的把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
)

三.对给定的基因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)

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

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

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

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

六.根据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] 82 83 84 85 86 87
# [1] 81 82 83 84 85 86
tmp[grep("refseq",tmp[,1]),]
##                        name                 description         page
## 82              refseq_mrna              RefSeq mRNA ID feature_page
## 83    refseq_mrna_predicted    RefSeq mRNA predicted ID feature_page
## 84             refseq_ncrna             RefSeq ncRNA ID feature_page
## 85   refseq_ncrna_predicted   RefSeq ncRNA predicted ID feature_page
## 86           refseq_peptide           RefSeq peptide ID feature_page
## 87 refseq_peptide_predicted RefSeq peptide predicted ID feature_page

很明显,是没有 refseq_dna 的,只有 refseq_mrna

然后我用了新的代码

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
## 1    NM_000546 IPR002117
## 2    NM_000546 IPR008967
## 3    NM_000546 IPR010991
## 4    NM_000546 IPR011615
## 5    NM_000546 IPR012346
## 6    NM_000546 IPR013872
## 7    NM_005359 IPR001132
## 8    NM_005359 IPR003619
## 9    NM_005359 IPR008984
## 10   NM_005359 IPR013019
## 11   NM_005359 IPR013790
## 12   NM_005359 IPR017855
##                                      interpro_description
## 1                            p53 tumour suppressor family
## 2              p53-like transcription factor, DNA-binding
## 3                             p53, tetramerisation domain
## 4                                 p53, DNA-binding domain
## 5  p53/RUNT-type transcription factor, DNA-binding domain
## 6                              p53 transactivation domain
## 7                               SMAD domain, Dwarfin-type
## 8                            MAD homology 1, Dwarfin-type
## 9                                         SMAD/FHA domain
## 10                                      MAD homology, MH1
## 11                                                Dwarfin
## 12                                       SMAD domain-like

就成功的完成了转换

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

entrez=c("673","7157","837")
getSequence(id = entrez, type="entrezgene",seqType="coding_gene_flank",upstream=100, mart=ensembl)
##                                                                                      coding_gene_flank
## 1 CCTCCGCCTCCGCCTCCGCCTCCGCCTCCCCCAGCTCTCCGCCTCCCTTCCCCCTCCCCGCCCGACAGCGGCCGCTCGGGCCCCGGCTCTCGGTTATAAG
## 2 CACGTTTCCGCCCTTTGCAATAAGGAAATACATAGTTTACTTTCATTTTTGACTCTGAGGCTCTTTCCAACGCTGTAAAAAAGGACAGAGGCTGTTCCCT
## 3 TCCTTCTCTGCAGGCCCAGGTGACCCAGGGTTGGAAGTGTCTCATGCTGGATCCCCACTTTTCCTCTTGCAGCAGCCAGACTGCCTTCCGGGTCACTGCC
##   entrezgene
## 1        673
## 2        837
## 3       7157

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

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

utr5 = getSequence(chromosome=3, start=185514033, end=185535839,
  type="entrezgene",seqType="5utr", mart=ensembl)
utr5
##                                                                                                                                             5utr
## 1                                                                                                                           Sequence unavailable
## 2 AGTCCCTAGGGAACTTCCTGTTGTCACCACACCTCTGAGTCGTCTGAGCTCACTGTGAGCAAAATCCCACAGTGGAAACTCTTAAGCCTCTGCGAAGTAAATCATTCTTGTGAATGTGACACACGATCTCTCCAGTTTCCAT
## 3                                                        TGAGCAAAATCCCACAGTGGAAACTCTTAAGCCTCTGCGAAGTAAATCATTCTTGTGAATGTGACACACGATCTCTCCAGTTTCCAT
## 4                                                                                                        ATTCTTGTGAATGTGACACACGATCTCTCCAGTTTCCAT
##   entrezgene
## 1     200879
## 2     200879
## 3     200879
## 4     200879
#这样就查询到了 chromosome=3, start=185514033, end=185535839,这个定位里面的所有utr5信息。
# 返回的urt5这个对象是一个data.frame,可以用exportFASTA()函数导出成fasta文件。
protein = getSequence(id=c(100, 5728),type="entrezgene",
  seqType="peptide", mart=ensembl)
protein
##                                                                                                                                                                                                                                                                                                                                                                                                                peptide
## 1                                                                                                                                                                                                                                                                                                                                         MAQTPAFDKPKVELHVHLDGSIKPETILYYGRRRGIALPANTAEGLLNVIGMDKPLTLPDFLAKFDYYMPAIARL*
## 2                                                                                                                                                                                                                                                                                                                                                                                                 Sequence unavailable
## 3                                                                                                                                                                                                                                                           ALLFHKMMFETIPMFSGGTCNPQFVVCQLKVKIYSSNSGPTRREDKFMYFEFPQPLPVCGDIKVEFFHKQNKMLKKDKMFHFWVNTFFIPGPEETSEKVENGSLCDQEIDSICSIERADNDKEYLVLTLTKNDLDKANKDKANRYFSPNFKVS*
## 4                                                                 MAQTPAFDKPKVELHVHLDGSIKPETILYYGRRRGIALPANTAEGLLNVIGMDKPLTLPDFLAKFDYYMPAIAGCREAIKRIAYEFVEMKAKEGVVYVEVRYSPHLLANSKVEPIPWNQAEGDLTPDEVVALVGQGLQEGERDFGVKARSILCCMRHQPNWSPKVVELCKKYQQQTVVAIDLAGDETIPGSSLLPGHVQAYQAVDILKTERLGHGYHTLEDQALYNRLRQENMHFEICPWSSYLTGAWKPDTEHAVIRLKNDQANYSLNTDDPLIFKSTLDTDYQMTKRDMGFTEEEFKRLNINAAKSSFLPEDEKRELLDLLYKAYGMPPSASAGQNL*
## 5                                                                                                                                             MAQTPAFDKPKVELHVHLDGSIKPETILYYGRRRGIALPANTAEGLLNVIGMDKPLTLPDFLAKFDYYMPAIAGCREAIKRIAYEFVEMKAKEGVVYVEVRYSPHLLANSKVEPIPWNQAEGDLTPDEVVALVGQGLQEGERDFGVKARSILCCMRHQPNWSPKVVELCKKYQQQTVVAIDLAGDETIPGSSLLPGHVQAYQEAVKSGIHRTVHAGEVGSAEVVKEAVDILKTERLGHGYHTLEDQALYNRLRQENMHFEAQK*
## 6                                                                                                                                                                                                                                                                                                                                                                                                 Sequence unavailable
## 7                                         MAQTPAFDKPKVELHVHLDGSIKPETILYYGRRRGIALPANTAEGLLNVIGMDKPLTLPDFLAKFDYYMPAIAGCREAIKRIAYEFVEMKAKEGVVYVEVRYSPHLLANSKVEPIPWNQAEGDLTPDEVVALVGQGLQEGERDFGVKARSILCCMRHQPNWSPKVVELCKKYQQQTVVAIDLAGDETIPGSSLLPGHVQAYQEAVKSGIHRTVHAGEVGSAEVVKEAVDILKTERLGHGYHTLEDQALYNRLRQENMHFEICPWSSYLTGAWKPDTEHAVIRLKNDQANYSLNTDDPLIFKSTLDTDYQMTKRDMGFTEEEFKRLNINAAKSSFLPEDEKRELLDLLYKAYGMPPSASAGQNL*
## 8 MTAIIKEIVSRNKRRYQEDGFDLDLTYIYPNIIAMGFPAERLEGVYRNNIDDVVRFLDSKHKNHYKIYNLCAERHYDTAKFNCRVAQYPFEDHNPPQLELIKPFCEDLDQWLSEDDNHVAAIHCKAGKGRTGVMICAYLLHRGKFLKAQEALDFYGEVRTRDKKGVTIPSQRRYVYYYSYLLKNHLDYRPVALLFHKMMFETIPMFSGGTCNPQFVVCQLKVKIYSSNSGPTRREDKFMYFEFPQPLPVCGDIKVEFFHKQNKMLKKDKMFHFWVNTFFIPGPEETSEKVENGSLCDQEIDSICSIERADNDKEYLVLTLTKNDLDKANKDKANRYFSPNFKVKLYFTKTVEEPSNPEASSSTSVTPDVSDNEPDHYRYSDTTDSDPENEPFDEDQHTQITKV*
##   entrezgene
## 1        100
## 2        100
## 3       5728
## 4        100
## 5        100
## 6       5728
## 7        100
## 8       5728
#这样就查到了这两个基因转录本的所有蛋白序列,100这个基因有5条序列,5728这个基因有四条序列。

getSequence函数返回的对象是一个data.frame,可以用exportFASTA()函数导出成fasta文件。

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

当然还有一些数据库的小技巧,第一个是参数 archive = TRUE,设置只用能获取的数据库

然后是设置特定选取hg19对应的信息。

ensembl = useMart("ensembl_mart_46", dataset="hsapiens_gene_ensembl", archive = TRUE)
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("ENSEMBL_MART_ENSEMBL",dataset="hsapiens_gene_ensembl")
head(keytypes(mart), n=3) 
## [1] "affy_hc_g110"        "affy_hg_focus"       "affy_hg_u133_plus_2"
length(keytypes(mart))
## [1] 334
length(columns(mart))
## [1] 1918
head(columns(mart), n=3) 
## [1] "3_utr_end"   "3_utr_end"   "3_utr_start"
k = keys(mart, keytype="chromosome_name") 
head(k, n=3) 
## [1] "1" "2" "3"
k = keys(mart, keytype="chromosome_name", pattern="LRG") 
head(k, n=3) 
## character(0)
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')
##   affy_hg_u133_plus_2 entrezgene
## 1           202763_at        836
## 2         209310_s_at        837
## 3           207500_at        838

而且这个包更新也比较频繁,所以大家如果看到比较老旧的教程,可能无法模仿,或者多年以后,你看到我这个教程,也会发现,重复不出来咯。

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2018-02-24

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏小樱的经验随笔

浅析Numpy.genfromtxt及File I/O讲解

Python 并没有提供数组功能,虽然列表 (list) 可以完成基本的数组功能,但它并不是真正的数组,而且在数据量较大时,使用列表的速度就会慢的让人难受。为此...

2894
来自专栏Vamei实验室

Django ORM模型:想说爱你不容易

作者:Vamei 出处:http://www.cnblogs.com/vamei 严禁转载。

1232
来自专栏灯塔大数据

每周学点大数据 | No.67 Hadoop 实践案例——记录去重

No.67 Hadoop 实践案例——记录去重 Mr. 王:现在我们看一个和 WordCount 很相似,在实际中应用也很多的例子——记录去重。 小可 :嗯,...

3198
来自专栏企鹅号快讯

LeetCode测试数据的爬虫

LeetCode的(包括付费)题目到处都有,可是测试数据怎么找呢?我设想了一种方法,来获得每道题的测试数据。 首先,对于权限不严格的在线评测系统,比如以前常做的...

3078
来自专栏王亚昌的专栏

合理使用const,慎用自运算

    项目最的出了几次运营事故,都是因为使用自乘、自加、自減运算,错改了非局部变量,导致将用户数据写溢出,最终只能进行回档处理。先给大家展示一下,漏出bug的...

711
来自专栏开源优测

编程入门的姿势-5月8日微信群语音分享

开头语 5月8日在微信群,语音分享了如何如何学习编程语言、并以python为例进行了分享相关经验,下面整理成文章共享给大家。 神马?还有微信群? 加入微信群正确...

3357
来自专栏大数据挖掘DT机器学习

Python]新手写爬虫全过程(已完成)

今天早上起来,第一件事情就是理一理今天该做的事情,瞬间get到任务,写一个只用python字符串内建函数的爬虫,定义为v1.0,开发中的版本号定义为v0.x。数...

3639
来自专栏C语言及其他语言

训练场题库中判题结果的详细解释

对于判题结果仅仅是大致的解释,仍不少同学感到迷惑,那今天我们就对这些结果一一详细解释并举例说明,让大家彻底觉悟! 等待 等待服务器正忙,请稍后查看运行并评判您的...

3045
来自专栏黑泽君的专栏

java多线程、集合和IO面试题_02

801
来自专栏Golang语言社区

Go语言是彻底的面向组合的并发语言

面向组合编程从AOP的Mixin,然后到Ruby的Traits,直至DCI设计,包括Scala的trait的组合设计,这些都有一个共同特点,组合特性是显式的,也...

3186

扫码关注云+社区