生信绘图生信生物信息

[富集分析] 4、clusterprofile富集分析--多组设

2021-08-15  本文已影响0人  小贝学生信

1、富集分析的背景知识 - 简书 (jianshu.com)
2、clusterprofile富集分析--上游分析 - 简书 (jianshu.com)
3、clusterprofile富集分析--下游可视化 - 简书 (jianshu.com)
4、clusterprofile富集分析--多组设计的富集及可视化 - 简书 (jianshu.com)

第一类: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)

第二类:复杂的分组设计

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)
上一篇下一篇

猜你喜欢

热点阅读