前言:
微博参与话题 #给你四年时间你也学不会生信#
主要参考：GEO数据挖掘小尝试：（三）利用clusterProfiler进行富集分析
Y叔开发的R包clusterProfiler 的确是最好用的，没有之一,可参看为Y叔疯狂打call

1、安装clusterProfiler

> source("http://bioconductor.org/biocLite.R")
> biocLite('clusterProfiler')

2、ID转换

对于没有转换的gene ID，clusterProfiler也提供了bitr方法进行转换ID：

# example:
> x <- c("GPX3",  "GLRX",   "LBP",   "CRYAB", "DEFB1", "HCLS1",   "SOD2",   "HSPA2", 
       "ORM1",  "IGFBP1", "PTHLH", "GPC3",  "IGFBP3","TOB1",    "MITF",   "NDRG1", 
       "NR1H4", "FGFR3",  "PVR",   "IL6",   "PTPRM", "ERBB2",   "NID2",   "LAMB1", 
       "COMP",  "PLS3",   "MCAM",  "SPP1",  "LAMC1", "COL4A2",  "COL4A1", "MYOC",  
       "ANXA4", "TFPI2",  "CST6",  "SLPI",  "TIMP2", "CPM",     "GGT1",   "NNMT",
       "MAL",   "EEF1A2", "HGD",   "TCN2",  "CDA",   "PCCA",    "CRYM",   "PDXK",  
       "STC1",  "WARS",  "HMOX1", "FXYD2", "RBP4",   "SLC6A12", "KDELR3", "ITM2B")
> head(eg)
  SYMBOL ENTREZID         ENSEMBL
1   GPX3     2878 ENSG00000211445
2   GLRX     2745 ENSG00000173221
3    LBP     3929 ENSG00000129988
4  CRYAB     1410 ENSG00000109846
5  DEFB1     1672 ENSG00000164825
6  DEFB1     1672 ENSG00000284881

3、GO、KEGG富集分析

OrgDb库
enrichGO默认gene type是entrezID，但其他OrgDb支持的类型（ENSEMBLE，SYMBOL等）都可以通过参数keyType指定。
gene的ID type不一样，富集的结果也会有稍微的差异。
原gene list是entrezID，直接通过bitr转换成ensembl和symbol，分别做enrichGO。
发现entrezedID可能对应多个ENSEMBL的。
entrezedID和SYMBOL是一一对应的

3.1 GO富集分析

> library(clusterProfiler)
> genelist <- egcounts$ENTREZID
> head(genelist)
[1] 2878 2745 3929 1410 1672 1672
> # 检查是否有重复
> genelist[duplicated(genelist)]
[1] 1672 6648 1428

由于clusterProfiler富集分析推荐的输入文件是Entrez ID，因此这里提取的是Entrez ID，接下来就可以进行富集分析了：

> go <- enrichGO(genelist, OrgDb = org.Hs.eg.db, ont='ALL',pAdjustMethod = 'BH',pvalueCutoff = 0.05, 
+                qvalueCutoff = 0.2,keyType = 'ENTREZID')
> head(go)
           ONTOLOGY         ID                          Description GeneRatio   BgRatio       pvalue
GO:0001503       BP GO:0001503                         ossification      9/56 371/17653 2.299994e-06
GO:0030198       BP GO:0030198    extracellular matrix organization      8/56 341/17653 1.131006e-05
GO:0045926       BP GO:0045926        negative regulation of growth      7/56 244/17653 1.151994e-05
GO:0043062       BP GO:0043062 extracellular structure organization      8/56 395/17653 3.250950e-05
GO:0003416       BP GO:0003416             endochondral bone growth      3/56  22/17653 4.461661e-05
GO:0001649       BP GO:0001649           osteoblast differentiation      6/56 205/17653 4.563530e-05
              p.adjust      qvalue                                        geneID Count
GO:0001503 0.003969790 0.003019045 5744/2719/3486/10140/2261/3569/6696/4653/6781     9
GO:0030198 0.006627806 0.005040480      22795/3912/1311/6696/3915/1284/1282/7077     8
GO:0045926 0.006627806 0.005040480             3929/1410/2719/2261/6696/978/5950     7
GO:0043062 0.013127755 0.009983723      22795/3912/1311/6696/3915/1284/1282/7077     8
GO:0003416 0.013127755 0.009983723                                2261/1311/6781     3
GO:0001649 0.013127755 0.009983723                5744/3486/10140/3569/6696/4653     6
> dim(go)
[1] 68 10
> dim(go[go$ONTOLOGY=='BP',])
[1] 46 10
> dim(go[go$ONTOLOGY=='CC',])
[1] 15 10
> dim(go[go$ONTOLOGY=='MF',])
[1]  7 10

> # 进行简单的可视化
> barplot(go,showCategory=20,drop=T)
> dotplot(go,showCategory=50)

Rplot20.png

Rplot21.png

3.2 KEGG通路富集
KEGG通路富集函数用法与GO富集分析方法类似：

> kegg <- enrichKEGG(genelist, organism = 'hsa', keyType = 'kegg', pvalueCutoff = 0.05,pAdjustMethod = 'BH', minGSSize = 10,maxGSSize = 500,qvalueCutoff = 0.2,use_internal_data = FALSE)
> head(kegg)
               ID                Description GeneRatio  BgRatio       pvalue     p.adjust
hsa04512 hsa04512   ECM-receptor interaction      6/35  82/7470 1.830700e-06 0.0002324989
hsa04151 hsa04151 PI3K-Akt signaling pathway      9/35 354/7470 2.568570e-05 0.0012944211
hsa04510 hsa04510             Focal adhesion      7/35 199/7470 3.057688e-05 0.0012944211
hsa05146 hsa05146                 Amoebiasis      5/35  96/7470 7.552538e-05 0.0023979308
hsa05222 hsa05222     Small cell lung cancer      4/35  93/7470 8.773963e-04 0.0222858657
hsa05134 hsa05134              Legionellosis      3/35  55/7470 2.092473e-03 0.0442906717
               qvalue                                       geneID Count
hsa04512 0.0001984864                3912/1311/6696/3915/1284/1282     6
hsa04151 0.0011050591 2261/3569/2064/3912/1311/6696/3915/1284/1282     9
hsa04510 0.0011050591           2064/3912/1311/6696/3915/1284/1282     7
hsa05146 0.0020471353                     3569/3912/3915/1284/1282     5
hsa05222 0.0190256458                          3912/3915/1284/1282     4
hsa05134 0.0378113484                               3306/3569/1917     3
> dim(kegg)
[1] 6 9
> # 简单可视化
> dotplot(kegg, showCategory=30)

Rplot22.png