clusterProfiler基因功能富集分析 +气泡图

1. 非模式生物

非模式生物需要提前从AnnotationHub抓取目标物种的OrgDb数据库，这个库里包括Uniport，ENSEMBL，ENTREZID，REFSEQ，GO...。拥有它，你就可以实现GENENAME，ENSEMBL，ENTREZID等多种格式的自由转换。不过这个功能一些在线网站也能做到，如果只有少量数据转换，DAVID/gProfiler还是挺香的。下面以油菜(Brassica_napus)为例介绍OrgDb：

#R script
setwd("C:\\Users\\Desktop\\Transcriptome\\class\\class3")

####非模式生物通过AnnotationHub在线检索并抓取OrgDb
BiocManager::install("AnnotationHub")
library(AnnotationHub)
hub<-AnnotationHub::AnnotationHub()
query(hub,"Brassica") #抓取包含Brassica字符串的物种OrgDb包，并获取其编号
bna <- hub[["AH75769"]] 
#玉米query(hub, "zea");maize <- hub[['AH55736']]

若你想深入了解OrgDb数据库，记住5个函数function():

  • 1 columns() 显示括号内对象拥有的数据名称, 每一类数据为key
  • 2 keys() 返回当前对象的keys
  • 2 keytypes() 字义理解keys的种类，用在select和mapIds中作参数转化
  • 3 select() 输入一组key，对应其column，可将其转化为目标keytype。输出的数据框同时包含转化前后数据。
  • 4 mapId() 输入输入一组key，对应其keytype，可将其转化为目标column中的key。输出的数据只转化后数据(即column)。
    这样的文字介绍还是抽象，不好理解。上机操作：
columns(bna)
# [1] "ACCNUM"  "ALIAS"  "CHR"  "ENTREZID" "EVIDENCE"   
# [6] "EVIDENCEALL" "GENENAME"   "GID"  "GO"  "GOALL"      
# [11] "ONTOLOGY"    "ONTOLOGYALL" "PMID"  "REFSEQ" "SYMBOL"     
# [16] "UNIGENE" 
keytypes(bna)
# [1] "ACCNUM"      "ALIAS"       "ENTREZID"    "EVIDENCE"    "EVIDENCEALL"
# [6] "GENENAME"    "GID"         "GO"          "GOALL"       "ONTOLOGY"   
# [11] "ONTOLOGYALL" "PMID"        "REFSEQ"      "SYMBOL"      "UNIGENE" 
head(keys(bna,keytype="SYMBOL"))
#[1] "BNAA01G00220D" "BNAA01G00960D" "BNAA01G01590D" "BNAA01G01860D"
#[5] "BNAA01G01890D" "BNAA01G02080D"
genename<-select(bna,keys=c("BNAA01G00220D","BNAA01G02080D"),column=c("GENENAME"),keytype = "SYMBOL")
head(genename)
#         SYMBOL                      GENENAME
# 1 BNAA01G00220D uncharacterized BNAA01G00220D
# 2 BNAA01G02080D uncharacterized BNAA01G02080D
entrezid<-mapIds(bna,keys=c("BNAA01G00220D","BNAA01G02080D"), column = c("ENTREZID"), keytype = "SYMBOL")
head(entrezid)
# BNAA01G00220D BNAA01G02080D 
# "106357927"   "106370211"

模式动物基因功能KEGG/GO注释(ClusterProfiler)

模式动物有现成的OrgDb包，安装加载目标物种的R package就可以了。

####模式动物(人为例)使用clusterProfiler进行GO/KEGG分析
#Bioconductor 安装 clusterProfiler
BiocManager::install("clusterProfiler")
BiocManager::install("org.Hs.eg.db")
install.packages("org.Hs.eg.db")
#加载clusterProfiler
library(clusterProfiler)
library(ggplot2)
library(org.Hs.eg.db)

hsa<-org.Hs.eg.db
keytypes(hsa)
head(keys(hsa, keytype="SYMBOL"))

#读入基因数据
gene<-read.csv("DEGs_identified.csv")
head(gene)
# Ensembl        padj log2FoldChange
# 1 ENSG00000001084 1.06000e-32       1.371851
# 2 ENSG00000002587 5.12000e-57       2.516021
# 3 ENSG00000003096 9.67315e-04      -1.088586
# 4 ENSG00000003137 6.01000e-18       1.170100
# 5 ENSG00000003400 8.16000e-07      -1.092483
# 6 ENSG00000004487 1.04000e-24       1.061878
dim(gene)
# [1] 934   3

#筛选|FC|>=2且FDR<=0.05的为显著差异基因
gene_ensembl<-na.omit(gene[(abs(gene$log2FoldChange)>= 1)&(gene$padj<= 0.05),1])
head(gene_ensembl);length(gene_ensembl)
#[1] ENSG00000001084 ENSG00000002587 ENSG00000003096 ENSG00000003137 ENSG00000003400 ENSG00000004487
# 934 Levels: ENSG00000001084 ENSG00000002587 ENSG00000003096 ENSG00000003137 ENSG00000003400 ... ENSG00000276600
# [1] 934
#通过ENSEMBL，获取ENTREZID，GENENAME
genelist<- bitr(gene_ensembl,fromType="ENSEMBL",toType=c("ENTREZID","GENENAME"),OrgDb=hsa)
#'select()' returned 1:many mapping between keys and columns
#Warning message:
#In bitr(gene_ensembl, fromType = "ENSEMBL", toType = c("ENTREZID",  :
#0.86% of input gene IDs are fail to map...

值得一提的是，除了bitr()还有个bitr_kegg()函数，也是转换用的，不过看名字就知道，这个函数服务于kegg。bitr_kegg中的转化参数有“Path”, “Module”, “ncbi-proteinid”, “ncbi-geneid”, “uniprot”, “kegg”。提供其中一项可以互相转化，举个栗子🌰:

kegg<-bitr_kegg(genelist$ENTREZID,fromType="ncbi-geneid",toType="kegg",organism='hsa')
head(kegg)
#  ncbi-geneid      kegg
# 1   100131827 100131827
# 2   100132074 100132074

Path<-bitr_kegg(genelist$ENTREZID,fromType="ncbi-geneid",toType="Path",organism='hsa')
head(Path)
#

人类KEGG库中对应编码hsa, 如果你不知道目标物种对应的缩写organism=?，请点击kegg_organism_list。而'kegg'和‘Path‘都是取自KEGG PATHWAY数据库，kegg表示KEGG PATHWAY数据库中基因命名方式，大多数采用ENTREZID，也有用其他方式，比如拟南芥。如果你不确定目标物种的kegg命名，可以使用ncbi-geneid(ENTREZID)转化为kegg。而Path则为目标基因对应的kegg pathway的编号。

在使用bitr_kegg时常会报错，比如

\color{red} {Error in KEGG_ convert("kegg", keyType, species) : ncbi-geneid is not supported for hsa ...}

如果你确定物种缩写没有错，keyType使用正确，那就不是你的问题，可能是网络的问题，多跑几遍🤷‍♂️。

差异基因功能GO分析

GO/kegg分析的在线工具很多，比如gProfiler。但是如果需要批量处理几组差异基因分析，在线工具还是挺费事的。你没时间研究R的话，用网页点点点也是错不了。

GO(GENE ONTOLOGY)对功能的注释分为三大类:分别是BP，MF，和CC。通过对ont的选择(ont='BP','MF','CC')就能分门别类地分析，但如果你想一次解决，直接令ont="ALL"。

下面进入ClusterProfiler的高光时刻，走起：

go.all<- enrichGO(gene=genelist$ENTREZID, 
                  OrgDb = org.Hs.eg.db, 
                  ont='ALL', #ont='BP'
                  pAdjustMethod = 'BH',
                  pvalueCutoff = 0.05, 
                  qvalueCutoff = 0.2,
                  keyType = 'ENTREZID')
go.all
# over-representation test
#...@organism    Homo sapiens 
#...@ontology    GOALL 
#...@keytype     ENTREZID 
#...@gene        chr [1:929] "2729" "9957" "90293" "56603" "843" "23028" "115703" "26224" "114881" "51330" "80853" "55200" ...
#...pvalues adjusted by 'BH' with cutoff <0.05 
#...627 enriched terms found
#'data.frame':   627 obs. of  10 variables:
# $ ONTOLOGY   : Factor w/ 2 levels "BP","MF": 1 1 1 1 1 1 1 1 1 1 ...
# $ ID         : chr  "GO:0001763" "GO:0061138" "GO:0045766" "GO:0048754" ...
# $ Description: chr  "morphogenesis of a branching structure" "morphogenesis of a branching epithelium" "positive regulation of angiogenesis" "branching morphogenesis of an epithelial tube" ...
# $ GeneRatio  : chr  "33/832" "30/832" "31/832" "26/832" ...
# $ BgRatio    : chr  "196/18670" "182/18670" "204/18670" "150/18670" ...
# $ pvalue     : num  4.02e-11 5.27e-10 2.13e-09 2.49e-09 2.89e-09 ...
# $ p.adjust   : num  2.10e-07 1.38e-06 3.03e-06 3.03e-06 3.03e-06 ...
# $ qvalue     : num  1.62e-07 1.06e-06 2.33e-06 2.33e-06 2.33e-06 ...
# $ geneID     : chr  "54903/6591/639/59277/2294/2138/374/7422/7475/1746/5228/53335/157506/6662/650/54567/10253/2247/1969/9353/170690/"| __truncated__ "54903/6591/59277/2294/2138/374/7422/7475/5228/157506/6662/650/54567/10253/2247/1969/9353/170690/133/10481/57216"| __truncated__ "9734/5743/23411/3162/3696/5054/7422/3552/2152/5228/3553/79924/694/2113/9314/9982/7057/1545/2247/133/7424/51738/"| __truncated__ "54903/2294/2138/374/7422/7475/5228/157506/6662/650/54567/10253/2247/1969/9353/170690/57216/64783/4824/2637/2303"| __truncated__ ...
# $ Count      : int  33 30 31 26 33 44 48 47 44 21 ...

#随后对富集结果进行总览，查看BP，CC，MF的个数
dim(go.all[go.all$ONTOLOGY=='BP',]);dim(go.all[go.all$ONTOLOGY=='CC',]);dim(go.all[go.all$ONTOLOGY=='MF',])
#保存结果
write.csv(go.all@result,'DEG_go.all.result.csv',row.names=F)

用enrichplot画个图看看：

dotplot(go.all,showCategory = 10, size = NULL,font.size = 10, title = "GO enrichment", split = "ONTOLOGY") + facet_grid(ONTOLOGY ~ ., scales = "free")
barplot(go.all,showCategory = 10, size = NULL,font.size = 10, title = "GO enrichment", split = "ONTOLOGY") + facet_grid(ONTOLOGY ~ ., scales = "free")

如果要将ont='BP'筛选出来，并且对P值设置更为严格阈值，有两种方式可以直接筛选:1.直接筛选；2.利用clusterProfiler.dplyr。

##1. 直接筛选，MF没有富集到terms，不作演示
enrich.BP<-go.all[go.all@result$pvalue<0.01&go.all$ONTOLOGY=='BP',]
enrich.CC<-go.all[go.all@result$pvalue<0.01&go.all$ONTOLOGY=='CC',]
dim(enrich.BP);dim(enrich.CC)

还可以安装个为cluster profiler服务的数据处理包clusterProfiler.dplyr。然后一条命令就解决了。

install.packages("devtools")
devtools::install_github("YuLab-SMU/clusterProfiler.dplyr")
BiocManager::install('clusterProfiler.dplyr')
library(clusterProfiler.dplyr)
enrich.BP<-filter(go.all,ONTOLOGY=='BP',p.adjust <.01, qvalue < 0.1)
enrich.BP@ontology<-"BP"
dim(enrich.BP)
BP.bar<-barplot(enrich.BP,showCategory=20);BP.dot<-dotplot(enrich.BP,showCategory=20)
BP.bar
BP.dot

通常得到的GOterms中有一些相似度较高的term，我们可以去除这些冗余terms，让结果更加简洁。

##去除冗余terms
go.filter<-simplify(enrich.BP,cutoff = 0.7,
                    by = "p.adjust",
                    select_fun = min)
dim(go.filter)
#将结果数据框中的ENTREZID转化为SYMBOL
y <-setReadable(go.all,OrgDb = org.Hs.eg.db,keyType="ENTREZID")
head(y)

接下来是KEGG，KEGG需要注意的是organism，keyType，use_internal_data设置。

enrich.kegg <- enrichKEGG(gene =genelist$ENTREZID,
                          organism ="hsa",
                          keyType = "kegg",
                          pvalueCutoff = 1,
                          pAdjustMethod = "BH",
                          minGSSize = 10,
                          maxGSSize = 500,
                          qvalueCutoff = 1,
                          use_internal_data =FALSE)
dim(enrich.kegg)

minGSSize：minimal size of genes annotated for testing，用于测试基因集的最小数目。
maxGSSize：maximal size of genes annotated for testing，用于测试的基因注释最大数目。
这两个参数常常被忽略。
use_internal_data : 是否使用内部数据，当use_internal_data =FALSE时，爬取在线KEGG PATHWAY数据库，特点就是实时更新，而内部数据就有滞后性了。

需要注意pvalue或p.adjust设置

sig.kegg<-filter(enrich.kegg,pvalue<.05,qvalue<0.2)
dim(sig.kegg)
kegg.bar<-barplot(sig.kegg,showCategory=20,color = "pvalue")
kegg.dot<-dotplot(sig.kegg,showCategory=20,color = "pvalue")
db<-plot_grid(kegg.bar,kegg.dot,ncol=2)
db
ggsave("kegg_bar_dot1.png",db,width=14,height=6)
kegg_bar_dot1.png
问题一：横坐标，气泡溢出
p1<-kegg.dot + xlim(NA,0.085) #数值>横坐标显示值
p1
db1<-plot_grid(kegg.bar,p1,ncol=2)
db1
ggsave("kegg_bar_dot2.png",db1,width=14,height=6)
kegg_bar_dot2.png
问题二：注释文字太长，截断成两行
library(stringr)
p2<- kegg.dot + scale_y_discrete(labels=function(y) stringr::str_wrap(y,width=36))
p2

#####问题三：改变颜色
p30 <-p2 + scale_color_continuous(low='purple',high='green')+ xlim(NA,0.084)
p30
p31 <-p2 + scale_color_gradientn(colours = rainbow(5))+ xlim(NA,0.084)
p31
p32 <-p2+scale_color_gradient(low="yellow", high="green")+ xlim(NA,0.084)
p32
db<-plot_grid(p30,p31,p32,ncol=3)
db
ggsave("kegg_bar_dot3.png",db,width=16,height=6)
kegg_bar_dot3.png
问题四：改变形状
p3 <-p31 + aes(shape=GeneRatio > 0.05) 
p3
问题五：气泡变大
p6<- p31 + scale_size(range=c(2,12))
p6
问题六：横坐标不是GeneRatio,而是richFactor或者FoldEnrichment
y1 <- mutate(sig.kegg, richFactor = Count / as.numeric(sub("/\\d+", "", BgRatio)))
#y2 <- mutate(sig.kegg, FoldEnrichment = parse_ratio(GeneRatio) / parse_ratio(BgRatio))
y1 

#若横轴是FoldEnrichment将y1替换成y2,
library(enrichplot)
library(forcats)
rf<-ggplot(y1, showCategory = 20, 
       aes(richFactor, fct_reorder(Description,richFactor))) + 
    geom_point(aes(color=pvalue, size = Count)) +
    scale_color_viridis_c(guide=guide_colorbar(reverse=TRUE)) +
    scale_size_continuous(range=c(2,10)) +
    #theme_minimal() + 
    xlab("Rich Factor") +
    ylab(NULL) + 
    ggtitle("Top 20 KEGG Enrichment")+
    theme_bw()
ggsave("kegg_bar_dot4.png",rf,width=7,height=5)
kegg_bar_dot4.png

如果你还想了解更多，点击下方链接，直接进入clusterProfiler官方解说
clusterProfiler 官方document
一本关于clusterProfiler的小册子 by Y叔

如果您有任何需要的组学服务，请关注伟寰生物公众号，扫码联系我们，谢谢支持！


伟寰生物.jpg
最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 206,214评论 6赞 481
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 88,307评论 2赞 382
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 152,543评论 0赞 341
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 55,221评论 1赞 279
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 64,224评论 5赞 371
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 49,007评论 1赞 284
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 38,313评论 3赞 399
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,956评论 0赞 259
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 43,441评论 1赞 300
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,925评论 2赞 323
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,018评论 1赞 333
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,685评论 4赞 322
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 39,234评论 3赞 307
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,240评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,464评论 1赞 261
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 45,467评论 2赞 352
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,762评论 2赞 345