GO富集结果聚类

生信喵

背景

富集分析得到太多GO term，其中好多还是相似的，怎样合理的合并？clusterProfiler有一个simplify函数，能给富集分析结果瘦身by removing redundancy of enriched GO terms，但有时瘦的不够多。例文列出了富集的 term，同时在旁边给出一个短语来概括性的描述这几个相似的term：

出自PMID: 28726821文章

方法探讨：

Extended Data Figure 5 | Candidate driver genes and pathways in MB subgroups. c, GO and pathway summary of recurrently mutated genes in MB. GO and pathway categories are grouped according to functional theme and the proportion of cases affected by individual pathway alterations are plotted per subgroup and across the series.

例文中关于这些相似term如何合并的描述只说grouped according to functional theme。因此，用Y叔的GOSemSim计算GO term之间的相似性，然后借助ggtree实现例文中相似GO germ的标注效果。

应用场景

多组GO富集分析结果放在一起对比；

富集分析得到太多相似的GO term，需要合并相似的，再结合背景知识概括出合适的pathway。

环境设置

使用国内镜像安装包

options("repos"= c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("GOSemSim", version = "3.8")
BiocManager::install("ggtree", version = "3.8")

加载包

library(plyr)
library(stringr)
library(ape)
library(GOSemSim)
library(ggtree)
library(scales)
library(cowplot)

Sys.setenv(LANGUAGE = "en") #显示英文报错信息
options(stringsAsFactors = FALSE) #禁止chr转成factor

输入文件

富集分析结果至少包含GO term ID（格式为GO:0046777）和一列统计量（例如pvalue、p.ajdust、FDR等等）。

此处用clusterProfier的enrichGO函数做的富集分析结果为例，其中有一列BgRatio，斜线左侧是全基因组范围内该term里基因的数量。例文以300为阈值，去掉了基因数量超过300的term。如果你的数据没有这一列，也可以不做这一步筛选。

#多组富集分析结果的文件名统一命名，便于导入
fnames <- Sys.glob("enrichGO*.csv")
fnames
#读取ID、Description、BgRatio、p.adjust四列
fdataset <- lapply(fnames, function(x){read.csv(x)[,c(2,3,5,7)]})
names(fdataset) <- fnames

纵向合并多组富集结果，用于给GO term分类

ego.all <- ldply(fdataset, data.frame)
ego.all$group <- unlist(strsplit(ego.all$.id, split = ".csv"))
head(ego.all)
dim(ego.all)

#用p.ajust<0.001来筛
ego.all <- ego.all[ego.all$p.adjust < 0.001,] 
dim(ego.all)

# GO term的合集
ego.ID <- unique(ego.all[,c(2:4)])
head(ego.ID)
dim(ego.ID)

##GO term筛选
#删掉超过300个基因的GO term
#如果不想删除，可以注释掉下面这段
ego.ID$Bg <- as.numeric(str_split_fixed(ego.ID$BgRatio, "/",2)[,1]) #提取GO term包含的基因数量
ego.ID <- ego.ID[ego.ID$Bg < 300,] 
dim(ego.ID)
#删掉少于100个基因的GO term
#具体怎样筛选，要根据你自己的数据灵活设置
ego.ID <- ego.ID[ego.ID$Bg > 100,] 
dim(ego.ID)
head(ego.ID)

横向合并多组富集结果，用于画热图

MyMerge <- function(x, y){
  df <- merge(x, y, by= "ID", all.x= TRUE, all.y= TRUE)
  return(df)
}
ego.m <- Reduce(MyMerge, fdataset)
head(ego.m)

#只保留GO term ID和各组的p.adjust
ego.m <- ego.m[,c(1,4,7,10)] #此处有三组，如果有更多组，要在后面继续加
#提取筛选过的GO term
ego.m <- merge(ego.ID[,1:2], ego.m, by= "ID", all.x= TRUE)
rownames(ego.m) <- ego.m$Description
ego.m$ID <- NULL
ego.m$Description <- NULL

#把列名改为组名
#colnames(ego.m)<- str_remove(fnames, ".csv")
colnames(ego.m) <- paste0("G", seq(1:length(fnames)))

#用全部GO term做富集分析，多组合并后会有缺失值
#解决方案是DIY富集分析
#不推荐的做法是用1补，或者换成其他你认为合理的数值
#ego.m[is.na(ego.m)] <- 1 
head(ego.m)

按GO term的相似性聚类

# 人
#hgGO <- godata('org.Hs.eg.db', ont="BP")
#save(hgGO, file="hgGO.rdata")
#(load("hgGO.rdata"))

# 小鼠
#这一步需要运行1min，可以把它保存到mmGO.rdata文件，便于重复使用
#mmGO <- godata('org.Mm.eg.db', ont="BP") ##或MF或CC，要跟你富集分析时的参数一致
#save(mmGO, file="mmGO.rdata")
(load("mmGO.rdata")) #导入之前保存的文件

#计算GO term之间的相似性
ego.sim <- mgoSim(ego.ID$ID, ego.ID$ID, semData=mmGO, measure="Wang", combine=NULL)
ego.sim[1:3, 1:3]

#用GO term作为行名、列名，便于查看和画图
rownames(ego.sim) <- ego.ID$Description
colnames(ego.sim) <- ego.ID$Description
ego.sim[1:3, 1:3]

#给GO term分类
tree <- nj(as.dist(1-ego.sim))
p <- ggtree(tree) + geom_tiplab() + #写GO term
  geom_text2(aes(subset=!isTip, label=node), hjust=-.3) + #写node编号
  coord_cartesian(xlim=c(-.1,1.3)) #左右两侧留出合适的空间
p

开始画图

**注：**根据实际情况，调整offset的值

#结合树的结构和背景知识，此处把GO term分为5类，记下每一类的node编号，写在node参数里
node <- c(21,22,24,27,28)
gtree <- groupClade(tree, .node=node)
#用不同颜色展示不同类GO term
pbase <- ggtree(gtree, 
                aes(color=group)) #每类用不同颜色画树枝

#给5类term分别总结出一个短语，作为分类名，标在树的旁边
fontsize <- 4 #字的大小
offset <- .9 #分类名向右移动到热图右侧，也可以设为0.3，让它显示在热图跟树之间
pnode <- pbase + 
  #如果不想显示每个GO term，就不运行这行，同时offset <- 0
  geom_tiplab(size=4, align=TRUE) + #GO germ对齐
  geom_cladelabel(node=node[1], align=TRUE, 
                  #文字颜色跟树枝颜色一致，也可以删掉下面这行，就会像例文那样全是黑色 
                  color = hue_pal()(length(node)+1)[2], 
                  fontsize = fontsize, offset=offset, label="pathway1") +
  geom_cladelabel(node=node[2], align=TRUE, color = hue_pal()(length(node)+1)[3], fontsize = fontsize, offset=offset, label="pathway2") +
  geom_cladelabel(node=node[3], align=TRUE, color = hue_pal()(length(node)+1)[4], fontsize = fontsize, offset=offset, label="pathway3") +
  geom_cladelabel(node=node[4], align=TRUE, color = hue_pal()(length(node)+1)[5], fontsize = fontsize, offset=offset, label="pathway4") +
  geom_cladelabel(node=node[5], align=TRUE, color = hue_pal()(length(node)+1)[6], fontsize = fontsize, offset=offset, label="pathway5") +
  #如果有更多分类，就继续往下粘贴
  
  coord_cartesian(xlim=c(-.1,1.5)) 
# 把p.ajust画在树形结构右侧
p2 <- gheatmap(pnode, ego.m, 
         offset=.7, #热图向右移动到合适的位置
         width=0.3, #heatmap格子的宽度
         colnames_angle=90, hjust=0, #组的名字竖着写
         low = "red", high = "white") 
attr(p2, "heatmap_mapping") <- mapping
p2

#保存到文件
ggsave("GOclustering.pdf", width = 12, height = 8)

最后生成的pdf文件是矢量图，可以在Illustrator或Inkscape中打开、编辑。

可以参考例文的布局，手动删掉树形结构或调整各部分的位置。

后记

因为ggplot2 4.0冲突，设计了https://github.com/YuLab-SMU/ggtree/issues/686，先把图拼起来。目前的树状图的label颜色丢失了，后续再想法子加上。

生信喵

文中所用数据可以关注公众号“生信喵实验柴”

发送关键词“20250508”获取