使用clusterProfiler进行基因集富分析两个函数就够了 2023-09-12

背景

做单细胞转录组分析时，通过找差异基因可以得到很多基因集，一方面我们需要看这些基因集的相对表达量是否显著上/下调，但我们往往更关注这些差异基因（DEGs）涉及的相关功能是否能对上注释好的细胞类型。因此我们需要进行功能验证，在没法进行湿实验的情况下，我们可以做的是就是基因富集分析了。根据选取的数据库不同，可以分为GO、KEGG和DO等等。clusterProfiler包已经非常方便，但为了更方便进行多种类型的富集分析，我根据官网教程最终整合成了两个函数，可以快速出图。

一些注意事项

1、org.*.eg.db系列包查询这个网址，人类的是org.Hs.eg.db，小鼠是org.Mm.eg.db。非模式物种参考这个网址自行构建.
2、KEGG数据库支持的物种使用search_kegg_organism('ece', by='kegg_code')查询，人类的是hsa，小鼠的是mmu。
3、KEGG第一次使用需要联网，设置use_internal_data=F，之后可以设置use_internal_data=T，更快进行分析。
4、treeplot在旧版本中并不支持，建议更新clusterProfiler到最新版本。
5、转换基因ID使用 bitr函数，例如test = bitr(gene, fromType="SYMBOL", toType=c("ENTREZID"), OrgDb=org.Hs.eg.db)，一般转换为ENTREZID。

基础依赖包

library(clusterProfiler)
library(org.Mm.eg.db)
library(org.Hs.eg.db)
library(ggplot2)
library(ggrepel)
library(ggpubr)
library(RColorBrewer)
library(stringr)
library(cowplot)
library(DOSE)
library(enrichplot)

单个基因集进行富集分析函数enrich_result

只有单个基因集的富集分析，输入为基因集向量，在官网使用enrich系列函数，这里则整合为一个函数，我们先看一下能输出的11张图，看看有没有你想要的。

1气泡图

2条形图1

3条形图2-计算了横轴qscore

4网络图1

5网络图2

6热图

7网络图3

8upsetplot图

9聚类树图1-使用原始值

10聚类树图2-使用平均值

11goplot-仅限GO富集分析

下面是主函数enrich_result的代码，
vgene是输入的基因集向量，
p.val=0.05是多重假设检验显著性阈值，
OrgDb='org.Hs.eg.db'是对应物种的org..eg.db系列包名称，
label='out'是输出文件前缀，
keyType='ENTREZID'是输入基因ID的类型，
colours = c('#336699','#66CC66','#FFCC33')是画图的色板
pAdjustMethod='BH'是矫正p值的方法，
fun= "GO"是进行富集分析的函数，可选GO、KEGG、DO、enricher等，
q.val=0.2是q值的阈值，
ont = "BP"是GO富集分析选择的类别，
showCategory = 10是展示通路的个数，
organism = "hsa"是KEGG分析的物种缩写，
use_internal_data=T是KEGG分析时是否使用内置数据（第一次跑需要联网），
minGSSize= 5和maxGSSize= 500是基因集大小的下限和上限，
categorySize="pvalue"是画图区分点大小的值，
foldChange=NULL是区分热图颜色深浅的表达量差异倍数向量，
node_label="all"是展示点的名称，可选只展示基因或者通路，
color_category='firebrick'是通路点的颜色，
color_gene='steelblue'是基因点的颜色，
interm=NULL是自定义基因集类型数据库的数据框，后面会细讲，
wid=18,hei=10是输出图片的宽和高，可以修改。

#画图函数
enrich_plot <- function(eobj,label='out',colours = c('#336699','#66CC66','#FFCC33'),
                       fun= "GO", showCategory = 10,
                       categorySize="pvalue",foldChange=NULL,node_label="all",color_category='firebrick',color_gene='steelblue',
                       cex_category=1.5,layout="kk",wid=18,hei=10) {
pdf(paste0(label,"_enrich_",fun,"_plot.pdf"),wid,hei)

ttl <- paste0(label,"_",fun)
gttl <-ggtitle(ttl)

p1 <- dotplot(eobj,showCategory = showCategory,title=ttl)+ scale_color_gradientn(values = seq(0,1,0.2),colours = colours)
p2 <- barplot(eobj, showCategory=showCategory, title=ttl)+ scale_fill_gradientn(values = seq(0,1,0.2),colours = colours)
p3 <- mutate(eobj, qscore = -log(p.adjust, base=10)) %>% barplot(x="qscore",showCategory=showCategory, title=ttl)+ scale_fill_gradientn(values = seq(0,1,0.2),colours = colours)
p4 <- cnetplot(eobj, categorySize=categorySize, foldChange=foldChange,node_label=node_label,color_category=color_category,color_gene=color_gene)+ gttl
p5 <- cnetplot(eobj, foldChange=foldChange, circular = TRUE, colorEdge = TRUE,node_label=node_label, color_category=color_category,color_gene=color_gene)+ gttl
p6 <- heatplot(eobj, foldChange=foldChange, showCategory=showCategory) + scale_color_gradientn(values = seq(0,1,0.2),colours = colours)+gttl

eobj1 <- pairwise_termsim(eobj)

p9 <- emapplot(eobj1, cex_category=cex_category,layout=layout) + gttl+ scale_color_gradientn(values = seq(0,1,0.2),colours = colours)
p10 <- upsetplot(eobj)+ gttl


print(p1)
print(p2)
print(p3)
print(p4)
print(p5)
print(p6)
try(print(p9))
print(p10)

try({
if (exists('treeplot')) {
p7 <- treeplot(eobj1)
p8 <- treeplot(eobj1, hclust_method = "average")
print(p7)
print(p8)
}
})
dev.off()
}
#主函数
enrich_result <- function(vgene,p.val=0.05,OrgDb='org.Hs.eg.db',label='out',
                          keyType='ENTREZID',colours = c('#336699','#66CC66','#FFCC33'), pAdjustMethod='BH',
                       fun= "GO", q.val=0.2, ont = "BP", showCategory = 10,organism = "hsa",use_internal_data=T,
                       minGSSize= 5,maxGSSize= 500,
                       categorySize="pvalue",foldChange=NULL,node_label="all",color_category='firebrick',color_gene='steelblue',interm=NULL,
                       wid=18,hei=10){

fun.use=paste0('enrich',fun)

if (fun=="GO"){
eobj <- enrichGO(gene         = vgene,
                OrgDb         = OrgDb,
                keyType       = keyType,
                ont           = ont,
                pAdjustMethod = pAdjustMethod,
                pvalueCutoff  = p.val,
                qvalueCutoff  = q.val,
                readable=TRUE)
p11 <- goplot(eobj)
png(paste0(label,"_GO_goplot.png"),1800,1000)
print(p11)
dev.off()
}

if (fun=="KEGG"){
kk <- enrichKEGG(gene         = vgene,
                 organism     = organism,
                 pvalueCutoff = p.val,
                 use_internal_data=use_internal_data)

eobj <- setReadable(kk,OrgDb=OrgDb,keyType=keyType)

}
if (fun=="DO"){
eobj <- enrichDO(gene       = vgene,
              ont           = "DO",
              pvalueCutoff  = p.val,
              pAdjustMethod = pAdjustMethod,
              minGSSize     = minGSSize,
              maxGSSize     = maxGSSize,
              qvalueCutoff  = q.val,
              readable      = TRUE)

}

if (fun=="NCG"){
eobj <- enrichNCG(gene      = vgene,
              pvalueCutoff  = p.val,
              pAdjustMethod = pAdjustMethod,
              minGSSize     = minGSSize,
              maxGSSize     = maxGSSize,
              qvalueCutoff  = q.val,
              readable      = TRUE)


}

if (fun=="DGN"){
eobj <- enrichDGN(gene      = vgene,
              pvalueCutoff  = p.val,
              pAdjustMethod = pAdjustMethod,
              minGSSize     = minGSSize,
              maxGSSize     = maxGSSize,
              qvalueCutoff  = q.val,
              readable      = TRUE)

}
if (fun=="enricher"){

x <- enricher(gene      = vgene,
              pvalueCutoff  = p.val,
              pAdjustMethod = pAdjustMethod,
              minGSSize     = minGSSize,
              maxGSSize     = maxGSSize,
              qvalueCutoff  = q.val,
              TERM2GENE = interm)

eobj <- setReadable(x,OrgDb=OrgDb,keyType=keyType)
}

saveRDS(eobj, paste0(label,"_",fun,"_enrich.rds"))
out=eobj@result
write.table(out,paste0(label,"_enrich_",fun,"List.xls"),row.names = FALSE,quote = FALSE,sep = "\t")
enrich_plot(eobj=eobj,label=label,colours = colours,
                       fun= fun, showCategory = showCategory,
                       categorySize=categorySize,foldChange=foldChange,node_label=node_label,color_category=color_category,color_gene=color_gene,
                       wid=wid,hei=hei)
return(eobj)
}

使用示例
加载数据

library(DOSE)
data(geneList, package="DOSE")
gene <- names(geneList)[abs(geneList) > 2]
head(gene)
#[1] "4312"  "8318"  "10874" "55143" "55388" "991"

绝大多数都可以使用默认参数，只需改变fun参数，最终生成三个文件，以GO为例，会生成out_enrich_GOList.xls、out_enrich_GO_plot.pdf和out_GO_enrich.rds三个文件，其中
out_enrich_GO_plot.pdf是输出的图片，
out_GO_enrich.rds是enrich对象，
out_enrich_GOList.xls则是可以直接查看的数据框文件。

#GO
ob1 <- enrich_result(gene,fun='GO',label='out')
#KEGG
ob1 <- enrich_result(gene,fun='KEGG',label='out')
#DO
ob1 <- enrich_result(gene,fun='DO',label='out')

多个基因集进行富集分析函数compare_result

多个基因集富集分析使用compareCluster函数，老规矩，先看能生成的4张图片。

1气泡图

2网络图1

3网络图2

4网络图3

compare_result和前面的enrich_result函数几乎是一样的，只是compare_result输入的是多个基因集的列表，而enrich_result的输入是单个基因集向量。相关参数，这里不再赘述。下面是compare_result的代码：

#作图函数
compare_plot <- function(eobj,label='out',colours = c('#336699','#66CC66','#FFCC33'),
                       fun= "GO", showCategory = 10,
                       categorySize="pvalue",foldChange=NULL,node_label="all",color_category='firebrick',color_gene='steelblue',
                       legend_n=2,inpie="count", cex_category=1.5, layout="kk",wid=18,hei=10) {
pdf(paste0(label,"_comparecluster_",fun,"_plot.pdf"),wid,hei)

ttl <- paste0(label,"_",fun)
gttl <-ggtitle(ttl)

p1 <- dotplot(eobj,showCategory = showCategory,title=ttl)+ scale_color_gradientn(values = seq(0,1,0.2),colours = colours)
p4 <- cnetplot(eobj, categorySize=categorySize, foldChange=foldChange,node_label=node_label,color_category=color_category,color_gene=color_gene)+ gttl
p5 <- cnetplot(eobj, foldChange=foldChange, circular = TRUE, colorEdge = TRUE,node_label=node_label, color_category=color_category,color_gene=color_gene)+ gttl

eobj1 <- pairwise_termsim(eobj)
p9 <- emapplot(eobj1, cex_category=cex_category,legend_n=legend_n,pie=inpie, layout=layout) + gttl+ scale_color_gradientn(values = seq(0,1,0.2),colours = colours)


print(p1)
print(p4)
print(p5)
try(print(p9))

dev.off()
}
#主函数
compare_result <- function(lgene,p.val=0.05,OrgDb='org.Hs.eg.db',label='out',
                          keyType='ENTREZID',colours = c('#336699','#66CC66','#FFCC33'), pAdjustMethod='BH',
                       fun= "GO", q.val=0.2, ont = "BP", showCategory = 10,organism = "hsa",use_internal_data=T,
                       categorySize="pvalue",foldChange=NULL,node_label="all",color_category='firebrick',color_gene='steelblue',
                       minGSSize= 5,maxGSSize= 500,interm=NULL,wid=18,hei=10){

fun.use=paste0('enrich',fun)

if (fun=="GO"){
eobj <- compareCluster(geneCluster   = lgene,
                            fun           = fun.use,
                            pvalueCutoff  = p.val,
                            pAdjustMethod = pAdjustMethod,
                            OrgDb = OrgDb,
                            ont = ont,
                            readable = TRUE)
}

if (fun=="KEGG"){
kk <- compareCluster(geneCluster = lgene,
                     fun = fun.use,
                     pvalueCutoff  = p.val,
                     pAdjustMethod = pAdjustMethod,
                     organism     = organism,
                     use_internal_data=use_internal_data
                     )

eobj <- setReadable(kk,OrgDb=OrgDb,keyType=keyType)

}
if (fun=="DO"){

eobj <- compareCluster(geneCluster = lgene,
              ont           = "DO",
              fun = fun.use,
              pvalueCutoff  = p.val,
              pAdjustMethod = pAdjustMethod,
              minGSSize     = minGSSize,
              maxGSSize     = maxGSSize,
              qvalueCutoff  = q.val,
              readable      = TRUE)

}

if (fun=="NCG"){
eobj <- compareCluster(geneCluster = lgene,
              fun = fun.use,
              pvalueCutoff  = p.val,
              pAdjustMethod = pAdjustMethod,
              minGSSize     = minGSSize,
              maxGSSize     = maxGSSize,
              qvalueCutoff  = q.val,
              readable      = TRUE)


}

if (fun=="DGN"){
eobj <- compareCluster(geneCluster = lgene,
              fun = fun.use,
              pvalueCutoff  = p.val,
              pAdjustMethod = pAdjustMethod,
              minGSSize     = minGSSize,
              maxGSSize     = maxGSSize,
              qvalueCutoff  = q.val,
              readable      = TRUE)

}
if (fun=="enricher"){
x <- compareCluster(geneCluster = lgene,
              fun = fun,
              pvalueCutoff  = p.val,
              pAdjustMethod = pAdjustMethod,
              minGSSize     = minGSSize,
              maxGSSize     = maxGSSize,
              qvalueCutoff  = q.val,
              TERM2GENE = interm)

eobj <- setReadable(x,OrgDb=OrgDb,keyType=keyType)
}

saveRDS(eobj, paste0(label,"_comparecluster_",fun,".rds"))

out=eobj@compareClusterResult
write.table(out,paste0(label,"_comparecluster_",fun,"List.xls"),row.names = FALSE,quote = FALSE,sep = "\t")

compare_plot(eobj=eobj,label=label,colours = colours,
                       fun= fun, showCategory = showCategory,
                       categorySize=categorySize,foldChange=foldChange,node_label=node_label,color_category=color_category,color_gene=color_gene,
                       wid=wid,hei=hei)
return(eobj)
}

使用示例，也是只需改fun参数
导入数据

> data(gcSample)
> str(gcSample)
List of 8
 $ X1: chr [1:216] "4597" "7111" "5266" "2175" ...
 $ X2: chr [1:805] "23450" "5160" "7126" "26118" ...
 $ X3: chr [1:392] "894" "7057" "22906" "3339" ...
 $ X4: chr [1:838] "5573" "7453" "5245" "23450" ...
 $ X5: chr [1:929] "5982" "7318" "6352" "2101" ...
 $ X6: chr [1:585] "5337" "9295" "4035" "811" ...
 $ X7: chr [1:582] "2621" "2665" "5690" "3608" ...
 $ X8: chr [1:237] "2665" "4735" "1327" "3192" ...

运行函数，也会生成三个文件

lgene <- gcSample
#GO
ob1 <- compare_result(lgene,fun='GO',label='test')
#KEGG
ob1 <- compare_result(lgene,fun='KEGG',label='test')
#DO
ob1 <- compare_result(lgene,fun='DO',label='test')

自定义基因集的涵义enricher函数

除了已有的数据框，可以自定义基因的涵义进行富集分析，例如可以自定义一个细胞类型的DEGs为一个小数据框然后进行富集分析，可以辅助进行细胞类型注释，这可以通过enricher函数，但前面的两个函数也已经包含了这个功能。
下面演示如何构建一个可用于富集分析的数据库:
首先在这个CellMarker下载细胞类型的markers列表，我下载的是Cell_marker_Human.xlsx，然后进行预处理，最终得到一个只有两列的tibble，第一列是基因注释信息，第二列是基因ID，相当于一个小的数据库。

library(readxl)
df1 <read_excel("Cell_marker_Human.xlsx")
df1 <- data.frame(df1)
cell_marker_data=df1
cell_marker_data$geneID <- cell_marker_data$GeneID
cells <- cell_marker_data %>%
    dplyr::select(cell_name, geneID) %>%
    dplyr::mutate(geneID = strsplit(geneID, ', ')) %>%
    tidyr::unnest()
head(cells)
# A tibble: 6 × 2
  cell_name       geneID
  <chr>           <chr>
1 Macrophage      10461
2 Macrophage      2215
3 Macrophage      4360
4 Macrophage      11326
5 Macrophage      9332
6 Brown adipocyte 2167

然后进行分析

x <- enricher(gene, TERM2GENE = cells)
x <- setReadable(x,OrgDb=OrgDb,keyType=keyType)

也可以使用上面两个函数

#单个基因集
ob1 <- enrich_result(gene,fun='enricher',interm=cells,label='term')
#多个基因集
ob1 <- compare_result(lgene,fun='enricher',interm=cells,label='term')

可以根据富集结果进行细胞类型注释。
以下面这个图为例，可以看到特定基因集合在对应细胞类型中的markers基因中富集。

细胞类型markers富集

总结与讨论

暂时没有，以后更新。

最后编辑于：2023.10.12 17:04:42

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