简单!目标基因怎么挑?这个小众分析,堪称高分文章标配!一刻钟包你学会!你一定要试试!
差异基因富集方法之Friends分析
大家好,我是阿琛。今天,我们来给大家介绍一种常见的拓展生信分析方法。在常规的⽣信文章中,我们往往能够得到⼀大堆的显著差异表达基因。然而,从⾥面挑选哪⼀个基因出来进行验证常常让我们感到困扰。通常,我们会从这些特定的蛋白功能信息出发,根据差异基因的条⽬进⾏富集分析,查找与其功能相似或相关的蛋白质,关联、比较、量化,看看是否富集在某个GOterm或者KEGG通路当中。没错,这就是我们常说的功能富集分析。
但是,如果富集到的某个我们感兴趣的通路中基因数目依然比较多,那么如何从这⼀堆基因中挑选最重要的那个就是⼀个问题。哪些基因会比较重要呢?是基因的差异改变程度吗?然而,基因的Friends分析对此进行了重新定义。假设一个基因与通路上的其它基因都存在相互作用的话,那么这个基因就可能⽐较重要,可能是所谓的hub-gene。
下面,我们来看下如何根据基因信息进行Friends分析。
1.1 R包安装与加载
#BiocManager::install("org.Hs.eg.db")
#BiocManager::install("GOSemSim")
#install.packages("ggplot2")
###R包加载
library(org.Hs.eg.db)
library(GOSemSim)
library(reshape2)
library(ggplot2)
这里,我们使用GOSemSim包进行差异蛋白的Friends分析;没错,这个包是由余叔在2010年发表于Bioinformatics杂志上,题目为GOSemSim: an R package for measuring semantic similarity among GO terms and gene products,并将该包储存于Bioconductor网站上。大家在后续使用这个R包进行后续分析时记得引用该文章。
1.2 数据读取
rt <-read.table("gene.txt", sep="\t", header = T,check.names = F)
str(rt)
接下来,我们将基因列表读取进来;通过str()函数,可以看到,其中包含了8个不同基因的相关信息。
rt$ENTREZID <- as.character(rt$ENTREZID)
head(rt)
同时,使用as.character()函数将数据框中ENTREZID号转换成字符串格式。
bp <- godata('org.Hs.eg.db', ont="BP", computeIC = FALSE)
cc <- godata('org.Hs.eg.db', ont="CC", computeIC = FALSE)
mf <- godata('org.Hs.eg.db', ont="MF", computeIC = FALSE)
首先,简单介绍一下,与大家熟悉的一样,常规的GO分析分为三种,分别为生物过程(biological process, BP),细胞组分(cellular component, CC),以及分子功能(molecular function, MF)。在此,我们同样从BP,CC和MF三个层面来进行分析。通过godata()函数来构建来构建相应物种三个不同方面的GOSemSim所需的GO注释数据。
simbp<- mgeneSim(rt$ENTREZID,
semData = bp,
measure = "Wang",
drop = NULL,
combine = "BMA")
simcc<- mgeneSim(rt$ENTREZID,
semData = cc,
measure = "Wang",
drop = NULL,
combine = "BMA")
simmf<- mgeneSim(rt$ENTREZID,
semData = mf,
measure = "Wang",
drop = NULL,
combine = "BMA")
随后,根据前面得到的注释数据,使用mgeneSim()函数,计算它们之间的语义相似度,最终得到它们之间的语义相似度和相应的GO语义相似度。
fsim <- (simmf * simcc * simbp)^(1/3)
基于分析得到的相似度结果,我们进一步计算这些基因在BP,CC和MF层面的几何平均值,得到最终的评分,以用于综合评估该基因在三个层面的功能。
colnames(fsim) = rt$SYMBOL
rownames(fsim) = rt$SYMBOL
由于前面的分析是基于ENTREZID号进行的。为了便于后续的分析,我们将基因名从ENTREZID改为Gene Symbol。
for (i in1:ncol(fsim)){
fsim[i,i] <- NA
}
dat <- melt(fsim)
dat <- dat[!is.na(dat$value),]
dat <- dat[,c(1,3)]
head(dat)
随后,进一步去除基因与本身基因之间的相关性,并使用melt()函数将宽格式数据转化成长格式数据。
dat.mean <- aggregate(value~Var1, dat, mean)
由于不同基因之间均可计算得到相似性评分,因此,我们使用aggregate()函数需要计算所有评分的平均数。
m <- dat.mean$value
names(m) <- dat.mean$Var1
#按平均值给基因名排序
dat$Var1 <- factor(dat$Var1,
levels=names(sort(m)))
str(dat)
同时,根据相似性评分的平均值,对其进行排序,并根据评分的高低,将基因名设置为因子(factor)格式。
ggplot(dat,
aes(x = Var1, y = value, fill = factor(Var1))) +
scale_fill_brewer(palette="Set2") +
geom_boxplot() +
coord_flip() +
xlab("") + ylab("") +
theme_bw()
最终,使用ggplot()函数对Friends分析结果进行可视化展示。结果显示,基因SIPA1在其中发挥重要功能的基因,可能是重要的基因。
好啦,今天对内容就分享到这里了~大家可以根据文介绍的内容,给自己到文章进一步增添一张美图。
回复“阿琛51”得到本次的代码和数据。
系列传送门
R语言小白入门课|一刻钟带你学会R数据转化
—
END—
撰文丨阿 琛
排版丨四金兄
值班 | 弘 毅
主编丨小雪球
最新评论
推荐文章
作者最新文章
你可能感兴趣的文章
Copyright Disclaimer: The copyright of contents (including texts, images, videos and audios) posted above belong to the User who shared or the third-party website which the User shared from. If you found your copyright have been infringed, please send a DMCA takedown notice to [email protected]. For more detail of the source, please click on the button "Read Original Post" below. For other communications, please send to [email protected].
版权声明:以上内容为用户推荐收藏至CareerEngine平台,其内容(含文字、图片、视频、音频等)及知识版权均属用户或用户转发自的第三方网站,如涉嫌侵权,请通知[email protected]进行信息删除。如需查看信息来源,请点击“查看原文”。如需洽谈其它事宜,请联系[email protected]。
版权声明:以上内容为用户推荐收藏至CareerEngine平台,其内容(含文字、图片、视频、音频等)及知识版权均属用户或用户转发自的第三方网站,如涉嫌侵权,请通知[email protected]进行信息删除。如需查看信息来源,请点击“查看原文”。如需洽谈其它事宜,请联系[email protected]。