使用clusterProfiler对非模式生物进行富集分析

最近，小编有很多同学问我，非模式生物如何做富集分析？

小编本身是做小麦的，也属于非模式生物的范畴。以前的话，非模式生物要用blast2go跑电子注释，而blast2go又需要使用MySQL，没有root权限的话非常麻烦。所以非模式生物如何做富集分析也困扰了小编很久，直到有一天，小编发现了Y叔的神包“ clusterProfiler ”！可以轻松做富集分析！

非模式生物的话，分为两种，一种是可以在AnnotationHub上在线抓取Org.Db的非模式生物，另一种是在AnnotationHub上没有Org.Db的生物。

下面我们先来讲讲可以在AnnotationHub上抓取到Org.Db的非模式生物如何做富集分析：

# 载入包
library("AnnotationHub")
library("biomaRt")
library("clusterProfiler")
library("topGO")
library("Rgraphviz")
library("pathview")
# 导入文件
data <- read.table(file,header = F,sep="t")
# 搜寻 OrgDb（比如我想查找物种 Theobroma cacao）
hub <- AnnotationHub::AnnotationHub()
query(hub,"Theobroma cacao")

我们会搜索到Org.Db“AH59191”。

# 抓取 OrgDb
Thecacao.OrgDb <- hub[["AH59191"]]
# 查看 OrgDb 里的Gene ID类型
columns(Thecacao.OrgDb)
# 查看Gene ID 类型
columns(Thecacao.OrgDb)
# ID转换（将我们的ID转换成与Org.Db一致的ID类型）
data <- as.character(data$V1)
data_id <- mapIds(x = Thecacao.OrgDb,keys = data,keytype = "SYMBOL",column = "ENTREZID")
# 去除NA值
na.omit(data_id)
# 富集（其他富集函数参见R包文档）
erich.go.BP <- enrichGO(gene=data,OrgDb = Thecacao.OrgDb,keyType = "SYMBOL",ont = "BP",pvalueCutoff = 0.01,qvalueCutoff = 0.05,readable = T)
# 绘制条形图
barplot(erich.go.BP)
# 绘制气泡图
dotplot(erich.go.BP)

以上部分，就是可以抓取到Org.Db的物种进行富集分析的步骤。

但是，对于抓取不到的物种，我们也有方式做富集分析。以小麦GO富集分析为例：

首先在IWGSC的官网上下载最新的注释文件：https://wheat-urgi.versailles.inra.fr/Seq-Repository/Annotations

得到注释文件后，我们就可以提取到每个基因对应的GO号，这时，我们可以准备三个文件。

第一个文件：所要富集基因的列表（格式如下）：

第二个文件，基因与GO号的对应关系（格式如下）：

小编之前整理的是csv文件，所以这里是逗号分隔，大家也可以用TAB分割，但是两列的顺序一定不要调换！

第三个文件，GO号与其对应注释的文件（格式如下）

有了这三个文件，我们就可以开始富集分析啦！

#载入包
library("clusterProfiler")
# 导入基因列表
gene <- read.csv(file,header = T,sep=",")
gene <- as.factor(gene$V1)
# 导入注释文件
term2gene <- read.csv(file,header=F,sep=",")
term2name <- read.csv(file,header=F,sep=",")
# 富集分析
x <- enricher(gene,TERM2GENE=term2gene,TERM2NAME=term2name,pvalueCutoff = 0.01, pAdjustMethod = "BH", qvalueCutoff = 0.05)
# 设置结果输出文件
ouf <- paste(out_file,sep ="t")
# 输出结果
write.csv(x,ouf)
# 绘制条形图
barplot(x)
# 绘制气泡图
dotplot(x)

有时，我们想做两列基因的富集结果比较时，可以用以下方式展现：

首先，我们要生成需要进行比较的两个基因列表（格式如下）：

# 导入文件
gene <- read.csv(file,header = T,sep=",")
# 富集分析
x <- compareCluster(gene, fun='enricher',TERM2GENE=term2gene,TERM2NAME=term2name,pvalueCutoff = 0.01, pAdjustMethod = "BH", qvalueCutoff = 0.05)
# 绘制气泡图
dotplot(x, showCategory=10,includeAll=TRUE)

这样，富集分析就完成啦！

当然，还有很多函数和细节，小编没有细讲，有兴趣的朋友可以查看“clusterProfiler”的文档：

https://www.bioconductor.org/packages/release/bioc/vignettes/clusterProfiler/inst/doc/clusterProfiler.html

参考资料：

https://www.bioconductor.org/packages/release/bioc/vignettes/clusterProfiler/inst/doc/clusterProfiler.html

https://guangchuangyu.github.io/cn/2017/07/clusterprofiler-maize/

https://blog.csdn.net/msw521sg/article/details/76167457

https://www.bilibili.com/video/av18311108/