使用topGO进行GO富集分析

作者: 生信修炼手册 | 来源:发表于2018-11-14 17:10 被阅读92次

    欢迎关注”生信修炼手册”!

    topGO是一个专门用于做GO富集分析的R包,它默认从GO.db中读取GO的分类和结构信息,结合富集分析的结果,它可以画出如下所示的GO有向无环图

    除了GO富集结果可视化这一特点,topGO还提供了多种富集分析的统计方法,示意如下

    甚至支持自定义统计算法和模型,当然,常规情况下我们使用经典的费舍尔精确检验就可以了。

    topGO的核心是构建一个topGOdata类型的对象,需要以下3种元素

    1. 基因组所有的基因

    2. 基因组所有基因对应的GO注释

    3. GO之间的分类和结构信息

    基因的列表在R中用向量存储,topGO会根据names属性来读取基因的标识符,而这个向量的内容是一个只包括01的因子,0表示这个基因为非差异基因,1表示这个基因为差异基因。示例如下

    > head(geneList)
    1  2  3  9 10
    1  0  1  0  1

    第一行为names属性,这里为entrez id, 第二行为向量的内容。

    基因的GO注释,可以从Bioconductor提供的注释包中获得,比如human就可以从org.Hs.eg.db包中得到,如果没有现成的注释包,也可以从文件中读取。

    所有GO的分类和结构信息是从GO.db这个R包中读取的,所以使用时,不需要考虑这个因素,但是需要注意,GO.db这个包的更新是没有GO官网更新的快的,这也是用这种R包的缺点。

    用human基因构建topGOdata的代码如下

    sampleGOdata <- new(
    "topGOdata",
    ontology = "BP",
    allGenes = geneList,
    nodeSize = 10,
    annot = annFUN.org,
    mapping = "org.Hs.eg.db",
    ID = "entrez")

    ontology指定GO的分类,GO的3大类别BP, CC, MF之间是独立,所以GO其实分为3个子数据库,做富集分析时,不同类别分开做。

    nodeSize用于对GO进行过滤,GO最少需要包含10个基因才会对该GO进行分析。

    annot指定基因对应的GO注释如何读取,annFUN.org代表从Bioconductor提供的org.xx.xx.db包中读取,mapping指定org包的名字,ID指定基因标识符的类型。

    对于那些没有现成的注释包的物种,可以从文件中读取所有基因的GO注释信息,文件内容示例如下

    121005    GO:0005576
    155158    GO:0005488
    160828    GO:0005488
    105778    GO:0016021, GO:0016020

    \t分隔的两列,第一列为基因ID,第二列为对应的GO注释,多个注释之间用逗号连接。这种文件可以通过readMapping函数直接读取,从文件读取注释信息,并构建topGOdata的代码如下

    geneID2GO <- readMappings("geneid2go.map")
    sampleGOdata <- new(
    "topGOdata",
    ontology = "BP",
    allGenes = geneList,
    annot = annFUN.gene2GO,
    gene2GO = geneID2GO)

    构建好topGOdata类型的对象后,可以直接进行富集分析,代码如下

    # 运行富集分析
    result <- runTest(
    sampleGOdata, 
    algorithm = "classic", 
    statistic = "fisher")
    
    # 提取显著富集的top100个GO
    sig.tab <- GenTable(
      sampleGOdata,
      Fis = result,
      topNodes = 100)

    富集分析的结果如下

    > head(sig.tab)
           GO.ID                                        Term Annotated Significant Expected     Fis
    1 GO:0006107              oxaloacetate metabolic process        12           3     0.20 0.00095
    2 GO:0065008            regulation of biological quality      3551          82    60.27 0.00139
    3 GO:0006116                              NADH oxidation         4           2     0.07 0.00168
    4 GO:0030837 negative regulation of actin filament po...        51           5     0.87 0.00169
    5 GO:0048872              homeostasis of number of cells       233          11     3.95 0.00218
    6 GO:0007010                   cytoskeleton organization      1199          34    20.35 0.00226

    构建富集GO的有向无环图有两种方式,第一种

    showSigOfNodes(
      sampleGOdata,
      score(result),
      firstSigNodes = 5,
      useInfo = 'all')

    fisrtSigNodes代表显示显著富集的top5个GO的层次结构,useInfo表示的图片中的节点上需要标注的信息,all表示全部信息,示例如下

    GO编号是都会标注的,第二行为GO的描述信息,对应的值为def, 第三行为富集分析的p值,对应的值pval, 第四行为位于该GO下的基因总数与差异基因的比例,对应的值为counts, all代表全部都标记,none代表只标记GO编号,如果只想要标记其中某一项内容,设置为对应的值即可。

    另外一种用法会在工作目录直接生成pdf文件,代码如下

    printGraph(
      sampleGOdata,
      result,
      firstSigNodes = 5,
      fn.prefix = "tGO",
      useInfo = "pval",
      pdfSW = TRUE)

    fn.prefix代表生成的PDF文件名的前缀, 更多的用法请参考官方文档。

    ·end·

    —如果喜欢,快分享给你的朋友们吧—

    扫描关注微信号,更多精彩内容等着你!

    相关文章

      网友评论

        本文标题:使用topGO进行GO富集分析

        本文链接:https://www.haomeiwen.com/subject/sqwjfqtx.html