【R】clusterProfiler的GO/KEGG富集分析用法小结

前言

关于clusterProfiler这个R包就不介绍了，网红教授宣传得很成功，功能也比较强大，主要是做GO和KEGG的功能富集及其可视化。简单总结下用法，以后用时可直接找来用。

首先考虑一个问题：clusterProfiler做GO和KEGG富集分析的注释信息来自哪里？

GO的注释信息来自Bioconductor，提供了19个物种的org类型的GO注释信息，如下表所示。Bioconductor中更多的注释包可参考http://www.bioconductor.org/packages/release/data/annotation/，很乱，大多数我都不知道干啥用的。

KEGG的注释信息clusterProfiler通过KEGG 数据库的API来获取，https://www.kegg.jp/kegg/rest/keggapi.html。

首先是一个物种所有基因对应的pathway注释文件，比如人的：http://rest.kegg.jp/link/hsa/pathway。
其次还需要pathway对应的描述信息，比如人的：
http://rest.kegg.jp/list/pathway/hsa。

关于KEGG数据库全部的物种及其简写（三个字母）如下列表：
https://www.genome.jp/kegg/catalog/org_list.html。

因此对于以上已有pathway注释的物种，只需要将物种简写输入给clusterProfiler， 它会通过联网自动获取该物种的pathway注释信息。

以上都是有物种信息的情况，那么对于无物种信息的项目怎么办？

GO可以通过读取外部的GO注释文件进行分析。关于基因的GO注释，interproscan、eggnog-mapper和blas2go等软件都可以做，不过输出格式有些不同。clusterProfiler需要导入的GO注释文件的格式如下：

需要包含以上三列信息，这3列信息任意顺序都可。

clusterProfiler包只针对含有OrgDb对象，如果是公共数据库中有该物种注释信息，只是未制作成org.db数据库（标准注释库），则可以不需要从头注释，只需手动制作org.db数据库类型，完成后直接使用即可，代码如下：

source("https://bioconductor.org/biocLite.R")
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("AnnotationHub") # 一个包含大量注释信息的数据库，里面有很多物种及来源于很多数据库的注释信息。
BiocManager::install("biomaRt")

library(AnnotationHub)
library(biomaRt)

hub <- AnnotationHub() #建立AnnotationHub对象（视人品，网不行加载不了）
# unique(hub$species) #查看AnonotationHub里面物种
hub$species[which(hub$species=="Solanum")] #看AnonotationHub里是否包含想要的物种
# Solanum是番茄的拉丁名
query(hub, "Solanum")  #查看该物种信息
hub[hub$species=="Solanum" & hub$rdataclass == "OrgDb"] #OrgDb属于rdataclass中，因此查看下该物种有没有OrgDb
Solanum.OrgDb <- hub[["AH59087"]]#AH59087是番茄对应的编号
#制作为标准注释库，就可和模式生物一样使用了

同样地，对于pathway数据库中没有的物种，也支持读取基因的pathway注释文件，然后进行分析，注释文件的格式如下：

以上三列信息的顺序也是任意的。

富集分析

通常用的富集分析有ORA、FCS和拓扑三种方法。ORA简单来说就是超几何检验或Fisher精确检验，大同小异，都符合超几何检验，这也是目前用的最多的方法，优劣不谈。FCS的代表就是GSEA，即基因集富集分析，优劣亦不谈。clusterProfiler提供了这两种富集分析方法。
1. ORA（Over-Representation Analysis）
GO富集参考代码：

#标准富集分析
ego <- enrichGO(
          gene  = gene$entrzID,
          keyType = "ENTREZID",
          universe = names(geneList), #背景基因集，可省
          OrgDb   = org.Hs.eg.db,
          ont     = "CC",
          pAdjustMethod = "BH",
          pvalueCutoff  = 0.01,
          qvalueCutoff  = 0.05,
          readable      = TRUE)

#通过导入外部注释文件富集分析
data <- read.table("go_annotation.txt",header = T,sep = "\t")
go2gene <- data[, c(2, 1)]
go2name <- data[, c(2, 3)]
x <- enricher(gene,TERM2GENE = go2gene,TERM2NAME = go2name)

gene差异基因对应的向量;
keyType指定基因ID的类型，默认为ENTREZID, 可参考keytypes(org.Hs.eg.db)类型 ；
OrgDb指定该物种对应的org包的名字；
ont代表GO的3大类别，BP, CC, MF，也可是全部ALL；
pAdjustMethod指定多重假设检验矫正的方法，有“ holm”, “hochberg”, “hommel”, “bonferroni”, “BH”, “BY”, “fdr”, “none”中的一种；
cufoff指定对应的阈值；
readable=TRUE代表将基因ID转换为gene symbol。

KEGG Pathway富集参考代码：

#标准富集分析
ego <- enrichKEGG(
          gene = gene,
          keyType = "kegg",
          organism  = 'hsa',
          pvalueCutoff  = 0.05,
          pAdjustMethod  = "BH",
          qvalueCutoff  = 0.05
)

#通过外部导入注释文件富集
data <- read.table("pathway_annotation.txt",header = T,sep = "\t")
go2gene <- data[, c(2, 1)]
go2name <- data[, c(2, 3)]
x <- enricher(gene,TERM2GENE = go2gene,TERM2NAME = go2name)

默认基因ID为kegg gene id，也可以是ncbi-geneid, ncbi-proteinid, uniprot等。
organism物种对应的三字母缩写，其他参数同GO富集。ID转换函数：

library(clusterProfiler)
bitr_kegg("1",fromType = "kegg",toType = 'ncbi-proteinid',organism='hsa')

library(org.Hs.eg.db)
keytypes(org.Hs.eg.db) #支持的ID类型
bitr(gene, fromType = "ENTREZID", toType = c("ENSEMBL", "SYMBOL"), OrgDb = org.Hs.eg.db)

#以上看出ID转换输入时，可以向量的形式，也可以单列基因名list导入，也可以是内置数据
gene <- c("AASDH","ABCB11","ADAM12","ADAMTS16","ADAMTS18")
gene  <-  data$V1 #字符串

data(geneList, package="DOSE") #富集分析的背景基因集
gene <- names(geneList)[abs(geneList) > 2]

2. GSEA（Gene Set Enrichment Analysis）
GO富集参考代码：

#标准富集分析
ego <- gseGO(
      geneList  = geneList,
      OrgDb  = org.Hs.eg.db,
      ont  = "CC",
      nPerm  = 1000,  #置换检验的置换次数
      minGSSize  = 100,
      maxGSSize  = 500,
      pvalueCutoff = 0.05,
      verbose  = FALSE)

#通过导入外部注释文件富集分析参考代码：
data <- read.table("go_annotation.txt",header = T,sep = "\t")
go2gene <- data[, c(2, 1)]
go2name <- data[, c(2, 3)]
x <- GSEA(gene,TERM2GENE = go2gene,TERM2NAME = go2name)

KEGG Pathway富集参考代码：

#标准富集分析
kk <- gseKEGG(
  geneList  = gene,
  keyType  = 'kegg',
  organism = 'hsa',
  nPerm  = 1000,
  minGSSize = 10,
  maxGSSize = 500,
  pvalueCutoff = 0.05,
  pAdjustMethod     = "BH"
)

#通过外部导入注释文件富集
data <- read.table("pathway_annotation.txt",header = T,sep = "\t")
go2gene <- data[, c(2, 1)]
go2name <- data[, c(2, 3)]
x <- GSEA(gene,TERM2GENE = go2gene,TERM2NAME = go2name)

可视化

1.GO富集分析结果可视化

#barplot
barplot(ego, showCategory = 10) #默认展示显著富集的top10个，即p.adjust最小的10个

#dotplot
dotplot(ego, showCategory = 10)

#DAG有向无环图
plotGOgraph(ego)  #矩形代表富集到的top10个GO terms, 颜色从黄色过滤到红色，对应p值从大到小。

#igraph布局的DAG
goplot(ego)

#GO terms关系网络图（通过差异基因关联）
emapplot(ego, showCategory = 30)

#GO term与差异基因关系网络图
cnetplot(ego, showCategory = 5)

2.Pathway富集分析结果可视化

#barplot
barplot(kk, showCategory = 10)

#dotplot
dotplot(kk, showCategory = 10)

#pathway关系网络图（通过差异基因关联）
emapplot(kk,  showCategory = 30)

#pathway与差异基因关系网络图
cnetplot(kk, showCategory = 5)

#pathway映射
browseKEGG(kk, "hsa04934") #在pathway通路图上标记富集到的基因，会链接到KEGG官网

Ref：
https://blog.csdn.net/weixin_43569478/article/details/83744242
https://blog.csdn.net/weixin_43569478/article/details/83744384
https://www.jianshu.com/p/065d38c28e2d
https://www.jianshu.com/p/47b5ea646932
https://www.cnblogs.com/yatouhetademao/p/8018252.html
https://zhuanlan.zhihu.com/p/35510434

原文地址：https://www.cnblogs.com/jessepeng/p/12159139.html