R包fastEnrich预开发一 -- 快速GO富集分析、自动化报告、优化气泡图
开发背景
- 生信分析中, 富集分析是非常常见的分析之一, 而且种类繁多, 至少包括通路富集分析、GSEA、单基因GSEA、GSVA、ssGSEA
- 而且基因ID的准备也需要转来转去, 很是麻烦
- 其次可视化的方式也非常多, 美观度参差不齐
- 最后对与富集结果的报告, 每次要复制粘贴翻译编辑很少麻烦
- 因此, 很久之前就准备开发个R包, 专门实现和可视化所有的富集分析相关的分析
- 虽然很多函数我都编写好了, 但要开发成R包且把我所有的想法都集成起来, 并编写好教学文档和函数文档, 还是非常费时间的
- 所以我就先一点点地开发吧, 先实现基础的GO富集分析和气泡图的美化
- 后续会慢慢开发, 【目前定价是50】, 微信 【直接转账】 给我即可
- 买过我的R包fastGEO V1.8.0 快速下游分析，多数据集批量下载处理，火山图优化都知道, R包的价格会随之更新的内容增加而上调了, 由20调到50再到100。
- 所以, 后续的R包也是一样, 越早买越便宜, 可永久免费获取最新版本！并且有交流群
# # 载入富集分析函数, 作者测试用的
# fun_dir = "W:/02_Study/R_build/build/05_fastEnrich/functions"
# invisible(sapply(lf("rf", fun_dir), source))
# 加载R包
library(fastEnrich)
========================================欢迎使用 fastEnrich (版本 1.0.0)- 作者: 生信摆渡- 微信: bioinfobaidu1- 帮助文档: help(package = fastEnrich)========================================
功能演示
run_GO - 快速GO富集分析
-
首先解决的就是ID自动转换问题
-
我们常用的是SYMBOL, 而GO要使用ENSEMBL, KEGG要使用ENTREZID 。。。
-
run_GO和run_KEGG可以根据输入基因的特征自动判断并转换ID, 无需再手动转换
-
需要注意的是目前只支持人和小鼠的自动富集分析, 因为我不做其他的物种
-
之后就是使用clusterProfiler进行富集分析了
-
然后解决的就是输出表格geneID这一列了, 原始表格依然是ENSEMBL或者ENTREZID, 然而我们还是需要SYMBOL进行快速查看, 所以就再给它转回去。。。
-
最后就是富集结果报告的自动生成, 报告富集到了多少通路, BP CC MF分别有多少，top5是哪些通路，并且进行翻译和文本拼接
-
这里使用limma的差异表达结果的差异基因进行测试
DEG_tb = read.csv("./test/01_GO/DEG_tb.csv", row.names = 1)
DEGs = rownames(DEG_tb)[DEG_tb$DEG != "Not_Change"]
hd(DEGs)
Type: character Dim: 3702 [1] "MMP1" "SPP1" "BPIFB1" "CP" "ND6"
ego = run_GO(DEGs, out_name = "./test/01_GO/01_GO")
Check packages ...Convert gene id ...Run GO ONTOLOGY: ALL ...Build report ...Finished! 727 terms enriched, BP:559 CC:93 MF:75
hd(ego, 5, 10)
Type: data.frame Dim: 727 × 10 ONTOLOGY ID Description GeneRatio BgRatio pvalue p.adjustGO:0002769 BP GO:0002769 natural killer cell inhibitory signaling pathway 34/3402 36/21288 3.281964e-25 2.102098e-21GO:0002765 BP GO:0002765 immune response-inhibiting signal transduction 35/3402 63/21288 6.761106e-13 2.165244e-09GO:0002767 BP GO:0002767 immune response-inhibiting cell surface receptor signaling pathway 34/3402 61/21288 1.268932e-12 2.709170e-09GO:0022904 BP GO:0022904 respiratory electron transport chain 55/3402 140/21288 2.577823e-11 4.127739e-08GO:0042773 BP GO:0042773 ATP synthesis coupled electron transport 48/3402 117/21288 7.858191e-11 8.388619e-08qvalueGO:0002769 1.719058e-21GO:0002765 1.770698e-09GO:0002767 2.215511e-09GO:0022904 3.375592e-08GO:0042773 6.860063e-08geneIDGO:0002769 KIR2DL1GO:0002765 CD33/KIR2DL1GO:0002767 KIR2DL1GO:0022904 ND6/DLD/NDUFA5/CCNB1/VPS54/NDUFS1/ETFDH/SLC25A12/PUM2/GPD2/SDHD/NDUFA7/NDUFS4/NDUFB4/COX4I2/COX7A2/COX4I1/NDUFB9/NDUFA12/GPD1/NDUFS5/NDUFB10/SDHB/NDUFA3/COX7A1/NDUFC2/DGUOK/NDUFA6/UQCR10/SDHC/NDUFS7/COX6B1/NDUFC2-KCTD14/DNAJC15/NDUFA2/COX6A1/COX8A/COX5BGO:0042773 ND6/DLD/NDUFA5/CCNB1/NDUFS1/SDHD/NDUFA7/NDUFS4/NDUFB4/COX4I2/COX7A2/COX4I1/NDUFB9/NDUFA12/NDUFS5/NDUFB10/SDHB/NDUFA3/COX7A1/NDUFC2/DGUOK/NDUFA6/UQCR10/SDHC/NDUFS7/COX6B1/NDUFC2-KCTD14/DNAJC15/NDUFA2/COX6A1/COX8A/COX5BCountGO:0002769 1GO:0002765 2GO:0002767 1GO:0022904 38GO:0042773 32
report_pathway - 自动化报告
- run_GO会自动调用report_pathway_GO, 而后者调用的是report_pathway
# 查看报告
readLines("./test/01_GO/01_GO_report.txt")
‘GO富集分析表明，共富集到727个通路。其中BP共富集到559个通路，富集的top5通路分别为：自然杀伤细胞抑制信号通路（natural killer cell inhibitory signaling pathway）、免疫反应抑制信号转导（immune response-inhibiting signal transduction）、免疫反应抑制剂细胞表面受体信号通路（immune response-inhibiting cell surface receptor signaling pathway）、呼吸电子传递链（respiratory electron transport chain）和ATP合成偶联电子传递（ATP synthesis coupled electron transport）。MF共富集到75个通路，富集的top5通路分别为：NAD（P）H脱氢酶（醌）活性（NAD§H dehydrogenase (quinone) activity）、NADH脱氢酶活性（NADH dehydrogenase activity）、NADH-脱氢酶（泛醌）活性（NADH dehydrogenase (ubiquinone) activity）、ATP水解活性（ATP hydrolysis activity）和NADH-脱氢酶（醌（NADH dehydrogenase (quinone) activity）。CC共富集到93个通路，富集的top5通路分别为：呼吸链复合体（respiratory chain complex）、线粒体呼吸体（mitochondrial respirasome）、线粒体呼吸链复合体I（mitochondrial respiratory chain complex I）、NADH脱氢酶复合体（NADH dehydrogenase complex）和呼吸链复合体I（respiratory chain complex I）。’
- 直接粘贴进分析报告就行了, 手动整理就问你得多久，逆天功能！
- 当然翻译是务必需要人工矫正的。。。
气泡图
-
默认每个ONTOLOGY绘制10个, 也可以选择绘制所有的, 当然这种情况是你已经对ego做了筛选, 不然一般会很多
-
自动化输出图片有个问题是图片宽度的设置, 这里也自动根据字符长度最长的通路进行自动调整宽度了, 当然如果我自动调整的不合适你可以再指定参数
w -
设置
out_name指定输出文件名, 不带后缀 -
由于有的通路名实在太长, 放在文章里实在不合适, 但又不想删掉, 解决策略就是省略后面一部分, 所以自动设置了最长显示字符宽度为60, 可以通过
max_omit调整和omit = FALSE进行关闭自动省略字符功能。 -
默认参数
set_image(10, 8) # set_image仅控制在jupyter里图片的宽度, 不影响Rstudio的显示和文件的输出
dotplot_GO(ego, out_name = "test/01_GO/02_GO_dotplot")
- 如果想少一些, 调整
ntop
set_image(10, 8) # set_image仅控制在jupyter里图片的宽度, 不影响Rstudio的显示和文件的输出
dotplot_GO(ego, out_name = "test/01_GO/02_GO_dotplot2", ntop = 5)
- 不进行字符省略
set_image(13, 8)
dotplot_GO(ego, out_name = "test/01_GO/02_GO_dotplot3", omit = FALSE)
- 虽然气泡按Gene ratio排序会更好看, 但是我们关注的还是p值
- 所以你也可以选择用P值进行排序
set_image(10, 8)
dotplot_GO(ego, out_name = "test/01_GO/02_GO_dotplot4", sort = "pvalue")
- 自定义气泡颜色, 颜色支持英文和十六进制
set_image(10, 8)
dotplot_GO(ego, out_name = "test/01_GO/02_GO_dotplot5", color_high = "orange", color_low = "#FBFFD4")
- 高亮感兴趣的通路, 并自定义颜色
highlight_paths = c("immune response-inhibiting signal transduction", "MHC class II protein complex", "NADH dehydrogenase complex", "MHC protein complex binding","mRNA splicing, via spliceosome")
dotplot_GO(ego, out_name = "test/01_GO/02_GO_dotplot5", color_high = "orange", color_low = "#FBFFD4", highlight = highlight_paths, highlight_color = "#FFA700")
- 还可以添加title, 不过这些基础的ggplot操作都可以后续自由发挥
- 如果有多个气泡图同时展示时, 比如不同的细胞类型、亚群、比较组别, 分别指定不同的颜色排列在一起会好看很多
dotplot_GO(ego, out_name = "test/01_GO/02_GO_dotplot6", color_high = "orange", color_low = "#FBFFD4", highlight = highlight_paths, highlight_color = "orange",title = "A vs B")
- 也可以手动挑选想要展示的通路, 这时候就不是展示前10了，而是给哪些展示哪些, 设置
show = "all"
# 三种ONTOLOGY分别取 15 10 5进行展示
neach = c(15, 10, 5)
ego2 = Reduce(rbind, lapply(1:length(neach), function(i) split(ego, f = ego$ONTOLOGY)[[i]][1:neach[i], ]))
hd(ego2)
Type: data.frame Dim: 30 × 10 ONTOLOGY ID Description GeneRatio BgRatioGO:0002769 BP GO:0002769 natural killer cell inhibitory signaling pathway 34/3402 36/21288GO:0002765 BP GO:0002765 immune response-inhibiting signal transduction 35/3402 63/21288GO:0002767 BP GO:0002767 immune response-inhibiting cell surface receptor signaling pathway 34/3402 61/21288GO:0022904 BP GO:0022904 respiratory electron transport chain 55/3402 140/21288GO:0042773 BP GO:0042773 ATP synthesis coupled electron transport 48/3402 117/21288
dotplot_GO(ego2, out_name = "test/01_GO/02_GO_dotplot7", color_high = "orange", color_low = "#FBFFD4", highlight = highlight_paths, highlight_color = "orange",show = "all",title = "A vs B")
软件安装
公共依赖包
-
依赖包自己安装一下好吧, 无非是那几种安装方法
-
install.packages、BiocManager::install、remotes::install_github
-
软件安装问题不免费答疑, 这是基本功, 也跟我的开发关系不大
-
ggplot2 (>= 3.5.2),
-
ggthemes (>= 5.1.0),
-
ggtext (>= 0.1.2),
-
grid (>= 4.4.1),
-
clusterProfiler (>= 4.12.3),
-
org.Hs.eg.db (>= 3.19.1),
-
rvest (>= 1.0.4),
自研依赖包
- 当然也依赖我自己开发的一些R包, 在这里提供了压缩包, 直接运行可快速安装
- fanyi2 (>= 1.0.0)
- fastR (>= 1.5.0)
# 安装 fastR
file_fastR = grep("^fastR.*.gz", dir(), value = TRUE)
install.packages(file_fastR, repos = NULL, upgrade = FALSE, force = TRUE)
# 安装 fanyi2
file_fanyi2 = grep("^fanyi2.*.gz", dir(), value = TRUE)
install.packages(file_fanyi2, repos = NULL, upgrade = FALSE, force = TRUE)