topGO 使用心得 转

genechip

topGO是Bioconductor的一个package，用于处理芯片数据。功能是“facilitate semi-automated enrichment analysis for Gene Ontology (GO) terms. The process consists of input of normalised gene expression measurements, gene-wise correlation or di fferential expression analysis, enrichment analysis of GO terms, interpretation and visualisation of the results.”我个人认为的核心过程在于差异表达显著基因的筛选以及GO term富集分析以及结果展示，给出的图挺拉风的。
资料是包在Bioconductor页给出的vignette，标题是Gene set enrichment analysis with topGO，October 27, 2010版。
前三章是功能介绍和快速演示，略过。
1. 首先是载入topGO和ALL包，利用Biocondutor人最爱的ALL包里的数据作为样例：
> library(topGO)
> library(ALL)
> data(ALL)
包导入时添加三个环境变量GOBPTerm、GOCCTerm和GOMFTerm，用ls()函数即可将其取出，保存为character vector，如：
> BPterms <- ls(GOBPTerm)
2. 导入后即可开始topGO分析。所有的分析均施用于topGOdata类对象，因此首先需要构建topGOdata对象。这类对象和BioC的很多对象一样，属于“大全”式的，构建好之后利用类方法大多两三下就能完成计算、给出结果。构建topGO这个大全，需要的数据包括：
a. 基因identifier，可以附上某种分数以便后面施用某种统计处理，分数可以是t检验的p值或者与某个表型的correlation等；
b. identifier和GO term间的map，如果是芯片数据的话BioC里有多种注释包，声明包的名称即可。至于我等蛋白界苦人，也能自己构建map，见下；
c. GO的层级结构，由GO.db提供，目前这个包只支持GO.db提供的结构：GOslim就再说了。
topGOdata对象构建函数的参数包括：
a. ontology，可指定要分析的GO term的类型，即BP、CC之类；
b. description：topGOdata的描述，可选；
c. allGenes：基因identifier的原始列表，和后面的geneSelectionFun联合作用，得出参与分析的基因，可以是numeric或factor。
d. geneSelectionFun：基因选择函数，如果前面allGenes是numeric的话就必须得指明此参数；
e. nodeSize：被认为富集的GO term辖下基因的最小数目（>=），默认为1。
f. annotationFun：基因identifier map到GO term的函数。
下面是蛋白人士的福音——自制用于将基因map到GO的注释库的解脱大法。样例文件在topGO/example下，geneid2go.map，方向为由基因id到GO。数据句法其实很简单，为gene_IDGO_ID1, GO_ID2, GO_ID3, ....。读map文件的方法是：
> readMappings(file = system.file("examples/geneid2go.map", package = "topGO"))
返回一个list。
还能倒着读，是个两步的过程，第一步见前，现在做第二步：
> inverseList(geneID2GO)
其中geneID2GO是刚才那一步的返回结果，这一步也返回一个list。强！
好了，矫情完了，开始真正干活儿了：
GOdata <- new("topGOdata",, allGenes = geneList, annot = annFUN.gene2GO, gene2GO = geneID2GO)
其中geneList是一个factor，必须named，即赋予了每个因子的名字，如：
> geneList <- factor(as.integer(geneNames %in% myInterestingGenes))
> names(geneList) <- geneNames
构建过程较为耗时，用时与所包含进去的基因的数目有关。完成后，输入GOdata即可得到该对象的概况信息。
4. topGOdata对象构建好后，即可利用这个包里的各种方法和函数做分析。首先可以做一些筛选工作以减小工作量、提高数据处理速度。不过暂时用不着，先跳过。
对于topGOdata对象，topGO包提供了多样的功能，比如修改对象的描述：
description(GOdata) <- paste(description(GOdata), "Object modified on:", format(Sys.time(), "%d %b %Y"), sep = " ")
利用函数descirption()修改描述，此外还可以查看。还可以查看对象里的基因identifier，如：
> genes(GOdata)
函数numGenes()可以查看对象包含的基因的数目。获取基因的分数需要用到函数geneScore()：
> selGenes <- sample(a, 10)
> gs <- geneScore(GOdata, whichGenes = selGenes)
gs为named vector，如果属性use.names为FALSE的话则vector为unnamed。
函数sigGenes()返回一个character vector，为各显著变化基因identifier。函数numSigGenes()则用于查看显著变化基因的数目。
还有函数updateGenes()可以修改topGOdata对象里所包含的基因，如果不想用new()重新构建，可以更换基因列表，在updateGenes()的参数中指定即可：
> GOdata <- updateGenes(GOdata, .geneList, topDiffGenes)
其中.geneList是另外构建的一个geneList对象。
想要看全部基因都对应上了哪些GO term，可用函数usedGO()得到一个character：
> ug <- usedGO(GOdata)
想看指定的几个GO term都对应上了多少基因，可以用：
> num.ann.genes <- countGenesInTerm(GOdata, sel.terms)
其中sel.terms是GO term构成的character vector，返回一个integer。
要是想知道具体是哪些基因对应到了这些GO term，可以：
> ann.genes <- genesInTerm(GOdata, sel.terms)
返回一个list。
此外还可以查看基因的分数：
> ann.score <- scoresInTerm(GOdata, sel.terms)
ann.score也是一个list，如果属性use.names设置为TRUE，则每个分数都会带上基因identier信息。
函数termStat()则返回一个dataframe，提供GO term的概况。每行为一个term，分别列出辖下基因数，其中属于显著变化的数目，以及预期的百分比（谁比谁还没闹明白……）。
5. 现在，终于可以做基因集富集分析（gene set enrichment analysis）啦！首先看看GSEA的三种方式：
a. 基于count，即仅要求输入一组基因，此种方式最为流行，可用Fisher's exact test、Hypegeometric
test和binomial test进行检验；
b. 基于基因的score或rank，可用Kolmogorov-Smirnov like tests（即05年那篇PNAS的GSEA文章所用方法），Gentleman's Category、t-test等方法；
c. 基于基因的表达，可从expression matrix直接分析，如Goeman's globaltest，以及GlobalAncova。
topGO提供两种分析方法，一种自由度更高但上手不易，本菜鸟还是跟着第二种走吧，较为用户友好但集成度较高。简单来说，就是用runTest()这个函数，要求三个主要的argument，一个是之前构建好的topGOdata类实例，第二个参数algorithm用于指定生成GO graph的方法，而参数statistic用于指定所用的检验方法，比如：
> resultFis <- runTest(GOdata, algorithm = "classic", statistic = "fisher")
> resultWeight <- runTest(GOdata, algorithm = "weight", statistic = "fisher")
> resultKS <- runTest(GOdata, algorithm = "classic", statistic = "ks")
> resultKS.elim <- runTest(GOdata, algorithm = "elim", statistic = "ks.elim")
此外algorithm还有weight01、weight和elim等，可用whichAlgorithms()查看；statistic还有t和ks等，可用whichTests()查看。不过有些algorithm和statistic的组合不能用，试试就晓得了。
6. runTest这一锤子买卖敲定后就能开始解读和展示结果了，得到的结果是topGOresult类的一个实例，其组成很简单，为对象的基本信息，以及各基因的分数（p值或者其他统计参数）。需要注意的是p值是未校正过的。想查看分数的话，可用函数score()：
> head(score(resultFis))
返回一个named numeric。如果指定了属性whichGO，则可以选择性返回一部分GO term所辖的基因的分数，如：
> pvalWeight <- score(resultWeight , whichGO = names(pvalFis))
返回resultFis和resultWeight共有的基因在resultWeight中的分数。有了这两组分数，可以做一些比较，比如关联分析：
> cor(pvalFis, pvalWeight)
或者作散点图：
> xyplot(pvalWeight ~ pvalFis)
通过函数geneData()，可以查看输入的基因的概况，包括有注释的基因、差异显著的基因、term最少节点数、辖制至少一个显著变化基因的term的数目。
数据输出方面，可以用函数GenTable()把所有分析结果打包存储和输出，构建这个结果的大杂烩是这样：
> GenTable(GOdata, classicFisher = resultFis, classicKS = resultKS, elimKS = resultKS.elim, orderBy = "elimKS", ranksOf = "classicFisher", topNodes = 10)
返回一个data frame。按列拷出然后自己加上三条线，就能看到有模有样的结果啦。
如果想查看单个term的分数分布，可以用函数showGroupDensity()，如：
> goID <- allRes[1, "GO.ID"] # allRes为上一个命令的返回结果，这行语句返回某个得到富集的GO term
> print(showGroupDensity(GOdata, goID, ranks = TRUE))
打印出一张图，显示输入的term所辖的基因的分数分布，及其不含有此条注释的基因的分数分布，以便比较这个term在这两类基因间的分数分布的差异，从而确认term是否得到了富集。
另外，函数printGene()用于得到某个term所辖的所有基因，并打印成一个表格，用法如下：
> gt <- printGenes(GOdata, whichTerms = goID, chip = affyLib, numChar = 40)
其中whichTerm参数可以指定多个term，此时则返回多个data fame；如果指定了参数file，则可保存为csv文件。需要注意的是目前只支持BioC提供的注释包，自制的不行。
7. 最后一项，结果的图形展示。我比较喜欢这个功能，呵呵。实现图形的打印输出的是两个函数，一个是showSigOfNodes()，将带有结构的富集term图打印到当前设备；另一个是printGraph()，可将图输出为PDF或者ps文件。例子：
> showSigOfNodes(GOdata, score(resultFis), firstSigNodes = 5, useInfo = "all")
> printGraph(GOdata, resultFis, firstSigNodes = 5, fn.prefix = "tGO", useInfo = "def", pdfSW = TRUE)
得到的图如http://topgo.bioinf.mpi-inf.mpg.de/topGO_grpah_logo.jpg中所示，得到富集的term为长方形，背景色越偏红色则富集的显著性越高。不爽的是节点所含的信息太多而字号太小，字体和背景色也不晓得能不能自主定制，如果能的话就太好了。
Win下装Rgraphviz有些问题，以我的水平暂时估计解决不了。不巧宿舍的ubuntu本本今天又抽疯了，饭后又不想回实验室。罢，晚上看KEGGgraph，明天来在lin下用topGO处理自己的数据吧。


转自http://marinzou.blogbus.com/logs/104671388.html