上一讲中刚教会大家如何进行差异基因的计算及可视化，然而大家拿着成千上百的差异基因直接调研、撰文显然是不现实的，因此需要通过富集分析去揭示这些基因群背后蕴含的 生物学意义 1. 视频及测试文件 1.1 视频教程 1.2 测试文件 如果你获取了当前教学文档，那么随之下载的便是测试文件，此外还可在推送中获取:手把手教你做单细胞测序数据分析
（七）:基因集富集分析 2. 代码 suppressWarnings ( suppressMessages ( if ( ! require (rvcheck))devtools :: install_version ( "rvcheck" , version = "0.1.8" , repos = "http://cran.us.r-project.org" ))) suppressWarnings ( suppressMessages ( if ( ! require (clusterProfiler))BiocManager :: install ( "clusterProfiler" ))) suppressWarnings ( suppressMessages ( if ( ! require (org.Mm.eg.db))BiocManager :: install ( "org.Mm.eg.db" ))) suppressWarnings ( suppressMessages ( if ( ! require (org.Hs.eg.db))BiocManager :: install ( "org.Hs.eg.db" ))) suppressWarnings ( suppressMessages ( library (dplyr))) 2.1 准备基因和绘图函数 # 读取marker列表 bcell.marker <- read.csv ( './bcell.marker.CSV' , row.names = 1 ) head (bcell.marker) # p_val avg_log2FC pct.1 pct.2 p_val_adj SYMBOL # GZMB 3.186146e-97 -5.507567 0.017 0.961 4.369480e-93 GZMB # NKG7 1.262672e-94 -6.262287 0.087 1.000 1.731628e-90 NKG7 # PRF1 1.835435e-93 -4.612813 0.029 0.948 2.517115e-89 PRF1 # GZMA 3.036581e-93 -4.466820 0.015 0.929 4.164367e-89 GZMA # CTSW 3.151498e-93 -4.203776 0.070 0.987 4.321965e-89 CTSW # CST7 9.767143e-92 -4.380031 0.038 0.948 1.339466e-87 CST7 # 筛选出B细胞marker中p_val_adj最显著(越小越显著)的100个基因用于下游分析 mygene <- bcell.marker %>% top_n ( n = - 100 , wt = p_val_adj) %>% rownames mygene[ 1 : 10 ] # [1] "GZMB" "NKG7" "PRF1" "GZMA" "CTSW" "CST7" "GNLY" "FGFBP2" # [9] "CD247" "FCGR3A" # 利用org.Hs.eg.db整合基因名、 ENTREZID、ENSEMBL gene.df <- bitr (mygene, fromType= "SYMBOL" , toType= c ( "ENTREZID" , "ENSEMBL" ), OrgDb = org.Hs.eg.db) 2.2 KEGG 与 GO 富集分析 # 使用在线数据,可能会受域名限制 # 设置clusterProfiler包的下载方法，保证在使用clusterProfiler进行数据库查询或注释数据下载时，可以自动选择合适的下载方式（如 curl、wget 或 internal） R.utils :: setOption ( "clusterProfiler.download.method" , "auto" ) # 进行kegg富集分析 ekegg <- enrichKEGG ( unique (gene.df $ ENTREZID), organism= 'hsa' , pvalueCutoff= 0.05 , pAdjustMethod= 'BH' , qvalueCutoff= 0.2 , minGSSize= 10 , maxGSSize= 500 , use_internal_data= F) # 把ENTREZID转换为SYMBOL ekegg <- setReadable (ekegg, 'org.Hs.eg.db' , 'ENTREZID' ) write.csv (ekegg @ result, '../result_enrich/kegg.result.csv' ) library ( "enrichplot" ) dotplot (ekegg, x = "GeneRatio" , color = "p.adjust" , size = "Count" , showCategory = 10 ) # 保存图片 png ( file = '../result_enrich/kegg.png' , width = 10 , height = 10 ) dotplot (ekegg, x = "GeneRatio" , color = "p.adjust" , size = "Count" , showCategory = 10 ) dev.off # png # 2 如果你在做KEGG富集分析时遇到这样的报错：那么你的 clusterProfiler 需要更新一下：# remove.packages('clusterProfiler') # devtools::install_github("YuLab-SMU/clusterProfiler") function、biological process、cell component三个类 # MF ggoMF <- enrichGO (gene.df $ SYMBOL, org.Hs.eg.db, keyType = "SYMBOL" , ont = "MF" , pvalueCutoff = 0.05 , pAdjustMethod = "BH" , qvalueCutoff = 0.2 , minGSSize = 10 , maxGSSize = 500 , readable = FALSE , pool = FALSE ) write.csv (ggoMF @ result, '../result_enrich/ggoMF.result.csv' ) dotplot (ggoMF, x = "GeneRatio" , color = "p.adjust" , size = "Count" , showCategory = 10 ) png ( file = '../result_enrich/ggoMF.png' , width = 10 , height = 10 ) dotplot (ggoMF, x = "GeneRatio" , color = "p.adjust" , size = "Count" , showCategory = 10 ) dev.off # png # 2 # CC ggoCC <- enrichGO (gene.df $ SYMBOL, org.Hs.eg.db, keyType = "SYMBOL" , ont = "CC" , pvalueCutoff = 0.05 , pAdjustMethod = "BH" , qvalueCutoff = 0.2 , minGSSize = 10 , maxGSSize = 500 , readable = FALSE , pool = FALSE ) write.csv (ggoCC @ result, '../result_enrich/GO.CC.result.csv' ) dotplot (ggoCC, x = "GeneRatio" , color = "p.adjust" , size = "Count" , showCategory = 10 ) png ( file = '../result_enrich/ggoCC.png' , width = 10 , height = 10 ) dotplot (ggoCC, x = "GeneRatio" , color = "p.adjust" , size = "Count" , showCategory = 10 ) dev.off # png # 2 ggoBP <- enrichGO (gene.df $ SYMBOL, org.Hs.eg.db, keyType = "SYMBOL" , ont = "BP" , pvalueCutoff = 0.05 , pAdjustMethod = "BH" , qvalueCutoff = 0.2 , minGSSize = 10 , maxGSSize = 500 , readable = FALSE , pool = FALSE ) write.csv (ggoBP @ result, '../result_enrich/GO.BP.result.csv' ) dotplot (ggoBP, x = "GeneRatio" , color = "p.adjust" , size = "Count" , showCategory = 10 ) png ( file = '../result_enrich/ggoBP.png' , width = 10 , height = 10 ) dotplot (ggoBP, x = "GeneRatio" , color = "p.adjust" , size = "Count" , showCategory = 10 ) dev.off # png # 2 2.3 GSEA分析 bcell.marker $ SYMBOL <- rownames (bcell.marker) genetable <- merge (gene.df,bcell.marker, by= 'SYMBOL' ) geneList <- genetable $ avg_log2FC names (geneList) <- genetable $ ENTREZID # 用于gsea分析的genelist本质上是一个名为ENTREZID的数值排序结果 geneList <- sort (geneList, decreasing = T) egseGO <- gseGO (geneList, ont = "MF" , org.Hs.eg.db, exponent = 1 , nPerm = 1000 , minGSSize = 10 , maxGSSize = 500 , pvalueCutoff = 1 , pAdjustMethod = "BH" , verbose = TRUE , seed = FALSE , by = "fgsea" ) BP kk2 <- gseKEGG ( geneList= geneList, organism = "hsa" , => human , mmu => mouse nPerm= 1000 , minGSSize = 120 , pvalueCutoff = 1 , verbose = FALSE ) library (patchwork) # 拼图神包，详细使用方法可以看：相关链接见专题页 for (i in 1 : length (egseGO $ Description)) { plot1 <- gseaplot (egseGO, geneSetID = i, title = egseGO $ Description[i]) ggsave ( filename = paste0 ( "../result_enrich/" ,egseGO $ Description[i], ".pdf" ), plot = plot1[[ 2 ]] / plot1[[ 1 ]], height = 10 , width = 10 ) } gseaplot (egseGO, geneSetID = 1 , title = egseGO $ Description[ 1 ]) 2.4 GSVA分析 一文了解GSVA原理 三行搞定GSVA分析 一文搞定GSVA及下游分析 单细胞中如何做GSVA