通过GO语义找到相似的基因

生信喵

背景

输入基因列表，用GO注释的相似性找它的好朋友，画出云雨图。

出自PMID: 23182705文章

应用场景

在⽣信分析中，⼀个常⻅的应⽤场景就是：从特定蛋⽩质的功能信息出发，查找与其功能相似或相关的蛋⽩质，并对这些蛋⽩质间的关联程度进⾏⽐较、量化。这是生物信息学中经常遇到的问题，GO语义相似度为这种分析提供了可能。

通常认为，如果两个基因产物的功能相似，那么他们在GO中注释的术语（term），在GOtree中所处的位置就⽐较相近，反映在语义相似度上，就是他们的语义相似度比较⾼。

在我们⽇常的分析实践中，经常能够拿到⼀大堆的差异表达基因，从⾥面挑选哪⼀个基因出来进行验证常常让我们感到困扰。通常我们会对差异基因的条⽬进⾏富集分析，看看是否富集在某个GOterm或者KEGG通路当中。这时候已经对结果进⾏了相当程度上的清晰化。

但是如果富集到的某个我们感兴趣的通路中基因数⽬依然⽐较多，那么如何从这⼀堆基因中挑选最重要的那个就是⼀个问题。哪些基因会比较重要呢？

• 第⼀个线索是基因的差异改变的程度⽐较大，但差异改变程度⼤并不一定代表重要。
• 第⼆个线索就是该基因的产物与通路上的其它基因产物都有互作的话。简⽽言之，该基因编码蛋⽩的“朋友”⽐较多的话，那么该基因就可能⽐较重要。

跳出差异基因的范畴，任何⼀个涉及到基因互作的场景，例如打质谱拉下来的⼀堆蛋⽩、酵⺟双杂交得到的⼀堆蛋⽩，都可以通过类似⽅法找到最重要的那个hub基因。

环境设置

使用国内镜像安装包

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("org.Hs.eg.db", version = "3.8")
BiocManager::install("GOSemSim", version = "3.8")
install.packages("ggplot2")

加载包

library(org.Hs.eg.db)
library(GOSemSim)
library(reshape2)
library(ggplot2)

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

输入文件

输入的数据为基因名，第一列是ENTREZ ID, 第二列是gene symbol

其中ENTREZ ID用来计算相似性，gene symbol用于画图

id_gsym <- read.csv("easy_input.csv")
id_gsym$ENTREZID <- as.character(id_gsym$ENTREZID)
head(id_gsym)

计算相似性

原文用Molecular Function和Cellular Component的几何平均值来衡量相似度。

#用godata()函数来构建相应物种的Molecular Function本体的GO DATA
mf <- godata('org.Hs.eg.db', ont="MF", computeIC = FALSE)
#用godata()函数来构建相应物种的Cellular Component本体的GO DATA
cc <- godata('org.Hs.eg.db', ont="CC", computeIC = FALSE)
#用godata()函数来构建相应物种的Biological Process本体的GO DATA
#bp <- godata('org.Hs.eg.db', ont="BP", computeIC = FALSE)

#用mgeneSim来计算MF本体，基因之间的语义相似度，结果为一个行列相同的矩阵
simmf <- mgeneSim(id_gsym$ENTREZID, semData = mf, measure = "Wang", drop = NULL, combine = "BMA")
#用mgeneSim来计算CC本体，基因之间的语义相似度，结果为一个行列相同的矩阵
simcc <- mgeneSim(id_gsym$ENTREZID, semData = cc, measure = "Wang", drop = NULL, combine = "BMA")
#用mgeneSim来计算BP本体，基因之间的语义相似度，结果为一个行列相同的矩阵
#simbp <- mgeneSim(id_gsym$ENTREZID, semData = bp, measure = "Wang", drop = NULL, combine = "BMA")

#计算基因在MF本体和CC本体下的几何平均值，一个打分值同时包括基因的分子功能和细胞定位两个信息
fsim <- sqrt(simmf * simcc)
#或者计算基因在MF、CC、BP本体下的几何平均值
#fsim <- (simmf * simcc * simbp)^(1/3)

声明！
由于GO注释在不断的完善，而当年Y叔发文章是2012年，本次测试是2025年，这几年GO数据库也在不断完善。所以用现在的GO注释算出来的结果画出来图和原文会有所不同，这不是原文的错误，只是对原文进行修正。这一次算出来的结果也更加的make sense，因为用抗FMNL1抗体拉下的蛋白中，理论上FMNL1应该是“朋友”最多的那个。

#如果想要完全重复Y叔文章里的图，请从这里导入当年计算的fsim。
#(load("fsim.rda"))

#将基因的名字由ENTREZID改为gene symbol，方便看懂。
colnames(fsim) = id_gsym$SYMBOL
rownames(fsim) = id_gsym$SYMBOL

#将基因自己和自己的相似度设为NA，方便接下来去掉。
for (i in 1:ncol(fsim)){
  fsim[i,i] <- NA
}

y <- melt(fsim) #把宽格式数据转化成长格式，其实就是把正方形矩阵转成三列
y <- y[!is.na(y$value),] #删掉带NA的行

# 把每两个基因之间的相似度保存到文件，只需要保存第一列基因名和第三列数值
write.csv(y[,c(1,3)], "very_easy_input.csv", row.names = F)

开始画图

画出原文的box plot

y <- read.csv("very_easy_input.csv")
head(y)

#计算每个基因跟其他基因相似度的平均值
y.mean <- aggregate(.~Var1,y,mean) 
m <- y.mean$value
names(m) <- y.mean$Var1
#按平均值给基因名排序，便于画图
y$Var1 <- factor(y$Var1, levels=names(sort(m)))

f <- function(y) {
  r <- quantile(y, probs = c(0.05, 0.25, 0.5, 0.75, 0.95))
  r[3] <- mean(y)
  names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
  r
}

p1 <- ggplot(y, aes(Var1, value, fill = factor(Var1))) + 
  scale_fill_brewer(palette="Set3") + #配色
  guides(fill=FALSE) + #不显示图例

  stat_summary(fun.data= f, geom='boxplot') + 
  geom_hline(aes(yintercept=0.75), linetype="dashed") + #画一条虚线

  coord_flip() + # x、y坐标轴互换
  xlab("") + ylab("") + 
  theme(axis.text.x = element_text(family = "Arial", size = 16, face = "bold"),
                          axis.text.y = element_text(family = "Arial", size = 16, face = "bold")) + 
  theme_bw() + 
  theme(panel.border=element_rect(size=1)) #边框粗细 
p1

# 保存到pdf文件
ggsave("friends_box.pdf")

补充：画成云雨图

有小伙伴想知道云雨图的输入数据格式，刚好用这套数据来做示例。

云雨图的输入只需要两列：

• 第一列是分组名称，此处是基因名，每个基因为一组
• 第二列是组内成员的数值，此处是每个基因跟其他基因相似度的值

#devtools::install_github("GuangchuangYu/gglayer")
require(gglayer)

y <- read.csv("very_easy_input.csv")
head(y)
unique(y$Var1)

#计算每个分组的平均值
y.mean <- aggregate(.~Var1,y,mean) 
m <- y.mean$value
names(m) <- y.mean$Var1
#按平均值给分组排序，便于画图
y$Var1 <- factor(y$Var1, levels=names(sort(m))) 
source('geom_flat_violin.R')
p2 <- ggplot(y, aes(Var1, value, fill = Var1)) +
  #scale_fill_brewer(palette="Set2") + #配色
  guides(fill=FALSE) +
  geom_flat_violin(position=position_nudge(x=.2)) +

  #分散不重叠的点图
  #geom_jitter(aes(color=Var1), width=.15) + guides(color=FALSE) +
  #堆叠的点图
  geom_dotplot(binaxis="y", stackdir="down", dotsize=.35) +

  geom_boxplot(width=.1, position=position_nudge(x=.1)) +
  geom_hline(aes(yintercept=0.75), linetype="dashed") + #画一条虚线

  coord_flip() + # x、y坐标轴互换
  xlab("") + ylab("") + 
  theme(axis.text.x = element_text(family = "Arial", size = 16, face = "bold"),
                          axis.text.y = element_text(family = "Arial", size = 16, face = "bold")) + 
  theme_bw() + 
  theme(panel.border=element_rect(size=1)) #边框粗细 

p2
ggsave("friends_raincloud.pdf")

组个图，对比两种展示效果

library(cowplot)
plot_grid(p1, p2, labels = c("A", "B"), align = "h")
ggsave("friends.pdf")

生信喵

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

发送关键词“20250501”获取