1、富集分析的背景知识 - 简书 (jianshu.com)
2、clusterprofile富集分析--上游分析 - 简书 (jianshu.com)
3、clusterprofile富集分析--下游可视化 - 简书 (jianshu.com)
4、clusterprofile富集分析--多组设计的富集及可视化 - 简书 (jianshu.com)
- clusterProfiler包的
compareCluster()函数提供了针对多组实验设计得到的若干差异基因列表的同时富集分析(ORA)结果，及友好的可视化。 - 该函数主要需要两个参数，一个是多个差异基因列表（可分为两类），另一个是富集分析method（One of "groupGO", "enrichGO", "enrichKEGG", "enrichDO" or "enrichPathway" .）
第一类：multiple gene lists
data(gcSample)
str(gcSample[5:8]) #取后4个差异基因列表做演示
# List of 4
# $ 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" ...
ck <- compareCluster(geneCluster = gcSample[5:8], fun = enrichKEGG)
ck <- setReadable(ck, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
table(ck@compareClusterResult$Cluster)
# X5 X6 X7 X8
# 10 1 15 19
t(ck@compareClusterResult[1,])
# 1
# Cluster "X5"
# ID "hsa04211"
# Description "Longevity regulating pathway"
# GeneRatio "15/446"
# BgRatio "89/8096"
# pvalue "9.100313e-05"
# p.adjust "0.02802896"
# qvalue "0.02461874"
# geneID "PRKAB1/FOXO1/IGF1R/PIK3CD/CREB1/ADCY9/IRS1/ADIPOQ/TSC1/IGF1/ATG5/AKT3/ADCY1/SOD2/TSC2"
# Count "15"
- 结果可视化
dotplot(ck)
第二类：复杂的分组设计
- 假如：分别用两种药物A、B进行干扰细胞系，经1h、3h后得到4组差异基因。想结合分组设计，观察这4组的富集分析结果的比较。
group_A_1h=gcSample[[2]]
group_A_3h=gcSample[[3]]
group_B_1h=gcSample[[4]]
group_B_3h=gcSample[[5]]
design=c(rep("A_1h",length(group_A_1h)),
rep("A_3h",length(group_A_3h)),
rep("B_1h",length(group_B_1h)),
rep("B_3h",length(group_B_3h)))
mydf = data.frame(Drug = stringr::str_split(design,"_",simplify = T)[,1],
Time = stringr::str_split(design,"_",simplify = T)[,2],
DEG = c(group_A_1h,group_A_3h,group_B_1h,group_B_3h))
formula_res <- compareCluster(DEG~Time+Drug, data=mydf, fun="enrichKEGG")
table(formula_res@compareClusterResult$Cluster)
# 1h.A 1h.B 3h.A 3h.B
# 3 18 3 10
t(formula_res@compareClusterResult[1,])
# 1
# Cluster "1h.A"
# Time "1h"
# Drug "A"
# ID "hsa04110"
# Description "Cell cycle"
# GeneRatio "18/384"
# BgRatio "124/8096"
# pvalue "1.983449e-05"
# p.adjust "0.00604952"
# qvalue "0.00565805"
# geneID "991/1869/890/1871/701/990/10926/9088/8317/9700/9134/1029/2810/699/11200/23594/8555/4173"
# Count "18"
- 结果可视化
dotplot(formula_res) +
dotplot(formula_res, x="Time") + facet_grid(~Drug)
