完整的生物信息学分析步骤往往会包含注释工作。在Bioconductor中,最方便的办法是使用注释包。注释资源除了以包的形式进行封装外,还可以通过诸如BiomaRt等工具获取在线的注释数据。使用在线资源为我们提供了更加及时以及丰富的注释资源。那么,什么是BiomaRt呢?如何理解和使用BiomaRt呢?
为了更好的理解和掌握biomaRt,我们可以先通过在线资源来了解一下它的原型biomart (http://www.biomart.org)。 biomart是为生物科研提供数据服务的免费软件,它为数据下载提供打包方案。它有许多成功的应用实例,比如欧洲生物信息学中心(The European Bioinformatics Institute ,EBI)维护的Ensembl数据库(http://www.ensembl.org/)就使用biomart提供数据批量下载服务, 还有COSMIC, Uniprot, HGNC, Gramene, Wormbase以及dbSNP等。
我们首先点击Ensembl主页上导航菜单中的BioMart链接可以进入下图所示的页面。我们可以通过页面下方的优酷链接查看视频教程。
这个页面是biomart提供的默认风格,布局分三个部分:主菜单,左侧导航条,右侧信息显示以及具体表单区。首先在页面左侧从上至下依次选择所需的数据源(dataset),过滤器(filters)以及数据组成(attributes)。
之后就可以点击主菜单中的结果(Results)按钮来查看结果了。我们可以看到,在Attributes中选中的每一项都会以列名的形式显示出来。在这一页中我们可以选择格式后点击GO按钮下载。
有了上面的介绍,我们就可以开始了解如何使用biomaRt软件包了。我们的任务是使用biomaRt实现基因名与Entrez Id及Ensemble ID之间的注释。来看代码:
> biocLite("biomaRt") #使用bioconnductor的biocLite安装biomaRt包 > library("biomaRt") #载入biomaRt包 > mart <- useMart("ensembl", "hsapiens_gene_ensembl") > entrez <- c("673","7157","837") > getBM(attributes=c("entrezgene","hgnc_symbol", "ensembl_gene_id", "affy_hg_u133_plus_2"), + filters = "entrezgene", + values = entrez, + mart = mart) entrezgene hgnc_symbol ensembl_gene_id affy_hg_u133_plus_2 1 673 BRAF ENSG00000157764 206044_s_at 2 673 BRAF ENSG00000157764 236402_at 3 673 BRAF ENSG00000157764 243829_at 4 7157 TP53 ENSG00000141510 5 7157 TP53 ENSG00000141510 211300_s_at 6 7157 TP53 ENSG00000141510 201746_at 7 837 CASP4 ENSG00000196954 209310_s_at 8 837 CASP4 ENSG00000196954 9 837 CASP4 ENSG00000196954 213596_at |
从上面的操作来看,使用biomaRt只需要两步,1,指定mart数据库,2,使用getBM获得注释。但是首先,我们如何知道有哪些服务器,以及这些服务器上哪些数据库呢?其次,我们如何获阳getBM中attributes,filters的正确设置呢?
关于第一个问题,我们可以使用biomaRt中的listMarts以及listDatasets两个函数来解决。
> marts <- listMarts(); head(marts) #查看当前可用的数据源 biomart version 1 ensembl ENSEMBL GENES 72 (SANGER UK) 2 snp ENSEMBL VARIATION 72 (SANGER UK) 3 functional_genomics ENSEMBL REGULATION 72 (SANGER UK) 4 vega VEGA 52 (SANGER UK) 5 fungi_mart_18 ENSEMBL FUNGI 18 (EBI UK) 6 fungi_variations_18 ENSEMBL FUNGI VARIATION 18 (EBI UK) > ensembl <- useMart("ensembl") #使用ensembl数据源 > datasets <- listDatasets(ensembl); datasets[1:10,] #查看ensembl中可用数据库 dataset description version 1 oanatinus_gene_ensembl Ornithorhynchus anatinus genes (OANA5) OANA5 2 tguttata_gene_ensembl Taeniopygia guttata genes (taeGut3.2.4) taeGut3.2.4 3 cporcellus_gene_ensembl Cavia porcellus genes (cavPor3) cavPor3 4 gaculeatus_gene_ensembl Gasterosteus aculeatus genes (BROADS1) BROADS1 5 lafricana_gene_ensembl Loxodonta africana genes (loxAfr3) loxAfr3 6 itridecemlineatus_gene_ensembl Ictidomys tridecemlineatus genes (spetri2) spetri2 7 mlucifugus_gene_ensembl Myotis lucifugus genes (myoLuc2) myoLuc2 8 hsapiens_gene_ensembl Homo sapiens genes (GRCh37.p11) GRCh37.p11 9 choffmanni_gene_ensembl Choloepus hoffmanni genes (choHof1) choHof1 10 csavignyi_gene_ensembl Ciona savignyi genes (CSAV2.0) CSAV2.0 |
对于第二个问题,我们使用biomaRt中的listFilters以及listAttributes两个函数来解决。
> mart <- useMart("ensembl", "hsapiens_gene_ensembl") > filters <- listFilters(mart); filters[grepl("entrez", filters[,1]),] name description 38 with_entrezgene with EntrezGene ID(s) 122 entrezgene EntrezGene ID(s) [e.g. 100287163] > attributes <- listAttributes(mart); attributes[grepl("^ensembl|hgnc", attributes[,1]), ] name description 1 ensembl_gene_id Ensembl Gene ID 2 ensembl_transcript_id Ensembl Transcript ID 3 ensembl_peptide_id Ensembl Protein ID 4 ensembl_exon_id Ensembl Exon ID 51 hgnc_id HGNC ID(s) 52 hgnc_symbol HGNC symbol 53 hgnc_transcript_name HGNC transcript name 134 ensembl_gene_id Ensembl Gene ID 135 ensembl_transcript_id Ensembl Transcript ID 136 ensembl_peptide_id Ensembl Protein ID 162 ensembl_exon_id Ensembl Exon ID 165 ensembl_gene_id Ensembl Gene ID 166 ensembl_transcript_id Ensembl Transcript ID 167 ensembl_peptide_id Ensembl Protein ID 175 ensembl_gene_id Ensembl Gene ID 176 ensembl_transcript_id Ensembl Transcript ID 177 ensembl_peptide_id Ensembl Protein ID 1616 ensembl_gene_id Ensembl Gene ID 1617 ensembl_transcript_id Ensembl Transcript ID 1618 ensembl_peptide_id Ensembl Protein ID 1691 ensembl_gene_id Ensembl Gene ID 1706 ensembl_transcript_id Ensembl Transcript ID 1707 ensembl_peptide_id Ensembl Protein ID 1715 ensembl_exon_id Ensembl Exon ID |
最后的问题是,biomaRt会被如何使用呢?我们做注释的时候,怎么就想到要使用biomaRt呢?因为在注释上,各种ID,symbol, name之间的转换都可以考虑使用biomaRt来做。更重要的是,biomaRt还会有很多SNP, alternative splicing, exon, intron, 5’utr, 3’utr等等信息。当然,只要能做也数据库并使用SQL访问的数据都可以使用biomaRt来获取。所以我们的思路可以更加发散一些。
很高兴看见你的博文,关于biomaRt我有个小问题。当我用biomaRt注释HuEx-1_0-st-v2的exon-level和gene-level的probesetid时发现,biomaRt似乎只能注释exon-level的probesetid,而对于gene-level的probesetid(transcription_cluster_ids)时却失败了,是因为biomaRt里没有它的数据吗?还是我操作失误?
mart <- useMart("ensembl", "hsapiens_gene_ensembl")
getBM(attributes=c("hgnc_symbol"),filters = "affy_huex_1_0_st_v2",values = "2319423", mart = mart)
只是举例一个value值的例子,不知道这种情况该怎么处理?
我也很久不做microarray了。我试了一下你的例子,的确是不行。对于基因水平,你需要去affymetrix网站上下载它的注释文件,然后自己parse一下。
大神我有一个蚊子的contigs的fasta文件,blast后是这个结果 都是蚊子的库
query_name query_length accession_number
1059N_ae_contig_1 676 XM_001650645
1059N_ae_contig_5 563 XM_001650988
1059N_ae_contig_6 2123 XM_001664234
我现在想对这些基因做pathway分析 但是显然这个genbank的accession id不给力,需要转成标准的名字 比如蚊子应该是Vectorbase ID。
请问我该怎么做?
我肯定不想一个个谷歌查手动改毕竟有万把个contigs。
此外我可以blast到人的蛋白质组上然后给每个基因附以基因名么?这样的话再去做pathway分析是不是就不客观了因为我的数据是蚊子的。
谢谢!
有几种解决方案,一个是使用类似http://www.ncbi.nlm.nih.gov/nuccore/XM_001650645这种的链接抓取网页,然后parse到Vectorbase ID进行分析,这个很多语言都可以做。
二个是象你说的,找同源蛋白,但是这种办法并不直接。
欧老师 我现在通过RNA-Seq的差分分析找到了一些基因 比如MX2 这个基因在人,牛,马,猪里都叫MX2 但是在斑马鱼叫MXE (我查ensembl的orthologues显示的)
现在我已经分别基于四种动物的基因组做了差分基因分析 得到了四个差分基因表(edgeR或者cuffdiff的输出表你们懂的) 想join成一个总表
范例:
基因 log2FC人 log2FC猪 log2FC马 log2FC鱼
MX2 1.5 2.6 1.6 -1.9
请问怎么做?
我现在的做法是 我先用Ensembl的Biomart把ensembl gene ID 转换为gene symbol 然后统一基于gene symbol去join或者grep
但我突然发现在鱼里面这个基因symbol叫MXE
但是如果用ensembl ID也不行 不同种类的ensembl注释也不一样 有没有什么统一的身份认证??
补上条,是不是不要用external gene name 而是用 hgnc symbol?
你自己建一个orthologous genes的映射表,给每一行一个key值,比如
key human pig horse zebrafish
1 MX2 MX2 MX2 MXE
然后你再把所有的差分基因用这个key值注释一下,最后以key值来merge就好了。
谢谢!那么问题来了。
(1)怎么建这个映射表?手动?我怎么知道哪些基因是orthologous 名称不同的?
比如
apol3 在牛只叫apol3,但是它的orthologous human gene有apol1,apol2,apol3,apol4。
(2)如果我要找人和牛共有的差分基因,而恰巧人的apol1和apol4也是差分基因。
那是不是sharing DE genes? 就是
人 牛 log2FC人 log2FC牛
apol1 apol3 1.2 3.2
apol3 apol3 3.4 3.2
apol4 apol3 1.4 3.2
无论问题1,2,都是一个多对多映射的问题,我想解决它并不难,你可以自己动手试试。在R当中,函数merge可以很好的解决这个问题。如果想让表格简洁一些,可以压缩表格,也可以使用mysql中的多表。反正方法多得是。自己动动手做吧。
欧老师,对于您说的“你自己建一个orthologous genes的映射表,给每一行一个key值,”我不知道如何自动去做。我现在的方法是
step1基于现有的牛的Ensembl cow gene ID ,然后用bioMart或者bioDBnet将这个比对到Homolog Ensembl human gene ID,
step2 再用这个Homolog Ensembl human gene ID去找它们的human gene name。
不知您看可否?
有很多在方可以下载到orthologous gene的表格,你只需要把这几个物种的merge到一起就可以了。你可以去http://www.ncbi.nlm.nih.gov/homologene试试运气。
欧老师我有个问题
我用GUI(http://www.biomart.org/biomart/martview)把牛基因 转化为Human Orthologs Ensembl ID,有一个基因ENSBTAG00000008021 转化后是 ENSG00000130303。
但是,如果我用biomaRt的R工具包
mart <- useMart("ensembl", "btaurus_gene_ensembl")
aa<-getBM(attributes=c("ensembl_gene_id","external_gene_name","hsapiens_homolog_ensembl_gene"),filters ="ensembl_gene_id",values = EnID[,1], mart = mart)
结果ENSBTAG00000008021 对应的Human Orthologs Ensembl ID 显示却是空格。。。为啥?难道GUI的数据库比R全?
不应该啊。会不会是biomart的指向数据库不同步的原因啊?
博主你好,如果bimaRt所有网站和数据库都没有这个物种,但是KEGG上有通路,NCBI有参考序列,比对之后怎么注释呢?比如Listeria monocytogenes单增李斯特菌
一般的思路就是找同源注释。你可以查查相关的homolog的数据库。
真高兴看到这篇文章,我是个菜鸟,想下载非编码RNA,或者是lncRNA该怎么选择参数呢,麻烦您抽空解答一下,万分感谢
试试