KeepNotes blog

Stay hungry, Stay Foolish.

0%

GO.db/topGO/goseq R包

Posted on Edited on In ,
Symbols count in article: 3.1k Reading time ≈ 3 mins.

GO.db

GO.db是一个用注释maps来描述Gene Ontology的一个R包，其在很多GO注释及富集的R包中被调用，应用广泛，每半年更新一次。主要用来描述各个Term之间的父子节点联系以及Term的信息。 至于使用的话，可以用names(GOMAPCOUNTS)来查看主要有几个函数，以BP为例：

  • 所有父节点信息 GOBPANCESTOR

  • 上一级父节点信息 GOBPPARENTS

  • 所有子节点信息 GOBPOFFSPRING

  • 下一级子节点信息 GOBPCHILDREN

  • 每个Term所有描述信息 GOTERM

使用方式个人总结下有以下几种：

  • 转化为list

          ancestor <- as.list(GOBPANCESTOR)

  • 转化为data.frame

          ancestor <- toTable(GOBPANCESTOR)

  • 直接查看某个GO id的信息

          ancestor <- GOBPANCESTOR$"GO:0001580"

用处：一些R包会调用GO.db获得某个GO id的描述信息，还有些会用GO.db（例如topGO）的每个GO之间的关系来绘制 directed acyclic graph (DAG 有向无环图)，而我使用这个R包是想将几个子节点的父节点也参与GO富集（通常的GO富集R包都不会有这个功能）。

topGO

topGO这个R包可以用来做GO富集，而我接触到这个R包则是因为它能将富集结果以有向无环图形式展示，这一点是其他大部分做GO富集的R包所有没有的，就算是有这个功能的R包一般也调用了topGO。

但从富集这个功能来说，topGO比其他的R包没有较大的优势，以及除了拥有不止一种的统计学方法以及算法来计算富集的P值，一般的步骤如下(一般用于有模式物种或者的有bioconductor注释包的物种):

sampleGOdata <- new("topGOdata",
    description = "Simple session", ontology = "BP",
    allGenes = geneList, geneSel = topDiffGenes,
    nodeSize = 10,
    annot = annFUN.db, affyLib = affyLib)

annFUN.db表示从bioconductor注释包中提取注释信息的函数（其实还有其他函数可以使用），affyLib则是对应注释包

resultFisher <- runTest(sampleGOdata, algorithm = "classic", statistic = "fisher")

然后就是使用经典的fisher检验进行富集分析

当然topGO这个R包并不是只能做有注释R包的物种，它提供了一种自定义的注释，其实就是输入你自己的GO注释，然后构建topGOdata对象，接着就跟上述模式物种富集方法类似了，从而可以使用其他功能。我们需要准备一个文件go_bp.txt（以BP为例），第一列是gene id，第二列则是gene对应的GO id，每个GO id以逗号分隔，然后使用readMappings读入，然后再输入差异基因列表geneid.txt

geneid2GO <- readMappings(file = "go_bp.txt")
genenames <- names(geneid2GO)
gene_data <- read.table(file = "geneid.txt")
gene_id <- gene_data[,1]
genelist <- factor(as.integer(genenames %in% gene_id))
names(genelist) <- genenames

GOdata <- new("topGOdata", allGenes = genelist, ontology = "BP", 
          annot = annFUN.gene2GO, gene2GO = geneid2GO)

这个也算一个不错的功能，但是这个功能有一些GO富集的R包也有。。。。。。

最后则是算是这个R包特有的将GO富集结果以DAG（有向无环图）展示

showSigOfNodes(GOdata, pvalue, firstSigNodes = 10, useInfo = 'all')
topGO.png

goseq

goseq也是一个做GO富集的R包（也做KEGG富集，但是KEGG的数据包已经在R里面好久没更新了，因为KEGG数据库收费了）。

goseq相比其他富集R包，有个特别的地方：其富集算法使用了Wallenius non-central hyper-geometric distribution，相对于普通的 Hyper-geometric distribution，此分布的特点是从某个类别中抽取个体的概率与从某个类别之外抽取一个个体的概率是不同的，这种概率的不同是通过对基因长度的偏好性进行估计得到的，从而其认为能更为准确地计算出 GOterm被差异基因富集的概率。

因此在做富集前，需要对gene长度进行计算，如果有参考基因组的话，直接选择这个R包所提供的基因组注释信息来提取gene长度

pwf=nullp(genes,"hg19","ensGene")

然后进行富集

GO.wall=goseq(pwf,"hg19","ensGene")

当然goseq也能像topGO一样接受自定义的注释信息，唯一的不同都是需要自行统计gene的长度，然后再导入，比如差异基因gene_id，自行提供的基因对应GO信息gene2GO

pwf=nullp(genes,bias.data=gene_lengths)

gene_lengths是每个基因的长度

然后进行自定义GO注释的富集分析

go <- goseq(pwf = pwf, id = pro_id, gene2cat = uniprotid2GO, method = "Hypergeometric", use_genes_without_cat = TRUE)

Summary

上述的3个R包，各有特点吧，如果没有特别要求的话，后两者都能满足GO富集的需求。GO.db包是很多GO分析包的基础，了解下也不错。

最近看了下clusterProfiler包代码，发现其也是调用了GO.db包进行一些GO信息的提取，并也调用了topGO包来展现富集分析的结果，而且也支持自定义的GO注释信息的输入来进行GO富集分析，所以说也蛮全面的，而且是使用起来也比较友好。

其实最好的说明还是查看三者的说明文档，里面都有详细的教程，看一下一般就懂了

http://www.bioconductor.org/packages/release/data/annotation/html/GO.db.html
http://www.bioconductor.org/packages/release/bioc/html/topGO.html
http://www.bioconductor.org/packages/release/bioc/html/goseq.html

参考文章：http://www.biotrainee.com/thread-558-1-1.html

本文出自于http://www.bioinfo-scrounger.com转载请注明出处

Related Posts