使用clusterProfiler进行基因集富分析两个函数就够了
2023-09-11 本文已影响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富集
总结与讨论
暂时没有，以后更新。
