富集分析的R包神器-ClusterProfiler使用

作者：MC学公卫

来源：百味科研芝士

 什么是GO富集分析？

GeneOntology富集分析是高通量数据分析的标配，不管是转录组、甲基化、ChIP-seq还是重测序，都会用到对一个或多个集合的基因进行功能富集分析。分析结果可以指示这个集合的基因具有什么样的功能偏好性，进而据此判断相应的生物学意义。 GO富集分析是先筛选差异基因，再判断差异基因在哪些注释的通路存在富集。 

进行富集分析

1. 安装r包clusterProfiler

在本地：

1. > source('http:///biocLite.R')

2. > biocLite(\\'clusterProfiler\\')

2.基因列表Entrez ID与FC值的准备

利用一些差异表达的基因作为输入的基因，之前Seurat包的 FindAllMarkers有一批这些基因，但是没有EntrezID，需要进行一些处理。总的来说，就是需要整理出需要作富集分析的基因列表 先查看一下：

spleen.markers <- FindAllMarkers(object = spleen, only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25)

1. print(x= head(x=spleen.markers,n = 10))

2. p_val avg_logFC pct.1 pct.2 p_val_adj cluster gene

3. MS4A1 2.410289e-189 1.592266 0.974 0.258 3.773307e-185 0 MS4A1

4. CD79A 1.245084e-177 1.485304 0.971 0.270 1.949180e-173 0 CD79A

5. HLA-DRA 1.792014e-151 1.370242 1.000 0.556 2.805397e-147 0 HLA-DRA

6. CD79B 2.280181e-149 1.265725 0.938 0.315 3.569624e-145 0 CD79B

7. CD74 3.217361e-144 1.256090 0.997 0.871 5.036778e-140 0 CD74

8. HLA-DPB1 4.238869e-137 1.154889 0.997 0.616 6.635949e-133 0 HLA-DPB1

9. HLA-DRB5 6.911780e-136 1.123850 0.990 0.437 1.082039e-131 0 HLA-DRB5

10. HLA-DQB1 4.004476e-134 1.121013 0.940 0.366 6.269007e-130 0 HLA-DQB1

11. HLA-DRB1 2.026996e-133 1.115709 0.992 0.452 3.173262e-129 0 HLA-DRB1

12. HLA-DPA1 3.492170e-130 1.099676 0.995 0.570 5.466992e-126 0 HLA-DPA1

spleen.markers一共有3,108个基因marker。

获取基因列表

genelist <- spleen.markers$gene
head(genelist)

1. [1] 'MS4A1' 'CD79A' 'HLA-DRA' 'CD79B' 'CD74' 'HLA-DPB1'

加载包，下载人类基因组注释的数据org.Hs.eg.db

library(clusterProfiler)
BiocInstaller::biocLite('org.Hs.eg.db')
library(org.Hs.eg.db)

安装与细节可在此网站查看Genome wide annotation for Human（含引用文献）：

转化为EntrezID

1. s.EntrezID <- bitr(genelist, fromType='SYMBOL', toType='ENTREZID', OrgDb='org.Hs.eg.db')

2. head(s.EntrezID)

    SYMBOL ENTREZID
1    MS4A1      931
2    CD79A      973
3  HLA-DRA     3122
4    CD79B      974
5     CD74      972
6 HLA-DPB1     3115

查看行数：

1. dim(s.EntrezID)

[1] 2043    2

原本是3,108个，现在剩下2043个基因，原因是去除了重复的个数。

3. 进行富集分析

接下来进行富集分析：

1. s.ego <- enrichGO(gene = s.EntrezID.df$ENTREZID, OrgDb = \\'org.Hs.eg.db\\', ont = \\'ALL\\', pAdjustMethod = \\'BH\\',pvalueCutoff = 0.05, qvalueCutoff = 0.2, keyType = \\'ENTREZID\\')

3. head(s.ego)

筛选的p值为0.05，q值为0.2

           ONTOLOGY         ID                                       Description GeneRatio
GO:0042119       BP GO:0042119                             neutrophil activation  211/1931
GO:0002283       BP GO:0002283 neutrophil activation involved in immune response  208/1931
GO:0043312       BP GO:0043312                          neutrophil degranulation  207/1931
GO:0002446       BP GO:0002446                      neutrophil mediated immunity  208/1931
GO:0006613       BP GO:0006613     cotranslational protein targeting to membrane   77/1931
GO:0045047       BP GO:0045047                           protein targeting to ER   78/1931
             BgRatio       pvalue     p.adjust       qvalue
GO:0042119 496/17381 2.402001e-74 1.352326e-70 9.264138e-71
GO:0002283 486/17381 7.667076e-74 2.158282e-70 1.478535e-70
GO:0043312 485/17381 3.218981e-73 6.040954e-70 4.138367e-70
GO:0002446 499/17381 2.372136e-71 3.338782e-68 2.287239e-68
GO:0006613  96/17381 5.653664e-56 6.366026e-53 4.361058e-53
GO:0045047 100/17381 5.659096e-55 5.310118e-52 3.637706e-52
geneID
GO:0042119 HVCN1/EEF2/SELL/IMPDH2/PTPN6/CTSH/CCL5/B2M/HSPA8/PTPRC/RAB27A/STOM/SURF4/HSP90AA1/BIN2/HSPA1A/HSPA6/HSPA1B/ATP6V0A1/CTSZ/RNASET2/HSP90AB1/SERPINA1/FCN1/CST3/CD68/CFP/FGL2/S100A12/FPR1/LYZ/CD14/CD36/MNDA/CFD/FCGR2A/MGST1/TLR2/CD93/C5AR1/TYROBP/FCER1G/LGALS3/
………………
………………
           Count
GO:0042119   211
GO:0002283   208
GO:0043312   207
GO:0002446   208
GO:0006613    77
GO:0045047    78

剩下1931个基因：

1. dim(s.ego)

2. [1] 1884 10

3. > dim(s.ego[s.ego$ONTOLOGY==\\'BP\\',])

4. [1] 1507 10

5. > dim(s.ego[s.ego$ONTOLOGY==\\'CC\\',])

6. [1] 215 10

7. > dim(s.ego[s.ego$ONTOLOGY==\\'MF\\',])

8. [1] 162 10

看来大部分的基因富集在BP中。 这里科普一下BP是指Biological Process，CC是指Cellular Component，MF是指Molecular Function。可以根据自己需要各取所需。

4. 数据可视化

可以作柱状图
s.barplot <- barplot(s.ego, showCategory = 20, drop = T)
s.barplot

1. s.barplot <- barplot(s.ego, showCategory = 25, drop = T)

2. s.barplot

3. 

可以作点图，作用与柱状图类似
dotplot(s.ego)

1. 也可以作富集网络

2. enrichMap(s.ego)

3. 

这里的GO较多，若要清晰展示排名前10的GO可加上参数n=10，即enrichMap(s.ego, n = 10)再进行绘图。

来源：简书