GO、KEGG富集分析（一）有参情况

2020-08-13 本文已影响0人大号在这里

一、GO分析的理论知识

what is Gene Ontology（GO）? 基因"本体论"
Gene Ontology中最基本的概念是 term。GO里面的每一个entry都有一个唯一的数字标记，形如 GO: nnnnnnn，还有一个term名，比如 "cell", "fibroblast growth factor receptor binding"，或者 "signal transduction"。每个term都属于一个ontology，总共有三个ontology，它们分别是molecular function, cellular component和biological process。

对基因的描述一般从三个层面进行：

Cellular component，CC 细胞成分
Biological process， BP 生物学过程
Molecular function，MF 分子功能

这三个层面具体是指：

Cellular component：解释的是基因存在在哪里，在细胞质还是在细胞核？如果存在细胞质那在哪个细胞器上？如果是在线粒体中那是存在线粒体膜上还是在线粒体的基质当中？这些信息都叫Cellular component。
Biological process：是在说明该基因参与了哪些生物学过程，比如，它参与了rRNA的加工或参与了DNA的复制，这些信息都叫Biological process
Molecular function：在讲该基因在分子层面的功能是什么？它是催化什么反应的？

得到GO注释

模式生物 ---> 有标准的注释数据库；
非模式生物 ---> 自己搜注释数据库(怎们搜后面具体介绍)，搜不到就用blast的办法解决。

做GO分析的思路：

gene expression distribution --->control VS treatment ---> DEG ---> GO enrichment analysis

也就是RNA-Seq先测出各组的基因表达分布：
通过比较 control 和 treatment 得到差异表达基因
再去做GO富集分析：
用找到的差异基因去做GO富集分析，希望能从这三方面找到和我们背景不一样的地方。

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

GO富集分析原理：
有一个term注释了100个差异表达基因参与了哪个过程，注释完之后（模式生物都有现成的注释包，不用我们自己注释），计算相对于背景它是否显著集中在某条通路、某一个细胞学定位、某一种生物学功能。

二、使用clusterProfiler进行GO/KEGG富集分析

clusterProfiler是一个功能强大的R包，同时支持GO和KEGG的富集分析，而且可视化功能非常的优秀，本章主要介绍利用这个R包来进行Gene Ontology的富集分析。
进行GO分析时，需要考虑的一个基础因素就是基因的GO注释信息从何处获取。Bioconductor上提供了以下19个物种的Org类型的包，包含了这些物种的GO注释信息

packages    organism
org.Ag.eg.db    Anopheles
org.At.tair.db  Arabidopsis
org.Bt.eg.db    Bovine
org.Ce.eg.db    Worm
org.Cf.eg.db    Canine
org.Dm.eg.db    Fly
org.Dr.eg.db    Zebrafish
org.EcK12.eg.db E coli strain K12
org.EcSakai.eg.db   E coli strain Sakai
org.Gg.eg.db    Chicken
org.Hs.eg.db    Human
org.Mm.eg.db    Mouse
org.Mmu.eg.db   Rhesus
org.Pf.plasmo.db    Malaria
org.Pt.eg.db    Chimp
org.Rn.eg.db    Rat
org.Sc.sgd.db   Yeast
org.Ss.eg.db    Pig
org.Xl.eg.db    Xenopus

对于以上19个物种，只需要安装对应的org包，clusterProfile就会自动从中获取GO注释信息，我们只需要差异基因的列表就可以了，使用起来非常方便。

clusterProfiler包的安装

#bioconductor 安装
#install.packages('BiocManager')  #需要首先安装 BiocManager，如果尚未安装请先执行该步
BiocManager::install('clusterProfiler')

1. GO富集分析

1.1 准备输入数据
待分析的数据就是一串基因名称了，可以是ensembl id、entrze id或者symbol id等类型都可以。把基因名称以一列的形式排开，放在一个文本文件中（例如命名“gene.txt”）。Excel中查看，就是如下示例这种样式。

1.2 加载参考物种的基因注释数据库

#（1）对于常见物种，例如人类，有些专门的 R 包数据库，例如人类参考基因组 hg19 的
library(org.Hs.eg.db)
#（2）对于不常见的物种，但却是存在参考基因组的情况
#通过 AnnotationHub 包索引基因组，例如期望找绵羊（Ovis aries）的注释库
library(AnnotationHub)
hub <- AnnotationHub()
query(hub, 'Ovis aries')  #输入绵羊（Ovis aries）的名称进行匹配
sheep <- hub[['AH72269']]  #返回了数据库编号 AH72269，就可以加载该库

1.3 GO富集分析
加载了注释库之后，读取基因列表文件，并使用clusterProfiler的内部函数enrichGO()即可完成GO富集分析。

library(clusterProfiler)
#读取基因列表文件中的基因名称
genes <- read.delim('gene.txt', header = TRUE, stringsAsFactors = FALSE)[[1]]

#GO富集分析
#对于加载的注释库的使用，以上述为例，就直接在 OrgDb 中指定人（org.Hs.eg.db）或绵羊（sheep）
enrich.go <- enrichGO(gene = genes,  #基因列表文件中的基因名称
    OrgDb = 'sheep',  #指定物种的基因数据库，示例物种是绵羊（sheep）
    keyType = 'ENTREZID',  #指定给定的基因名称类型，例如这里以 entrze id 为例
    ont = 'ALL',  #可选 BP、MF、CC，也可以指定 ALL 同时计算 3 者
    pAdjustMethod = 'fdr',  #指定 p 值校正方法
    pvalueCutoff = 0.05,  #指定 p 值阈值，不显著的值将不显示在结果中
    qvalueCutoff = 0.2,  #指定 q 值阈值，不显著的值将不显示在结果中
    readable = FALSE)

#例如上述指定 ALL 同时计算 BP、MF、CC，这里将它们作个拆分后输出
BP <- enrich.go[enrich.go$ONTOLOGY=='BP', ]
CC <- enrich.go[enrich.go$ONTOLOGY=='CC', ]
MF <- enrich.go[enrich.go$ONTOLOGY=='MF', ]

write.table(as.data.frame(BP), 'go.BP.txt', sep = '\t', row.names = FALSE, quote = FALSE)
write.table(as.data.frame(CC), 'go.CC.txt', sep = '\t', row.names = FALSE, quote = FALSE)
write.table(as.data.frame(MF), 'go.MF.txt', sep = '\t', row.names = FALSE, quote = FALSE)

GO富集分析结果表格

2. KEGG富集分析

读取基因列表文件，并使用clusterProfiler的内部函数enrichKEGG()即可完成KEGG富集分析。

library(clusterProfiler)

#读取基因列表文件中的基因名称，注意这里只能用 entrze id
genes <- read.delim('gene.txt', header = TRUE, stringsAsFactors = FALSE)[[1]]

#每次打开R计算时，它会自动连接kegg官网获得最近的物种注释信息，因此数据库一定都是最新的
kegg <- enrichKEGG(
    gene = genes,  #基因列表文件中的基因名称
    keyType = 'kegg',  #kegg 富集
    organism = 'oas',  #例如，oas 代表绵羊，其它物种更改这行即可
    pAdjustMethod = 'fdr',  #指定 p 值校正方法
    pvalueCutoff = 0.05,  #指定 p 值阈值，不显著的值将不显示在结果中
    qvalueCutoff = 0.2,  #指定 q 值阈值，不显著的值将不显示在结果中
)
#输出结果
write.table(kegg, 'kegg.txt', sep = '\t', quote = FALSE, row.names = FALSE)

3. clusterProfiler的作图功能

此外，clusterProfiler中也额外提供了一系列的可视化方案用于展示本次富集分析结果，具有极大的便利。

#clusterProfiler 包里的一些默认作图方法，例如
barplot(ego, showCategory = 10)#散点图
dotplot(ego, showCategory = 10)#气泡图
plotGOgraph(ego)#无环图1
goplot(ego)#无环图2
emapplot(ego, showCategory = 30)#对于富集到的GO terms之间的基因重叠关系进行展示，如果两个GO terms系的差异基因存在重叠，说明这两个节点存在overlap关系，在图中用线条连接起来
cnetplot(ego, showCategory = 5)#对于基因和富集的GO terms之间的对应关系进行展示，如果一个基因位于一个GO Terms下，则将该基因与GO连线

参考：
https://www.jianshu.com/p/47b5ea646932?utm_source=desktop&utm_medium=timeline