上一期我们一起学习了如何差异表达基因的注释，接下来的工作就是统计分析和可视化。

对于差异表达分析，最主要的两种统计分析就是GO的富集分析和KEGG通路的富集分析。

今天我们一起学习GO分析。

GO分析即基因本体论(Gene Onotology)分析，包括GO注释和GO富集分析。

GO是基因本体联合会(Gene Ontology Consortium)所建立的数据库，它由一组预先定义好的术语(GOterm)组成，这组术语对基因和蛋白质功能进行限定和描述，适用于各种物种，并能随着研究不断深入而更新。

GO中的每个术语都有唯一的一个GOID，GO的所有术语由有向无环图(DAG)来相互联系，术语之间通常有三种关系:"is_ a"、"part_ of" 和"regulates"。

在R语言中GO分析主要是由Bioconductor的GOstats包实现的。

我们这一期的代码内容是在上一期的内容的基础上运行的

Bioconductor基因芯片数据分析之注释

#加载所需R包

BiocManager::install("GOstats")library(GOstats)

#提取HG_U95Av2芯片中所有探针组对应的EntrezID，注意保证uniq

entrezUniverse <- unique(unlist(mget(rownames(eset),hgu95av2ENTREZID)))

#提取所有差异表达基因及其对应的EntrezID,注意保证uniq

entrezSelected <- unique(dif[!is.na(dif$EntrezID),"EntrezID"])

#设置GO富集分析的所有参数

params <- new("GOHyperGParams",geneIds=entrezSelected, universeGeneIds =entrezUniverse,annotation = affydb,ontology = "BP",pvalueCutoff = 0.001,conditional = FALSE,testDirection = "over")

#对所有的GOterm根据params参数做超几何检验。

hgOver <- hyperGTest(params)

#生成所有GOterm的检验结果报表

bp <- summary(hgOver)

#同时生成所有GOterm的检验结果文件，每个GOterm都有指向官方网站的链接，可以获得其详细信息

htmlReport(hgOver, file='ALL_ go.html')

#显示结果的前六行

head(bp)

从变量"bp"的结果可以看到每个显著性富集的GOterm含有六列信息（不包括行号）：

第1列是GO term的ID，该ID对应的内容在后面列出，如"GO:0022900"对应后面的"respiratory electron transport chain"；第2列"Pvalue"是超几何检验的P值；第3列"OddsRatio"是超几何分布中的比值比；第4列"ExpCount"是根据超几何分布，差异表达基因中期望属于这个GO term的基因数量；第5列"Count"是差异表达基因中实际属于这个GO term的基因数量；第6列"Size"是总基因中属于这个GO term的基因数量。

以"GO:0022904"为例，此次分析的总基因数量为8804，差异表达基因数量为113个基因中有75个基因（即"Size"）属于"GO:0022904"，如果从8804个基因中随机抽取113个基因，那么113个基因中期望属于"GO:0022904"的基因数量应该是2.25 （即"ExpCount"），而实际上是12个（即"Count"），根据这个情况，计算出来的P值应该是1.506871e-10（远远小于 0.01），因此可以说差异显著基因在"GO:0022904"上是显著富集的。

我们还通过函数htmlReport输出一个HTML的报告文件，它在前面六列信息的基础上，多加了一列GO term的描述，并且链接到GO的官方网站，便于读者进一步查看相关信息。