差异分析得到的结果注释一文就够

通过前面的讲解，我们顺利的了解了GEO数据库以及如何下载其数据，得到我们想要的表达矩阵，也学会了两个常用的套路分析得到的表达矩阵，就是GSEA分析和差异分析。

历史目录： 解读GEO数据存放规律及下载，一文就够解读SRA数据库规律一文就够从GEO数据库下载得到表达矩阵一文就够 GSEA分析一文就够（单机版+R语言版）根据分组信息做差异分析- 这个一文不够的

但是差异分析通过自定义的阈值挑选了有统计学显著的基因列表后我们其实是需要对它们进行注释才能了解其功能，最常见的就是GO/KEGG数据库注释咯，当然也可以使用Reactome和Msigdb数据库来进行注释。而最常见的注释方法就是超几何分布检验了咯。

超几何分布检验原理

其实我在生信菜鸟团博客里面讲的非常清楚了这个统计学原理：用超几何分布检验做富集分析 http://www.bio-info-trainee.com/1225.html

直接上代码： https://github.com/jmzeng1314/humanid/blob/master/R/hyperGtest_jimmy.R

维基百科的解释是：超几何分布是统计学上一种离散概率分布。它描述了由有限个物件中抽出n个物件，成功抽出指定种类的物件的个数（不归还）。

例如在有N个样本，其中m个是不及格的。超几何分布描述了在该N个样本中抽出n个，其中k个是不及格的概率

若n=1，超几何分布还原为伯努利分布。

若N接近∞，超几何分布可视为二项分布。

作为离散概率分布的超几何分布尤其指在抽样试验时抽出的样品不再放回去的分布情况。在一个容器中一共有N个球，其中M个黑球,(N-M)个红球，通过下面的超几何分布公式可以计算出，从容器中抽出的n个球中（抽出的球不放回去）有k个黑球的概率是多少：

例如，容器中一共10个球，其中6个黑色，4个白色，一共抽5次（抽出的球不放回去），在这5个球中有3个黑球的概率是：

但是我们要算的不是恰好有3个黑球的概率，而是我们的基因富集问题，那么在R里面如何实现呢？

phyper(62-1, 1998, 5260-1998, 131, lower.tail=FALSE)

下面是我自己的理解：

超几何分布很简单，球分成黑白两色，数量已知，那么你随机抽有限个球，应该抽多少白球的问题！

公式就是exp_count=n*M/N

然后你实际上抽了多少白球，就可以计算一个概率值！

换算成通路的富集概念就是，总共有多少基因（这个地方值得注意，主流认为只考虑那些在KEGG等数据库注释的背景基因），你的通路有多少基因，你的通路被抽中了多少基因（在差异基因里面属于你的通路的基因），这样的数据就足够算出上面表格里面所有的数据啦！

原理讲完了，就该讲如何做了：https://github.com/jmzeng1314/humanid/blob/master/R/batch_enrichment.R

需要自己搜索什么是GO/KEGG/BIOCARTA/REACTOME等数据库

http://www.cnblogs.com/emanlee/archive/2011/08/02/2125314.html

虽然懂了原理可以让我们更方便的理解结果，但实际数据分析过程中我们通常不会自己写代码完成这个超几何分布检验，因为有现成的，而且非常好用的轮子。

强烈推荐Y叔的包clusterProfiler

首先需要理解下面的 geneList和 gene这两个数据集。然后，理解 GO/KEGG/REACTOME/MSIGDB 这4个数据库结构，及对应的生物学一样。接着，理解超几何分布建议，GSEA这两个算法。最后把下面的代码跑一遍即可。(因为Y叔的包更新太频繁，我只能保证下面的代码，今天是有效的，所以赶快试一下吧)

library(clusterProfiler)

### get the universal genes and sDEG 

data(geneList, package="DOSE")
gene <- names(geneList)[abs(geneList) > 2]
gene.df <- bitr(gene, fromType = "ENTREZID",
                toType = c("ENSEMBL", "SYMBOL"),
                OrgDb = org.Hs.eg.db)
head(gene.df)

### Using MSigDB gene set collections
gmtfile <- system.file("extdata", "c5.cc.v5.0.entrez.gmt", package="clusterProfiler")
c5 <- read.gmt(gmtfile)
egmt2 <- GSEA(geneList, TERM2GENE=c5, verbose=FALSE)
head(egmt2)[,1:6]
gseaplot(egmt2, geneSetID = "EXTRACELLULAR_REGION")




## GO analysis 

ego <- enrichGO(gene          = gene,
                universe      = names(geneList),
                OrgDb         = org.Hs.eg.db,
                ont           = "CC",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.01,
                qvalueCutoff  = 0.05,
                readable      = TRUE)
head(ego)[,1:6]
ego2 <- gseGO(geneList     = geneList,
              OrgDb        = org.Hs.eg.db,
              ont          = "CC",
              nPerm        = 1000,
              minGSSize    = 100,
              maxGSSize    = 500,
              pvalueCutoff = 0.05,
              verbose      = FALSE)
head(ego2)[,1:6]
## KEGG pathway analysis
kk <- enrichKEGG(gene         = gene,
                 organism     = 'hsa',
                 pvalueCutoff = 0.05)
head(kk)[,1:6]
kk2 <- gseKEGG(geneList     = geneList,
               organism     = 'hsa',
               nPerm        = 1000,
               minGSSize    = 120,
               pvalueCutoff = 0.05,
               verbose      = FALSE)
head(kk2)[,1:6]
gseaplot(kk2, geneSetID = "hsa04145")

### Reactome pathway analysis
library(ReactomePA)
pp <- enrichPathway(gene         = gene,
                 organism     = 'human',
                 pvalueCutoff = 0.05)
head(pp)[,1:6]
pp2 <- gsePathway(geneList     = geneList,
               organism     = 'human',
               nPerm        = 1000,
               minGSSize    = 120,
               pvalueCutoff = 0.05,
               verbose      = FALSE)
head(pp2)[,1:6]

### Disease analysis
library(DOSE)
dd <- enrichDO(gene         = gene)
head(dd)[,1:6]
dd2 <- gseDO(geneList     = geneList,
                  nPerm        = 1000,
                  minGSSize    = 120,
                  pvalueCutoff = 0.05,
                  verbose      = FALSE)
head(dd2)[,1:6]

下游的GO/KEGG注释一般是得到如下表格：

也可以有一些简单的可视化展现：