关于富集分析这件事，我有话想说…

那些年，

我们曾经见过的富集化分析结果，

你确定都见过吗？

Hi，大家好，我是晨曦，今天是晨曦碎碎念的第六期，这期推文的灵感来源于以下和学员的对话：

A同学：晨曦，富集分析有多种展现方式，我该如何选择呢？

晨曦：可视化的本质是数据，挑选自己觉得好看同时杂志社要求的格式即可

A同学：那么曦曦，能否整理出一篇笔记，告诉我们富集分析究竟有多少种展现形式？

晨曦：好的，关注下期推文哦~

本篇推文基于以上问题诞生，将聚焦在富集分析究竟有多少种展现形式，书写这篇推文的同时也是晨曦对于富集分析这一方面知识的又一次梳理和总结，那么，我们就开始吧~

干货预警！！预计阅读时间10min，建议配上音乐一起食用哦~

晨曦碎碎念系列传送门

想白嫖单细胞生信文章？这五大源头数据库，是你发文章的源泉！高频预警！你一定要收藏！

盘活国自然的新思路！你研究的热点真的是热点吗？大数据帮你定位！

好家伙！90%以上审稿人都会问到的问题，今天帮你解决！就是这么齐齐整整！

没想到！生信分组还有这个大坑！你被坑过吗？！

富集分析的相关概念

这里并不准备用太多官方的语言来给大家进行解释，大家只需要知道，我们筛选出来的差异基因（DEGs），需要给其明确一个功能，也就是说给其一个具体用途，DEGs究竟是干什么的、行使了什么功能，那么我们就需要进行富集分析，而我们富集分析最常用的R包就是——clusterProfifiler包，最近已经升级为4.0全新版本，所以本篇推文也可以算作是对clusterProfifiler包进行一个全面的梳理

Ps:毫不夸张的说，掌握了这个R包，基本上整个富集分析就已经掌握了，所以性价比还是十分高的

富集分析结果剖析

在进行可视化的前提，我们需要对我们富集分析后的结果有一个了解

了解输出数据即可以帮助我们更好的选择可视化的参数和方式，同时通过了解输出数据，我们也可以对富集分析的结果用其它R包，诸如ggplot2包、igraph包进行更加多元的可视化展现，那么我们接下来首先来获得一个富集分析的结果

#准备工作（加载R包+获得输入数据）library(clusterProfiler)library(ggplot2)library(org.Hs.eg.db)library(enrichplot)library(ggupset)data(geneList,package = "DOSE")#geneList所包含的信息其实就是两个变量，一个变量是Gene ID，一个变量是表达量gene <- names(geneList[abs(geneList)>2])#对gene的表达进行筛选head(gene)#[1] "4312" "8318" "10874" "55143" "55388" "991"

晨曦解读

clusterProfiler包支持的输入数据建议是Entrez gene ID，而不建议适用Gene symbol，而且富集分析我们只需要一个信息，就是Gene ID

#GO富集分析的两种形式#方法一ggo <- groupGO(gene = gene, OrgDb = org.Hs.eg.db, ont = "CC", level = 3, readable = TRUE)#gene参数就是我们的输入数据-Gene ID#OrgDb参数为我们富集的物种参考信息（目前来说支持二十个物种）#ont参数为我们GO富集的方向，包括MF、CC、BP以及all

#level参数为富集水平，可以理解为：级别1提供了最高的列表覆盖率和最少的术语特异性。随着每增加一级的覆盖范围减少，而特异性增加，因此第5级提供的覆盖范围最小，特异性最高，简单来说就是数值越大富集到的功能特异性越高，当然默认也是可以的，默认数值为2

#readable参数为True的时候会在结果部分把Gene ID转换为Gene symbol#方法二ego <- enrichGO(gene = gene, universe = names(geneList), OrgDb = org.Hs.eg.db, ont = "CC", pAdjustMethod = "BH", pvalueCutoff = 0.01, qvalueCutoff = 0.05, readable = TRUE)

晨曦解读

经常有小伙伴问可不可以富集到小鼠的GO或者KEGG，经过上面的解析答案当然是可以的，只需要切换OrgDb参数即可，下面展示了部分支持物种

接下来，我们对方法一和方法二富集的结果来查看一下，其实也不用区分的很清晰，方法二富集得到的信息更多，且可以对结果进行筛选

#方法一的结果展示## ID Description Count GeneRatio geneID## GO:0000133 GO:0000133 polarisome 0 0/207 ## GO:0000408 GO:0000408 EKC/KEOPS complex 0 0/207 ## GO:0000417 GO:0000417 HIR complex 0 0/207 ## GO:0000444 GO:0000444 MIS12/MIND type complex 0 0/207 ## GO:0000808 GO:0000808 origin recognition complex 0 0/207 ## GO:0000930 GO:0000930 gamma-tubulin complex 0 0/207#方法二的结果展示## ID Description GeneRatio## GO:0005819 GO:0005819 spindle 26/201## GO:0000779 GO:0000779 condensed chromosome, centromeric region 16/201## GO:0072686 GO:0072686 mitotic spindle 17/201## GO:0000775 GO:0000775 chromosome, centromeric region 18/201## GO:0098687 GO:0098687 chromosomal region 23/201## GO:0000776 GO:0000776 kinetochore 15/201## BgRatio pvalue p.adjust qvalue## GO:0005819 306/11853 1.072029e-11 3.151766e-09 2.888837e-09## GO:0000779 114/11853 7.709944e-11 8.659125e-09 7.936756e-09## GO:0072686 133/11853 8.835841e-11 8.659125e-09 7.936756e-09## GO:0000775 158/11853 1.684987e-10 1.179661e-08 1.081250e-08## GO:0098687 272/11853 2.006225e-10 1.179661e-08 1.081250e-08## GO:0000776 106/11853 2.733425e-10 1.339378e-08 1.227644e-08## geneID

## GO:0005819 CDCA8/CDC20/KIF23/CENPE/ASPM/DLGAP5/SKA1/NUSAP1/TPX2/TACC3/NEK2/CDK1/MAD2L1/KIF18A/BIRC5/KIF11/TRAT1/TTK/AURKB/PRC1/KIFC1/KIF18B/KIF20A/AURKA/CCNB1/KIF4A

## GO:0000779 CENPE/NDC80/HJURP/SKA1/NEK2/CENPM/CENPN/ERCC6L/MAD2L1/KIF18A/CDT1/BIRC5/TTK/NCAPG/AURKB/CCNB1## GO:0072686 KIF23/CENPE/ASPM/SKA1/NUSAP1/TPX2/TACC3/CDK1/MAD2L1/KIF18A/KIF11/TRAT1/AURKB/PRC1/KIFC1/KIF18B/AURKA## GO:0000775 CDCA8/CENPE/NDC80/TOP2A/HJURP/SKA1/NEK2/CENPM/CENPN/ERCC6L/MAD2L1/KIF18A/CDT1/BIRC5/TTK/NCAPG/AURKB/CCNB1

## GO:0098687 CDCA8/CENPE/NDC80/TOP2A/HJURP/SKA1/NEK2/CENPM/RAD51AP1/CENPN/CDK1/ERCC6L/MAD2L1/KIF18A/CDT1/BIRC5/EZH2/TTK/NCAPG/AURKB/CHEK1/CCNB1/MCM5

## GO:0000776 CENPE/NDC80/HJURP/SKA1/NEK2/CENPM/CENPN/ERCC6L/MAD2L1/KIF18A/CDT1/BIRC5/TTK/AURKB/CCNB1## Count## GO:0005819 26## GO:0000779 16## GO:0072686 17## GO:0000775 18## GO:0098687 23## GO:0000776 15

晨曦解读

很直观的就可以看到，方法二富集得到的信息要比方法一富集得到的信息多，对于GO分析，我们也使用方法二的次数要比较多，上面都是对结果的一个微观的展示，我们直接来看一下我们的环境中得到的富集分析结果，以下剖析的富集分析结果我们都默认用方法二得到的结果展示

可以很清晰的看到，我们富集得到的结果是一个S4对象的形式，也就是可以简单来说数据类型是一个列表的形式，然后我们也可以很清楚的看到各个数据的类型是什么，下面我们分别提取富集分析结果中的两类信息来看一下

#提取结果head(ego@result,2)

从图中可以看到result总共有9列

第一列是ID,也就是富集通路的编号；

第二列是Description，也就是富集通路的名称；

第三列是GeneRatio,也就是要富集的基因中在对应通路中的比例；

第4列是BgRation,也就是对应的通路基因在全基因组注释中的比例；

第5,6,7列都是统计检验的结果；

第8列是geneID，也就是富集到基因的名字，它的格式是以斜线隔开的；

第9列是Count，也就是富集到的基因数目

#提取Genesymboldata.frame(symbol=head(ego@gene2Symbol,2))

晨曦解读

很清楚的看到我们得到了一个Gene ID和Gene symbol的对应数据框，也就是说我们也可以通过富集分析来转化我们的Gene ID，当然了，Gene ID的转换在clusterProfifiler包中也有专门的函数可以做到，放在后面的富集分析技巧中讲解~

至此，我们就讲解完了GO富集分析的结果内容，其实KEGG富集分析的结果和这个大同小异，掌握一个即可类推到另一个，那么我们下面进行富集分析可视化的展示

富集分析可视化展示

Bar Plot

barplot(ego, showCategory=10)

Dot plot

dotplot(ego, showCategory=10) + ggtitle("dotplot for ORA")

Gene-Concept Network

barplot()和dotplot()都只显示最重要的或选择的术语，而用户可能想知道哪些基因参与了这些重要术语。为了考虑一个基因可能属于多个注释类别的潜在生物复杂性，并提供数字变化的信息，可以尝试使用下面这个可视化来展示相关性

## convert gene ID to Symboledox <- setReadable(ego, 'org.Hs.eg.db', 'ENTREZID')

#上一步转换如果富集的时候readable参数为TRUE，其实这一步就没有必要了，也就是说如果我们在上一步readable参数选择FALSE，这里我们也可以通过额外的函数完成Gene ID的转换

p1 <- cnetplot(edox, foldChange=geneList)## categorySize can be scaled by 'pvalue' or 'geneNum'p2 <- cnetplot(edox, categorySize="pvalue", foldChange=geneList)p3 <- cnetplot(edox, foldChange=geneList, circular = TRUE, colorEdge = TRUE) #这一步需要提供足够大的画板，否则会报错无法出图cowplot::plot_grid(p1, p2, p3, ncol=3, labels=LETTERS[1:3], rel_widths=c(.8, .8, 1.2))#获取函数帮助文档可以获得更多的参数调节，默认参数依旧颜值抗打

当然这个图因为看上去比较复杂，所以我们仍旧可以DIY一下：

#我们想要着重凸显功能，让Gene symbol消失p1 <- cnetplot(edox, node_label="category", cex_label_category = 1.2) #让功能消失，我们想凸显出Gene symbolp2 <- cnetplot(edox, node_label="gene", cex_label_gene = 0.8) #全部存在，但是代码简单点p3 <- cnetplot(edox, node_label="all") #想要颜色更加自定义p4 <- cnetplot(edox, node_label="none", color_category='firebrick', color_gene='steelblue')

Heatmap

p1 <- heatplot(edox, showCategory=5)p2 <- heatplot(edox, foldChange=geneList, showCategory=5)cowplot::plot_grid(p1, p2, ncol=1, labels=LETTERS[1:2])

Tree plot

如果你对于富集分析的结果觉得太散乱，需要一个专门的分类，那么下面这个绘图功能将会帮助到你。

edox2 <- pairwise_termsim(edox)#这一步骤可以简单理解计算了GO结果之间彼此的相似度，其实就是在我们的富集结果中添加了这个信息，然后我们才可以绘制下面这个图形

p1 <- treeplot(edox2)p2 <- treeplot(edox2, hclust_method = "average")aplot::plot_list(p1, p2, tag_levels='A')#这里建议最好不要把图拼在一起，会因为图展不开而导致分类术语无法显示完全#当然也可以通过调节nWords和fontsize参数对字体进行限制，两种思路

#p1和p2本质上是聚集方法不同，也就是GO语义聚集因为何种算法聚集在一起中的算法不同，可以通过hclust_method参数设置，默认为ward.D

晨曦解读

treeplot函数将把树切割成几个子树(由nCluster参数指定(默认为5))，并使用高频词对子树进行标签。这将降低丰富结果的复杂性，提高用户的解释能力，所以对于算法的选择，可以参考帮助文档或者直接默认即可

注意：这里如果出现报错除了画板大小的问题，还有就是R包版本的问题，这里我们使用的R包版本如下：

Enrichment Map

富集图将丰富的术语组织成一个网络，网络的边缘连接重叠的基因集。这样，相互重叠的基因集往往聚在一起，便于识别功能模块。

ego <- pairwise_termsim(ego)p1 <- emapplot(ego)p2 <- emapplot(ego, cex_category=1.5)#调节点的大小p3 <- emapplot(ego, layout="kk")#调节整体布局p4 <- emapplot(ego, cex_category=1.5,layout="kk")cowplot::plot_grid(p1, p2, p3, p4, ncol=2, labels=LETTERS[1:4])#直接运行容易报错的原因是画板太小，建议分步展开结果#下面只展现p4的结果

UpSet Plot

可以使基因和基因组之间的复杂关联可视化。它强调不同基因组之间的基因重叠。

upsetplot(edox)

ridgeline plot for expression distribution of GSEA result

当然上面的可视化大部分都是用来展示GO和KEGG的

Ps：其实我们常用的富集分析GO和KEGG，这个名称是按照注释来源的数据库来命名的，除此之外，我们如果更换注释数据库，包括MeSH、Reactome、DOSE等这类数据库也会得到不同的注释结果，好在Y叔把以上不同注释来源的数据库写成了不同的函数，我们直接选择即可

当然目前来看富集分析的算法就是两种即Over Representation Analysis（ORA）和Gene Set Enrichment Analysi（GSEA）

也就是说GO和KEGG这两个注释来源也同样支持不同算法，这里大家需要做出区分，不要混淆了

我们接下来看一下GSEA算法下的可视化展示

#首先进行GSEA分析ego3 <- gseGO(geneList = geneList, OrgDb = org.Hs.eg.db, ont = "CC", minGSSize = 100, maxGSSize = 500, pvalueCutoff = 0.05, verbose = FALSE)#Beware that only gene Set size in [minGSSize, maxGSSize] will be tested.#也就是说如果我们想要我们GSEA富集分析的结果更加广泛，我们可以调节上面这两个参数#可视化ridgeplot(ego3,showCategory=20)

#注意这里默认显示的是p值最小的，因为我们富集分析得到的结果也是按照p值有小到大进行排序的

晨曦解读

山脊图将显示GSEA富集类别的核心富集基因的表达分布。它帮助用户解释向上/向下的路径。

running score and preranked list of GSEA result

当然我们这里也可以展示我们最常见的GSEA展现的形式

gseaplot2(ego3, geneSetID = 1:3, pvalue_table = TRUE, color = c("#E495A5", "#86B875", "#7DB0DD"), ES_geom = "dot")

至此，富集分析基本上所有的可视化展现形式就给大家介绍到这里，当然了我们已经了解了富集分析结果的存储形式，我们完全可以提取富集分析的结果然后用外在R包去碰撞出无限可能，在此也感谢Y叔可以为我们提供如此优秀的R包，不管是配色还是可视化的展现形式都十分的优秀！

你以为结束了？

哦不，并没有，下面我们将介绍一下基于clusterProfiler包的一些操作技巧，以及如何让我们富集分析更加自由的技巧

富集分析技巧讲解

其实这块主要是clusterProfiler包升级后所带来的更丰富的功能，因为开放了更多接口，所以富集分析结果可以直接与dplyr和ggplot2对接，对图形可视化非常友好，且Y叔考虑到了下游分析，所以直接给对象整体赋予了数据框的性能，也就是说可以直接从富集分析的结果中直接调用行、列等信息

注意：并不是说富集分析的结果是数据框，而是可以参考数据框的提取子集、列等操作，其实更像是列表那样，加上双取子集符号即可以完成提取

提问：如何筛选满足条件的富集分析结果？

filter(ego, p.adjust < .05, qvalue < 0.2)

提问：如何按照富集比例（GeneRatio）降序排列富集结果？

mutate(ego, geneRatio = parse_ratio(GeneRatio)) %>% arrange(desc(geneRatio))

提问：如何剔除掉富集分析的一些结果（变量）或者提取某些内容？

#剔除select(ego, -geneID) %>% head#提取ego[1:2,c("ID", "Description", "pvalue", "p.adjust")]head(ego$Description)ego[["GO:0072686"]]

提问：除了常见的GeneRation和BgRatio，还有没别的富集分析评价指标？

#rich factor: the ratio of input genes (eg,DEGs) that are annotated in a term to all genes that annotated in this term.#某一术语注释的输入基因与用该术语注释的所有基因的比率y <- mutate(ego, richFactor = Count / as.numeric(sub("/\\d+", "", BgRatio)))#然后我们还可以绘制可视化展示library(ggplot2)library(forcats)library(enrichplot)ggplot(y, showCategory = 20, aes(richFactor, fct_reorder(Description, richFactor))) + geom_segment(aes(xend=0, yend = Description)) + geom_point(aes(color=p.adjust, size = Count)) + scale_color_viridis_c(guide=guide_colorbar(reverse=TRUE)) + scale_size_continuous(range=c(2, 10)) + theme_minimal() + xlab("rich factor") + ylab(NULL) + ggtitle("Enriched Disease Ontology")

提问：如何再次筛选富集分析的结果？

#方法一down_ego <- simplify(ego, cutoff=0.7, by="p.adjust", select_fun=min)#方法二函数 dropGO 可以在由 enrichGO 得到的结果中移除特定的 GO term 或 GO level.

提问：如何调整我们富集分析的结果顺序？

#准备工作library(forcats)library(ggplot2)library(ggstance)library(enrichplot)library(dplyr)#因为新版的富集分析R包支持直接使用哈德雷大神的R包体系，所以我们可以直接使用#按照富集分数的大小排列由大到小排列y <- arrange(ego3, abs(NES)) %>% group_by(sign(NES))%>% slice(1:5)#绘图ggplot(y, aes(NES, fct_reorder(Description, NES), fill=qvalues), showCategory=10) + geom_col(orientation='y') + scale_fill_continuous(low='red', high='blue', guide=guide_colorbar(reverse=TRUE)) + theme_minimal() + ylab(NULL)