GO富集结果聚类

生信喵 发表于 2025-11-29 23:51:17

本帖最后由生信喵于 2025-12-8 16:18 编辑

<h1>背景</h1>
<p>富集分析得到太多GO term，其中好多还是相似的，怎样合理的合并？clusterProfiler有一个simplify函数，能给富集分析结果瘦身by removing redundancy of enriched GO terms，但有时瘦的不够多。例文列出了富集的 term，同时在旁边给出一个短语来概括性的描述这几个相似的term：<br />
<img src="https://tcapi.voiceclouds.cn:8444/uploads/2025/11/29/mikcns8u8nvagnpfco8.png" alt="" /></p>
<p>出自PMID: 28726821文章</p>
<p><strong>方法探讨：</strong></p>
<p>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 <strong>grouped according to functional theme</strong> and the proportion of cases affected by individual pathway alterations are plotted per subgroup and across the series.</p>
<p>例文中关于这些相似term如何合并的描述只说grouped according to functional theme。因此，用Y叔的GOSemSim计算GO term之间的相似性，然后借助ggtree实现例文中相似GO germ的标注效果。<br />
<img src="https://tcapi.voiceclouds.cn:8444/uploads/2025/11/29/mikcq1q7ccizom1t6i8.png" alt="" /></p>
<h1>应用场景</h1>
<p>多组GO富集分析结果放在一起对比；</p>
<p>富集分析得到太多相似的GO term，需要合并相似的，再结合背景知识概括出合适的pathway。</p>
<h1>环境设置</h1>
<p>使用国内镜像安装包</p>
<pre><code class="language-r">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")
</code></pre>
<p>加载包</p>
<pre><code class="language-{r}">library(plyr)
library(stringr)
library(ape)
library(GOSemSim)
library(ggtree)
library(scales)
library(cowplot)

Sys.setenv(LANGUAGE = "en") #显示英文报错信息
options(stringsAsFactors = FALSE) #禁止chr转成factor
</code></pre>
<h1>输入文件</h1>
<p>富集分析结果至少包含GO term ID（格式为GO:0046777）和一列统计量（例如pvalue、p.ajdust、FDR等等）。</p>
<p>此处用clusterProfier的enrichGO函数做的富集分析结果为例，其中有一列BgRatio，斜线左侧是全基因组范围内该term里基因的数量。例文以300为阈值，去掉了基因数量超过300的term。如果你的数据没有这一列，也可以不做这一步筛选。</p>
<pre><code class="language-{r}">#多组富集分析结果的文件名统一命名，便于导入
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
</code></pre>
<h1>纵向合并多组富集结果，用于给GO term分类</h1>
<pre><code class="language-{r}">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
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
dim(ego.ID)
#删掉少于100个基因的GO term
#具体怎样筛选，要根据你自己的数据灵活设置
ego.ID <- ego.ID
dim(ego.ID)
head(ego.ID)
</code></pre>
<h1>横向合并多组富集结果，用于画热图</h1>
<pre><code class="language-{r}">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 <- 1
head(ego.m)
</code></pre>
<h1>按GO term的相似性聚类</h1>
<pre><code class="language-{r,"># 人
#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

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

#给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
</code></pre>
<p><img src="https://tcapi.voiceclouds.cn:8444/uploads/2025/11/29/mikdorycbod1ruchoaa.png" alt="" /></p>
<h1>开始画图</h1>
<p>**注：**根据实际情况，调整offset的值</p>
<pre><code class="language-{r,">#结合树的结构和背景知识，此处把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, align=TRUE,
               #文字颜色跟树枝颜色一致，也可以删掉下面这行，就会像例文那样全是黑色
               color = hue_pal()(length(node)+1),
               fontsize = fontsize, offset=offset, label="pathway1") +
geom_cladelabel(node=node, align=TRUE, color = hue_pal()(length(node)+1), fontsize = fontsize, offset=offset, label="pathway2") +
geom_cladelabel(node=node, align=TRUE, color = hue_pal()(length(node)+1), fontsize = fontsize, offset=offset, label="pathway3") +
geom_cladelabel(node=node, align=TRUE, color = hue_pal()(length(node)+1), fontsize = fontsize, offset=offset, label="pathway4") +
geom_cladelabel(node=node, align=TRUE, color = hue_pal()(length(node)+1), 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, "mapping") <- NULL
p2

#保存到文件
ggsave("GOclustering.pdf", width = 12, height = 8)
</code></pre>
<p><img src="https://tcapi.voiceclouds.cn:8444/uploads/2025/11/29/mikgtycy0ro35xk3lt5.png" alt="" /><br />
最后生成的pdf文件是矢量图，可以在Illustrator或Inkscape中打开、编辑。</p>
<p>可以参考例文的布局，手动删掉树形结构或调整各部分的位置。</p>
<h1>后记</h1>
<p>因为ggplot2 4.0冲突，设计了https://github.com/YuLab-SMU/ggtree/issues/686，先把图拼起来。目前的树状图的label颜色丢失了，后续再想法子加上。</p>

生信喵 发表于 2025-11-29 23:53:12

<p>文中所用数据可以关注公众号“生信喵实验柴”<br />
<img src="https://tcapi.voiceclouds.cn:8444/uploads/2025/09/07/mf9nu72s90ofyseu9hl.png" alt="" /><br />
发送关键词“20250508”获取</p>

页: [1]

bioinfoer's Archiver

GO富集结果聚类