查看: 456|回复: 6

[技术前沿] 报错不怕 计算不慌 美图流畅 关于功能富集,它一直很可以

[复制链接]

管理员

Rank: 15Rank: 15Rank: 15Rank: 15Rank: 15

主题
666
注册时间
2020.6.16
在线时间
386 小时

发表于 2022.6.13 10:02:10 | 显示全部楼层 |阅读模式
前几期我们学习了基于clusterProfiler的功能富集分析,包括KEGG、GO及GSEA富集和可视化,详情戳:

《KEGG富集分析及可视化,一把子拿捏!》
《GO富集分析及可视化,一把子拿捏!》
《快速拿捏KEGG/GO/Reactome/Do/MSigDB的GSEA富集分析!》

今天是基于这三期内容的补充:为什么使用enrichKEGG/GO等富集函数时会出现报错?如何解决?如何将富集结果中的ENTREZID直接转换为symbol?如何计算富集因子(RichFactor)和富集倍数(FoldEnrichment)?可视化补充(富集条形图、富集气泡图、富集棒棒糖图、富集网络图、富集热图、富集upset图、文章趋势图)

由于是补充,不涉及具体的富集流程,想学习富集流程的可以戳前三期(下文实操所使用全部数据也源自前三期)。下面开始今天的内容!

  1. #所需R包载入:
  2. library(dplyr)
  3. library(org.Hs.eg.db)
  4. library(clusterProfiler)
  5. library(ggplot2)
  6. library(RColorBrewer)
  7. library(enrichplot)
  8. library(DOSE)
  9. library(ggupset)
复制代码

一、富集时报错如何处理?

之前小编是没出现过报错问题的,不过本周将RStudio更新了一下,意外发现跑富集流程时竟然error了:


想起近期在生信技能树刚好看到了相关文章,报错原因是嵌套的两个下载函数(utils::download.file和downloader::download)是使用默认协议“auto”才能正常访问https://rest.kegg.jp/link/hsa/pathway进行下载,现在下载协议变为了"libcurl",就会访问失败。如果大家出现了报错,可以使用如下代码更改下载协议:

  1. install.packages("R.utils")
  2. R.utils::setOption( "clusterProfiler.download.method","auto")
复制代码

然后就可以正常进行富集分析啦!

二、entrezID如何转symbol?

在使用enrichGO函数进行GO富集时,可以直接设置参数readable = TRUE将entrez转为symbol:


但在进行KEGG和GSEA富集时不可以,我们拿到的富集结果(以KEGG为例)是entrezID:


转化非常简单,使用setReadable函数即可:

  1. #将ENTREZID重转为symbol:
  2. KEGG_diff2 <- setReadable(KEGG_diff,
  3. OrgDb = org.Hs.eg.db,
  4. keyType="ENTREZID")
  5. #转换前:
  6. KEGG_diff@result$geneID[1]
  7. #转换后:
  8. KEGG_diff2@result$geneID[1]
复制代码


三、如何计算RichFactor(富集因子)和FoldEnrichment(富集倍数)?

富集得到的结果表格中是没有这两列的,但我们绘制富集柱形图或者富集气泡图时,经常会需要用到,计算方法如下:

富集因子=差异基因中富集到该通路的基因数(count)/背景基因中富集到该通路的基因数(BgRatio分子)
富集倍数=GeneRatio/BgRatio

下面直接计算并添加到富集结果表(result):

  1. #计算Rich Factor(富集因子):
  2. KEGG_diff2 <- mutate(KEGG_diff2,
  3. RichFactor = Count / as.numeric(sub("/\\d+", "", BgRatio)))

  4. #计算Fold Enrichment(富集倍数):
  5. KEGG_diff2 <- mutate(KEGG_diff2, FoldEnrichment = parse_ratio(GeneRatio) / parse_ratio(BgRatio))
  6. KEGG_diff2@result$RichFactor[1:6]
  7. KEGG_diff2@result$FoldEnrichment[1:6]
复制代码


现在就得到了一个完整的富集结果表,作图时需要哪一列直接调用即可,也可以输出到本地:

  1. write.csv(KEGG_diff2@result,"./KEGG_diff.csv",row.names = TRUE)
复制代码


四、基于clusterProfiler的快速可视化补充

《KEGG富集分析及可视化,一把子拿捏!》可视化部分,我们学习了如何调出感兴趣KEGG通路、使用enrichplot包进行富集条形图和气泡图的快速可视化,以及使用ggplot2进行更佳个性化图表绘制。

《GO富集分析及可视化,一把子拿捏!》可视化部分,除了富集条形图和气泡图(包含ggplot2绘图小技巧),还补充了GO有向无环图及富集网络图。

《快速拿捏KEGG/GO/Reactome/Do/MSigDB的GSEA富集分析!》可视化部分,绘制了气泡图、山峦图及GSEA的enrichmentScore折线图等。

除了上述高频使用的图表,下面在给大家补充几种,可以按需使用。

#1.富集棒棒糖图(条形图+气泡图的组合)

  1. mytheme <- theme(axis.title = element_text(size = 13),
  2. axis.text = element_text(size = 11),
  3. plot.title = element_text(size = 14,
  4. hjust = 0.5,
  5. face = "bold"),
  6. legend.title = element_text(size = 13),
  7. legend.text = element_text(size = 11))

  8. ggplot(GO_CC_diff2,
  9. showCategory = 20,
  10. aes(x = RichFactor,
  11. y = reorder(Description, RichFactor))) + #将Description按照RichFactor进行排序
  12. geom_segment(aes(xend=0, yend = Description),size = 1.1,color = "grey") +
  13. geom_point(aes(color=-log10(pvalue), size = Count)) +
  14. scale_color_distiller(palette = "Spectral",direction = -1)+
  15. scale_size_continuous(range=c(2, 8)) + #气泡大小范围
  16. theme_minimal() +
  17. labs(x = "Rich factor",
  18. y = NULL,
  19. title = "GO CC terms Enrichment") + mytheme
复制代码


#2.富集基因网络图

用于展示感兴趣的通路所富集到的差异基因有哪些及上下调情况(logFC值),还能够看出哪些基因同时出现在多个通路(重叠关系)。

  1. ##先制作差异基因的genelist(需要差异分析结果表格中的log2FoldChange值):
  2. genelist <- DESeq2$log2FoldChange
  3. names(genelist) <- rownames(DESeq2)
  4. length(genelist)#全部基因,16345个

  5. ##筛选差异基因:
  6. genelist <- genelist[names(genelist) %in% diff]
  7. length(genelist)#3404个差异基因
  8. head(genelist)
复制代码


  1. cnetplot(KEGG_diff2,
  2. showCategory = 5,#显示top5 pathway
  3. categorySize = "pvalue",
  4. foldChange = genelist)
复制代码


  1. #自定义展示感兴趣的通路:
  2. interest <- c("Butanoate metabolism","Histidine metabolism","Proximal tubule bicarbonate reclamation","Propanoate metabolism","Primary bile acid biosynthesis")
  3. cnetplot(KEGG_diff2,
  4. showCategory = interest,
  5. categorySize="pvalue",
  6. foldChange=genelist)
  7. #可以用layout参数调整网络图布局,可选'star', 'circle', 'gem', 'dh', 'graphopt', 'grid', 'mds', 'randomly', 'fr', 'kk', 'drl' or 'lgl'.
复制代码


  1. #另一种形式:
  2. interest2 <- c("Butanoate metabolism","Histidine metabolism","Propanoate metabolism")

  3. cnetplot(KEGG_diff2,
  4. showCategory = interest2,
  5. foldChange=genelist,
  6. circular = TRUE,
  7. colorEdge = TRUE)
复制代码


网络图中pathway或基因名如果出现重叠,可以导出矢量图后期用AI、PPT等软件进行微调。

#3.富集通路网络图

主要用于展示不同节点/通路间的基因重叠关系。如果两条通路内富集到的差异基因存在重叠,则会用线条进行连接。节点的大小表示富集到基因的数量。

  1. GO_CC_diff2_edo <- pairwise_termsim(GO_CC_diff2)
  2. emapplot(GO_CC_diff2_edo,
  3. cex_category=1, #节点大小
  4. showCategory = 40) #展示terms的数量
复制代码


不仅是富集网络图,富集条形图、富集气泡图等等都可以自定义展现感兴趣的通路,而非仅展示top。

#4.富集热图

用热图的形式展示感兴趣通路所富集到的差异基因及上下调情况。

  1. heatplot(KEGG_diff2,
  2. foldChange = genelist,
  3. showCategory = interest)
复制代码


#5.富集upset韦恩图

主要用于直观展示不同通路间富集到的差异基因数量和不同通路间基因的重叠关系及重叠基因的数量。

  1. upsetplot(KEGG_diff2)
复制代码


#6.相关发表文章数趋势图

  1. #根据PubMed Central查询结果绘制出版文章趋势图:
  2. install.packages("europepmc")
  3. library(europepmc)

  4. terms <- KEGG_diff2$Description[1:5] #取感兴趣的pathway(这里取top5)
  5. p <- pmcplot(terms, 2010:2021,proportion=FALSE)#在PubMed Central中查询2010-2021年跟着5个pathway有关的出版物文章数量
  6. p
复制代码


感觉很适合用于汇报,或者查询感兴趣的pathway是否有研究价值。

更详细说明都可以在clusterProfiler包网页说明中Part II第15小节去查阅。今天的分享就到这里!

参考资料

1.https://mp.weixin.qq.com/s/3MBYjm0bNvN6kxairREljA
‍2.https://yulab-smu.top/biomedical-knowledge-mining-book/index.html



本文作者:基迪奥-喵酱

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
3
注册时间
2021.6.22
在线时间
54 小时

发表于 2022.6.13 18:59:05 | 显示全部楼层
周五啦!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
0
注册时间
2020.11.21
在线时间
101 小时

发表于 2022.6.14 12:50:55 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
3
注册时间
2021.6.22
在线时间
54 小时

发表于 2022.6.14 22:37:29 | 显示全部楼层
再次路过
周五啦!
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2018.8.16
在线时间
38 小时

发表于 2022.6.15 08:35:26 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
0
注册时间
2020.11.21
在线时间
101 小时

发表于 2022.6.16 14:00:57 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
0
注册时间
2020.11.21
在线时间
101 小时

发表于 7 天前 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

快速回复 返回顶部 返回列表