主讲老师第二十九课：GO分析
得到差异基因以后，常规的操作自然是进行功能的富集分析，关于功能分析，主要分为GO 分析，KEGG分析，以及 GSEA 分析，首先，我们来看一下 GO 分析（GO.zip）GO 是基因本体论联合会建立的一个数据库，旨在建立一个适用于各种物种的、对基因和蛋白功能进行限定和描述的，并能随着研究不断深入而更新的语义词汇标准。GO 注释分为三大类：
分子生物学功能（Molecular Function，MF）生物学过程（Biological Process，BP）细胞学组分（Cellular Components，CC）通过这三个功能大类，对一个基因的功能进行多方面的限定和描述，GO 功能分类是在某一功能层次上统计蛋白或者基因的数目或组成，此外也有研究都挑选一些 Term，而后统计直接对应到该 Term 的基因或蛋白数，对于GO 分析，除了 DAVID 等在线分析网站外，我们主要介绍一下如何使用 R 来进行相应的富集分析。
在分析之前，我们需要将差异分析所得到的结果进行导入。
与前面可视化部分一样，这里还是使用 Deseq2 包分析得到的结果
1.添加 ENTREZID 列
在得到的差异分析结果，基因的名字主要以 Gene symbol 的形式，做 GO 分析是不能直接用基因名的，必须得先转化成 entre id。
简单介绍一下几种常用的 ID：
Ensemble id ：一般以 ENSG 开头，后边跟 11 位数字Entrez id ：通常为纯数字。如 TP53 基因：7157
Symbol id ：为我们常在文献中报道的基因名称
Refseq id ：NCBI 提供的参考序列数据库，可以是 NG、NM、NP 开头，代表基因，转录本和蛋白质。
首先，读取需要的 R 包，使用 clusterProfiler 包的 bitr()函数进行基因 id 的转换，此外 org.Hs.eg.db是人类基因组注释的 R 包这里，与前面介绍的方法不一样，介绍一个简单的函数，bitr()函数。来看一下转换函数的相关参数：
geneID ：输入待转换的 geneID； fromType ：输入的 ID 类型； toType ：希望输出的 ID 类型；
OrgDb ：注释对象的信息，此处用的是人类注释信息
根据提示可以看到，在转换过程中，有 9.75%的基因未能匹配上，其中 3389 个基因成功进行了转换。
注意，关于参数 fromType 和 toType 的转换类型，这里是进行了限定的，必须是以下的几种：
接着，我们使用 dplyr 包中的 inner_join()函数来将两个结果取交集，其中 DEG 中的 Gene和 genelist中的 SYMBOL 是同一类信息。
最终，经过合并，最后得到 3389 个基因的信息。
随后，选取差异分析列表中显著上调和显著下调基因的 entre id，形成一个新的向量，用于后续的分析。
2.GO 分析
2.1 GO 分析
其中，参数 OrgDb 指定该物种对应的 org 包的名字，ont 代表 GO 的 3 大类别，BP, CC, MF，也可以是全部的 all，cutoff 指定对应的阈值，pAdjustMethod 指定多重假设检验矫正的方法，readable=TRUE代表将基因 ID 转换为 gene symbol。此外，GO 分析也可以在 BP(生物过程），MF(分子功能），CC(细胞组分）三个角度进行注释
2.2 细胞组分
2.3 生物过程
2.4 分子功能
最终，分别得到了全部的 all，以及三种不同生物学功能所对应的结果，分别将其进行输出保存。
3.可视化
接下来，我们介绍几个常见的函数用于结果的可视化展示：
barplot() #富集柱形图dotplot() #富集气泡图cnetplot() #网络图展示富集功能和基因的包含关系emapplot() #网络图展示各富集功能之间共有基因关系heatplot() #热图展示富集功能和基因的包含关系
3.1 柱状图
3.2 气泡图
3.3 分类展示
将三个 GO 层次分开来展示。
3.4 网络图
对差异基因按 logFC 进行排序。
3.5 网络图 2
3.6 富集分布图
3.7 热图
3.8 upset 图
关于 GO 分析及其可视化结果，先介绍这几种方法。
那些年，我们曾经见过的富集化分析结果，你确定都见过吗？
Hi，大家好，我是主讲老师，今天是主讲老师碎碎念的第六期，这期推文的灵感来源于以下
和学员的对话：
A 同学：主讲老师，富集分析有多种展现方式，我该如何选择呢？
主讲老师：可视化的本质是数据，挑选自己觉得好看同时杂志社要求的格式即可A 同学：那么曦曦，能否整理出一篇笔记，告诉我们富集分析究竟有多少种展现形式？
主讲老师：好的，关注下期推文哦~
本篇推文基于以上问题诞生，将聚焦在富集分析究竟有多少种展现形式，书写这篇推文的同时也是主讲老师对于富集分析这一方面知识的又一次梳理和总结，那么，我们就开始吧~干货预警！！预计阅读时间 10min，建议配上音乐一起食用哦~
富集分析的相关概念
这里并不准备用太多官方的语言来给大家进行解释，大家只需要知道，我们筛选出来的差异基因（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:0000779CENPE/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:0000776CENPE/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 富集分析的结果和这个大同小异，掌握一个即可类推到另一个，那么我们下面进行富集分析可视化的展示
富集分析可视化展示
1.Bar Plot
barplot(ego, showCategory=10)
2.Dot plot
dotplot(ego, showCategory=10) + ggtitle("dotplot for ORA")
3.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')
5.Heatmap
p1 <- heatplot(edox, showCategory=5)p2 <- heatplot(edox, foldChange=geneList, showCategory=5)cowplot::plot_grid(p1, p2, ncol=1, labels=LETTERS[1:2])
6.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 包版
本如下：
7.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 的结果
8.UpSet Plot
可以使基因和基因组之间的复杂关联可视化。它强调不同基因组之间的基因重叠。
upsetplot(edox)
9.ridgeline plot for expression distribution of GSEA result
当然上面的可视化大部分都是用来展示 GO 和 KEGG 的Ps：其实我们常用的富集分析 GO 和 KEGG，这个名称是按照注释来源的数据库来命名的，除此之外，我们如果更换注释数据库，包括 MeSH、Reactome、DOSE 等这类数据库也会得到不同的注释结果，好在 Y 叔把以上不同注释来源的数据库写成了不同的函数，我们直接选择即可当然目前来看富集分析的算法就是两种即 Over Representation Analysis（ORA）和 Gene SetEnrichment 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 富集类别的核心富集基因的表达分布。它帮助用户解释向上/向下的路径。
10.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 ina 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)至此，本篇推文到这里就结束啦，零零散散写了也有接近 1w5 左右了，能看到这里的小伙伴大大点赞哦~主讲老师碎碎念本质上来说是主讲老师的学习笔记，所以各位小伙伴有感兴趣的内容也欢迎在评论区留言哦~最后，如果大家用这款 R 包进行富集分析别忘了引用下面的文献哦clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. *TheInnovation*. 2021, 2(3):100141. doi: 10.1016/j.xinn.2021.100141我是主讲老师，我们下期再见~Ps：回复“主讲老师碎碎念 06”可以获得本篇推文的代码以及参考文献哦~