背景
做單細(xì)胞轉(zhuǎn)錄組分析時叁巨，通過找差異基因可以得到很多基因集记盒，一方面我們需要看這些基因集的相對表達(dá)量是否顯著上/下調(diào)，但我們往往更關(guān)注這些差異基因（DEGs）涉及的相關(guān)功能是否能對上注釋好的細(xì)胞類型菊碟。因此我們需要進(jìn)行功能驗證撼短，在沒法進(jìn)行濕實驗的情況下筛严，我們可以做的是就是基因富集分析了。根據(jù)選取的數(shù)據(jù)庫不同忘巧，可以分為GO恒界、KEGG和DO等等。clusterProfiler包已經(jīng)非常方便砚嘴，但為了更方便進(jìn)行多種類型的富集分析十酣，我根據(jù)官網(wǎng)教程最終整合成了兩個函數(shù)，可以快速出圖际长。
一些注意事項
- 1耸采、org.*.eg.db系列包查詢這個網(wǎng)址，人類的是org.Hs.eg.db工育，小鼠是org.Mm.eg.db虾宇。非模式物種參考這個網(wǎng)址自行構(gòu)建.
- 2、KEGG數(shù)據(jù)庫支持的物種使用search_kegg_organism('ece', by='kegg_code')查詢如绸，人類的是hsa嘱朽，小鼠的是mmu。
- 3怔接、KEGG第一次使用需要聯(lián)網(wǎng)搪泳，設(shè)置use_internal_data=F，之后可以設(shè)置use_internal_data=T扼脐，更快進(jìn)行分析岸军。
- 4、treeplot在舊版本中并不支持，建議更新clusterProfiler到最新版本艰赞。
- 5佣谐、轉(zhuǎn)換基因ID使用 bitr函數(shù)，例如test = bitr(gene, fromType="SYMBOL", toType=c("ENTREZID"), OrgDb=org.Hs.eg.db)方妖，一般轉(zhuǎn)換為ENTREZID台谍。
基礎(chǔ)依賴包
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)
單個基因集進(jìn)行富集分析函數(shù)enrich_result
只有單個基因集的富集分析，輸入為基因集向量吁断，在官網(wǎng)使用enrich系列函數(shù)趁蕊，這里則整合為一個函數(shù)，我們先看一下能輸出的11張圖仔役，看看有沒有你想要的掷伙。
1氣泡圖
2條形圖1
3條形圖2-計算了橫軸qscore
4網(wǎng)絡(luò)圖1
5網(wǎng)絡(luò)圖2
6熱圖
7網(wǎng)絡(luò)圖3
8upsetplot圖
9聚類樹圖1-使用原始值
10聚類樹圖2-使用平均值
11goplot-僅限GO富集分析
下面是主函數(shù)enrich_result的代碼，
vgene是輸入的基因集向量又兵，
p.val=0.05是多重假設(shè)檢驗顯著性閾值任柜，
OrgDb='org.Hs.eg.db'是對應(yīng)物種的org..eg.db系列包名稱，
label='out'是輸出文件前綴沛厨，
keyType='ENTREZID'是輸入基因ID的類型宙地，
colours = c('#336699','#66CC66','#FFCC33')是畫圖的色板
pAdjustMethod='BH'是矯正p值的方法，
fun= "GO"是進(jìn)行富集分析的函數(shù)逆皮，可選GO宅粥、KEGG、DO电谣、enricher等秽梅，
q.val=0.2是q值的閾值，
ont = "BP"是GO富集分析選擇的類別剿牺，
showCategory = 10是展示通路的個數(shù)企垦，
organism = "hsa"是KEGG分析的物種縮寫，
use_internal_data=T是KEGG分析時是否使用內(nèi)置數(shù)據(jù)（第一次跑需要聯(lián)網(wǎng)）晒来，
minGSSize= 5和maxGSSize= 500是基因集大小的下限和上限钞诡，
categorySize="pvalue"是畫圖區(qū)分點大小的值，
foldChange=NULL是區(qū)分熱圖顏色深淺的表達(dá)量差異倍數(shù)向量湃崩，
node_label="all"是展示點的名稱荧降，可選只展示基因或者通路，
color_category='firebrick'是通路點的顏色竹习，
color_gene='steelblue'是基因點的顏色誊抛，
interm=NULL是自定義基因集類型數(shù)據(jù)庫的數(shù)據(jù)框，后面會細(xì)講整陌，
wid=18,hei=10是輸出圖片的寬和高，可以修改。
#畫圖函數(shù)
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()
}
#主函數(shù)
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)
}
使用示例
加載數(shù)據(jù)
library(DOSE)
data(geneList, package="DOSE")
gene <- names(geneList)[abs(geneList) > 2]
head(gene)
#[1] "4312" "8318" "10874" "55143" "55388" "991"
絕大多數(shù)都可以使用默認(rèn)參數(shù)泌辫，只需改變fun參數(shù)随夸，最終生成三個文件，以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則是可以直接查看的數(shù)據(jù)框文件。
#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')
多個基因集進(jìn)行富集分析函數(shù)compare_result
多個基因集富集分析使用compareCluster函數(shù)墨礁，老規(guī)矩幢竹，先看能生成的4張圖片。
1氣泡圖
2網(wǎng)絡(luò)圖1
3網(wǎng)絡(luò)圖2
4網(wǎng)絡(luò)圖3
compare_result和前面的enrich_result函數(shù)幾乎是一樣的恩静，只是compare_result輸入的是多個基因集的列表焕毫，而enrich_result的輸入是單個基因集向量。相關(guān)參數(shù)驶乾，這里不再贅述邑飒。下面是compare_result的代碼：
#作圖函數(shù)
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()
}
#主函數(shù)
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參數(shù)
導(dǎo)入數(shù)據(jù)
> 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" ...
運行函數(shù)级乐，也會生成三個文件
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函數(shù)
除了已有的數(shù)據(jù)框疙咸，可以自定義基因的涵義進(jìn)行富集分析，例如可以自定義一個細(xì)胞類型的DEGs為一個小數(shù)據(jù)框然后進(jìn)行富集分析风科，可以輔助進(jìn)行細(xì)胞類型注釋罕扎，這可以通過enricher函數(shù)，但前面的兩個函數(shù)也已經(jīng)包含了這個功能丐重。
下面演示如何構(gòu)建一個可用于富集分析的數(shù)據(jù)庫:
首先在這個CellMarker下載細(xì)胞類型的markers列表腔召，我下載的是Cell_marker_Human.xlsx，然后進(jìn)行預(yù)處理扮惦，最終得到一個只有兩列的tibble臀蛛，第一列是基因注釋信息，第二列是基因ID崖蜜，相當(dāng)于一個小的數(shù)據(jù)庫浊仆。
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
然后進(jìn)行分析
x <- enricher(gene, TERM2GENE = cells)
x <- setReadable(x,OrgDb=OrgDb,keyType=keyType)
也可以使用上面兩個函數(shù)
#單個基因集
ob1 <- enrich_result(gene,fun='enricher',interm=cells,label='term')
#多個基因集
ob1 <- compare_result(lgene,fun='enricher',interm=cells,label='term')
可以根據(jù)富集結(jié)果進(jìn)行細(xì)胞類型注釋。
以下面這個圖為例豫领，可以看到特定基因集合在對應(yīng)細(xì)胞類型中的markers基因中富集抡柿。
細(xì)胞類型markers富集
總結(jié)與討論
暫時沒有，以后更新等恐。
