使⽤clusterProfiler进⾏GO、KEGG富集分析（有参情况）
倒垃圾
寻找差异表达的基因并识别它们的功能，是我们进⾏RNA测序的最主要⽬的。很明显，这些差江苏灌云
异的基因必然与功能改变密切相关，例如，⽐较患病个体与正常个体的组织表达谱，不难想到
这些显著失调的基因参与了⽣物学过程、信号通路等，导致了疾病的发⽣。
前⾯已经讲了如何使⽤DESeq2、edgeR基于转录组测序获得的基因表达值鉴定差异表达基因。
那么，后续如何继续通过⽣信分析的⽅法，探索差异表达的基因发挥了怎样的功能，参与了哪
些调控通路呢？
我们平时看RNA-q相关的⽂献时，⽂章中在鉴定了差异表达的基因后，⼤都会在随后承接⼏
句关于这些失调基因所涉及通路的描述。例如，讨论这些差异基因主要映射到哪些GO或KEGG
分类条⽬中，以说明基因表达的改变会导致哪些调控途径原有功能失调，进⽽与表型联系起
来。通常称这种分析为GO、KEGG富集分析。
⽬前，能够进⾏GO、KEGG富集分析的⼯具有很多，不同的⼯具之间在算法、数据库组成上略
巨大的英语有不同，因此结果也可能⼤相径庭。本篇就先以R包clusterProfiler的⽅法为例，展⽰如何基于给
定的基因列表分析它们的GO、KEGG功能。
clusterProfiler包的安装
对于clusterProfiler的安装也很简单，⼀般情况下，直接通过Bioconductor安装clusterProfiler就
可以了。
#bioconductor 安装
#install.packages('BiocManager') #需要⾸先安装 BiocManager，如果尚未安装请先执⾏该步
BiocManager::install('clusterProfiler')
clusterProfiler的GO富集分析（有参向）
⾸先来看GO富集分析。
注：这⾥均对于有参考基因组的情况⽽⾔的，⽆参分析暂不涉及。
1、准备输⼊数据
待分析的数据就是⼀串基因名称了，可以是enmbl id、entrze id或者symbol id等类型都可
以。把基因名称以⼀列的形式排开，放在⼀个⽂本⽂件中（例如命名“”）。Excel中查
看，就是如下⽰例这种样式。全国交通日
学习党章的意义
输⼊数据，待分析的基因名称
常见交通标志2、加载参考物种的基因注释数据库
对于有参考基因组物种的分析，⾸先需要指定该物种的基因数据库。
存在两组情况，⼀种是常见物种，另⼀种是不常见物种。
#（1）对于常见物种，例如⼈类，有些专门的 R 包数据库，例如⼈类参考基因组 hg19 的
library(db)
#（2）对于不常见的物种，但却是存在参考基因组的情况
#通过 AnnotationHub 包索引基因组，例如期望找绵⽺（Ovis aries）的注释库
library(AnnotationHub)
hub <- AnnotationHub()
query(hub, 'Ovis aries') #输⼊绵⽺（Ovis aries）的名称进⾏匹配
sheep <- hub[['AH72269']] #返回了数据库编号 AH72269，就可以加载该库
3、GO富集分析
加载了注释库之后，读取基因列表⽂件，并使⽤clusterProfiler的内部函数enrichGO()即可完成GO富集分析。
library(clusterProfiler)
#读取基因列表⽂件中的基因名称
genes <- read.delim('', header = TRUE, stringsAsFactors = FALSE)[[1]]
#GO富集分析
#对于加载的注释库的使⽤，以上述为例，就直接在 OrgDb 中指定⼈（db）或绵⽺（sheep）
< <- enrichGO(gene = genes, #基因列表⽂件中的基因名称
OrgDb = 'sheep', #指定物种的基因数据库，⽰例物种是绵⽺（sheep）
keyType = 'ENTREZID', #指定给定的基因名称类型，例如这⾥以 entrze id 为例
ont = 'ALL', #可选 BP、MF、CC，也可以指定 ALL 同时计算 3 者
pAdjustMethod = 'fdr', #指定 p 值校正⽅法
pvalueCutoff = 0.05, #指定 p 值阈值，不显著的值将不显⽰在结果中
qvalueCutoff = 0.2, #指定 q 值阈值，不显著的值将不显⽰在结果中
readable = FALSE)
#例如上述指定 ALL 同时计算 BP、MF、CC，这⾥将它们作个拆分后输出
BP <- [$ONTOLOGY=='BP', ]
CC <- [$ONTOLOGY=='CC', ]
MF <- [$ONTOLOGY=='MF', ]
write.table(as.data.frame(BP), '', p = '\t', row.names = FALSE, quote = FALSE)
write.table(as.data.frame(CC), '', p = '\t', row.names = FALSE, quote = FALSE)
write.table(as.data.frame(MF), '', p = '\t', row.names = FALSE, quote = FALSE)
GO富集分析结果表格
对于各列内容：
ONTOLOGY，GO分类BP（⽣物学过程）、CC（细胞组分）或MF（分⼦功能）；
ID和Description，富集到的GO id和描述；
GeneRatio和BgRatio，分别为富集到该GO条⽬中的基因数⽬/给定基因的总数⽬，以及该条⽬中背景基因总数⽬/该物种所有已知的GO功能基因数⽬；
pvalue、p.adjust和qvalue，p值、校正后p值和q值信息；
geneID和Count，富集到该GO条⽬中的基因名称（分析中使⽤的entrze id，故这⾥也显⽰的entrze id）和数⽬。
注：如期望显⽰其它类型的基因id，如通俗的symbol id等类型，除了更改为使⽤symbol id的基因名称做分析外，还可以通过基因名称转换的⽅式对entrze id和symbol id作个匹配转换。
clusterProfiler的KEGG富集分析（有参向）然后是KEGG富集分析。
同样地，这⾥均对于KEGG数据库中已经收录的物种⽽⾔的，⽆参分析暂不涉及。
1、准备输⼊数据
相⽐上述GO富集，clusterProfiler的KEGG富集分析⽅法特殊，它⽆需加载本地注释库，⾃动使⽤KEGG的在线数据库进⾏注释，因此给定的基因名称只能识别entrze id。把entrze id的基因名称以⼀列的形式排开，放在⼀个⽂本⽂件中（例如命名“”）。Excel中查看，就是如下⽰例这种样式。
输⼊数据，待分析的基因名称
2、KEGG富集分析
读取基因列表⽂件，并使⽤clusterProfiler的内部函数enrichKEGG()即可完成KEGG富集分析。library(clusterProfiler)
#读取基因列表⽂件中的基因名称，注意这⾥只能⽤ entrze id
genes <- read.delim('', header = TRUE, stringsAsFactors = FALSE)[[1]]
#每次打开R计算时，它会⾃动连接kegg官⽹获得最近的物种注释信息，因此数据库⼀定都是最新的
kegg <- enrichKEGG(
gene = genes, #基因列表⽂件中的基因名称
一颗纽扣keyType = 'kegg', #kegg 富集
organism = 'oas', #例如，oas 代表绵⽺，其它物种更改这⾏即可
pAdjustMethod = 'fdr', #指定 p 值校正⽅法
pvalueCutoff = 0.05, #指定 p 值阈值，不显著的值将不显⽰在结果中
qvalueCutoff = 0.2, #指定 q 值阈值，不显著的值将不显⽰在结果中
)
#输出结果
write.table(kegg, '', p = '\t', quote = FALSE, row.names = FALSE)
KEGG富集分析结果表格
对于各列内容：
ID和Description，富集到的KEGG id和描述；
GeneRatio和BgRatio，分别为富集到该KEGG条⽬中的基因数⽬/给定基因的总数⽬，以及该条⽬中背景基因总数⽬/该物种所有已知的KEGG功能基因数⽬；
pvalue、p.adjust和qvalue，p值、校正后p值和q值信息；
geneID和Count，富集到该KEGG条⽬中的基因名称（分析中使⽤的entrze id，故这⾥也显⽰的entrze id）和数⽬。
注：如期望显⽰其它类型的基因id，如通俗的symbol id等类型，由于该分析中只能输⼊entrze id，因此可以通过基因名称转换的⽅式对entrze id和symbol id作个匹配转换。clusterProfiler附带的作图功能
此外，clusterProfiler中也额外提供了⼀系列的可视化⽅案⽤于展⽰本次富集分析结果，具有极⼤的便利。
#clusterProfiler 包⾥的⼀些默认作图⽅法，例如
barplot(kegg) #富集柱形图
dotplot(kegg) #富集⽓泡图
cnetplot(kegg) #⽹络图展⽰富集功能和基因的包含关系
emapplot(kegg) #⽹络图展⽰各富集功能之间共有基因关系
heatplot(kegg) #热图展⽰富集功能和基因的包含关系
左边小肚子疼
⼀些富集柱形图、⽓泡图、⽹络图等
