还在为基因通路富集担心你的发际线么？ClusterProfiler通路富集-让你的发际线无忧！

单个基因水平上能反映的生物学信息有限，很多时候要进行通路富集分析，来从系统水平上反映出一组基因与哪些生物学通路相关。

今天我们就来谈谈Y叔开发的ClusterProfiler包做通路富集时的应用场景和详细步骤（以Homo sapiens为例）。

Step1.文件准备（如图1）。

1. 将所有文件存放于本地路径中（如D:/Bioming）。

2. 感兴趣的基因集合文件，每个基因为一行。（感兴趣的基因集合可以是差异表达基因、差异甲基化基因、突变基因集合等）。文件格式如图2。

3. 若感兴趣的基因集合是基于特定的panel芯片得到的，而并非全基因组数据，则需要准备背景基因集合文件，文件中包含了panel中所覆盖的所有基因，格式同图2。

4. 若需要将两组感兴趣的基因集合进行通路比较（如原发灶突变基因与转移灶突变基因所富集通路比较等），需要准备另一组文件感兴趣的基因集合2。

图1

图2

Step2.R语言执行通路富集分析。

进入到R界面，若未安装过ClusterProfiler包，用如下命令安装：

if (!requireNamespace("BiocManager", quietly = TRUE))

install.packages("BiocManager")

BiocManager::install("clusterProfiler", version = "3.8")

加载ClusterProfiler包：

library(org.Hs.eg.db)

library(clusterProfiler)

分析开始：

setwd("D:\\ Bioming") #设置默认路径，即存放文件的本地路径，需要注意的是在windows系统中路径层级的符号为“\”，输入到R中时需要改成双反斜杠”\\”或正斜杠”/”。

CaseGene=read.table("Case.txt",header=F)[,1] #读取感兴趣基因集合文件，第一列是基因名。

CaseGeneSet=bitr(CaseGene, 'SYMBOL', "ENTREZID", "org.Hs.eg.db")[, "ENTREZID"] #将Gene名字转化为基因ID，若输入的是基因SYMBOL需要用该命令转化成基因ID，若输入的是基因ID，则忽略该命令。

Case_KEGG <- enrichKEGG(CaseGeneSet, keyType = "kegg",pvalueCutoff=0.05,qvalueCutoff=0.05,pAdjustMethod = "BH", minGSSize = 5, maxGSSize = 500,organism = "hsa", use_internal_data=FALSE) #KEGG富集分析

###Parameter解析

#CaseGeneSet: 感兴趣基因集合

#keyType: 设置类型，可选择的有"kegg", 'ncbi-geneid', 'ncib-proteinid' and 'uniprot'

#pvalueCutoff：设置p值的阈值

#qvalueCutoff：设置q值阈值（既校正后的阈值）

# pAdjustMethod：对p值进行校正的方法，可选方法有"holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none"

#minGSSize: 设置最小基因个数

#maxGSSize：设置最大基因个数

#orgnaism：物种选择，可参考：http://www.genome.jp/kegg/catalog/org_list.html

# use_internal_data：设置注释参考的数据库，TRUE—参考最新的线上KEGG数据库，FALSE—参考KEGG.db数据库。

运行完上述命令后得到的Case_KEGG结果如下图。

接下来进行富集结果存储。

运行命令：

write.csv(summary(Case_KEGG),"KEGG-enrich.csv",row.names =F)#导出Case_KEGG结果至默认路径下。

Step3.可视化通路富集结果。

运行命令：

barplot(Case_KEGG, showCategory=15,title="Barplot of Case KEGG Enrichment")#绘制条形图

运行命令：

dotplot(Case_KEGG,title=" Dotplot of Case KEGG Enrichment ")#绘制气泡图

运行命令：

enrichMap(Case_KEGG, layout=igraph::layout.kamada.kawai)#网络图，每个节点是一个富集到的pathway，若pathways之间有重叠的感兴趣基因，则自动将这两个通路用线连接。默认画top30个富集到的pathways, 节点大小对应该pathway下富集到的感兴趣基因个数，节点的颜色对应p.adjust的值，从小到大，对应蓝色到红色。

到此，已经可以完整呈现出感兴趣基因的通路富集分析了。那么如果这些感兴趣的基因并非在全基因组范围内找到的，而且来源于特定的芯片，如何科学地做出通路富集分析呢？

很简单，只需要将芯片所覆盖的所有基因作为背景基因，存放于文件中(背景基因集合)。

运行命令：

BackGround=read.table("BackGroundGene.txt",header=F)[,1]#读取背景基因，第一列是基因名。

BackGroundSet=bitr(BackGround, 'SYMBOL', "ENTREZID", "org.Hs.eg.db")[, "ENTREZID"]#将Gene名字转化为基因ID，若输入的是基因ID，则忽略该命令。

Case_KEGG <- enrichKEGG(CaseGeneSet, keyType = "kegg",pvalueCutoff=0.05,qvalueCutoff=0.05,pAdjustMethod = "BH", minGSSize = 5, maxGSSize = 500,organism = "hsa", use_internal_data=FALSE, universe=BackGroundSet)# 只需要在Step2的enrichKEGG这一步设置universe参数。

假如你想比较两组感兴趣基因集合的通路异同又应如何操作呢？需准备另一组感兴趣的基因集合（如图1中的感兴趣基因集合2，格式如图2）。导入第二组感兴趣的基因集合，并进行通路富集分析。

运行命令：

ControlGene=read.table("Control.txt",header=F)[,1]

ControlGeneSet=bitr(ControlGene, 'SYMBOL', "ENTREZID", "org.Hs.eg.db")[, "ENTREZID"]

Control_KEGG <- enrichKEGG(ControlGeneSet, keyType = "kegg",pvalueCutoff=0.05,qvalueCutoff=0.05,pAdjustMethod = "BH", minGSSize = 5, maxGSSize = 500,organism = "hsa")

cp = list(Case=CaseGeneSet, Control=ControlGeneSet) # 合并两个数据集，并转换为列表

compare <- compareCluster(cp, fun="enrichKEGG", organism="hsa", pvalueCutoff=0.05,pAdjustMethod = "BH",qvalueCutoff = 0.05)#比较两组通路富集结果。

dotplot(compare, showCategory=15,font.size = 9,includeAll=TRUE)#画两组集合的通路比较图

今天的基因通路富集就到此结束了，不知道你学会了么？