使用BiomaRt获得在线注释信息 17

完整的生物信息学分析步骤往往会包含注释工作。在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

这个页面是biomart提供的默认风格,布局分三个部分:主菜单,左侧导航条,右侧信息显示以及具体表单区。首先在页面左侧从上至下依次选择所需的数据源(dataset),过滤器(filters)以及数据组成(attributes)。
Filters

Attributes

之后就可以点击主菜单中的结果(Results)按钮来查看结果了。我们可以看到,在Attributes中选中的每一项都会以列名的形式显示出来。在这一页中我们可以选择格式后点击GO按钮下载。

Results

有了上面的介绍,我们就可以开始了解如何使用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来获取。所以我们的思路可以更加发散一些。

17 thoughts on “使用BiomaRt获得在线注释信息

  1. Reply ljy 10月 24,2014 9:03 上午

    很高兴看见你的博文,关于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值的例子,不知道这种情况该怎么处理?

    • Reply admin 10月 24,2014 9:14 上午

      我也很久不做microarray了。我试了一下你的例子,的确是不行。对于基因水平,你需要去affymetrix网站上下载它的注释文件,然后自己parse一下。

  2. Reply super 10月 27,2014 12:49 下午

    大神我有一个蚊子的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分析是不是就不客观了因为我的数据是蚊子的。
    谢谢!

    • Reply admin 10月 30,2014 1:57 下午

      有几种解决方案,一个是使用类似http://www.ncbi.nlm.nih.gov/nuccore/XM_001650645这种的链接抓取网页,然后parse到Vectorbase ID进行分析,这个很多语言都可以做。
      二个是象你说的,找同源蛋白,但是这种办法并不直接。

  3. Reply super 6月 10,2015 5:37 上午

    欧老师 我现在通过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注释也不一样 有没有什么统一的身份认证??

    • Reply super 6月 10,2015 8:55 上午

      补上条,是不是不要用external gene name 而是用 hgnc symbol?

    • Reply admin 6月 10,2015 2:07 下午

      你自己建一个orthologous genes的映射表,给每一行一个key值,比如
      key human pig horse zebrafish
      1 MX2 MX2 MX2 MXE
      然后你再把所有的差分基因用这个key值注释一下,最后以key值来merge就好了。

      • Reply super 6月 10,2015 4:18 下午

        谢谢!那么问题来了。
        (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

        • Reply admin 6月 12,2015 7:56 上午

          无论问题1,2,都是一个多对多映射的问题,我想解决它并不难,你可以自己动手试试。在R当中,函数merge可以很好的解决这个问题。如果想让表格简洁一些,可以压缩表格,也可以使用mysql中的多表。反正方法多得是。自己动动手做吧。

      • Reply super 6月 11,2015 5:18 上午

        欧老师,对于您说的“你自己建一个orthologous genes的映射表,给每一行一个key值,”我不知道如何自动去做。我现在的方法是
        step1基于现有的牛的Ensembl cow gene ID ,然后用bioMart或者bioDBnet将这个比对到Homolog Ensembl human gene ID,
        step2 再用这个Homolog Ensembl human gene ID去找它们的human gene name。
        不知您看可否?

        • Reply admin 6月 12,2015 8:04 上午

          有很多在方可以下载到orthologous gene的表格,你只需要把这几个物种的merge到一起就可以了。你可以去http://www.ncbi.nlm.nih.gov/homologene试试运气。

  4. Reply super 6月 29,2015 11:07 上午

    欧老师我有个问题
    我用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全?

  5. Reply kaji331 7月 6,2015 2:12 上午

    博主你好,如果bimaRt所有网站和数据库都没有这个物种,但是KEGG上有通路,NCBI有参考序列,比对之后怎么注释呢?比如Listeria monocytogenes单增李斯特菌

  6. Reply liuxyanbio 12月 1,2015 11:52 上午

    真高兴看到这篇文章,我是个菜鸟,想下载非编码RNA,或者是lncRNA该怎么选择参数呢,麻烦您抽空解答一下,万分感谢

    • Reply admin 12月 1,2015 12:11 下午

      试试

      library("biomaRt")
      
      mart < - useMart("ensembl", "hsapiens_gene_ensembl")
      
      attributes <- listAttributes(mart)
      
      attributes[grepl("gene type", attributes[,2], ignore.case=T), ]
      
      x <- getBM(c("gene_biotype", "ensembl_gene_id"), mart=mart)
      
      levels(factor(x$gene_biotype))
      
      lincRNA <- x[x$gene_biotype=="lincRNA",]
      

Leave a Reply

  

  

  

%d 博主赞过: