clusterprofiler,小鼠clusterProfiler：GO KEGG

clusterprofiler,小鼠clusterProfiler：GO KEGG详细介绍

本文目录一览： GEO数据挖掘小尝试：（三）利用clusterProfiler进行富集分析

安装clusterProfiler：

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

可以看到，这里转换ID的对应文件来源于"org.Hs.eg.db"这个包。

在开始富集分析之前先看看GO和KEGG富集分析的方法以及参数：

导入数据，这是一个整合数据，在这里我们要用到的只是entrez ID列和最后一列（logFC）：

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

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

我们继续使用上面的数据进行KEGG富集分析：

这里使用clusterProfiler里面的 GSEA 函数进行GSEA富集分析，并与使用超几何分布富集（ enricher 函数）的结果进行简单比较， enricher 函数与 GSEA 函数用法基本相同，因此这里只给出 GSEA 的用法及参数。

在进行富集分析之前需要对数据做一个预处理——排序。

这里使用的是broad GSEA提供的gene sets 来提供TERM2GENE：

万事俱备，只欠东风。现在可以开始分析了，先进行超几何分布的富集分析：

再做GSEA富集分析，在此之前需要对输入gene list做一下处理，包括三步：

输入文件准备好了尽可以进行GSEA富集分析了：

不要问我为什么要pvalueCutoff设置0.8，因为一直调大到0.8才富集到结果。。。可见数据应该是有问题的，但这里作为一个实践就不管那么多了。
而且，clusterProfiler还支持GSEA的GO、KEGG富集。

顺便再利用上面处理好的glist进行一下KEGG富集到的某一条通路的可视化：

上面的命令会在当前目录生成3个文件：一个原始KEGG通路图片，一个标注了上下调基因的，最后一个文本文件则是一些KEGG通路信息。

参考：
clusProfiler命令参考手册
GSEA的分析汇总
Statistical analysis and visualization of functional profiles for genes and gene clusters
GSEA分析是个什么鬼？(上)
GSEA是个什么鬼?(下)

clusterProfiler导入失败解决办法

最近在安装光叔的clusterProfiler包时，遇到了以下报错：

即可以正常安装，但不能导入。

显示错误：package or namespace load failed for 'clusterProfiler': 'namespace:rvcheck'没有出口'get_fun_from_pkg'这个对象，随后电脑右侧弹出360杀毒软件的提示，把某dll文件列入了风险清单，准备删除。于是我思考，可能是由于杀毒软件的原因，于是我重装了Rstudio(version 1.4.1717)和R(version 4.0.5)，关掉了360杀毒和安全卫士，利用github接口安装包：

导入成功了！

R包 clusterProfiler： GO和KEGG富集结果显示基因symbol

注：1）MF和CC方法同BP,将BP改为MF,CC即可。
2）可视化中，showCategory为显示的item数，scale_y_discrete则调节label过长的情况，让图片看起来
更美观。
3)检查结果，可见geneID展示为gene symbol。

（1）在enrichGO函数中，设置readable = TRUE;
(2)用setReadable函数，对GO或者KEGG结果进行转化即可。

R语言：clusterProfiler进行GO富集分析和Gene_ID转换

ID转换用到的是 bitr() 函数，bitr()的使用方法：
org.Hs.eg.db包含有多种gene_name的类型
keytypes() ：keytypes(x)，查看注释包中可以使用的类型 columns() ：类似于keytypes()，针对org.Hs.eg.db两个函数返回值一致 select() ：select(x, keys, columns, keytype, ...) eg.
函数enrichGO()进行GO富集分析，enrichGO()的使用方法：
举例：

快捷查找KEGG里的通路和基因

1.快捷查找ID对应的description，知道通路对应的编号是多少。
2.找出某一个/几个通路里的全部基因，用来做单独的下游分析。

如果是要做KEGG的富集分析，clusterProfiler可以搞定： https://bioconductor.org/packages/release/bioc/vignettes/clusterProfiler/inst/doc/clusterProfiler.html

想看kegg通路图的话，用R包pathview来看，看函数的帮助文档就行。

不批量找的话，直接网站搜最简单 https://www.genome.jp/kegg/kegg2.html

需要找全部的对应关系，基于前面讲的msigdbr可以完成： https://www.jianshu.com/p/0098baf2df46

msigdb里面本来就包括了kegg,而且挺齐全的,ID,description,基因，全都有啦。

在org.Hs.eg.db包里有：

看起来像一堆密码？这个列表，名字是通路的id，只是省略了hsa，内容是基因的entrizid。

举个栗子，提取hsa03030里的基因,并且转换成symbol。

网页在线做富集分析

? ?? 用clusterProfiler做其实够用了，网页在线做没多大必要 ，能够起起对照作用吧。网页在线做，后续挑出自己想要的模块，会方便一些，但自己用包做手动挑也还好。记录这篇，是因为网页做了几次，想记录下来留个印象。

? ??做富集分析，对于输出结果，不能只看前多少个有显著富集的term，应该要看符合自己设的padjust阈值的所有term，看完再挑选自己想要展示的term成图。这样做能避免得出的结论不全面，对于事先没有预想的term或者是事先预想的term不全面这些情况有帮助。

????clusterProfiler做富集，可以对冗余 GOterm 去冗余，如：merge_go <- clusterProfiler::simplify(go, cutoff=0.7, by="p.adjust", select_fun=min)。要不要去冗余看去冗余前后展示出的效果是不是自己想要的，来决定。

????动态GO富集分析：富集结果统计图、条形图、气泡图，可动态调整? ??

????https://www.omicshare.com/tools/home/report/goenrich.html ? ?

????GO富集分析高级版：富集结果统计图、条形图、气泡图、富集圈图、富集差异气泡图、有向无环网络图，不可动态调整 ? ? ????https://www.omicshare.com/tools/Home/Soft/gogseasenior

????动态KEGG富集分析：富集结果统计图、条形图、气泡图，可动态调整? ??????

????https://www.omicshare.com/tools/home/report/koenrich.html

????KEGG富集分析高级版： 富集结果统计图、条形图、气泡图、富集圈图、富集差异气泡图、kegg网络图，不可动态调整?? ? ???? https://www.omicshare.com/tools/Home/Soft/pathwaygseasenior

????GO总共有三个ontology（本体），分别描述基因的 分子功能（molecular function）、细胞组分（cellular component）、参与的生物过程（biological process） 。GO的基本单位是term（词条、节点），每个term都对应一个属性。

???? x轴是Rich Factor，表示目的基因富集到该通路的基因数目与背景基因富集到该通路的基因数目的比值，所以比值越大，富集到该通路的基因数目越多；

????y轴是富集出来的通路名称，一般根据P-value或Q-value从小到大排列顺序（最上面是最小的），挑选富集通路前20或30的通路来绘图；点的大小表示Gene数目，点越大，表示富集到该通路的基因越多；点的颜色渐变最为重要，代表P值的高低，-log10(Pvalue)越大，P值越小，表示该通路越显著。

???? P值是在进行富集分析时利用超几何检验计算出来的结果。Q值是计算得到的P值进一步经过多重检验校正后的值。所以一般情况下Q值比P值的检验更严格。 这时候也可以分两种情况，第一种是经过Q值的检验筛选后得到了与实验设计相关的通路，那此时就可以用Q值来绘图。第二种是经过Q值检验没有得到满意的结果，那可以换P值检验，得到与实验相关的通路。

????P-value是正常数值显示还是取10的对数的负值来显示，主要看P-value的大小，如果P-value非常小时，就可以对数据归一化处理；如果数据较大，但又小于0.05时，可以不用对数据做处理。

? ?? 圆圈的大小代表基因的数目，圆圈的颜色代表P-value，也就是说 Rich Factor越大，P-value越小，gene count圈越大，这事就越可信。

三个Ontology（C, F, P）会分别展示。以生物过程（biological process）为例子，如下表：

第一列为GO term的ID，点击GO ID，可显示这个GO term包含的所有基因：

????再点击这个GO ID，就可以链接到 http://amigo.geneontology.org 官网，可以查看GO的具体信息。

????第二列为GO term的功能描述；

???? 第三列：数字为目的基因中富集到这个GO term的基因数，out 括号里数字为目的基因在 BP/MF/CC 里的总数，括号内百分比等于Generatio。

????第四列：数字为背景基因中富集到这个GO term的基因数，All括号里数字为背景基因在 BP/MF/CC 里的总数，百分比为Bgratio

????第五列：P value，即第三列的百分比与第四列的百分比相比，是否有显著差异。将小于0.05的P value标红显示； 这些GO term是按照P value从小到大排列的，方便找差异富集结果。即最上面的 term 为在目的基因中富集最显著的GO term，说明目的基因显著富集于这个功能。

????包含了pathway名称、 目的基因中注释到该pathway的基因数目以及占总目的基因数目的比例 、 所有背景基因中注释到该pathway的基因数目以及占总背景基因数目的比例 、P值、Q值、pathway 的ID（ko号）。点击pathway名称可以查看该pathway包含的基因ID，继续点pathway名称可以链接到KEGG官网上pathway相应的通路图

????存放的就是每个pathway的map图和相应的KEGG官网链接。

????Pathway的B级分类基因注释数目的统计图：纵坐标黑色字体为A级分类名称，彩色字体为B级分类名称。横坐标表示注释到相应B级pathway的基因数目。

? ? 显著富集pathway去前多少个根据 富集结果表（out.htm）再调整，用p值比用q值条件要松，颜色代表p值，柱状长短（气泡大小）代表富集到此通路基因的数量，横坐标RichFactor代表目的基因中位于该pathway条目的基因数目与背景基因中位于该pathway条目的基因总数的比值，比值越大，富集程度则越大。

参考：

https://www.omicshare.com/forum/thread-826-1-1.html

https://www.omicshare.com/forum/thread-6821-1-1.html

https://www.omicshare.com/forum/thread-6822-1-1.html

https://www.omicshare.com/forum/forum.php?mod=viewthread&tid=7295&highlight=KEGG

https://zhuanlan.zhihu.com/p/35065777

https://www.jieandze1314.com/post/cnposts/164/

R安装clusterProfiler 包遇到无法安装依赖包stringi

install.packages("BiocManager")

BiocManager::install("clusterProfiler")

stringi 报错，本地安装

wget https://github.com/gagolews/stringi/archive/master.zip -O stringi.zip

unzip stringi.zip

sed -i '/\/icu..\/data/d' stringi-master/.Rbuildignore

R CMD build stringi-master

R CMD INSTALL?stringi_1.7.5.9001.tar.gz

BiocManager::install("clusterProfiler")

如何将clusterProfiler富集结果中entrezID转化为基因名

注意：使用前提是你的富集结果对象类型必须是：'enrichResult' , 'gseaResult' or 'compareClusterResult'其中的一个，也就是用clusterProfiler里的富集函数直接给出的结果

注意：使用前提是你的富集结果对象类型必须是：'enrichResult' , 'gseaResult' or 'compareClusterResult'其中的一个，也就是用clusterProfiler里的富集函数直接给出的结果

小鼠clusterProfiler：GO KEGG

##################################################富集分析

library(clusterProfiler)

#人种 org.Hs.eg.db 小鼠 org.Mm.eg.db?

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

? install.packages("BiocManager")

BiocManager::install("org.Mm.eg.db")

library(org.Mm.eg.db)

##转换ID

trans <- bitr(gene, fromType="SYMBOL",

? ? ? ? ? ? ? toType="ENTREZID", OrgDb="org.Mm.eg.db")

#######################GO

ego <- enrichGO(gene? ? ? ? ? = trans$ENTREZID,

? ? ? ? ? ? ? ? OrgDb? ? ? ? ?= org.Mm.eg.db,

? ? ? ? ? ? ? ? ont? ? ? ? ? ?= "ALL",

? ? ? ? ? ? ? ? pAdjustMethod = "BH",

? ? ? ? ? ? ? ? pvalueCutoff? = 0.05,

? ? ? ? ? ? ? ? qvalueCutoff? = 1,

? ? ? ? ? ? ? ? readable? ? ? = TRUE)

head(ego@result)

result<-ego@result

write.csv(result,file = 'go.result.csv',quote = F,row.names = T,col.names = T)

##########################KEGG

#organism supported organism listed in ' http://www.genome.jp/kegg/catalog/org_list.html '

#hsa 人；mmu 小鼠

kk <- enrichKEGG(gene? ? ? ? ?= trans$ENTREZID,

? ? ? ? ? ? ? ? ?organism? ? ?= 'mmu',

? ? ? ? ? ? ? ? ?pvalueCutoff = 0.05,

? ? ? ? ? ? ? ? ?qvalueCutoff = 1)

tmp <- data.frame(kk)

kkx <- setReadable(kk, 'org.Mm.eg.db', 'ENTREZID')

kk<- data.frame(kkx)

write.csv(kk,file = 'kegg.result.csv',quote = F,row.names = T,col.names = T)

蛋白ID转基因ID

将Ensembl 中的蛋白ID转化成基因ID,可以通过clusterProfiler这个包。
如以大鼠的基因与蛋白转化为例;
安装clusterProfiler与大鼠org.Rn.eg.db，如果是人的注释包为org.Hs.eg.db,小鼠的注释包为org.Mm.eg.db.

查看可以转化的ID:

[1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS"
[6] "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME"
[11] "GO" "GOALL" "IPI" "ONTOLOGY" "ONTOLOGYALL"
[16] "PATH" "PFAM" "PMID" "PROSITE" "REFSEQ"
[21] "SYMBOL" "UNIGENE" "UNIPROT"

将蛋白ID转为基因ID: