【R】clusterProfiler的GO/KEGG富集分析用法小结 – 小xuo生 – 博客园

【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/,很乱,大多数我都不知道干啥用的。

packages organism
org.Ag.eg.db Anopheles
org.At.tair.db Arabidopsis
org.Bt.eg.db Bovine
org.Ce.eg.db Worm
org.Cf.eg.db Canine
org.Dm.eg.db Fly
org.Dr.eg.db Zebrafish
org.EcK12.eg.db E coli strain K12
org.EcSakai.eg.db E coli strain Sakai
org.Gg.eg.db Chicken
org.Hs.eg.db Human
org.Mm.eg.db Mouse
org.Mmu.eg.db Rhesus
org.Pf.plasmo.db Malaria
org.Pt.eg.db Chimp
org.Rn.eg.db Rat
org.Sc.sgd.db Yeast
org.Ss.eg.db Pig
org.Xl.eg.db Xenopus

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注释文件的格式如下:

GeneID GO GO_Description
1 GO:0005819 spindle
2 GO:0072686 mitotic spindle
3 GO:0000776 kinetochore

需要包含以上三列信息,这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注释文件,然后进行分析,注释文件的格式如下:

GeneID Pathway Path_Description
1 ko:00001 spindle
2 ko:00002 mitotic spindle
3 ko:00003 kinetochore

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

富集分析

通常用的富集分析有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

点击数:0

发表评论

邮箱地址不会被公开。 必填项已用*标注