RNA-Seq(8):差异基因KEGG/GO富集分析

2022-01-14 本文已影响0人 Z_bioinfo

1.什么是功能富集分析

高通量的数据的分析，可以让我们得到很多候选的结果。但是如果只是把结果这样的平铺开的话，反正不利于我们去发现事情的本质。所以为了更情况的看清楚这些基因的功能，我们就使用了富集分析。我们可以把富集分析理解为在把很零零碎碎的东西，通过一个整体来反应出来，类似于从微观到宏观的变化

功能富集分析是将基因或者蛋白列表分成多个部分，即将一堆基因进行分类，而这里的分类标准往往是按照基因的功能来限定的。换句话说，就是把一个基因列表中，具有相似功能的基因放到一起，并和生物学表型关联起来。

利用富集分析，我们就可以把很多看着杂乱的差异基因总结出一个比较整体反应事件发生的概述性的句子。例如：TP53信号通路和胃癌的发生有关。而不是说BAX、BID、ABL1、ATM、BCL2、BOK、CDKN1A这7个基因和胃癌的发生有关系。

评估差异基因主要影响的生物学功能和通路。

2 什么是GO和KEGG

1.GO（找差异表达基因主要富集在哪些GO term（分子到生物过程，分三类：molecular function、cellular component、biological process）中，评估「差异表达基因与哪些生物学功能显著相关」，对生物学功能起上调还是下调作用。）

GO注释（GO annotations）库，它主要是为GO terms提供注释，也就是描述这个GO terms有什么功能（例如某些基因的产物是什么，是蛋白质，还是非编码RNA，还是大分子等）,GO注释分为三大类，
Cellular component，CC 细胞成分，细胞成分解释的是基因存在在哪里，在细胞质还是在细胞核？如果存在,那么在细胞质那在哪个细胞器上？如果是在线粒体中那是存在线粒体膜上还是在线粒体的基质当中？这些信息都叫Cellular component。

Biological process， BP 生物学过程，生物学过程是在说明该基因参与了哪些生物学过程，比如，它参与了rRNA的加工或参与了DNA的复制，这些信息都叫Biological process。用的比较多的是GO:BP

Molecular function，MF 分子功能，分子功能在讲该基因在分子层面的功能是什么？它是催化什么反应的？

通过这三个功能大类，对一个基因的功能进行多方面的限定和描述。
对基因的描述一般从三个层面进行：

做GO分析的思路：

control VS treatment ---> DEG ---> GO enrichment analysis

用找到的差异基因去做GO富集分析，希望能从这三方面找到和我们背景不一样的地方。比如，在疾病研究的时候，进行药物治疗之后某些基因的表达量明显的发生了变化，拿这些基因去做GO分析发现在Biological process过程当中集中在RNA修饰上，然后在此基础上继续进行挖掘。这个例子就是想启示大家拿到差异表达基因DEG只是一个开始，接下来就应该去做GO注释，之后需要进行一个分析看这些注释主要集中在哪个地方。假如我们有100个差异表达基因其中有99个都集中在细胞核里，那我们通过GO分析就得到了一个显著的分布。

代码部分

---------------拿什么数据来富集---------------
得到的差异表达基因的名字就可以，也就是说不需要其他的值
---------------用什么工具富集--------------
Y叔的clusterProfiler

通过第6步得到了251个差异基因，现在需要拿这些基因id名就可以去做GO分析，看这些注释主要集中在哪个地方

查看一下第6步得到的差异基因
> head(diff_gene,3)
DataFrame with 3 rows and 9 columns
  ensembl_gene_id  baseMean log2FoldChange     lfcSE      stat       pvalue
      <character> <numeric>      <numeric> <numeric> <numeric>    <numeric>
1 ENSG00000152583  1007.391        4.60333  0.211291   21.7867 3.10467e-105
2 ENSG00000179094   805.846        3.20082  0.197830   16.1796  7.02025e-59
3 ENSG00000125148  3587.633        2.18717  0.144117   15.1764  5.07003e-52
          padj external_gene_name            description
     <numeric>        <character>            <character>
1 5.49991e-101            SPARCL1 SPARC like 1 [Source..
2  4.14546e-55               PER1 period circadian reg..
3  1.79631e-48               MT2A metallothionein 2A [..
library(clusterProfiler) #用来做富集分析
library(DOSE)
library(topGO)#画GO图用的
library(pathview)  #看KEGG pathway的
library(org.Hs.eg.db)#这个包里存有人的注释文件
library(ggplot2)
GO富集分析
----------方法一：用差异基因的ensembl_gene_id（身份证号）进行分析的------
BP层面上的富集分析：
go_bp<-enrichGO(gene = diff_gene$ensembl_gene_id,OrgDb  = org.Hs.eg.db,keyType    = 'ENSEMBL', ont  = "BP", pAdjustMethod = "BH",pvalueCutoff = 0.01, qvalueCutoff = 0.05)
CC层面上的富集分析：
go_cc<-enrichGO(gene  = diff_gene$ensembl_gene_id,OrgDb  = org.Hs.eg.db,keyType    = 'ENSEMBL', ont = "CC", pAdjustMethod = "BH",pvalueCutoff = 0.01, qvalueCutoff = 0.05)
----------方法二：如果不想用差异基因的ensembl_gene_id（身份证号）进行分析，想通过基因的名字进行功能分析，请看以下代码------
先进行id到名字转换
library(clusterProfiler)
library(DOSE)
library(org.Hs.eg.db)
ensembl_gene_id = diff_gene[,1]
gene.name <-bitr(ensembl_gene_id, fromType = "ENSEMBL", 
              toType = c("SYMBOL","ENTREZID"),
              OrgDb = org.Hs.eg.db,drop =  FALSE )
BP层面上的富集分析：

go_bp<-enrichGO(gene = gene.name$SYMBOL,OrgDb  = org.Hs.eg.db,keyType    = 'SYMBOL', ont  = "BP", pAdjustMethod = "BH",pvalueCutoff = 0.01, qvalueCutoff = 0.05)
把结果导出保存
write.csv(go_bp@result,"go_bp.csv") 
CC层面上的富集分析：
go_cc<-enrichGO(gene  = gene.name$SYMBOL,OrgDb  = org.Hs.eg.db,keyType    = 'SYMBOL', ont = "CC", pAdjustMethod = "BH",pvalueCutoff = 0.01, qvalueCutoff = 0.05)

可视化
##分析完成后，作图，自己选择喜欢的一款,一般GO分析画这两个图就可以了
做成气泡图形式
dotplot(go_bp)
dotplot(go_cc)
做成条形图形式
> barplot(go_bp,showCategory = 12,title="The GO_BP enrichment analysis of all DEGs ", drop=TRUE)

气泡图：横坐标是GeneRatio，意思是说输入进去的基因，它每个term（纵坐标）占整体基因的百分之多少。圆圈的大小代表基因的多少，图中给出了最大的圆圈代表60个基因，圆圈的颜色代表P-value，也就是说P-value越小gene count圈越大，这事就越可信。

image.png

可以看到基因描述，与原文中找到的基因对比,选择CRISPLD2为新型的类固醇反应基因.我们分析了两项公开发表的基因表达微阵列研究的数据，这些研究测量了GCs对人类ASM细胞的影响，以确定这些先前的研究是否支持我们的CRISPLD2差异表达结果。尽管在这些研究中，CRISPLD2不属于差异表达最大的基因之一，但用GC处理的ASM细胞与基线条件的所有比较表明，CRISPLD2具有显著调整的P值。具体而言，GSE34313研究发现，在用DEX处理ASM细胞后4小时和24小时，CRISPLD2均有差异表达，而GSE13168研究发现，当用GC(即氟替卡松)处理ASM细胞时，CRISPLD2的差异表达比用促炎细胞因子(即EGF和IL1b)刺激细胞时最强

image.png

2.KEGG（找差异表达基因主要显著影响了哪些「生化代谢途径和信号转导途径」。）

KEGG，除了对基因本身功能的注释，我们也知道基因会参与人体的各个通路，基于人体通路而形成的数据库就是通路相关的数据库。而KEGG就是通路相关的数据库的一种,大多数听说过KEGG的人都会把它当做一个基因通路（Pathway）的数据库，其实人家的功能远不止于此。KEGG是一个整合了基因组、化学和系统功能信息的综合数据库。KEGG下属4个大类和17和子数据库，而其中有一个数据库叫做KEGG Pathway，专门存储不同物种中基因通路的信息，也是用的最多的一个，所以，久而久之，KEGG就被大家当做是一个通路数据库了.

KEGG pathway 分析和上面介绍的GO分析是一样的只是把enrichGO()函数改成 enrichKEGG()

> kk <- enrichKEGG(gene = gene.name$ENTREZID,organism = 'hsa',  pvalueCutoff = 0.1)
对前10条通路可视化
> ggplot(head(kk@result,10),aes(x=GeneRatio,y=Description))+
labs(title="The KEGG enrichment analysis of all DEGs")+
geom_point(aes(size = `Count`, color = `pvalue`), shape = 16)+  scale_size(range = c(1, 10))+scale_color_gradient(low = "#cf0000",high = "#f5c6c6")
或者条形图
barplot(kk,showCategory = 25, title="The KEGG enrichment analysis of all DEGs")+
scale_size(range=c(2, 12))+
scale_x_discrete(labels=function(kk) str_wrap(kk,width = 25))
或者点图
dotplot(kk,showCategory = 25, title="The KEGG enrichment analysis of all DEGs")+
scale_size(range=c(2, 12))+
scale_x_discrete(labels=function(kk) str_wrap(kk,width = 25))

分别与苯丙氨酸代谢，寿命调节途径-多个物种，昼夜节律夹带，安非他明成瘾信号通路相关

image.png
或者使用DAVID工具分析使用DAVID进行GO/KEGG分析 - 简书 (jianshu.com)

GO分析：enrichGO() —---—> KEGG pathway分析：enrichKEGG()

3. GO、KEGG和富集分析有什么关系呢？

通过上面的解释，我们知道，其实GO和KEGG是两个数据库，里面有每个基因相关的功能信息，而富集分析就是一个把这些功能进行进行整合计算的算法。

GO和KEGG是基础，而富集是过程，最后得到的结果就是整合后的宏观的结果。