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

作者: 黄甫一 | 来源:发表于2023-09-11 16:30 被阅读0次

背景

做单细胞转录组分析时，通过找差异基因可以得到很多基因集，一方面我们需要看这些基因集的相对表达量是否充足，但我们往往更关注这些差异基因（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富集

总结与讨论

暂时没有，以后更新。

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

背景

一些注意事项

基础依赖包

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

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

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

总结与讨论

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读