如何手动做GO分析（包含p值，p.adj的计算）

最新推荐文章于 2024-07-03 18:00:31 发布

qq_26309777

最新推荐文章于 2024-07-03 18:00:31 发布

阅读量1.8w

点赞数 6

文章标签：大数据编程语言生物学

本文链接：https://blog.csdn.net/qq_26309777/article/details/106818657

版权

本文深入探讨RNA-seq技术在生物学研究中的应用，特别是通过R语言进行GO分析的过程。从理解GO分析原理出发，讲解如何统计基因组中GO术语的频数，以及如何计算差异表达基因在特定GO术语上的显著性。此外，还介绍了p值和调整后的p值（padj）的计算方法，以及如何判断差异基因在GO术语上的富集情况。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

引言

不论你是主攻基因的上下游关系，还是想鉴定在外界刺激下生物体的响应基因，亦或是仅仅想凑一个文章毕业，只要是涉及到活细胞，大都会涉及到RNA-seq。

RNA-seq的测序技术本身已经很成熟，同样一批处理过的样品，在不发生意外事故的前提下，不论送哪个公司，测序结果都大同小异。在这个前提成立的基础上，文章的层次主要由分析手法的高低决定。在一些高手手中，除了一些生物公司就会给做的QC之外，还要进行GO 分析，KEGG分析，PPI网络和代谢网络分析，因果推理，还会结合实验设计的目的和已有的文献，进行合理的推论，等等。在下面的内容中，我们就以R语言为工具，来看一看如何一步一步地做GO分析。

GO分析的网站有很多，笔者见过并且使用过的就不下十个。动物植物通用的有：DAVID，PATHER；植物特有：（AgriGo ，PlantGSEA）。然而，理想很美好，现实却不一定。莎翁说过：1000个人眼中的1000个哈姆雷特。同一批差异基因提交到不同的网站后生成的结果差别巨大，有些甚至彼此冲突。到底哪一种才是对的呢？

在回答这个问题之前，我们先看看GO分析的原理：在基因组中，每一个基因往往参与了1到多个GO（biological process, molecular function, cellular component）中，以人类P53基因为例（该基因至少参与到34个molecular function和125个biological process中）；而每个GO又被1到多个基因所参与，以ATP binding（GO:0005524）为例，有上百个基因参与其中，包括：Acute AML1、Protein translocase subunit SecA 、ATPase、RNA-directed RNA polymerase等等。我们假设我们的RNA-seq结果找到了k个上调差异表达基因，正确GO的方法应该是这样：

对于每一个GO，先后统计出

①统计全基因组上每个GOterm出现的频数，定义为m;

②统计全基因组上基因数量，定义为n;

③统计本次GO分析中提交差异基因的数量，定义为k，这里就是上面的k;

④统计在k个上调差异基因中的GO频数，定义为o。

在R语言中，可以用以下代码来实现p值的计算：

p[i,1]=phyper(o-1, k, n-k, m, lower.tail=FALSE)

p值计算出来之后，还要计算p-adjusted，即校正后的p值（qvalue=padj=FDR=Corrected p-Value=p-adjusted），是对p值进行了多重假设检验，能更好地控制假阳性率。

在R语言中，可以用

padj=p.adjust(p, method = "BH", n = length(p))

来实现。

对于这个GO，倘若p.adj<0.05。那就说明差异基因在这个GO term上显著富集。